server-scripts/qhtcp-workflow/apps/matlab/easy/NCfitImCFparforFailGbl2.m

%% CALLED BY par4GblFnc8c.m %%
function [par4scanselIntensStd,par4scanselTimesStd,par4scanTimesELr,par4scanIntensELr,par4scanCFparameters,par4scanCFdate,outC,outCstd]= ...
  NCfitImCFparforFailGbl2(parMat,times, values, timeOffsets, fileSuffix, AUCfinalTime, ~, spotAreas, printResultsDir, ~,~, sols, ~)  %,scan)

  % Preallocation for parfor loop
  st(1,1:size(times,2))=1111;
  resMat(1,1:27)=0; 
  resMatStd=resMat;
  outC=zeros(384,27);  
  outCstd=zeros(384,27);  
  for m=1:384
    pa{m}=st; 
    par4scanCFparameters{m}=parMat; 
    par4scanCFdate{m}=datestr((now),31); 
  end
  par4scanselTimesStd=pa;   
  par4scanselIntensStd=pa;
  par4scanTimesELr=pa;
  par4scanIntensELr=pa;
  par4resMat=zeros(384,27);
  par4resMatStd=zeros(384,27);

  % Spot (cultures) loop
  for ii=1:384   % startSpot:numCultures
    ii; % db print out the culture number 
    timepts=[];
    currValues=[];
    currSpotAreas=[];
    currSpotArea=[];
    dataMatrix=[];
    selTimesStd=[]; % 191024 parfor
    selIntensStd=[]; % 191024 parfor
    FiltTimesELr=[]; % 191024 parfor
    NormIntensELr=[]; % 191024 parfor
    
    % add offset...1 offset PER PLATE
    timepts=times + timeOffsets;  % (floor((ii-1)/arrayFormat) + 1);    
    currValues=values(ii,:); % change values(spotNum,:);
    % get spot areas for this culture
    currSpotArea=spotAreas(:,ii);
      
    % just use the area at the last time point
    % currSpotArea=currSpotAreas(1);
  
    % Preallocate to accomodate parfor loop 
    resMatStd=zeros(1,27);
    resMat=zeros(1,27);
    currNormIntens=currValues/currSpotArea;
    tmpx=find(currNormIntens>5); %15jh % 2.3);
    validSpot=true;
    if(isempty(tmpx) || length(tmpx)<3)
        validSpot=false;
        normIntens=currNormIntens;
        filterTimes=timepts; %filterTimes; %currTimes;
        selTimesStd=timepts;
        selIntensStd=currNormIntens;
        FiltTimesELr=timepts;
        NormIntensELr=currNormIntens;
    else
      % NCfilImCF.m
      % Preallocate incase something bails in NCscurImCFparfor
      resMatStd=zeros(1,27);
      resMat=zeros(1,27);
      hold off;
      dataMatrix=[];
      K=0;r=0;l=0;Klow=0;Kup=0;rlow=0;rup=0;llow=0;lup=0;AUC=0;MSR=0;rsquare=0; 
      bl=0;
      Tpt1=0;numFitTpts=0;thresGT2=0;minTime=0;fitbl=0; % diagnostic outputs only 
      
      timepts=timepts; % timepts=currTimes; parfor
      normIntens=currNormIntens;

      dataMatrix=[]; % arfor move clear from NCfitImCF...m
      loIntensThres=parMat(4);
      stdLoIntLim=parMat(5);
          
      % Basic filtering
      % [loInten Thres,stdbased Trim before Scurve start, and dropout detection]    
      if(max(normIntens) > 2.29)
        threshold=loIntensThres; %1.9; %Increase this value to reduce low data point (flag=2)
      else
        threshold=0;
      end

      dropThreshold=-0.0001*max(normIntens);

      % Initialize dataMatrix
      dataMatrix(1,:)=timepts;
      dataMatrix(2,:)=normIntens;
      dataMatrix(3,:)=ones;
      dataMatrix(4,:)=normIntens;
          
      % Determine a mean Intensity Index point and assoc'd TimePt    
      a=min(normIntens(normIntens>=0));   %(find(normIntens>=0)));
      b=max(normIntens(normIntens>=0));   %(find(normIntens>=0)));
      c=0.5*(b-a);
      d=b-c;
      meanIntIndPt=find(normIntens>d,1);
      meanInt=normIntens(meanIntIndPt);
      meanTime=times(meanIntIndPt);

      % NCLoIntstdTrim
      % NCLoSstdTrim.m called by NCfilImCF and NCfil.m
      flg1=0;
      loScurvLim=stdLoIntLim;
      loStimeN=1;
      last2n=1;
      stdDev=[];
      nrmIntens0=normIntens;
      for n=1:meanIntIndPt
        if nrmIntens0(n)<=0
          nrmIntens0(n)=0;
        end    
        if(nrmIntens0(n)<threshold)
          if (loStimeN-2)>0
            dataMatrix(3,1:(n-2))=2; % add to lowIntensity cull flags, the pre S cull data
          else 
            dataMatrix(3,1:n)=2;
          end
          dataMatrix(3,1:(n-2))=2;
          last2n=n; 
        end
        if n<(length(nrmIntens0)-3)
          x=nrmIntens0(n:(n+3));
          stdDev(n)=std(x);
          if (stdDev(n)<loScurvLim && flg1~=1)
            loStime=timepts(n); 
            loStimeN=n;
          end
          if stdDev(n)>6
            flg1=1;
          end
        end
      end

      % TODO repetitive code
      if (loStimeN-2)>0
        dataMatrix(3,1:(loStimeN-2))=2; % add to lowIntensity cull flags, the pre S cull data
      else 
        dataMatrix(3,1:(loStimeN-2))=2
      end

      qcutoff=2;
      qind=find(normIntens>2);  %,:,'first');
      if ~isempty(qind(3))
        qcutoff=qind(3);
      end
      [minInt,I]=min(normIntens(2:qcutoff));
      bl=minInt;
      minTime=dataMatrix(1,I);

      if (length(qind)>5)&&I>1
        dataMatrix(3,1:(I-1))=5;
      end

      tGT2Flg=0;
      for n=1:length(normIntens)
        dataMatrix(4,n)=normIntens(n)-bl; 
        if n>I && dataMatrix(4,n)>=2 && tGT2Flg==0
          thresGT2=n;thresGT2Tm=dataMatrix(1,n);tGT2Flg=1;
        end   
      end

      resMat(18)=bl;
      resMatStd(18)=bl;   
      resMatStd(20)=minTime;
      resMat(20)=minTime;

      % DropOut cull section (single drop points)
      DropOutStartPt=length(normIntens);   
      for n=1:length(normIntens) 
        if(n>1)
          if(((normIntens(n)- normIntens(n-1))< dropThreshold)) && ...
            (n > max(meanIntIndPt,thresGT2) )
            dataMatrix(3,n)=6;
          end
        end
      end

      % TODO should/could this be removed as recreated in%NCscurImCF_3parfor.m
      % Post Stdtest cull for low intensities inclusion of additional low value points 
      % selTimes=[--,--] %don't know size before as it is a filtered output
      tmpIndx=0;
      for n=1:length(normIntens) 
        if (dataMatrix(3,n)==1)
          tmpIndx=tmpIndx+1;
          selTimes(tmpIndx)=dataMatrix(1,n); % selTimes(nn)=dataMatrix(1,n);
          selIntens(tmpIndx)=dataMatrix(4,n); %s elIntens(nn)=dataMatrix(4,n);
        
        end
      end

      selTimes=selTimes';
      selIntens=selIntens';
      filtNormIntens=normIntens;
      dataMatrix0=dataMatrix;
      filterTimes=timepts; % parfor
      
      %NCscurImCF
      %NCscurImCF_1
      %NCscurImCF_2
      %NCscurImCF_3
      %NCscurImCF_3parfor
      %NCscurImCF_3parfor(dataMatrix, AUCfinalTime)
      %dataMatrix  %debugging parfor gbl  ok 85.7145; 126.4579,6, 124.5264  37tPt
      %adsfj  %debugging parfor gbl

      [resMatStd, resMat, selTimesStd, selIntensStd, FiltTimesELr, NormIntensELr] =...
        NCscurImCF_3parfor(dataMatrix0, AUCfinalTime, currSpotArea, sols, bl, minTime); 

    end

    par4scanselTimesStd{ii}=selTimesStd  %timepts'; %filterTimes';
    par4scanselIntensStd{ii}=selIntensStd;  %normIntens';
    par4scanTimesELr{ii}=FiltTimesELr; % 19_1021 preserve for CurveDisplay and EZview
    par4scanIntensELr{ii}=NormIntensELr; % 19_1021 preserve for CurveDisplay and EZview

    outC(ii,:)=resMat; %{ii, par4resMat};
    outCstd(ii,:)=resMatStd; %{ii, par4resMatStd};
  end  

  % To accomodate parfor 
  % Copy par4scan thru global p4 functions inside of parfor loop --then outside to par4Gbl_Main8b.m
  fileExt='.txt';
  filePrefix='FitResultsComplete_';
  fileNamePlate=[filePrefix fileSuffix fileExt];
  fileName=fullfile(printResultsDir, fileNamePlate);  %[printResultsDir fileNamePlate];
  fid=fopen(fileName,'w');
  fprintf(fid, 'Fit Results Complete\n');
  %fprintf(fid, 'Num.\tAUC\tMSR\tK\tr\tl\tR-squared\tK-lower\tK-upper\tr-lower\tr-upper\tl-upper\tl-lower\tArea\tLastInten\tSpineMaxRateTimePt\tLastFitTimePt\n');
  fclose(fid);
end