function [results,biftfig]=biftfull(input1,Full,NoGui,biftfig,BiftParams)
% C***********************************************************************
%
% results=biftfull(input1,Full,NoGui,biftfig,BiftParams)
%
% Input parameter input1 is optional and will be either a file name
% or a matrix containing the data (x,y,Delta(y))
%
% Full is a parameter which indicates how much of the routine should be
% executed
%    Full=0 is that only a standard downhill simplex to find the best
%    solution should be done
%    Full=1 is that a full Bayesian approach should be executed
%    Comment to full...another version of this program exists which simply
%    takes
%
% NoGui is an optional input parameter 1 if the routine should run without GUI's,
% with standard parameters, 0 if it should run with GUI features.
%
% biftfig is an optional input parameter that puts the GUI information into
% the figure with number biftfig
%
% BiftParams is an optional input parameter (a structure) containing
% the non-default values required to run the Bift routine. Presently there
% is not GUI input for these parameters...but there will be

%
% This program calculates the IFT of SAS data
% using a Bayesian a approach developed
% by Dr. S. Hansen, of the Danish Agricultural University.
%
% It has the advantage that the ususally required Dmax and Relaxation
% parameters are found using statistics of the data
%
% Details in Hansen, S. (2000) J.Appl.Chryst 33, 1415-1421
%
% It is a matlab revision of a fortran program provided kindly by
% Dr. S. Hansen.
%
% C***********************************************************************
if nargin<4 
    biftfig=0;
end
if nargin == 5 && ~isempty(BiftParams),
    %if parameters are passed to the file then we assume that also a value for
    %Full, NoGui is passed along
    
    % For explanations see the "else"
    nmax=BiftParams.nmax;
    ndist=BiftParams.ndist;
    ncalcmax=BiftParams.ncalcmax;
    alphastart=BiftParams.alphastart
    alphaband=BiftParams.alphaband;
    alpharange=BiftParams.alpharange;
    dmaxrange=BiftParams.dmaxrange;
    ndiam=BiftParams.ndiam;
    nalpha=BiftParams.nalpha;
    
    
    searchparam.maxit=BiftParams.searchparam.maxit;
    searchparam.minit=BiftParams.searchparam.minit;
    searchparam.omegaorig=BiftParams.searchparam.omegaorig;
    searchparam.omegareduction=BiftParams.searchparam.omegareduction;
    searchparam.omegamin=BiftParams.searchparam.omegamin;
    searchparam.dotsptol=BiftParams.searchparam.dotsptol;
    
    searchparam.errorcorr=1;
    searchparam.errorcorrthresh=3;
    searchparam.targeterr=0.9;
    
    verboseparam.on=BiftParams.verboseparam.on;
    verboseparam.viewstep=BiftParams.verboseparam.viewstep;
else
    default_bift_params
end

if nargin==1
    Full=1; % will do full analysis
    NoGui=0; % will have gui interfaces
end
if nargin==2
    NoGui=0; % will have GUI interfaces
end

if nargin>0
    if isstr(input1)
        %if it's a string we assume it is a filename
        aname=input1;
        fid_aname=fopen(aname,'r');
        %determine whether it is space or comma separated
        firstline=fscanf(fid_aname,'%s',120);
        fclose(fid_aname);
        fid_aname=fopen(aname,'r');
        if  isempty(findstr(',',firstline)) % not commaseparated
            IDATA=fscanf(fid_aname,'%g %g %g ',[3 inf]); % a more general input routine should be established in the future
        else
            IDATA=fscanf(fid_aname,'%g, %g, %g ',[3 inf]);
        end
        
        xtemp=IDATA(1,:);
        ytemp=IDATA(2,:);
        sdtemp=IDATA(3,:);
        %should get an errormessage in here if the file is not found
    elseif size(input1,2)==3
        %if it's not a string then it should by x,y,DeltaY data in row format
        
        aname='saxsgui';
        xtemp=input1(:,1)';
        ytemp=input1(:,2)';
        sdtemp=input1(:,3)';
    elseif size(input1,1)==3 %maybe it has been transposed
        xtemp=input1(1,:);
        ytemp=input1(2,:);
        sdtemp=input1(3,:);
        aname='saxsgui';
    else
        %errormessage should be put in
        disp('We need to have either a filename(fullpath) or a data-set with errorbars (X,Y,DeltaY)')
    end
end

if (nargin==0)
    Full=1; % will do full analysis
    NoGui=0; % will have gui interfaces
    % if no input params are given than the routine requires entering a file with data
    [aname,adir]=uigetfile('*.*','Name of inputfile');
    fid_aname=fopen([adir,aname],'r');
    IDATA=fscanf(fid_aname,'%g %g %g ',[3 inf]); % this should be changed to accept either comma or space separated variables
    xtemp=IDATA(1,:);
    ytemp=IDATA(2,:);
    sdtemp=IDATA(3,:);
end
fclose('all');

ytemp=max(ytemp,0); % negative Intensity values are turned into zeroes
sdtemp=sdtemp.*(sdtemp>0)+max(sdtemp)*(1-(sdtemp>0)); % negative or zero standard deviations are given maximum errorbars

