function model = Fitting(handles)
% fit specs is a handle structure bearing details of the fit

disp('This program fits a simple binomial model to a given set of data,')
disp('and tests whether the resulting fit is adequate.')
disp('Author Jenny Read; this version May 2008.')
disp('Any problems, email j.c.a.read@ncl.ac.uk')
disp('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&')
disp(' ')


% NB this is a daft kludge. Unless I use NLLF_gaussian outside an if-loop, it doesn't get compiled for standalone mode, and i don't seem to be able to get around this with a pragma
%negloglikelihood = NLLF_gaussian(1,data);

% Open file for writing results of fitting 
resultsfile = [handles.workdir '\results.' handles.fitspecs.datafile];
fid=fopen(resultsfile,'a');
if fid==-1
    errordlg(['Unexpected error: Cannot open the file ' resultsfile ' for writing results.'])
    confirmandexit
end
datafile = [handles.workdir '\' handles.fitspecs.datafile] ;
fprintf(fid,'%s %s %s\n','Fitting data from ' ,datafile, '********************************************');
fprintf(fid,'%s\n','Fit specification:');

% Read off fitting params from handles.fitspecs.
n = handles.fitspecs.n;
p = handles.fitspecs.p;
q = handles.fitspecs.q;
qsigmapercentage = handles.fitspecs.qs;
QVtype = handles.fitspecs.qvtype;
noise = handles.fitspecs.noise;
pstim = handles.fitspecs.pcond;
offset = handles.fitspecs.offset;
glitchrate = handles.fitspecs.glitchrate;
glitchmin = handles.fitspecs.glitchmin;
glitchmax = handles.fitspecs.glitchmax;
fittype = handles.fitspecs.fittype;
data = handles.fitspecs.data;
% NB some of my functions need data to be a row not a column vector. Next few lines ensure this:
[tmp1,tmp2] = size(data);
if tmp1>tmp2, data = data', end
options = optimset('MaxFunEvals',handles.fitspecs.maxnevals);

disp('******************************************************************************')

% Record fit constraints to results file
if isnan(n)
    fprintf(fid,'%s\n','n is unconstrained');
else
    fprintf(fid,'%s %d\n','n = ',n);
end
if isnan(p)
    fprintf(fid,'%s\n','p is unconstrained');
else
    fprintf(fid,'%s %d\n','p = ',p);
end
if isnan(q)
    fprintf(fid,'%s\n','q is unconstrained');
else
    fprintf(fid,'%s %d\n','q = ',q);
end
if isnan(qsigmapercentage)
    fprintf(fid,'%s\n','sigma_Q is unconstrained');
else
    fprintf(fid,'%s %d %s\n','sigma_Q = ',qsigmapercentage,'% of q');
end
if isnan(QVtype)
    fprintf(fid,'%s\n','QV type is unconstrained');
elseif QVtype==1
    fprintf(fid,'%s\n','QV type is constrained to Type I');
elseif QVtype==2
    fprintf(fid,'%s\n','QV type is constrained to flat (~ mixed Type I/II)');
end
if isnan(noise)
    fprintf(fid,'%s\n','noise is unconstrained');
else
    fprintf(fid,'%s %d\n','noise = ',noise);
end
if isnan(pstim)
    fprintf(fid,'%s\n','pcond is unconstrained');
else
    fprintf(fid,'%s %d\n','pcond = ',pstim);
end
fprintf(fid,'Glitch rate is constrained to %f\n',glitchrate);
fprintf(fid,'Fit seeks to maximise %s\n',handles.fitspecs.maximand);
 
% If n is free, try all n between 1 and the maximum
if isnan(n)
    n_totry = [1:handles.fitspecs.maxn];
else
    n_totry = n;
end
% If QVtype is free, try both Type I and flat.
if isnan(QVtype)
    QVtype_totry = [1 2];
else
    QVtype_totry = QVtype;
end


datamean=mean(data);
dataSD=std(data);

initialparams=0;
params=0;
set(handles.text37,'visible','on','string','Current values');
set(handles.text38,'visible','on','string','Current values');
set(handles.NLLFlabel,'visible','on','string',sprintf('%s for current best fit = ',handles.fitspecs.maximand));
bestNLLF_ever = inf;
fprintf(fid,'%10s %10s %10s %10s %10s %10s %10s %10s %10s %10s %30s\n'...
    ,'n','p','q','qsig','qsig as %','QVtype','noise','pstim','offset','glitchrate',[handles.fitspecs.maximand ' (being maximised)']);
for thisQVtype = QVtype_totry
    for thisn = n_totry
        % Keep track of which parameters are to be varied by the optimiser.
        % n and QVtype are not (they are discrete)
        givenparams = [thisn p q qsigmapercentage thisQVtype noise pstim offset glitchrate fittype glitchmin glitchmax];
        % All the NaNs in the vector given params represent values which are to be varied.
        tobevaried = find(isnan(givenparams));
        thisNLLF=inf;
        % Now optimise
        % Run it several times with different initial conditions
        for initial=1:handles.fitspecs.nstartingpoints
            % First find good initial conditions - keep picking them at random until LLF is non-zero
            NLLF=inf;
            
            while isinf(NLLF)
                clear initialparams params
                % params is a vector containing [n p q qsigmapercent QVtype noise pstim offset glitchrate fittype]
                % Go through params vector, randomly initialising the unconstrained variables
                [newrandomparams,changedparams] = RandomiseNaN(givenparams,datamean,dataSD);
                % See whether this gives a finite NLLF
                NLLF = NLLF_optall(changedparams,data,newrandomparams,tobevaried);
            end
            [optimisedparams,NLLF,exitstatus] = fminsearch('NLLF_optall',changedparams,options,data,newrandomparams,tobevaried);
            % Now reconstruct the whole parameter vector
            params = newrandomparams;
            params(tobevaried) = optimisedparams;
            if (exitstatus<=0)
                disp('Fit failed')
                NLLF = 9999999999999999;
            end
            if (NLLF<thisNLLF) % if it is the best of this n and QVtype.........
                thisNLLF = NLLF;
                bestparams_this=params;
            end
            
            fprintf('Seeking best fit ... %s currently %f\n', handles.fitspecs.maximand, -NLLF )
            % Update figure to let user know what is happening
            if isfield(handles,'nbest')
                set(handles.nbest,'string',num2str(thisn),'visible','on');
                set(handles.pbest,'string',num2str(params(2),'%5.3f'),'visible','on');
                set(handles.qbest,'string',num2str(params(3),'%5.3f'),'visible','on');
                set(handles.qsbest,'string',num2str(params(4),'%3.0f'),'visible','on');
                if thisQVtype == 1
                set(handles.qvbest,'string','I','visible','on');
                elseif thisQVtype == 2
                    set(handles.qvbest,'string','flat','visible','on');
                end
                set(handles.noisebest,'string',num2str(abs(params(6)),'%5.3f'),'visible','on');
                set(handles.pcondbest,'string',num2str(params(7),'%5.3f'),'visible','on');
                set(handles.offsetbest,'string',num2str(params(8),'%5.3f'),'visible','on');
                set(handles.text23,'string',num2str(-NLLF),'visible','on');
            end
            drawnow        
        end  % choose a new set of initial conditions
        if (thisNLLF<bestNLLF_ever) % if it is the best ever encountered......
            bestNLLF_ever = thisNLLF;
            bestparams_ever=bestparams_this;
        end
        
        fprintf(fid,'%10d %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f %10.2f %30.6f\n',...
            bestparams_this(1),bestparams_this(2),bestparams_this(3),bestparams_this(4)*bestparams_this(3)/100,...
            bestparams_this(4),bestparams_this(5),...
            bestparams_this(6),bestparams_this(7),bestparams_this(8),glitchrate,-thisNLLF);
        % Cut short the search over n if it looks like we're getting steadily worse:
        if bestparams_ever(1) + 4 < bestparams_this(1)
            thisn=handles.fitspecs.maxn;
            break
        end
    end % next n 
end % next QVtype
model.n = bestparams_ever(1);
model.p = bestparams_ever(2);
model.q = bestparams_ever(3);
model.qsigma = bestparams_ever(4)*bestparams_ever(3)/100;
model.QVtype = bestparams_ever(5);
model.noise = bestparams_ever(6);
model.pstim = bestparams_ever(7);
model.offset = bestparams_ever(8);
model.glitchrate = glitchrate;
model.glitchmin = glitchmin;
model.glitchmax = glitchmax;
model.fittype = fittype;

CheckParamsAreValid(model);

fprintf(fid,'%s\n','********************************** Best model encountered:');
fprintf(fid,'%10d %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f %10.2f %30.6f\n',...
    model.n,model.p,model.q,model.qsigma,bestparams_ever(4),model.QVtype,...
    model.noise,model.pstim,model.offset,model.glitchrate,-bestNLLF_ever);
fprintf(fid,'\n');
fclose(fid)


disp('Best model found within the required constraints : ')
disp(['Negative log likelihood = -log(L) = ' num2str(bestNLLF_ever)])


