dynamicSPECT/matlab/fitFromClusters.m

function [globalPar, cPars, xa, U] = ...
    fitFromClusters(path, cax, cm, data, nclass, realizationId)
    %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=3;
    nglob=3;
    
    %correct plotFunctions and evaluateCenterDiff
    
    c1=size(cax,1);
    A=reshape(data,[],c1);

    %find most common responses
    %nclass=10;
    [centers,U]=fcm(A,nclass);

    %data(ix,iy,iz,:)=A((iz-1)*i1*i2+(iy-1)*i1+ix,:)

    cm2=ones(size(centers,1),1)*cm';
    
    %w=sqrt(cm2.*centers)./sum(cm2.*centers,2);
    w=ones(size(centers,1),1)*sqrt(cm)';
    
    pars=cax(1:2)';
    
    fun=@(x)evaluateCenterDiff(cax,centers,w,pars,x(1:nglob),...
        reshape(x(nglob+1:nglob+npar*nclass),[],nclass)');
    %cax are the delays
    %c
    
    valueGeneratorSettings.type='poisson';
    valueGeneratorSettings.parameters=[25];%mu
    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;
    
    
    for j=1:nclass
        
        valueGeneratorSettings(ivar).type='loggaus';
        valueGeneratorSettings(ivar).parameters=[-3,1];
        valueGeneratorSettings(ivar).u=[0,Inf];
        ivar=ivar+1;
        
        valueGeneratorSettings(ivar).type='flat';
        valueGeneratorSettings(ivar).parameters=[0,1];
        valueGeneratorSettings(ivar).u=[0,1];
        ivar=ivar+1;
    
        valueGeneratorSettings(ivar).type='loggaus';
        valueGeneratorSettings(ivar).parameters=[-7,2];
        valueGeneratorSettings(ivar).u=[0,Inf];
        ivar=ivar+1;
    end
    
    %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];
        
        
    [x0,ul,ub]=generateInitialValues(valueGeneratorSettings);
    options = optimoptions(@lsqnonlin);
    options.Display='off'; 
    %[options.OptimalityTolerance,options.FunctionTolerance,options.StepTolerance];
    %holds best estimate
    x=x0;
    fmin=1e30;
    n=20;
    n1=20;
    imin=-1;
    for i=1:n
    
        %find best fit for nclass responses. 
        %total number of points that are fitted is 20*nclass, total number of
        %parameters estimated is 4+3*nclass
        [x1,resnorm,residual,exitflag,output,lambda,jacobian]  = lsqnonlin(fun,x0,ul,ub,options);
        
        %this gets tau and A
        
        %now fit every curve individually with the known input function
        for j=1:nclass
            cay=centers(j,:)';
            w=ones(size(cay));
            
            fun1=@(par)evaluateSingleDiff(cax,cay,w,pars,x1,par);
            
            qpar=zeros(1,3);
            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
            x1(nglob+(j-1)*3+1:nglob+j*3)=qpar;
            fprintf('Center: %d (%f %f %f)/(%f %f %f)\n',j,...
                par0(1),par0(2),par0(3),...
                par1(1),par1(2),par1(3));
            %residual
        
        end
        %writematrix(ci,fullfile(path,sprintf('CI%d.txt',i)),'Delimiter','tab');
        %output
        chi2=fun(x1);
        chi2=chi2'*chi2;
        
        if chi2<fmin
            fmin=chi2;
            x=x1;
            imin=i;
            %fci=ci;
        end
            
        %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
    ofile=fullfile(path,sprintf('fitParFinal_%d_%d.txt',nclass,realizationId));
    try 
      writematrix([fmin x],ofile,'Delimiter','tab');
    catch ME
      dlmwrite(ofile,[fmin x],'\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 cax(1:2)'];
    
    cPars=reshape(x(nglob+1:nglob+npar*nclass),[],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;
    
    plotFunctions(path,globalPar,cPars,cax,centers);
    
    
    dt=size(data);
    
    k1v=U'*cPars(:,1);
    xa(:,:,:,1)=reshape(k1v,dt(1:3));

    vbfv=U'*cPars(:,2);
    xa(:,:,:,2)=reshape(vbfv,dt(1:3));
    
    k2v=U'*cPars(:,3);
    xa(:,:,:,3)=reshape(k2v,dt(1:3));
end