(************************************************************************)
(*			    Theranostics_laser treatment analysis		         *)
(*	    an algorithm for quantification of the treatment effect       *)
(*																	 *)
(* 					         Rok Podlipec					       *)
(*							 january, 2021					 	  *)
(*																	 *)
(************************************************************************)


$HistoryLength=0;
(RefLink[Needs, paclet:ref/Needs])["ErrorBarPlots`"];
SetDirectory[NotebookDirectory[]];
ClearAll;

(******Tested modeled functions to best quantify the extent of laser treatment effect induced on tissue********)
(*The algorithm is useful for the treatment evaluation and laser mode adaptation during theranostics operation*)

theranosticsFunc2[descriptors_,B_,C_,D_,E_,func6Const_]:=Module[{part1,part2,part3,func},

(*descriptor1 ... FLIM parameter ... tau2/a1_init       ... column 5*)
(*descriptor2 ... FLIM parameter ... tau2/a1_afterTreat ... column 7*)
(*descriptor3 ... FLIM parameter ... I0_init            ... column 1*)
(*descriptor4 ... AF parameter   ... I_afterTreat       ... column 2*)
(*descriptor5 ... AF parameter   ... ΔI0_init           ... column 3*)
(*descriptor6 ... AF paraemter   ... ΔI_afterTreat      ... column 4*)
(*descriptor7 ... FLIM parameter ... w_init             ... column 6*)
(*descriptor8 ... FLIM parameter ... w_afterTreat       ... column 8*)

part1=Abs[(descriptors[[2]]-descriptors[[1]])]/(B*descriptors[[1]]);

(*1st option*)
(*part2=Power[Abs[(descriptors[[6]]/descriptors[[4]]-descriptors[[5]]/descriptors[[3]])],1]*C;
part3=Power[(descriptors[[8]]-descriptors[[7]])/descriptors[[7]],1]*D ;*)

(*2nd option*)
(*Which[Abs[(descriptors[[6]]/descriptors[[4]]-descriptors[[5]]/descriptors[[3]])]<0.07,
	part2=Abs[(descriptors[[6]]/descriptors[[4]]-descriptors[[5]]/descriptors[[3]])]*C,
	0.07<Abs[(descriptors[[6]]/descriptors[[4]]-descriptors[[5]]/descriptors[[3]])]<0.2,
	part2=0.7,
	0.2<Abs[(descriptors[[6]]/descriptors[[4]]-descriptors[[5]]/descriptors[[3]])],
	part2=1
	];
Which[Abs[(descriptors[[8]]-descriptors[[7]])/descriptors[[7]]]<0.25,
	part3=Abs[(descriptors[[8]]-descriptors[[7]])/descriptors[[7]]]*D,
	0.25<Abs[(descriptors[[8]]-descriptors[[7]])/descriptors[[7]]]<1,
	part3=0.5,
	1<Abs[(descriptors[[8]]-descriptors[[7]])/descriptors[[7]]],
	part3=1];
func=part1+part2+part3;*)

(*3rd option*)
(*part2=Exp[-Power[Abs[(descriptors[[6]]/descriptors[[4]]-descriptors[[5]]/descriptors[[3]])]/0.05,2]];
part3=Exp[-Power[Abs[(descriptors[[8]]-descriptors[[7]])/descriptors[[7]]/0.5],2]];
func=part1/(part2*part3);*)

(*4th option*)
(*intensity*)
(*Which[
(descriptors[[6]]/descriptors[[4]]-descriptors[[5]]/descriptors[[3]])<0,
(*part2=0,*)
part2=(descriptors[[6]]/descriptors[[4]]-descriptors[[5]]/descriptors[[3]])*E,
0<(descriptors[[6]]/descriptors[[4]]-descriptors[[5]]/descriptors[[3]])<0.06,
part2=(descriptors[[6]]/descriptors[[4]]-descriptors[[5]]/descriptors[[3]])*C,
0.06<(descriptors[[6]]/descriptors[[4]]-descriptors[[5]]/descriptors[[3]]),
part2=Power[(descriptors[[6]]/descriptors[[4]]-descriptors[[5]]/descriptors[[3]])-0.06,1/4]];

(*FLIM descriptor width*)
Which[Abs[(descriptors[[8]]-descriptors[[7]])/descriptors[[7]]]<0.2,
	part3=Abs[(descriptors[[8]]-descriptors[[7]])/descriptors[[7]]]*D,
	0.2<Abs[(descriptors[[8]]-descriptors[[7]])/descriptors[[7]]],
	part3=Power[Abs[(descriptors[[8]]-descriptors[[7]])/descriptors[[7]]]-0.2,1/4]];
func=part1+part2+part3;
If[func<0, func=0];*)

(*5th option*)
part2=(1-Exp[-(Abs[(descriptors[[6]]/descriptors[[4]]-descriptors[[5]]/descriptors[[3]])]/0.1)^3])*func6Const[[1]];
part3=( 1-Exp[-(Abs[(descriptors[[8]]-descriptors[[7]])/descriptors[[7]]]/0.4)^3])*func6Const[[2]];

func=part1+part2+part3;

If[func<0, func=0];

Return[func]

];


