function [meanError, meanErrorDRF, meanErrorDumb, meanErrorVeryDumb, ...
            meanError1, meanErrorDRF1, meanErrorDumb1, meanErrorVeryDumb1] = ...
        trainData(trainOrTest, displayFlag, dimensions, scales, trainingType, modelType, smoothingCoeff, DRFFlag, ...
        specialFilename)
% this trains the parameters
% Gaussian model, to train using Laplacian model uncomment the liprog lines,
% linprog requires significant amount of RAM.
%  and also see the trainDataRelative_someotherversion.m, etc.


% Version 0.1, Feb 2006.
% This code is copyrighted by Ashutosh Saxena, Stanford University
% It is available for non-commercial use only. For commericial use, please
% contact the authors.
% Any publication or report resulting from Use of this code, should cite:
% Learning Depth from Single Monocular Images, Ashutosh Saxena, Sung H. Chung, Andrew Y. Ng, NIPS 2005.
% 3-D Depth Reconstruction from a Single Still Image, Ashutosh Saxena, Sung H. Chung, Andrew Y. Ng, To appear in IJCV 2007.
% For more info and dataset, visit: http://ai.stanford.edu/~asaxena/learningdepth/
% For a better version of the algorithm, also visit: 
% http://ai.stanford.edu/~asaxena/reconstruction3d/



global nGridRow nGridCol maxDistance minDistance
global featureDirectory
global nDim nLaw nStatistics
global nThetaStep

upsilon = 1e-22;

if nargin < 9
    specialFilename = '';
end
configurationFile;
if nargin < 5
    trainingType = 'gaussian';
end
trainingType = lower(trainingType);
if nargin < 6
    modelType = 'absolute';
end
if nargin < 7
    smoothingCoeff = [1 .01 .01];
end
if nargin < 8
    %if strcmp(trainingType, 'gaussian') && strcmp(modelType, 'absolute')
    %    DRFFlag = true;
    %else
        DRFFlag = false;
        %end
end
if nargin < 4 || scales(1) == 0
    scales = [1 2 5 -1];
end
if nargin < 3 %|| dimensions == -1
    if sum(scales == -1) <= 0          % column feature is not appended with neighbors
        dimensions = 1:((nDim*5)*length(scales));
    else
        dimensions = 1:((nDim*5)*(length(scales)-1)+nDim);
    end
end
if nargin < 2
    displayFlag = false;
end
if nargin < 1
    trainOrTest = 1;        % 2 for test
end

if length(dimensions) == 1 && dimensions == -1
    %H2 only
    dimensions = [];
    dimensionsH1 = [];
    for d = 0:(5*sum(scales>0)-1)
        dimensions = [dimensions nDim*d+(1:nLaw) nDim*d+2*nLaw+(1:nLaw)];
        dimensionsH1 = [dimensionsH1 nDim*d+2*nLaw+(1:nLaw)];
    end
    if sum(scales == -1) > 0
        flagColumn = 1;
        dimensions = [dimensions 5*nDim*sum(scales>0)+(1:4*nLaw)];
        dimensionsH1 = [dimensionsH1 5*nDim*sum(scales>0)+(1:4*nLaw)];
    else
        flagColumn = 0;
    end
    if sum(scales == -2) > 0
        dimensions = [dimensions 5*nDim*sum(scales>0)+flagColumn*4*nLaw+(1:4*nLaw)];
        dimensionsH1 = [dimensionsH1 5*nDim*sum(scales>0)+(1:4*nLaw)];
    end
elseif length(dimensions) == 1 && dimensions == -2
    dimensions = [1:15 46:60 91:105 136:150 181:195];
elseif length(dimensions) == 1 && dimensions == -3
    dimensions = [31:45 76:90 121:135 166:180 211:225];
end

totalImages = 425;%325;    
imageNStep = 25;

nAppendedDim = length(dimensions);
if sum(scales == -1) <= 0
    nFeatures = nAppendedDim / length(scales) / 5;
else
    nFeatures = (nAppendedDim-nDim) / (length(scales)-1) / 5;
end

tmp = 1:totalImages;        %[1:6];%4 6:7];
if trainOrTest == 1
    dataSetsNumber = tmp(find(mod(tmp,3) > 0));
elseif trainOrTest == 2
    dataSetsNumber = tmp(find(mod(tmp,3) == 0));
end
nImages = length(dataSetsNumber);

%weights = zeros(2*nAppendedDim+1,nGridRow);
%absoluteVarianceWeights = zeros(2*nAppendedDim+1,nGridRow);       %1 for bias

