The last version flattened data in order [op1fx1 op1fx2 op2fx1 op2fx2...]. This version flattens data in order [op1fx1 op2fx2 op3fx1 ...]. These two version give same results. For both it bothers me that Chi2 converged is larger than initial guess. I thought reordering of flattened data would change something but it didn't. I guess that makes sense.

Author: hayakawa16

inverse-solutions/r-of-fx-multi-op-prop-inversion.py

@@ -45,15 +45,17 @@
 opsMeasured = Tissue(chromophoresMeasuredData, scattererMeasuredData, "", n=1.4).GetOpticalProperties(wavelengths)
 # Create measurements using Nurbs-based white Monte Carlo forward solver
 measurementForwardSolver = NurbsForwardSolver()
-# len(measuredROfFx)=16 (#wavelengths)x(#fxs) organized [op1fx1, op1fx2, op2fx1, op2fx2, ...]
-measuredROfFx= measurementForwardSolver.ROfFx(opsMeasured, fxs)
+# len(measuredROfFx)=(#fxs)x(#wavelengths) flattened
+measuredROfFx=np.concatenate(
+        [np.array(measurementForwardSolver.ROfFx(opsMeasured, fxs[0]), dtype=float), 
+         np.array(measurementForwardSolver.ROfFx(opsMeasured, fxs[1]), dtype=float)])
 # Create a forward solver as a model function for inversion
 forwardSolverForInversion = PointSourceSDAForwardSolver()
 
 # Declare local forward reflectance function that computes reflectance 
 # from chromophores and scatterer values
 # valuesSought = [Hb, HbO2, A, b]
