Variable name inconsistency corrected. Trying to add plot to r-of-fx-... script but having trouble converting IScatterer.GetMusp(wavelength) to an array that python understands.

Author: hayakawa16

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

@@ -46,7 +46,7 @@
 # Create measurements using Nurbs-based white Monte Carlo forward solver
 measurementForwardSolver = NurbsForwardSolver()
 # len(measuredROfFx)=(#fxs)x(#wavelengths) flattened
-measuredROfFx=np.concatenate(
+rOfFxMeasured=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
@@ -92,7 +92,7 @@ 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=np.concatenate(
+rOfFxInitialGuess=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
@@ -101,7 +101,7 @@ def residual(valuesSought, wavelengths, fxs, measuredROfFx, forwardSolver):
    residual,
    initialGuess, 
    ftol=1e-9, xtol=1e-9, max_nfev=10000, # max_nfev needs to be integer
-   args=(wavelengths, fxs, measuredROfFx, forwardSolverForInversion), 
+   args=(wavelengths, fxs, rOfFxMeasured, forwardSolverForInversion), 
    method='lm')
 # Calculate final reflectance from model at fit values
 chromophoresFit = Array.CreateInstance(IChromophoreAbsorber, 2)
@@ -110,25 +110,44 @@ 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=np.concatenate(
+rOfFxFit=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
+# plot Reflectance: flattened so have to separate
 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]
+meas= [m for m in rOfFxMeasured]
 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]
+ig = [i for i in rOfFxInitialGuess]
 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]
+conv = [f for f in rOfFxFit]
+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)
+# plot Mus': flattened so have to separate
+chart = go.Figure()
+xLabel = "wavelength [nm]"
+yLabel = "us'(wavelength)"
+wvs = [w for w in wavelengths]
+# plot measured data first fx first
+meas= [m for m in scattererMeasuredData]
+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 [scattererInitialGuess]]
+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 [scattererFit]]
 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)

inverse-solutions/r-of-rho-multi-op-prop-inversion-using-csharp-solve.py

@@ -45,7 +45,7 @@
 opsMeasured = Tissue(chromophoresMeasuredData, scatterer, "", n=1.4).GetOpticalProperties(wavelengths)
 # Create measurements using Nurbs-based white Monte Carlo forward solver
 measurementForwardSolver = NurbsForwardSolver()
-measuredROfRho = measurementForwardSolver.ROfRho(opsMeasured, rho)
+rOfRhoMeasured = measurementForwardSolver.ROfRho(opsMeasured, rho)
 # Create a forward solver as a model function for inversion
 forwardSolverForInversion = PointSourceSDAForwardSolver()
 