if nargin<4 && NoGui==0
    biftfig=figure;
    plot(xtemp,ytemp)
    xlabel('q (1/A)')
    ylabel('intensity (a.u.)')
    set(gcf,'Name','Data for Bayesian IFT')
end

%now enter init parameters
if ~NoGui
    h=figstatus('First choose q-range and p(r) parameters',biftfig);
end
[qmax,qmin,ntot0]= get_init_bift_params(xtemp,NoGui);


%taking out data that's below qmin and above qmax
gooddata=(xtemp>=qmin & xtemp<=qmax);
x=xtemp(gooddata);
y=ytemp(gooddata);
sd=sdtemp(gooddata);
sigy=sdtemp(gooddata); %since sd changes (is squared) later on we use sigy to plot the errorbars on data later on
mtot=max(size(x));
sd=sd.^2;

%End of Data input

% This next section tries to come up with a good center to start the
% alpha, dmax variation from. First we estimate Rg from a guinieir
% expression, then we use this to set a fixed Dmax and let do a large alpha
% variation. The best result is then used as a starting point for full
% simplex optimization.


%first guess followed by fitting the first 3rd of the data to a Gunier
%function
p=[50,0.3];
[bx,by]=guiniercalc(p,x,y,sigy);
p(2)=p(2)*sum(y)/sum(by); %normalizing to get a better multiplication factor (P(2))
halflength=floor(length(x)/3);
p=fminsearch(@guinierchi2,p,optimset('MaxIter',10000,'MaxFunEvals',10000),x(1:halflength),y(1:halflength),sd(1:halflength));
disp ('1st guinier guess')
rg=p(1); % this is the first estimate of Rg

%trying out a very broad range of alphas to get an approximate starting
%value for alpha
if ~NoGui
    delete(h)
    h=figstatus(['Calculating rough estimate of alpha-parameter'],biftfig);
end
[alphacen,dmax]= guess_rough_bift_params(rg,alphaband,ntot0,searchparam,x,y,sd);


%now doing a simplex to obtain a better maximum
%finding optimum convergence criteria has proven difficult so therefore we
%only use a small number of iterations and leave it at that.
if ~NoGui
    delete(h)
    h=figstatus(['Rough Estimate: alpha=',num2str(alphacen),' Now calculating fine estimate alpha-parameter'],biftfig);
end

[dmaxcen,alphacen]=guess_fine_bift_params(dmax,alphacen,ntot0,searchparam,x,y,sd);

disp(['Found better estimate of best hyper-parameters, alpha = ',num2str(alphacen),' Dmax = ',num2str(dmaxcen)])

% Error correction: for many datasets there are artifacts in the data that
% are not adequately adressed through the provide errorbars on each point
% These could results from inaccurate flat-field correction, lines in the
% images, splicing of data etc.
% In any case the following allows for correction of this insufficient
% error and is activitated if the parameter

if searchparam.errorcorr==1
    if ~NoGui
        delete(h)
        h=figstatus(['Checking if errorbars are OK'],biftfig);
    else
        h=0;
    end
    
    firstresults=biftcalc(aname,x,y,sd,sigy,NoGui,1,h,nmax,alphacen,alphacen,1,1,dmaxcen,1,ntot0,mtot,searchparam,verboseparam);
    %now correcting sigy errors
    errorcorrfact=std((firstresults.yin-firstresults.yfit)./firstresults.sigmayin)/searchparam.targeterr;
    if errorcorrfact>searchparam.errorcorrthresh
        errorcorrfact=1/searchparam.targeterr*2;
    end
    sigy=max(sigy,sigy.*errorcorrfact);
    sd=sigy.*sigy;
    if ~NoGui
        h=figstatus(['Error checked: correctionfact=', num2str(errorcorrfact),'.  Now re-calculating fine estimate alpha-parameter'],biftfig);
    end
    [dmaxcen,alphacen]=guess_fine_bift_params(dmaxcen,alphacen,ntot0,searchparam,x,y,sd); %redo calculation
end


disp(['Found better estimate of best hyper-parameters, alpha = ',num2str(alphacen),' Dmax = ',num2str(dmaxcen)])
if ~NoGui
    delete(h)
    h=figstatus(['Found better estimate of best hyper-parameters, alpha = ',num2str(alphacen),' Dmax = ',num2str(dmaxcen)],biftfig);
end

