Format NCscurImCF_3parfor.m

workflow/templates/easy/NCfitImCFparforFailGbl2.m

@@ -218,11 +218,8 @@ function [par4scanselIntensStd,par4scanselTimesStd,par4scanTimesELr,par4scanInte
     outCstd(ii,:)=resMatStd; %{ii, par4resMatStd};
   end  
 
-  %*********************************19_1001***********************************
+  % To accomodate parfor 
-  %To accomodate parfor copy par4scan thru global p4 functions inside of
+  % Copy par4scan thru global p4 functions inside of parfor loop --then outside to par4Gbl_Main8b.m
-  %parfor loop --then outside to par4Gbl_Main8b.m
-
-  %**************************************************************************
   fileExt='.txt';
   filePrefix='FitResultsComplete_';
   fileNamePlate=[filePrefix fileSuffix fileExt];

workflow/templates/easy/NCscurImCF_3parfor.m

@@ -1,25 +1,18 @@
 %% CALLED BY NCfitImCFparforFailGbl2.m %%
 function [resMatStd, resMat, selTimesStd, selIntensStd, FiltTimesELr, NormIntensELr] =...
   NCscurImCF_3parfor(dataMatrix, AUCfinalTime, currSpotArea, sols, bl, minTime)
-%Major revision for Early-Late data cuts to improve accuracof 'r'. Removed legacy iterative method.
-%Significant Modification for Parfor  
-%***************************************************************
 
-%##########################################################################
- %*******************************************  New Stage 1***************************************************************
   % Preallocate
   resMatStd=zeros(1,27);
   resMat=zeros(1,27);
   % Set internal variables sent to matlab fit function 
   me=200;
   meL=750;
-mi=25;  %50
+  mi=25;
   miL=250;
-%***********************************
   rmsStg1=0;
   rmsStg1I(1)=0;
   slps=1;
-
   filterTimes=[];
   normIntens=[];
   nn=1;
@@ -33,28 +26,17 @@ for n=1:size(dataMatrix,2)
       nn=nn+1;
     end
   end
-
-%------------------------------------------------------------------
-
-%++++++++++++++++++++++++++++++++++++++++
-
   filterTimes=filterTimes';  
   selTimesStd=filterTimes; 
   normIntens=normIntens';  
   selIntensStd=normIntens;
-
-%normIntens  %debugging parfor gbl  200330 good values
-%afgj
-
   lastTptUsed=1;
   lastIntensUsed=1;
   thresGT2TmStd=0;  
   try
-    
     lastTptUsed=max(filterTimes);
     lastIntensUsed=normIntens(length(normIntens));
     lastIntensUsedStd=lastIntensUsed;
-    
     lastTptUsedStd=lastTptUsed;
     Tpt1Std=filterTimes(1);    
     numFitTptsStd=nnz((normIntens(:)>=0)==1);
@@ -67,12 +49,12 @@ selIntensStd=normIntens;
     numTptsGT2Std=nnz((normIntens(:)>=2)==1); % nnz(filterTimes(find(filterTimes>=thresGT2Tm)));
     K_Guess=max(normIntens);
     numTimePts=length(filterTimes);
-    opts =  fitoptions('Method','Nonlinear','Robust','On',...
+    opts=fitoptions('Method','Nonlinear','Robust','On','DiffMinChange',1.0E-11,'DiffMaxChange',0.001,...
-        'DiffMinChange',1.0E-11,'DiffMaxChange',0.001,...
+      'MaxFunEvals',me, 'MaxIter', mi, 'TolFun', 1.0E-12, 'TolX', 1.0E-10, 'Lower', [K_Guess*0.5,0,0],...
-        'MaxFunEvals',me, 'MaxIter', mi, 'TolFun', 1.0E-12, 'TolX', 1.0E-10, 'Lower', [K_Guess*0.5,0,0], 'StartPoint', [K_Guess,filterTimes(floor(numTimePts/2)),0.30], 'Upper', [K_Guess*2.0,max(filterTimes),1.0],'Display','off');
+      'StartPoint', [K_Guess,filterTimes(floor(numTimePts/2)),0.30], 'Upper', [K_Guess*2.0,max(filterTimes),1.0],'Display','off');
     ftype=fittype('K / (1 +  exp(-r* (t - l )))','independent','t','dependent',['K','r','l'],'options',opts);
     
-    % carry out the curve fitting process
+    % Carry out the curve fitting process
     [fitObject, errObj]=fit(filterTimes,normIntens,ftype);
     coeffsArray=coeffvalues(fitObject);
     rmsStg1=errObj.rsquare;
@@ -82,17 +64,15 @@ selIntensStd=normIntens;
     l=coeffsArray(2); sDat(slps,3)=coeffsArray(2); % lag time
     r=coeffsArray(3); sDat(slps,4)=coeffsArray(3); % rateS
 
-    % integrate (from first to last time point)
+    % Integrate (from first to last time point)
     numVals=size(filterTimes);
     numVals=numVals(1);
     t_begin=0;
     t_end=AUCfinalTime;
     AUC=(K/r*log(1+exp(-r*(t_end-l)))-K/r*log(exp(-r*(t_end-l)))) - (K/r*log(1+exp(-r*(t_begin-l)))-K/r*log(exp(-r*(t_begin-l))));
     MSR=r;
-
     rsquare=errObj.rsquare;  
     confObj=confint(fitObject,0.9); % get the 90% confidence 
-    
     NANcond=0; stdNANcond=0;  % stdNANcond added to relay not to attempt ELr as there is no curve to find critical point
     confObj_filtered=confObj; 
     Klow=confObj(1,1); sDat(slps,5)=confObj(1,1);
@@ -104,34 +84,22 @@ selIntensStd=normIntens;
     if(isnan(Klow)||isnan(Kup)||isnan(llow)||isnan(lup)||isnan(rlow)||isnan(rup))
       NANcond=1; stdNANcond=1; % stdNANcond added to relay not to attempt ELr as there is no curve to find critical point
     end  
-    
+  catch ME
-    %rup  %debugging parfor gbl 200330
-    %Klow
-    %asdfjj114
-% {   
-catch %ME
     % if no data is given, return zeros
     AUC=0;MSR=0;K=0;r=0;l=0;rsquare=0;Klow=0;Kup=0;
     rlow=0;rup=0;lup=0;llow=0; 
     NANcond=1; stdNANcond=1;  %stdNANcond added to relay not to attempt ELr as there is no curve to find critical point
-end   %end Try
+  end
-%}
+
   if (exist('K','var')&& exist('r','var') && exist('l','var'))
     t=(0:1:200);   
     Growth=K ./ (1 +  exp(-r.* (t - l )));
-    fitblStd= min(Growth);   %jh diag
+    fitblStd=min(Growth);
   end
    
   cutTm(1:4)=1000;  %-1 means cuts not used or NA
-%{    
-l  %debugging parfor
-K
-r
-Klow
-k
-%}    
-%***Preserve for ResultsStd+++++++++++++++++++++++++++++
   
+  % Preserve for ResultsStd
   resMatStd(1)=AUC;
   resMatStd(2)=MSR;
   resMatStd(3)=K;
@@ -147,18 +115,12 @@ resMatStd(12)= lup;
   resMatStd(13)=currSpotArea;
   resMatStd(14)=lastIntensUsedStd; % filtNormIntens(length(filtNormIntens));
 
-%spline fit unneccessary and removed;therefore No maxRateTime assoc'd w/spline fit
-%try
-%resMatStd(15)= maxRateTime;
-%catch
   maxRateTime=0; %[];  %Std shows []; ELr shows 0;  %parfor
   resMatStd(15)=0; %maxRateTimestdMeth;
-%end
 
-
+  resMatStd(16)=lastTptUsedStd;
-resMatStd(16)= lastTptUsedStd;   %filterTimes(length(filterTimes));
   if isempty(Tpt1Std)
-        Tpt1Std= 777;  %0.000002;  %
+    Tpt1Std=777;
   end
   resMatStd(17)=Tpt1Std;
   resMatStd(18)=bl; %perform in the filter section of NCfitImCFparfor
@@ -167,13 +129,11 @@ resMatStd(20)= minTime; %Not affected by changes made in NCscur...for refined 'r
   resMatStd(21)=thresGT2TmStd;
   resMatStd(22)=numFitTptsStd;
   resMatStd(23)=numTptsGT2Std;
-
   resMatStd(24)=999; %The Standard method has no cuts .:.no cutTm  
   resMatStd(25)=999;
   resMatStd(26)=999;
   resMatStd(27)=999;
 
-
   %if SwitchEvsEL==3  %Remove 'SwitchEvsEL==...' temporary SWITCH when Hartman decides what he wants
   %*********************************************************************************
   %ELr         New Experimental data through L+deltaS Logistic fit for 'Improved r'  Fitting
@@ -199,17 +159,13 @@ rsTmStd= LL-tc1;   %%riseTime (first critical point to L)
     deltS=rsTmStd/2;
     tc1Early=tc1-deltS;    %AKA- tc1AdjTm %2*tc1 -LL
     L_Late=LL+deltS;
-    
-           
     tc1EdatPt=find(filterTimes>(tc1Early),1,'first');
     cutTm(1)=filterTimes(2);
     cutDatNum(1)=2;
     cutTm(2)=tc1Early;
     cutDatNum(2)=tc1EdatPt-1;
-    
     L_LDatPt=find(filterTimes< L_Late,1,'last');
     tc2LdatPt=find(filterTimes< tc2+rsTmStd,1,'last');
-    
     cutTm(3)=L_Late;
     cutDatNum(3)=L_LDatPt;
     
@@ -222,19 +178,16 @@ tms=[];
     tms(1:L_LDatPt-tc1EdatPt+2)=(stdTimes(L_LDatPt));
     tms(2:end)=stdTimes(tc1EdatPt:L_LDatPt);
     tms(1)=stdTimes(1);
-%-----------------------------------------------
+
-%Include/Keep late data that define K *********
+    % Include/Keep late data that define K
     if length(tmpIntens(tc2LdatPt:end))> 4
       KlastInts=tmpIntens(tc2LdatPt:end);
       KlastTms=stdTimes(tc2LdatPt:end);
       lengthKlast=length(tmpIntens(tc2LdatPt:end));
       ints(end:(end+ lengthKlast-1))=KlastInts;
       tms(end:(end+ lengthKlast-1 ))=KlastTms;
-   
       cutTm(4)=tc2+rsTmStd;
       cutDatNum(4)=tc2LdatPt-1;
-   
-   
     else
       lengthKlast=length(tmpIntens(tc2LdatPt-1:end)); 
       if lengthKlast>1
@@ -243,26 +196,18 @@ else
         ints(end:(end+ lengthKlast-1 ))=KlastInts;
         tms(end:(end+ lengthKlast-1 ))=KlastTms;
       end
-
       cutTm(4)=stdTimes(tc2LdatPt-1);
       cutDatNum(4)=tc2LdatPt-2;   %length(stdTimes(end-(lengthKlast-1):end)); 
-   
     end
-    
-
-%************************************************
     Ints=[];
     Tms=[];
     Ints=ints';
     Tms=tms';
-
     try
       filterTimes=Tms; filterTimes4=Tms; 
       normIntens=Ints;   normIntens4=Ints;
 
-%----------------------------------------------------------------------------------------------
+      % Classic symmetric logistic curve fit setup restated as COMMENTS for reference convenience
-    
-   %classic symetric logistic curve fit setup restated as COMMENTS for reference convenience----------------------------
       % opts=fitoptions is the same as for Std and so is redundant     
       % opts=fitoptions('Method','Nonlinear','Robust','On',...
       % 'DiffMinChange',1.0E-11,'DiffMaxChange',0.001,...
@@ -280,29 +225,25 @@ try
         GrowthELr = K ./ (1 +  exp(-r.* (t - l )));
         fitblELr= min(GrowthELr);   %jh diag
       end
-   
     catch ME
       % if no data is given, return zeros
       AUC = 0;MSR = 0;K = 0;r = 0;l = 0;rsquare = 0;Klow = 0;Kup = 0;
       rlow = 0;rup = 0;lup = 0;llow = 0; %normIntens=[];
+    end
+  end
 
-end   %end Try
-
-end  %if stdNANcond=0   
- %++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
   % Update values if r is better(higher) with removal of early data
   try
     if r3>r && stdNANcond==0 
       r=r3; sDat(slps,4)=sDat(slps,4); % rateS 
       K=coeffsArray(1); sDat(slps,2)=coeffsArray(1);  % Carrying Capacity
       l=coeffsArray(2); sDat(slps,3)=coeffsArray(2); % lag time
-     
       coeffsArray=coeffvalues(fitObject);
       rmsStg1=errObj.rsquare;
       rmsStg1I(slps)=errObj.rsquare;
       sDat(slps,1)=errObj.rsquare; 
     
-     %jh diagnostics*************    
+      % JH diagnostics   
       numFitTpts=nnz((normIntens(:)>=0)==1);
       thresGT2=find(((normIntens(:)>2)==1), 1);
       thresGT2Tm=filterTimes(thresGT2); 
@@ -311,12 +252,9 @@ try
     
       AUC=(K/r*log(1+exp(-r*(t_end-l)))-K/r*log(exp(-r*(t_end-l)))) - (K/r*log(1+exp(-r*(t_begin-l)))-K/r*log(exp(-r*(t_begin-l))));
       MSR=r3;
-    %***************************
-    
       rsquare=errObj.rsquare;  
       confObj=confint(fitObject,0.9); % get the 90% confidence 
       NANcond=0;
-
       confObj_filtered=confObj; 
       Klow=confObj(1,1); sDat(slps,5)=confObj(1,1);
       Kup=confObj(2,1); sDat(slps,6)=confObj(2,1);
@@ -327,7 +265,6 @@ try
       if(isnan(Klow)||isnan(Kup)||isnan(llow)||isnan(lup)||isnan(rlow)||isnan(rup))
         NANcond=1;    
       end
-    
       filterTimes=Tms; 
       normIntens=Ints;   
       resMat(17)=.00002;  
@@ -337,9 +274,7 @@ try
     else % r is better than r3 so use the Std data in the ELr result sheet
       filterTimes=selTimesStd; 
       normIntens=selIntensStd;
-     
       lastTptUsed=lastTptUsedStd; % Reinstall Std values for jh diags
-     
       Tpt1=filterTimes(1);
       try
         if isempty(Tpt1)
@@ -353,24 +288,17 @@ resMat(17)= Tpt1;
       numFitTpts=numFitTptsStd;
       numTptsGT2=numTptsGT2Std;
       thresGT2Tm=thresGT2TmStd;
-     
+      cutTm(1:4)=1000;  % 1 means cuts not used or NA
-     cutTm(1:4)= 1000;  %-1 means cuts not used or NA
-     
       resMat(18)=bl;  % only applicable to Std curve Fit; ELr superceeds and makes meaningless
       resMat(19)=fitblStd; % only applicable to Std curve Fit; ELr superceeds and makes meaningless
       resMat(20)=minTime; % only applicable to Std curve Fit; ELr superceeds and makes meaningless
     end  % if r3>r1
-    
-   %rup
-   %asdf352
-% {   
   catch ME
     % if no data is given, return zeros
     AUC=0;MSR=0;K=0;r=0;l=0;rsquare=0;Klow=0;Kup=0;
     rlow=0;rup=0;lup=0;llow=0; % normIntens=[];
+  end
    
-end   %end Try
-%}    
   resMat(1)=AUC;
   resMat(2)=MSR;
   resMat(3)=K;
@@ -385,14 +313,10 @@ resMat(11)= llow;
   resMat(12)=lup;
   resMat(13)=currSpotArea;
   resMat(14)=lastIntensUsed;  %filtNormIntens(length(filtNormIntens));
-
   % spline fit unneccessary and removed therfor no max spline rate time->set 0
   maxRateTime=0; % ELr will show 0; Std will show [] 
   resMat(15)=maxRateTime;
-
   resMat(16)=lastTptUsed; % filterTimes(length(filterTimes));
-
-
   try % if Std fit used no cuts .:. no cutTm 
     resMat(24)=cutTm(1);
     resMat(25)=cutTm(2);
@@ -407,10 +331,6 @@ end
 
   FiltTimesELr=filterTimes;
   NormIntensELr=normIntens;
-
-%**********************************************************************************
-%########################################################################## 
-
   lastTptUsed=max(filterTimes);
   lastIntensUsed=normIntens(length(normIntens));  
   
@@ -419,9 +339,10 @@ NormIntensELr= normIntens;
     Growth=K ./ (1 +  exp(-r.* (t - l )));
     fitbl=min(Growth);  % jh diag
   end 
-  %} 
+  try % jh diag
-  try
+    if isempty(thresGT2Tm)
- if isempty(thresGT2Tm),thresGT2Tm=0;end  %jh diag
+      thresGT2Tm=0
+    end
   catch
     thresGT2Tm=0;
     numTptsGT2=0;
@@ -430,8 +351,6 @@ NormIntensELr= normIntens;
   resMat(21)=thresGT2Tm;
   resMat(22)=numFitTpts;
   resMat(23)=numTptsGT2;
- 
+end
- 
-end  %function end