ClearAll;

(*Import data*)

(*theranostics parameters*)
startRow=34;
endRow=167;
(*parameters from fitted histogram distributions*)
theranosticsParam=Flatten[Import["Dataset_Theranostics laser parameters and descriptors for treatment effect quantif.xlsx",{"Data",2}][[startRow;;endRow,2;;9]],0];

(*Print[{"wtf",theranosticsParam[[2,1]]}];
emptyRow2=Select[theranosticsParam[[1;;2]],Length@#≠0&];
emptyRow=DeleteCases[theranosticsParam[[1;;2]],{}];*)

(*delete all rows where we don't have the data of parameters for therapy function*)
theranosticsParam=Replace[theranosticsParam,""->0,{2}];
emptyRow=Position[theranosticsParam,{0,0,0,0,0,0,0,0}];
Print[emptyRow];
theranosticsParam=Delete[theranosticsParam,emptyRow];

(**********************************)
(*IMPORT RAW DATA from DATASET file*)
(***********************************)

(*peak irradiance [10^10 W/cm2]*)
theranosticsParam2=Delete[Flatten[Import["Dataset_Theranostics laser parameters and descriptors for treatment effect quantif.xlsx",{"Data",1}][[startRow;;endRow,12]],0],emptyRow];
(*peak irradiance_pop [10^10 W/cm2]*)
theranosticsParam2Pop=Delete[Flatten[Import["Dataset_Theranostics laser parameters and descriptors for treatment effect quantif.xlsx",{"Data",1}][[startRow;;endRow,13]],0],emptyRow];
(*energy per pixel [uJ]*)
theranosticsParam3=Delete[Flatten[Import["Dataset_Theranostics laser parameters and descriptors for treatment effect quantif.xlsx",{"Data",1}][[startRow;;endRow,18]],0],emptyRow];
(*experiment done more times on one site and done single time: more - 3; single- 1*)
theranosticsParamExpApproach=Delete[Flatten[Import["Dataset_Theranostics laser parameters and descriptors for treatment effect quantif.xlsx",{"Data",1}][[startRow;;endRow,21]],0],emptyRow];
(*qualitative therapy effect 0-5*)
theranosticsParam4=Delete[Flatten[Import["Dataset_Theranostics laser parameters and descriptors for treatment effect quantif.xlsx",{"Data",1}][[startRow;;endRow,5]],0],emptyRow];

(*data1=Flatten[Import["Data.txt","Table"][[1;;All(*, {,}*)]],0];*)

theranosticsParamLength=Length[theranosticsParam];
Print[theranosticsParamLength];

descriptorsOrder1={5,7,1,2,6,8};
descriptorsOrder2={5,7,1,2,3,4,6,8};
theranosticsParamDescriptors1=Table[theranosticsParam[[#,i]],{i,descriptorsOrder1}]&/@Range[theranosticsParamLength];
Print[{"theranosticsParam",theranosticsParamDescriptors1[[1]]}];
theranosticsParamDescriptors2=Table[theranosticsParam[[#,i]],{i,descriptorsOrder2}]&/@Range[theranosticsParamLength];

a=0.25;
b=0.25;
c=0.1;
d=0.1;
e=0.1;
func6Const={1,1};

therapyEffect2=theranosticsFunc2[theranosticsParamDescriptors2[[#]],b,c,d,e,func6Const]&/@Range[theranosticsParamLength];
Print[{"therapyEffect2",therapyEffect2}];


(*Representation of energy per pixel..the range is approx 1 uJ - 1000 uJ per pixel*)
(*we might use logarithnic scale using Hue function*)

theranosticsParam3Log=N[Log10[theranosticsParam3[[#]]]]&/@Range[theranosticsParamLength];


(*represenation of results - therapy effect vs peak irradiance, using measured descriptors in real-time*)
therapyEffect={therapyEffect1,therapyEffect2};
therapyEffectVSIrradiance=Transpose[{theranosticsParam2,therapyEffect[[#]],theranosticsParam3Log/3}]&/@Range[2];
Print[therapyEffectVSIrradiance[[2]]];
ListPlot[therapyEffectVSIrradiance[[#,All,1;;2]],PlotRange->All]&/@Range[2]
ListPlot[Table[Style[{theranosticsParam2[[i]],therapyEffect[[#,i]]},Hue[theranosticsParam3Log[[i]]/3+0.15]],{i,1,theranosticsParamLength}],PlotRange->All]&/@Range[2]

(*representation of data with error bars which represent "Std Error" of descriptor2*)
(*possibility is ErrorBarFunction*)
(*ErrorListPlot[Table[Style[{theranosticsParam2[[i]],therapyEffect[[i]]},Hue[theranosticsParam3Log[[i]]/3+0.15]],{i,1,theranosticsParamLength}],PlotRange->All]*)

(**)
Graphics[Table[Style[Line[{{theranosticsParam2[[i]],therapyEffect[[#,i]]-theranosticsParamDescriptors1[[i,6]]/10},{theranosticsParam2[[i]],therapyEffect[[#,i]]+theranosticsParamDescriptors1[[i,6]]/10}}],Hue[theranosticsParam3Log[[i]]/3+0.15]],{i,1,theranosticsParamLength}],PlotRange->All]&/@Range[2]

(*graphic represetation of function based quantitive therapy effect vs laser peak power - irradiance*)
(*linewidth depends ont he experiment type, more therapies on the same site, wider the line*)

(*returns Hue[H,S,V], needed for dose per pixel!*)
Color[{intensity_,parameterValue_}] := Hue[.8(1-parameterValue^1.7(*1.3*))^1.7(*1.5*),1.,1.]; 

errorBar=40;
Graphics[Table[Style[Line[{{theranosticsParam2Pop[[i]],therapyEffect[[2,i]]-theranosticsParamDescriptors1[[i,6]]/errorBar},{theranosticsParam2Pop[[i]],therapyEffect[[2,i]]+theranosticsParamDescriptors1[[i,6]]/errorBar}}],Hue[.8(1-(Abs[theranosticsParam3Log[[i]]/3])^1.7)^1.1,1.,1.],Thickness[theranosticsParamExpApproach[[i]]*0.005]],{i,1,theranosticsParamLength}],PlotRange->All,Frame->True,FrameLabel->{"j [10^10 W/cm^2]","Therapy effect"},LabelStyle->Directive[14,Bold]]

(*legend color scheme*)
Show[Table[Graphics[{RGBColor[i,0,1-i],Line[{{0,i},{0.2,i}}]}],{i,0,1,0.001}]]

rowsToInclude=Drop[Table[i,{i,1,theranosticsParamLength}],{59,61}];
Graphics[Table[Style[{RGBColor[Abs[theranosticsParam3Log[[i]]/2.5],0,1-Abs[theranosticsParam3Log[[i]]/2.5],0.75](*ColorData["FallColors"][Abs[theranosticsParam3Log[[i]]/3]]*),Line[{{theranosticsParam2Pop[[i]],therapyEffect[[2,i]]-theranosticsParamDescriptors1[[i,6]]/errorBar},{theranosticsParam2Pop[[i]],therapyEffect[[2,i]]+theranosticsParamDescriptors1[[i,6]]/errorBar}}]},Thickness[theranosticsParamExpApproach[[i]]*0.02]],{i,rowsToInclude}],PlotRange->All,Frame->True,FrameLabel->{"j [(10^10) W/(cm^2)]","Therapy effect"},LabelStyle->Directive[18,Bold]]

rowsToInclude=Drop[Table[i,{i,1,theranosticsParamLength}],{59,61}];
Graphics[Table[Style[{Hue[0.6+0.2*(Abs[theranosticsParam3Log[[i]]/3.1])](*ColorData["FallColors"][Abs[theranosticsParam3Log[[i]]/3]]*),Line[{{theranosticsParam2Pop[[i]],therapyEffect[[2,i]]-theranosticsParamDescriptors1[[i,6]]/errorBar},{theranosticsParam2Pop[[i]],therapyEffect[[2,i]]+theranosticsParamDescriptors1[[i,6]]/errorBar}}]},Thickness[theranosticsParamExpApproach[[i]]*0.02]],{i,rowsToInclude}],PlotRange->All,Frame->True,FrameLabel->{"j [(10^10) W/(cm^2)]","Therapy effect"},LabelStyle->Directive[18,Bold]]


Graphics[Table[Style[{ColorData["FallColors"][Abs[theranosticsParam3Log[[i]]/3.1]],Line[{{theranosticsParam2Pop[[i]],therapyEffect[[2,i]]-theranosticsParamDescriptors1[[i,6]]/errorBar},{theranosticsParam2Pop[[i]],therapyEffect[[2,i]]+theranosticsParamDescriptors1[[i,6]]/errorBar}}]},Thickness[theranosticsParamExpApproach[[i]]*0.02]],{i,rowsToInclude}],PlotRange->All,AspectRatio->0.6,Frame->True,FrameLabel->{"j [(10^10) W/(cm^2)]","Therapy effect"},LabelStyle->Directive[22,Bold]]

Graphics[
Table[Style[{ColorData["FallColors"][Abs[theranosticsParam3Log[[i]]/3.1]](*ColorData["Rainbow"][Abs[theranosticsParam3Log[[i]]/3]]*),Line[{{theranosticsParam2Pop[[i]],therapyEffect[[2,i]]-theranosticsParamDescriptors1[[i,6]]/errorBar},{theranosticsParam2Pop[[i]],therapyEffect[[2,i]]+theranosticsParamDescriptors1[[i,6]]/errorBar}}]},Thickness[theranosticsParamExpApproach[[i]]*0.15]],{i,{59,60,61,70,71,72}}],PlotRange->{{4.2,5},{0,3}},Frame->True,FrameTicks->{{{0,0.5,1,1.5,2,2.5,3,3.5},None},{{4.2,4.6,5},None}},(*ImageSize->{40,100},*)FrameLabel->{"j [(10^10) W/(cm^2)]","Therapy effect"},LabelStyle->Directive[22,Bold]]

(*********************************************************************************************************)
(*representation of data in other way..x-axis (irradiance), y - axis (Log(dose)), color is therapy efffect*)
(*********************************************************************************************************)
Graphics[Table[Style[{ColorData["FallColors"][Abs[therapyEffect[[2,i]]/1.1]],Line[{{theranosticsParam2Pop[[i]],theranosticsParam3Log[[i]]-theranosticsParam3Log[[i]]/errorBar},{theranosticsParam2Pop[[i]],theranosticsParam3Log[[i]]+theranosticsParam3Log[[i]]/errorBar}}]},Thickness[theranosticsParamExpApproach[[i]]*0.008]],{i,rowsToInclude}],PlotRange->All,Frame->True,FrameLabel->{"j [(10^10) W/(cm^2)]","(Subscript[Log, 10])(Dose) [μJ/pixel]"},LabelStyle->Directive[18,Bold]]



(*Graphics[Table[{RGBColor[i/5,0,0],Line[{{i,i-0.1},{i,i+0.1}}]},{i,1,5}]]
Graphics[Table[{ColorData["BlueGreenYellow"][i/5],Line[{{i,i-0.1},{i,i+0.1}}]},{i,1,5}]]*)

(*compariosn of qualitative and quatitative therapy effect*)
qualitativeVSquantitativeTherapyEffect=Transpose[{theranosticsParam4,therapyEffect[[#]]}]&/@Range[2];
Print[{"qualitativeVSquantitativeTherapyEffect2",qualitativeVSquantitativeTherapyEffect[[2]]}];

ListPlot[qualitativeVSquantitativeTherapyEffect[[#]],PlotRange->All]&/@Range[2]