for row = 1:nThetaStep:nGridRow        %neglecting first row as  features are heavily copied
    for n=1:imageNStep:(totalImages-imageNStep+1)
        for callRow = 1:nThetaStep
            rowDepthMaps(callRow,:,:,n-1+(1:imageNStep)) = makeAppendedDepths(min(nGridRow,row+callRow-1),n,scales);
        end
    end

    rowDepthMaps = rowDepthMaps(:,:,:,dataSetsNumber);
    actualDepthMap(min(nGridRow,row-1+(1:nThetaStep)),:,:) = permute( rowDepthMaps(:,:,1,:), [1 2 4 3]);  % the actual depth map
    
    depthVector = double(reshape( permute(rowDepthMaps(:,:,1,:), [1 2 4 3]), [nThetaStep*nGridCol*nImages 1]) );
    clear rowDepthMaps;
    
    %=======Loading feature vectors===========
    fV = zeros(nThetaStep,nGridCol,5*sum(scales>0)*nDim+nLaw*4*sum(scales<0),totalImages);
    for n=1:imageNStep:(totalImages-imageNStep+1)
        for callRow = 1:nThetaStep
            fV(callRow,:,:,n-1+(1:imageNStep)) = makeAppendedFeatures(min(nGridRow,row+callRow-1), n, scales);%, dimensions);   
        end
    end
    disp('Feature Vector loaded......');
    
    fV = fV(:,:,:,dataSetsNumber);
    %permuting to take "dimensions" in 2, and nImages*nCols into 1
    fV = reshape( permute( fV, [1 2 4 3]), [nThetaStep*nGridCol*nImages size(fV,3)]);
    dumbFeatureVector = ones(size(fV,1),1);
    
    if strcmp(trainingType, 'gaussian') == 1
        featureVector = [ones(size(fV,1),1), fV(:,dimensions) ];%, fV(:,[1:nLaw:length(dimensions)]).^2 ];%, ...
    elseif strcmp(trainingType, 'laplacian') == 1