if ~Full
    alpha=alphacen;
    dmax=dmaxcen;
    verboseparam.on=0;
    
    mtot=max(size(x));
    ntot=ntot0;
    dfx=dmax/ntot;
    
    [M,xf]=bift_prior(ntot,y,dfx);
    [A]=bift_trans(x,xf);
    
    ysum=(y./sd)*A;
    
    onetemp=1+0*(1:ntot);
    SDMAT=sd'*onetemp;
    B=A'*(A./SDMAT);
    clear SDMAT
    clear onetemp
    
    %This calculates the special determinant matrix.
    Bsize=max(size(B));
    %creating matrix with "next-to-diagonal"=1 and diagonal=0
    uoffset=circshift(eye(Bsize,Bsize),[1,0])+circshift(eye(size(B)),[0,1]);
    uoffset(1,Bsize)=0;
    uoffset(Bsize,1)=0;
    
    %this calculates the intensity close to 0 for normalization
    fm=(A*M')';
    summ=sum((y(1:3)./sd(1:3)))/sum(fm(1:3)./sd(1:3));
    
    M=summ.*M;
    fprior=1.0005*M;
    
    sigma=(1:ntot)*0+1;
    evmax=-1.e10;
    searchparam.omega=searchparam.omegaorig;
    verboseparam.on=0;
    
    [f,c,s,dotsp,fm]=bift_search(alpha,dmax,fprior,M,A,B,ysum,searchparam,verboseparam,x,y,sd,xf);
    alpha=alpha/(2.);
    evidence=bift_calc_evidence(alpha,s,c,B,uoffset,mtot);
    
    ftot(:,1)=f*dfx;
    ftot(ntot+1:nmax-4,nof)=0.;
    ftot(nmax-1,1)=dmax;
    ftot(nmax-2,1)=c;
    ftot(nmax-3,1)=alpha;
    
    dot(1)=dotsp;
    prob(1)=evidence;
    evimatrix(1,1:4)=[1,dmax,alpha,evidence];
    
    
    results.name=aname;
    results.xin=x;
    results.yin=y;
    results.sigmayin=sigy;
    results.yfit=fm;
    results.r=xf;
    results.pr=f;
    results.sigmapr=0;
    results.evix=0;
    results.eviy=p(1);
    results.eviz=p(2);
    results.ndmax=ndiam;
    results.rg=rg;
    results.sigmarg=sigmar;
    results.Io=sumtemp*const;
    results.sigmaIo=sqrt(esum);
    results.DmaxMean=diammean;
    results.SigmaDmaxMean=sigmad;
    
    
    
else  %do full calculation
    if ~NoGui
        delete(h)
        h=figstatus(['Full Bayesian Analysis'],biftfig);
    else
        h=0;
    end
    
    % First define the desired alpha and dmax range
    dmax=dmaxcen/dmaxrange;
    dmaxx=dmaxcen*dmaxrange *dmaxrange;
    dmax1=1.5*pi/x(1);
    
    if(dmax1 < dmaxx) %ensuring dmax is not larger than data warrants
        dmaxx=dmax1;
        dmax=min(dmaxx/2,dmax);
    end
    
    alphamin=alphacen/alpharange;
    alphamax=alphacen*alpharange;
    
    dstep=(dmaxx-dmax)/(ndiam-1);
    
    % in a later version we may do a small scanning of the terrain
    % surrounding the maximum, in order to determine which alpharange and
    % dmax range would be appropriate for now.
    
    ncalc=inf;
    while ncalc>ncalcmax
        % entering d_range
        if ~NoGui %GUI version will give possibility to change alpha-Dmax space
            
            prompt={'Enter minimum value for Dmax (A):','Enter maximum value for Dmax (A):',....
                'Enter number of D-values:'};
            def={num2str(dmax),num2str(dmaxx),num2str(ndiam)};
            dlgTitle='Input for D-max calculation';
            lineNo=1;
            qanswer=inputdlg(prompt,dlgTitle,lineNo,def);
            
            dmax=max(str2num(qanswer{1}),0);
            dmaxx=str2num(qanswer{2});
            ndiam=floor(max(str2num(qanswer{3}),1));
            dstep=(dmaxx-dmax)/(ndiam-1);
            
            % entering alpha_range
            prompt={'Enter minimum value for alpha (A):','Enter maximum value for Alpha (A):','Number of Alphapoints:'};
            def={num2str(alphamin),num2str(alphamax),num2str(nalpha)};
            dlgTitle='Input for Alpha calculation';
            lineNo=1;
            qanswer=inputdlg(prompt,dlgTitle,lineNo,def);
            
            alphamin=min(str2num(qanswer{1}),str2num(qanswer{2}));
            alphamax=max(str2num(qanswer{1}),str2num(qanswer{2}));
            nalpha=floor(max(str2num(qanswer{3}),1));
            
        end
        
        dmax0=dmax;
        
        %C     check number of estimates
        
        ncalc=ndiam*nalpha;
        disp(['Total number of calculations = ',num2str(ncalc)])
        
        if ~NoGui
            if(ncalc > ncalcmax)
                disp('That is too many calculations. It will take too long')
                disp('Please change the parameters')
            end
        else
            ncalc=0; % if NoGui then we just accept the number of calculations
        end
    end % end of alpha,dmax span definition
    if NoGui
        biftfig=0;
    end
    results=biftcalc(aname,x,y,sd,sigy,NoGui,biftfig,h,nmax,alphamin,alphamax,nalpha,ndiam,dmax0,dstep,ntot0,mtot,searchparam,verboseparam);
    end %end of full bayesian analysis