@@ -77,53 +77,53 @@ def CalculateReflectanceFuncVsWavelengthFromChromophoreConcentration(
 chromophoresInitialGuess[2] = ChromophoreAbsorber(ChromophoreType.H2O, initialGuess[2])
 # Compose tissue for initial guess data to obtain OPs
 opsInitialGuess = Tissue(chromophoresInitialGuess, scatterer, "", n=1.4).GetOpticalProperties(wavelengths)
-initialGuessROfRho = forwardSolverForInversion.ROfRho(opsInitialGuess, rho)
+rOfRhoInitialGuess = forwardSolverForInversion.ROfRho(opsInitialGuess, rho)
 # Run the levenberg-marquardt inversion
 optimizer = MPFitLevenbergMarquardtOptimizer()
 initialGuessCopy = Array.CreateInstance(float, 3)
 initialGuessCopy = initialGuess
 parametersToFit = Array.CreateInstance(bool, 3)
 parametersToFit = [True, True, True]
-measuredDataWeight = Array.CreateInstance(float, len(measuredROfRho))
+measuredDataWeight = Array.CreateInstance(float, len(rOfRhoMeasured))
 measuredDataWeight = [1] * len(measuredDataWeight)
 params = Array.CreateInstance(Object, 3)
 params = [wavelengths, rho, scatterer]
 
 # try calling our LM
-fit = optimizer.Solve(initialGuessCopy, parametersToFit, measuredROfRho, measuredDataWeight, 
+fit = optimizer.Solve(initialGuessCopy, parametersToFit, rOfRhoMeasured, measuredDataWeight, 
    forward_func, params)
 # Calculate final reflectance from model at fit values
 chromophoresFit = Array.CreateInstance(IChromophoreAbsorber, 3)
 chromophoresFit[0] = ChromophoreAbsorber(ChromophoreType.HbO2, fit[0])
 chromophoresFit[1] = ChromophoreAbsorber(ChromophoreType.Hb, fit[1])
 chromophoresFit[2] = ChromophoreAbsorber(ChromophoreType.H2O, fit[2])
 opsFit = Tissue(chromophoresFit, scatterer, "", n=1.4).GetOpticalProperties(wavelengths)
-fitROfRho= forwardSolverForInversion.ROfRho(opsFit, rho)
+rOfRhoFit= forwardSolverForInversion.ROfRho(opsFit, rho)
 # plot the results using Plotly
 xLabel = "wavelengths [nm]"
 yLabel = "R(wavelength) [mm-2]"
 wvs = [w for w in wavelengths]
 # plot measured data
-meas = [m for m in measuredROfRho]
+meas = [m for m in rOfRhoMeasured]
 chart = go.Figure()
 chart.add_trace(go.Scatter(x=wvs, y=meas, mode='markers', name='measured data'))
 # plot initial guess data
-ig = [i for i in initialGuessROfRho]
+ig = [i for i in rOfRhoInitialGuess]
 chart.add_trace(go.Scatter(x=wvs, y=ig, mode='markers', name='initial guess'))
 # plot fit
-conv = [f for f in fitROfRho]
+conv = [f for f in rOfRhoFit]
 chart.add_trace(go.Scatter(x=wvs, y=conv, mode='lines', name='converged'))
 chart.update_layout( title="ROfRho (inverse solution for chromophore concentrations, multiple wavelengths, single rho)", xaxis_title=xLabel, yaxis_title=yLabel)
 chart.show(renderer="browser")
 # output results
 print("Meas =    [%5.3f %5.3f %5.3f]" % (measuredData[0], measuredData[1], measuredData[2]))
 print("IG   =    [%5.3f %5.3f %5.3f] Chi2=%5.3e" % (
                 initialGuess[0], initialGuess[1], initialGuess[2],
-                np.dot(np.subtract(measuredROfRho,initialGuessROfRho),
-                       np.subtract(measuredROfRho,initialGuessROfRho))))
+                np.dot(np.subtract(rOfRhoMeasured,rOfRhoInitialGuess),
+                       np.subtract(rOfRhoMeasured,rOfRhoInitialGuess))))
 print("Conv =    [%5.3f %5.3f %5.3f] Chi2=%5.3e" % (fit[0], fit[1], fit[2],
-                np.dot(np.subtract(measuredROfRho,fitROfRho),
-                       np.subtract(measuredROfRho,fitROfRho))))
+                np.dot(np.subtract(rOfRhoMeasured,rOfRhoFit),
+                       np.subtract(rOfRhoMeasured,rOfRhoFit))))
 print("error =   [%5.3f %5.3f %5.3f]%%" % (
                 100*abs((measuredData[0]-fit[0])/measuredData[0]),
                 100*abs((measuredData[1]-fit[1])/measuredData[1]),

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

@@ -44,7 +44,7 @@
 opsMeasured = Tissue(chromophoresMeasuredData, scatterer, "", n=1.4).GetOpticalProperties(wavelengths)
 # Create measurements using Nurbs-based white Monte Carlo forward solver
 measurementForwardSolver = NurbsForwardSolver()
-measuredROfRho= measurementForwardSolver.ROfRho(opsMeasured, rho)
+rOfRhoMeasured = measurementForwardSolver.ROfRho(opsMeasured, rho)
 # Create a forward solver as a model function for inversion
 forwardSolverForInversion = PointSourceSDAForwardSolver()
 
@@ -84,46 +84,46 @@ def residual(chromophoreConcentration, wavelengths, rho, scatterer, measuredROfR
 chromophoresInitialGuess[2] = ChromophoreAbsorber(ChromophoreType.H2O, initialGuess[2])
 # Compose tissue for initial guess data to obtain OPs
 opsInitialGuess = Tissue(chromophoresInitialGuess, scatterer, "", n=1.4).GetOpticalProperties(wavelengths)
-initialGuessROfRho = forwardSolverForInversion.ROfRho(opsInitialGuess, rho) 
+rOfRhoInitialGuess= forwardSolverForInversion.ROfRho(opsInitialGuess, rho) 
 # Run the levenberg-marquardt inversion
 from scipy.optimize import least_squares
 fit = least_squares(
    residual,
    initialGuess, 
    ftol=1e-9, xtol=1e-9, max_nfev=10000, # max_nfev needs to be integer
-   args=(wavelengths, rho, scatterer, measuredROfRho, forwardSolverForInversion), 
+   args=(wavelengths, rho, scatterer, rOfRhoMeasured, forwardSolverForInversion), 
    method='lm')
 # Calculate final reflectance from model at fit values
 chromophoresFit = Array.CreateInstance(IChromophoreAbsorber, 3)
 chromophoresFit[0] = ChromophoreAbsorber(ChromophoreType.HbO2, fit.x[0])
 chromophoresFit[1] = ChromophoreAbsorber(ChromophoreType.Hb, fit.x[1])
 chromophoresFit[2] = ChromophoreAbsorber(ChromophoreType.H2O, fit.x[2])
 opsFit = Tissue(chromophoresFit, scatterer, "", n=1.4).GetOpticalProperties(wavelengths)
-fitROfRho= forwardSolverForInversion.ROfRho(opsFit, rho)
+rOfRhoFit= forwardSolverForInversion.ROfRho(opsFit, rho)
 # plot the results using Plotly
 xLabel = "wavelength [nm]"
 yLabel = "R(wavelength)"
 wvs = [w for w in wavelengths]
 # plot measured data
-meas = [m for m in measuredROfRho]
+meas = [m for m in rOfRhoMeasured]
 chart = go.Figure()
 chart.add_trace(go.Scatter(x=wvs, y=meas, mode='markers', name='measured data'))
 # plot initial guess data
-ig = [i for i in initialGuessROfRho]
+ig = [i for i in rOfRhoInitialGuess]
 chart.add_trace(go.Scatter(x=wvs, y=ig, mode='markers', name='initial guess'))
 # plot fit
-conv = [f for f in fitROfRho]
+conv = [f for f in rOfRhoFit]
 chart.add_trace(go.Scatter(x=wvs, y=conv, mode='lines', name='converged'))
 chart.update_layout( title="ROfRho (inverse solution for chromophore concentrations, multiple wavelengths, single rho)", xaxis_title=xLabel, yaxis_title=yLabel)
 chart.show(renderer="browser")
 # output results
 print("Meas =    [%5.3f %5.3f %5.3f]" % (measuredData[0], measuredData[1], measuredData[2]))
 print("IG   =    [%5.3f %5.3f %5.3f] Chi2=%5.3e" % (initialGuess[0], initialGuess[1], initialGuess[2],
-                np.dot(np.subtract(measuredROfRho,initialGuessROfRho),
-                       np.subtract(measuredROfRho,initialGuessROfRho))))
+                np.dot(np.subtract(rOfRhoMeasured,rOfRhoInitialGuess),
+                       np.subtract(rOfRhoMeasured,rOfRhoInitialGuess))))
 print("Conv =    [%5.3f %5.3f %5.3f] Chi2=%5.3e" % (fit.x[0], fit.x[1], fit.x[2],
-                np.dot(np.subtract(measuredROfRho,fitROfRho),
-                       np.subtract(measuredROfRho,fitROfRho))))
+                np.dot(np.subtract(rOfRhoMeasured,rOfRhoFit),
+                       np.subtract(rOfRhoMeasured,rOfRhoFit))))
 print("error =   [%5.3f %5.3f %5.3f]%%" % (
                 100*abs((measuredData[0]-fit.x[0])/measuredData[0]),
                 100*abs((measuredData[1]-fit.x[1])/measuredData[1]),