%        featureVector = fV(:,dimensionsH1);%, fV(:,[1:nLaw:length(dimensions)]).^2 ];%, ...
        featureVector = fV(:,dimensions);%, fV(:,[1:nLaw:length(dimensions)]).^2 ];%, ...
        if trainOrTest == 2
            featureVector = [ ones(size(featureVector,1),1), featureVector];
        end
    end
    % the first one is most effective, so introducing square of the first
    % one (L3'*L3) for  every neighbor and scale
    
    clear fV;
   
    if trainOrTest == 1
            dumbWeights(:,row) = pinv(dumbFeatureVector) * depthVector;
            dumbDepthVector(:,row) = dumbFeatureVector * dumbWeights(:,row);
            disp([num2str([row nGridRow]) '  <<<<< Training with ' num2str(nThetaStep) ' rows * ' num2str(nImages) ' images * ' ...
                    num2str(nGridCol) ' cols = ' num2str(nThetaStep*nImages*nGridCol) ' samples; in ' ...
                    num2str( size(featureVector,2)) ' dimensions....']);
            
            if strcmp(trainingType, 'gaussian') == 1
                invFeatureVector = pinv(featureVector);
                %weights(:,row) = featureVector \ depthVector;
                weights(:,row) = invFeatureVector * depthVector;
                predictedDepthVector(:,row) = featureVector * weights(:,row);
                absoluteVarianceVector = (depthVector - predictedDepthVector(:,row)) .^ 2;
            
                absoluteVarianceWeights(:,row) = invFeatureVector * absoluteVarianceVector;
            elseif strcmp(trainingType, 'laplacian') == 1
                disp('Robust fitting using Huber.......');
                %invFeatureVector = pinv(featureVector);
                weights(:,row) = robustfit(featureVector, depthVector, 'huber');
                featureVector = [ ones(size(featureVector,1),1), featureVector];
                %weights(:,row) = invFeatureVector * depthVector;
                
                % nConstraints = size(featureVector,1);
                % nVariables = size(featureVector,2);
                %tmp = linprog([repmat(0,nVariables,1); repmat(1,nConstraints,1)], ...
                %              [sparse(featureVector), -speye(nConstraints); -sparse(featureVector), -speye(nConstraints)], ...
                %              [depthVector; -depthVector] );
                %weights(:,row) = tmp(1:nVariables);
         
                predictedDepthVector(:,row) = featureVector * weights(:,row);
                absoluteVarianceVector = abs(depthVector - predictedDepthVector(:,row));
            
                absoluteVarianceWeights(:,row) = featureVector \ absoluteVarianceVector;
                
            end
            
            actualAbsoluteVariance(row) = mean( absoluteVarianceVector );
            predictedAbsoluteVarianceVector(:,row) = featureVector * absoluteVarianceWeights(:,row);
            actualVarianceinVariance(row) = mean( (absoluteVarianceVector-predictedAbsoluteVarianceVector(:,row)).^2 );
            %warning on;
            clear invFeatureVector;
            
            residualAbsError(row) = sum(absoluteVarianceVector) / sum( (depthVector - mean(depthVector) ).^2 );
            residualAbsVarianceError(row) = sum( (absoluteVarianceVector-predictedAbsoluteVarianceVector(:,row)).^2) / ...
                        sum( (absoluteVarianceVector-mean( absoluteVarianceVector)).^2);
    
            disp( num2str([residualAbsError(row) residualAbsVarianceError(row)] ));
                    
            if displayFlag
                disp('Absolute depths......');
                figure(200),
                subplot(3,1,1); hist( depthVector, [0:.05:4.5]);
                subplot(3,1,2); hist( depthVector - featureVector * weights(:,row),[-4.5:.05:4.5]);% ./ mean((depthVector - mean(depthVector) ).^2), 100);
                subplot(3,1,3), hist( featureVector * weights(:,row), [0:.05:4.5]);
                mean( absoluteVarianceVector) % / sum( depthVector.^2)
                mean( depthVector )
            end
            

            %elseif strcmp(trainingType, 'laplacian')
            %disp('Starting linear program...');
            %size(featureVector);
            %nAppendedDim;
            %tmp = linprog(  [repmat(0,nAppendedDim+1,1); repmat(1,nImages*nGridCol,1)], ...
            %                [sparse(featureVector), -speye(nImages*nGridCol); -sparse(featureVector), -speye(nImages*nGridCol)], ...
            %                [depthVector; -depthVector] );
            %weights(:,row) = tmp(1:(nAppendedDim+1));
            %
            %tmp = linprog(  [repmat(0,1,1); repmat(1,nImages*nGridCol,1)], ...
            %                [sparse(dumbFeatureVector), -speye(nImages*nGridCol); -sparse(dumbFeatureVector), -speye(nImages*nGridCol)], ...
            %                [depthVector; -depthVector] );
            %dumbWeights(:,row) = tmp(1:(nAppendedDim+1));
            %end
        
            %if DRFFlag
            %varianceWeightsX(:,row) = quadprog(diffXFeatureVector'*diffXFeatureVector, -diffXFeatureVector'*varianceVectorX, ...
            %                            -speye(size(diffXFeatureVector,2)), zeros(size(diffXFeatureVector,2),1) );
            %varianceWeightsY(:,row) = quadprog(diffYFeatureVector'*diffYFeatureVector, -diffYFeatureVector'*varianceVectorY, ...
            %                            -speye(size(diffYFeatureVector,2)), zeros(size(diffYFeatureVector,2),1) );
    elseif trainOrTest == 2
        if mod(row,5) == 0
            disp([num2str([row nGridRow]) '  >>>>>>>> Testing on ' num2str( length(dataSetsNumber) ) ' images, for ' num2str(nGridRow) ' rows x ' ...
                num2str(nGridCol) ' columns depthmap resolution; in ' num2str( size(featureVector,2)) ' dimensions....']);
        end
        load(['weightsAbsoluteDepths_' trainingType '.mat']);
        predictedDepthVector(:,row) = featureVector * weights(:,row);
        predictedAbsoluteVarianceVector(:,row) = featureVector * absoluteVarianceWeights(:,row);
        dumbDepthVector(:,row) = dumbFeatureVector * dumbWeights(:,row);
        clear featureVector;
        actualAbsoluteVariance = 0;     residualAbsError = 0;
        actualVarianceinVariance = 0;   residualAbsVarianceError = 0;
        
        varianceVector = ( depthVector - predictedDepthVector(:,row) ).^2;
       
        residualAbsError(row) = sum(varianceVector) / ...
                                        sum( (depthVector - mean(depthVector) + upsilon).^2 );
        disp( num2str(residualAbsError(row) ));
    end

    
    if mod(row-1,10) == 0
        if trainOrTest == 1
            save(['weightsAbsoluteDepths_' trainingType '.mat'], 'weights', 'dumbWeights', 'absoluteVarianceWeights');
        end

        filename = ['predictedAbsoluteDepthVectorsEarly' num2str(trainOrTest) '_' trainingType '.mat'];
        save(filename, 'row', 'predictedDepthVector', 'actualDepthMap', 'dumbDepthVector', ...
                'actualAbsoluteVariance', 'residualAbsError', 'predictedAbsoluteVarianceVector', ...
                'actualVarianceinVariance', 'residualAbsVarianceError');
    end
end

if trainOrTest == 1
    save(['weightsAbsoluteDepths_' num2str(trainingType) '.mat'], 'weights', 'dumbWeights', 'absoluteVarianceWeights');
end

filename = ['predictedAbsoluteDepthVectorsEarly' num2str(trainOrTest) '_' trainingType '.mat'];
save(filename, 'predictedDepthVector', 'actualDepthMap', 'dumbDepthVector', ...
                'actualAbsoluteVariance', 'residualAbsError', 'predictedAbsoluteVarianceVector', ...
                'actualVarianceinVariance', 'residualAbsVarianceError');