-def CalculateReflectanceVsWavelengthFromChromophoreConcentration(
+def CalculateReflectanceVsWavelengthFromChromophoreConcAndScatterer(
     valuesSought, wavelengths, fxs, forwardSolver):
 # Create a forward solve model function to solve inverse
    forwardSolverForInversion = PointSourceSDAForwardSolver()
@@ -66,19 +68,20 @@ def CalculateReflectanceVsWavelengthFromChromophoreConcentration(
    # Compose local tissue to obtain optical properties
    opsLocal = Tissue(chromophoresLocal, scattererLocal, "", n=1.4).GetOpticalProperties(wavelengths)
    # Compute reflectance for local absorbers
-   modelDataLocal = Array.CreateInstance(float, len(wavelengths) * len(fxs))
-   modelDataLocal = forwardSolver.ROfFx(opsLocal, fxs) 
-   modelDataForReturn= Array.CreateInstance(float, len(wavelengths) * len(fxs))
-   for i in range(0, len(wavelengths) * len(fxs)):
+   modelDataLocal=np.concatenate(
+        [np.array(forwardSolver.ROfFx(opsLocal, fxs[0]), dtype=float), 
+         np.array(forwardSolver.ROfFx(opsLocal, fxs[1]), dtype=float)])
+   modelDataForReturn=Array.CreateInstance(float, len(wavelengths)*len(fxs))
+   for i in range(0, len(fxs)*len(wavelengths)-1):
         modelDataForReturn[i] = modelDataLocal[i]
    return modelDataForReturn
 
 # func for residual
 def residual(valuesSought, wavelengths, fxs, measuredROfFx, forwardSolver):
-   predictedROfFx= CalculateReflectanceVsWavelengthFromChromophoreConcentration(
+   predictedROfFx= CalculateReflectanceVsWavelengthFromChromophoreConcAndScatterer(
       valuesSought, wavelengths, fxs, forwardSolver) 
-   difference = Array.CreateInstance(float,len(wavelengths) * len(fxs))
-   for i in range(0,len(wavelengths) * len(fxs)):
+   difference = Array.CreateInstance(float,len(wavelengths)*len(fxs))
+   for i in range(0, len(fxs)*len(wavelengths)-1):
        difference[i] = predictedROfFx[i] - measuredROfFx[i]
    return difference
 
@@ -90,7 +93,9 @@ def residual(valuesSought, wavelengths, fxs, measuredROfFx, forwardSolver):
 scattererInitialGuess = PowerLawScatterer(initialGuess[2], initialGuess[3])
 # Compose tissue for initial guess data to obtain OPs
 opsInitialGuess = Tissue(chromophoresInitialGuess, scattererInitialGuess, "", n=1.4).GetOpticalProperties(wavelengths)
-initialGuessROfFx = forwardSolverForInversion.ROfFx(opsInitialGuess, fxs) 
+initialGuessROfFx=np.concatenate(
+        [np.array(forwardSolverForInversion.ROfFx(opsInitialGuess, fxs[0]), dtype=float), 
+         np.array(forwardSolverForInversion.ROfFx(opsInitialGuess, fxs[1]), dtype=float)])
 # Run the levenberg-marquardt inversion
 from scipy.optimize import least_squares
 fit = least_squares(
@@ -106,24 +111,27 @@ def residual(valuesSought, wavelengths, fxs, measuredROfFx, forwardSolver):
 scattererFit = PowerLawScatterer(fit.x[2], fit.x[3])
 # Compose tissue for fit data to obtain OPs
 opsFit= Tissue(chromophoresFit, scattererFit, "", n=1.4).GetOpticalProperties(wavelengths)
-fitROfFx = forwardSolverForInversion.ROfFx(opsFit, fxs) 
-# plot the results using Plotly data organization [op1fx1, op1fx2, op2fx1, op2,fx2,...]
+fitROfFx=np.concatenate(
+        [np.array(forwardSolverForInversion.ROfFx(opsFit, fxs[0]), dtype=float), 
+         np.array(forwardSolverForInversion.ROfFx(opsFit, fxs[1]), dtype=float)])
+# plot the results using Plotly results flattened so have to unflat
 chart = go.Figure()
 xLabel = "wavelength [nm]"
 yLabel = "R(wavelength)"
 wvs = [w for w in wavelengths]
 # plot measured data first fx first
 meas= [m for m in measuredROfFx]
-chart.add_trace(go.Scatter(x=wvs, y=meas[::2], mode='markers', name='measured data: fx1'))
-chart.add_trace(go.Scatter(x=wvs, y=meas[1:len(meas):2], mode='markers', name='measured data: fx2'))
+midpoint=len(meas) // 2
+chart.add_trace(go.Scatter(x=wvs, y=meas[:midpoint], mode='markers', name='measured data: fx1'))
+chart.add_trace(go.Scatter(x=wvs, y=meas[midpoint:], mode='markers', name='measured data: fx2'))
 # plot initial guess data
 ig = [i for i in initialGuessROfFx]
-chart.add_trace(go.Scatter(x=wvs, y=ig[::2], mode='markers', name='initial guess: fx1'))
-chart.add_trace(go.Scatter(x=wvs, y=ig[1:len(ig):2], mode='markers', name='initial guess: fx2'))
+chart.add_trace(go.Scatter(x=wvs, y=ig[:midpoint], mode='markers', name='initial guess: fx1'))
+chart.add_trace(go.Scatter(x=wvs, y=ig[midpoint:], mode='markers', name='initial guess: fx2'))
 # plot fit: need to organize by fx
 conv = [f for f in fitROfFx]
-chart.add_trace(go.Scatter(x=wvs, y=conv[::2], mode='lines', name='converged: fx1'))
-chart.add_trace(go.Scatter(x=wvs, y=conv[1:len(conv):2], mode='lines', name='converged: fx2'))
+chart.add_trace(go.Scatter(x=wvs, y=conv[:midpoint], mode='lines', name='converged: fx1'))
+chart.add_trace(go.Scatter(x=wvs, y=conv[midpoint:], mode='lines', name='converged: fx2'))
 chart.update_layout( title="ROfFx (inverse solution for chromophore concentrations, multiple wavelengths, multiple fx)", xaxis_title=xLabel, yaxis_title=yLabel)
 
 chart.show(renderer="browser")