function model= dbnFit(X, numhid, y, varargin)
%fit a DBN to bianry data in X

%INPUTS: 
%X              ... data. should be binary, or in [0,1] interpreted as
%               ... probabilities
%numhid         ... list of numbers of hidden units
%y              ... List of discrete labels

%OUTPUTS:
%model          ... A cell array containing models from all RBM's

%varargin may contain options for the RBM's of this DBN, in row one by one
%for example:
%dbnFit(X, [500,400], opt1, opt2) uses opt1 for 500 and opt2 for 400
%dbnFit(X, [500,400], opt1) uses opt1 only for 500, and defaults for 400

numopts=length(varargin);
H=length(numhid);
model=cell(H,1);
cant_clases = length(unique(y));

if H>=2
    
    %train the first RBM on data
    if(numopts>=1)
        model{1}= rbmBB(X, numhid(1),varargin{1});
    else
        model{1}= rbmBB(X, numhid(1));
    end
    
%     %extra-debug
%     visualize(model{1}.W);
%     drawnow;
    
    %train all other RBM's on top of each other
    for ii=2:H
        if(numopts>=ii)
            model{ii}=rbmBB(model{ii-1}.top, numhid(ii), varargin{ii});
        else
            model{ii}=rbmBB(model{ii-1}.top, numhid(ii));
        end
        
%         %extra-debug
%         visualize(model{ii}.W);
%         drawnow;
        
    end

%     %the last RBM has access to labels too
%     if(numopts>=H)
%         model{H}= rbmFit(model{H-1}.top, numhid(end), y, varargin{H});
%     else
%         model{H}= rbmFit(model{H-1}.top, numhid(end), y);
%     end
    model{H}.Wc = 0.1*randn(cant_clases, size(model{H}.W,2) );
    model{H}.cc = 0.1*randn(1, cant_clases);
    model{H}.labels = 1:cant_clases;

    
    %fine-tuning of the DBN
    
    targets = zeros(length(y), cant_clases);
    for ii = 1:cant_clases
        targets(y == ii, ii) = 1;
    end
    
    for epoch = 1:200
    
        %%%%%%%%%%%%%%% PERFORM CONJUGATE GRADIENT WITH 3 LINESEARCHES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        max_iter=3;

        if epoch<6  % First update top-level weights holding other weights fixed. 
%             N = size(data,1);
%             XX = [data ones(N,1)];
%             w1probs = 1./(1 + exp(-XX*w1)); w1probs = [w1probs  ones(N,1)];
%             w2probs = 1./(1 + exp(-w1probs*w2)); w2probs = [w2probs ones(N,1)];
%             w3probs = 1./(1 + exp(-w2probs*w3)); %w3probs = [w3probs  ones(N,1)];

            testdata = X;
            for ii=1:length(model)
                testdata= rbmVtoH(model{ii}, testdata);
            end
            
            VV = colvec([ model{H}.Wc'; model{H}.cc ]);
            Dim = [ numhid(end); cant_clases];
            XX = minimize(VV,'CG_CLASSIFY_INIT',max_iter,Dim,testdata,targets);
            aux = reshape(XX, numhid(end)+1, cant_clases);
            model{H}.Wc = aux(1:end-1,:)';
            model{H}.cc = aux(end,:);

        else
            
%             VV = [w1(:)' w2(:)' w3(:)' w_class(:)']';
%             Dim = [l1; l2; l3; l4; l5];
            VV = [];
            for ii=1:H
                VV = [VV; colvec([model{ii}.W; model{ii}.b ]) ];
            end
            VV = [VV; colvec([model{ii}.Wc'; model{ii}.cc ])];
            
            Dim = [ size(X,2); colvec(numhid); cant_clases];

            XX = minimize(VV,'CG_CLASSIFY',max_iter,Dim,X,targets);

            xxx = 0;
            for ii=1:length(model)
                aux = reshape(XX(xxx+1:xxx+(Dim(ii)+1)*Dim(ii+1)),Dim(ii)+1,Dim(ii+1));
                model{ii}.W = aux(1:end-1,:);
                model{ii}.b = aux(end,:);
                xxx = xxx+(Dim(ii)+1)*Dim(ii+1);
            end
            aux = reshape(XX(xxx+1:xxx+(Dim(ii+1)+1)*Dim(ii+2)),Dim(ii+1)+1,Dim(ii+2));
            model{ii}.Wc = aux(1:end-1,:)';
            model{ii}.cc = aux(end,:);
            
%             w1 = reshape(X(1:(l1+1)*l2),l1+1,l2);
%             xxx = (l1+1)*l2;
%             w2 = reshape(X(xxx+1:xxx+(l2+1)*l3),l2+1,l3);
%             xxx = xxx+(l2+1)*l3;
%             w3 = reshape(X(xxx+1:xxx+(l3+1)*l4),l3+1,l4);
%             xxx = xxx+(l3+1)*l4;
%             w_class = reshape(X(xxx+1:xxx+(l4+1)*l5),l4+1,l5);

        end
        %%%%%%%%%%%%%%% END OF CONJUGATE GRADIENT WITH 3 LINESEARCHES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        prediction = dbnPredict(model, X);
        errors = sum(prediction~=y);
        fprintf(1, 'Errores %d\n', errors);
        
    end
    
    
%     %extra-debug
%     visualize(model{H}.W);
%     drawnow;
    
else
    
    %numhid is only a single layer... but we should work anyway
    if (numopts>=1)
        model{1}= rbmFit(X, numhid(1), y, varargin{1});
    else
        model{1}= rbmFit(X, numhid(1), y);
    end
end