dynamicSPECT/matlab/fitCenters.m

function [globalPar, cPars] = ...
    fitCenters(path, patientID, realizationId, cax, cm, centers)
    %path is the core directory
    %cax are the mid intervals in images
    %cm are the interval durations (1,n)
    %data are the images
    %nclass is number of typical responses sought
    npar=1;
    nglob=2;
    nclass=size(centers,1);
    
    cm2=ones(nclass,1)*cm';
    
    %w=sqrt(cm2.*centers)./sum(cm2.*centers,2);
    w=ones(nclass,1)*sqrt(cm)';
    %make sure that the function that is most similar to input function has
    %a high weight on maximum value, which comes as the second reading
    %abandonded, because decay time gets shortened and fitting of the slow
    %rise is given up in favor of meeting the high point through k1 value
    %[~,im]=max(centers(:,2));
    %w(im,2)=w(im,2)+0.1*(sqrt(max(cm))-w(im,2));
    
    
    %cax are the delays
    %c
    valueGeneratorSettings.type='poisson';
    valueGeneratorSettings.parameters=[25];
    valueGeneratorSettings.u=[0,Inf];
        
    ivar=2;
    valueGeneratorSettings(ivar).type='loggaus';
    valueGeneratorSettings(ivar).parameters=[1.5,1];%log10(mean),log10(std)
    valueGeneratorSettings(ivar).u=[0,Inf];
    ivar=ivar+1;
    
    valueGeneratorSettings(ivar).type='loggaus';
    valueGeneratorSettings(ivar).parameters=[-1,1];%log10(mean),log10(std)
    valueGeneratorSettings(ivar).u=[0,Inf];
    ivar=ivar+1;

    valueGeneratorSettings(ivar).type='gaus';
    valueGeneratorSettings(ivar).parameters=[10,2];%log10(mean),log10(std)
    valueGeneratorSettings(ivar).u=[-Inf,Inf];
    ivar=ivar+1;

    %for individual fits
    valueGeneratorSettings1.type='loggaus';
    valueGeneratorSettings1.parameters=[-3,1];
    valueGeneratorSettings1.u=[0,Inf];
    ivar=2;
        
    valueGeneratorSettings1(ivar).type='flat';
    valueGeneratorSettings1(ivar).parameters=[0,1];
    valueGeneratorSettings1(ivar).u=[0,1];
    ivar=ivar+1;
    
    valueGeneratorSettings1(ivar).type='loggaus';
    valueGeneratorSettings1(ivar).parameters=[-7,2];
    valueGeneratorSettings1(ivar).u=[0,Inf];
    ivar=ivar+1;
    
    valueGeneratorSettings1(ivar).type='gaus';
    valueGeneratorSettings1(ivar).parameters=[10,5];
    valueGeneratorSettings1(ivar).u=[0,Inf];
        
        
    options = optimoptions(@lsqnonlin);
    options.Display='off'; 
    %[options.OptimalityTolerance,options.FunctionTolerance,options.StepTolerance];
    %holds best estimate
    fmin=1e30;
    n=20;
    n1=20;
    imin=-1;
    nFitIndividual=4;
    qfit=zeros(1,nFitIndividual*nclass);
    qfitmin=qfit;
    [x,ul,ub]=generateInitialValues(valueGeneratorSettings);
    for i=1:n
    
        %first estimate input functions
        %take function with maximum response as IVF surrogate
        [m, mi]=max(centers);
        [~, mj]=max(m);
        cay=centers(mi(mj),:)';
        w=ones(size(cay));
        
        fun=@(x)evaluateSingleDiffZero(cax,cay,w,x);
        [x0,ul,ub]=generateInitialValues(valueGeneratorSettings);
        
        [x1,resnorm,residual,exitflag,output,lambda,jacobian]  = lsqnonlin(fun,x0,ul,ub,options);
        fprintf('A %f\n',x1(1));
        
        while x1(1)<1
            [x0,ul,ub]=generateInitialValues(valueGeneratorSettings);
            [x1,resnorm,residual,exitflag,output,lambda,jacobian]  = lsqnonlin(fun,x0,ul,ub,options);
            fprintf('A %f\n',x1(1));
        end
        %this gets A, tau and alpha
        %A shouldn't be trivial
        
        fprintf('[%d] Global fit: (%f,%f,%f %f) (%f,%f,%f %f)\n',i,x0(1),x0(2),x0(3),x0(4),...
            x1(1),x1(2),x1(3),x1(4));
        
        chiSum=0;
        %now fit every curve individually with the known input function
        for j=1:nclass
            cay=centers(j,:)';
            w=ones(size(cay));
            w=sqrt(cm);
            
            fun1=@(par)evaluateSingleDiff(cax,cay,w,x1,par);
            
            qpar=zeros(1,4);
            qmin=1e30;
            for k=1:n1
                [par0,ul1,ub1]=generateInitialValues(valueGeneratorSettings1);
                [par1,resnorm,residual,exitflag,output,lambda,jacobian]  = ...
                    lsqnonlin(fun1,par0,ul1,ub1,options);
                fchi2=sqrt(residual'*residual);
                if fchi2<qmin
                    qmin=fchi2;
                    qpar=par1;
                end
            end
            qfit((j-1)*nFitIndividual+1:j*nFitIndividual)=qpar;
            fprintf('Center: %d (%f %f %f %f)/(%f %f %f %f)\n',j,...
                par0(1),par0(2),par0(3),par0(4),...
                par1(1),par1(2),par1(3),par1(4));
            %residual
            chiSum=chiSum+qmin;
        end
        %writematrix(ci,fullfile(path,sprintf('CI%d.txt',i)),'Delimiter','tab');
        %output
        chi2=chiSum;
        
        if chi2<fmin
            fmin=chi2;
            x=x1;
            imin=i;
            qfitmin=qfit;
            %fci=ci;
        end
        fprintf('chi2: %d,%d\n',size(chi2,1),size(chi2,2));    
        %writematrix([chi2 x1],fullfile(path,sprintf('fitPar%d.txt',i)),'Delimiter','tab');
        fprintf('Done %d %f (%.2f,%.2f,%.2f)/(%.2f,%.2f,%.2f)\n',...
            i,chi2,x0(1),x0(2),x0(3),x1(1),x1(2),x1(3)); 
        [x0,ul,ub]=generateInitialValues(valueGeneratorSettings);
        
        
    end
    
    %this is the result that will be used by k1 calculator
    code=sprintf('%s_%d_%d',patientID,nclass,realizationId);
    ofile=fullfile(path,patientID,sprintf('%s_fitParFinal.txt',code));
    try
      writematrix([fmin x qfitmin],ofile,'Delimiter','tab');
    catch ME
      dlmwrite(ofile,[fmin x qfitmin],'\t')
    end
    %ofile=fullfile(path,sprintf('CIFinal_%d_%d.txt',nclass,realizationId));
    %writematrix(fci,ofile,'Delimiter','tab');
    disp(sprintf('Found min at %d',imin)); 
    
    globalPar=[x(1:3) 0 0 0];
    
    cPars=reshape(qfitmin,[],nclass)';
    %reset k1 for those where BVF is nearly one
    %cPars(cPars(:,2)>0.4,1)=0;
    %reset k1 for those where k2 is abnormally large
    %cPars(cPars(:,3)>1,1)=0;
    %opath=fullfile(path,patientID);
    %plotFunctions(opath,code,globalPar,cPars,cax,centers);
    
    
end