72 Implement and test script-style use of the VTS

[dcuccia]

Initial translated samples that demo "script-style" use the Vts, like in the Matlab scripts, as a foundation for consumption in .Net Polyglot and Jupyter Notebooks. For now, all 12 Monte Carlo demos are implemented and 11/18 of the Forward Solver demos are implemented.

Commented here to see how much more we ultimately want to do. IMO those could be added now, or at some time in the future. It is highly likely that these will need to evolve anyway once we implement the notebooks, as we'll want to keep some level of parity on how the data is being sliced/diced and plotted. So, I suggest we consider having additional changes to these in the core library tracked in a separately Item, such as "Enhance and align script demos to those used in notebooks."

Note: It's unclear how/when we should include the Vts.Scripting.Tests in a CI process - it should probably be part of an integration test step, but we don't have that layer of granularity (yet - it's something I'd like to help with). Like the Matlab scripts they currently open up new windows that need to be closed. For now, just knowing that these libraries compile, and that the script tests pass when run ad-hoc is ok with me.

[lmalenfant]

@dcuccia I notice that you are still pushing commits to this PR, I don't want to have to review it multiple times, could you let us know when it is finished and ready for review?

[dcuccia]

Sorry - will do. When you mentioned that you'd probably take a few nights to get to it, I thought I'd squeeze a few more in the initial batch. I'm now close enough with the solvers that it's worth pushing through the rest of them. I'll let you know when it's ready again.

[lmalenfant]

Thanks! Sounds good, I pulled the code locally so I could review it better and I will make some general comments to start with but leave the overall review until you have completed the solvers.

[lmalenfant]

@dcuccia I converted this to a draft, I think that will stop the email notifications until it is ready. I merged my branch with the NSubstitute changes so you will need to merge in master (The button below "Update branch" should be able to do that).

When you are ready for review, you can click the "Ready for Review" button.

If I'm telling you things you already know, please disregard :)

[dcuccia]

Hi @lmalenfant, so sorry for the spam. I actually just finished the demos, so I'm marking it as ready for review now. In addition to a general code review, it would be great if @hayakawa16 could run some or all of the demos to spot-check my work. There was a lot of copy-paste, and I tried to go back and review, but it's a lot of stuff!

[dcuccia]

I'd add one more thing: there's a lot to explore with Plotly, and we'll get a chance to dial in the visualizations, but hopefully (beyond a tweak here/there) that won't inhibit us from closing this initial effort. There is one thing in particular I haven't figured out, which is how to keep colorbars from overlapping. That, and how to do better aspect ratio layout for multiple charts when combined. My guess is that as we go to Python demos, we'll learn what we really want, and then have an opportunity to back-port those styles, etc to the C# side.

[lmalenfant]

@dcuccia the majority of my comments are coming from Resharper suggestions, I think I remember you said you had Resharper so would you like to go though and fix these before our review?

If you do not have Resharper, I would be happy to make these changes for you.

[dcuccia]

Thanks @lmalenfant! I'll go through and make the changes, no problem.

[dcuccia]

Changes applied, except where noted above in conversations.

[lmalenfant]

@dcuccia I am still seeing issues in the files like making variables global, typos and unnecessary string interpolation, I'm not sure if these are caught by Resharper or SonarLint (which is another good plugin that we use). Would you be able to take another pass?

At any time, I am more than willing to dig in, fixing is easier than pointing out all the issues.

[dcuccia]

The constant and interpolation comments are overkill and not important for this use case. Will update the typos.

[dcuccia]

Installing ReSharper to see what else might be work fixing, and will update when I'm done.

[dcuccia]

Ok, all done. @hayakawa16 let me know when you might have a chance to glance through the samples for correctness (labels, axes, etc). No hurry!

[hayakawa16]

Hi, I just pulled latest. I've been following along with the code cleanup but haven't run demos in a while. I set Vts.Scripting to the startup project, then ran with green right arrow. Only one demo ran and I think it is ForwardSolvers Demo18. This looks like the x-axis label should be wavelength, not rho, and likewise y-axis should be log(R(lambda)). How can I run individual demo class so I know what I'm running and so I can inform you as to if I see error in axes labels?

[dcuccia]

You can change the name of the demo to run in the Program.cs file of the script project, or change the constructor to call the RunAll... Function for forward solvers, mc, short course, and then step through them all sequentially

[hayakawa16]

Thanks. I have the following questions: 1) Demo04FluenceOfRhoAndZAndFt: The title is "modulation" rather than "fluence" and there are 2 extra right ")" 2) Demo05PhdOfRhoAndZ: I don't really see "banana" even with a very major log scale. If you pull up WPF GUI and set same settings, a very different plot is shown. First, the isotropic point source under rho=0 is not seen, it looks more like a linear decay of that source under rho=0. and there is no banana halo, is this due to the color scheme? 3) Demo07PhdOfRhoAndZTwoLayer: I see the "banana" here! Although on this and Demo06FluenceOfRhoAndZTwoLayer there is a lighter line at about z=0.7mm, which is not the top layer thickness=5mm. I brought up the WPF GUI and don't see that lighter line. I wonder what is causing it?
4) Demo08AbsorbedEnergyOfRhoAndZ: I don't think the units are correct. The "mm-3" is correct, however is "mW-1" suppose to mean inverse milli-Watts? If so, I don't think this is correct. I would believe "W" watts in numerator, not sure where the "mW" is coming from. 5) Demo10ROfRhoAndTimeMulti: the MATLAB script defines a single rho and this script defines 19, not a problem but on the 2 plots it is hard to see any difference between all of the cases, possibly because it is on log scale? Possibly fewer rho, and/or smaller max t, and/or regular scale rather than log? Also on the units there is a "s-1" which I think should be "ns-1". 6) Demo11ROfFxAndTime and Demo12ROfFxSingle: The former has the y-axis label with absolute value signs, however the latter does not have these signs. In both cases, the plot is of amplitude correct? 7) Demo16ROfFxMultiPowerLaw: interesting plot, at first didn't know what "A" was but then looked at code. Nice! 8) Demo18ROfRhoMultiOpPropInversion: x-axis should be "wavelength" and y-axis "log(R(lambda))". 9) Demo21ROfRhoAndFtTwoLayerMultiOpProp: same comment about units with respect to time, y-axis units should be "ns-1" instead of "s-1"

Great work! A lot of effort here. Thank you. I will work on the demos in MonteCarlo and ShortCourse next.

[dcuccia]

Thanks Carole! I will take a look at all of these, might take me a day or so, depending on how my Sunday goes.

[hayakawa16]

I should mention that my suggested modifications are just suggestions. For example on 8. above I think you've been using "wavelength" as the axis label, so the y-axis of "log(R(wavelength))" makes more sense. I think if it is consistent across plots then it would be good.

[hayakawa16]

For the MonteCarlo folder, all looks terrific! For the ShortCourse folder: 1) Demo01APhotonCountWithFluence: you mentioned this already about the overlapping colorbar. Otherwise, these are great plots! 2) Demo01BAnalogVsDiscreteWithFluence: this puts out 8 tabs. I happen to know these plots well from the short course so I know what they are. This might be a big ask, but at the top where is says "log(Phi(rho,z))", could this be extended to list which random walk process is being used and how big N is? Something like "log(Phi(rho,z)) Analog N=100" or something like that? 3) Demo02AnalogVsContinuousWithReflectance: on the top plot, the y-axis label should just be "log(R(rho)) [mm-2])" and then a length with red line = analog, blue line=CAW somewhere. On bottom plot, I would prefer "Analog" instead of "AAW" because "AW" stands for "absorption weighting" and analog has no absorption weighting.

Otherwise, looking good!

[lmalenfant]

@dcuccia this is looking good, none of my recent comments are critical to fix or even need to be fixed if that was your intention. Mainly pointing them out because I noticed them.

There are some of the local variable names that are a little ambiguous to me (like tIdx0p01), not sure if they could be named to be more clear what they are. Maybe they would be clear to someone who knows the subject better.

[dcuccia]

Carole, I've tried to answer/address your questions inline below

Thanks. I have the following questions: Demo04FluenceOfRhoAndZAndFt: The title is "modulation" rather than "fluence" and there are 2 extra right ")"

Yes, three plots: 1) Fluence at 0GHz, 2) Fluence at 1GHz, 3) the Modulation - the ratio of 1ghz to DC. That make sense?

Demo05PhdOfRhoAndZ: I don't really see "banana" even with a very major log scale. If you pull up WPF GUI and set same settings, a very different plot is shown. First, the isotropic point source under rho=0 is not seen, it looks more like a linear decay of that source under rho=0. and there is no banana halo, is this due to the color scheme?

No, I'm not sure if it's just the colormap. I changed it to Hot, and got this:

But..AHA! It was a typo and computing FluenceOfFxAndZ :) That would do it! Works now:

Thanks a lot for catching this.

Demo07PhdOfRhoAndZTwoLayer: I see the "banana" here! Although on this and Demo06FluenceOfRhoAndZTwoLayer there is a lighter line at about z=0.7mm, which is not the top layer thickness=5mm. I brought up the WPF GUI and don't see that lighter line. I wonder what is causing it?

I don't know either. It's pretty weird, but I don't have a good guess after beard scratching for a little bit.

Demo08AbsorbedEnergyOfRhoAndZ: I don't think the units are correct. The "mm-3" is correct, however is "mW-1" suppose to mean inverse milli-Watts? If so, I don't think this is correct. I would believe "W" watts in numerator, not sure where the "mW" is coming from.

I mulled this one for a bit and thought it would need to be per unit light dose or power put into to the tissue, thinking the input would be in Watts or Joules. But, I guess it's just a "normalized" mm-3 (per depth, per surface area). Fluence is a value normalized to unit input - I guess that input could be a power (W) or energy dose (J)? I always thought of it as a fluence rate - so W would be the right units for a non-normalized rate of Energy Absorption. I'll keep it normalized for sure. Would be interesting to have your and the team's input on that, though.

Demo10ROfRhoAndTimeMulti: the MATLAB script defines a single rho and this script defines 19, not a problem but on the 2 plots it is hard to see any difference between all of the cases, possibly because it is on log scale? Possibly fewer rho, and/or smaller max t, and/or regular scale rather than log? Also on the units there is a "s-1" which I think should be "ns-1".

Good point, thanks. Might have just been copy-paste from the previous example to try to keep all rhos. I guess if you plot log-log, it looks like this:

...but I went ahead and simplified and just put on a linear axis: Look ok?

Demo11ROfFxAndTime and Demo12ROfFxSingle: The former has the y-axis label with absolute value signs, however the latter does not have these signs. In both cases, the plot is of amplitude correct?

Thanks for finding the inconsistency. Removed the absolute value signs - think that would probably just confuse people anyway (we can have phase shifts in the SFD, but we don't ever discuss those in the current demo, and probably only useful to use magnitude in that future case where we're differentiating from spatial phase).

Demo16ROfFxMultiPowerLaw: interesting plot, at first didn't know what "A" was but then looked at code. Nice!

Thanks :)

Demo18ROfRhoMultiOpPropInversion: x-axis should be "wavelength" and y-axis "log(R(lambda))".

Thanks!

Demo21ROfRhoAndFtTwoLayerMultiOpProp: same comment about units with respect to time, y-axis units should be "ns-1" instead of "s-1"

Thanks, fixed that in 3 files

[dcuccia]

@lmalenfant

There are some of the local variable names that are a little ambiguous to me (like tIdx0p01), not sure if they could be named to be more clear what they are. Maybe they would be clear to someone who knows the subject better.

Maybe I'm too close to it, but that's 0p123 is a relatively common way we computational people get around decimals not working well in variable names. Defer to @hayakawa16 ...

[dcuccia]

@hayakawa16

For the ShortCourse folder:

Demo01APhotonCountWithFluence: you mentioned this already about the overlapping colorbar. Otherwise, these are great plots!

Yeah, sorry - I'll spend some time on a bug branch trying to fix these. Too confusing for them to be separate, I thought.

Demo01BAnalogVsDiscreteWithFluence: this puts out 8 tabs. I happen to know these plots well from the short course so I know what they are. This might be a big ask, but at the top where is says "log(Phi(rho,z))", could this be extended to list which random walk process is being used and how big N is? Something like "log(Phi(rho,z)) Analog N=100" or something like that?

Done!

Demo02AnalogVsContinuousWithReflectance: on the top plot, the y-axis label should just be "log(R(rho)) [mm-2])" and then a length with red line = analog, blue line=CAW somewhere. On bottom plot, I would prefer "Analog" instead of "AAW" because "AW" stands for "absorption weighting" and analog has no absorption weighting.

Thanks, fixed! (though, isn't Analog just a binary "all at once" weighting?)

[hayakawa16]

Regarding forward solver demos: thanks for addressing all of my comments! Just a few responses. 1) I did not realize that Demo04FluenceOfRhoAndZAndFt put out 3 plots. Its hard to know what tabs are new when running in succession. Now that I see the 3 plots, I see why the third is titled "Modulation". Does everyone know that it is the ratio of the prior 2 plots? 2) glad my comments helped to debug :) 3) I'm going to see if I can debug the line issue with Demo07PhdOfRhoAndZTwoLayer and Demo06FluenceOfRhoAndZTwoLayer 4) Steve Jacques chapter 5 in the purple book says: Total energy absorbed [J/cm3] can be computed when the energy [J] of the light source is given. When talking about Monte Carlo solutions, he says: It is convenient to consider the Monte Carlo simulation as yielding the fractional density matrix of incident light absorbed, A [1/cm3], in response to one unit of delivered power or energy. I think all of our models assume a normalized source. 5) Demo10ROfRhoAndTimeMulti looks good.

Regarding ShortCourse folder: thanks for all of the fixes, they look great! Analog is indeed a binary estimator and so no weight modification is performed to the photon due to absorption. That is why I didn't like the name "AW" (absorption weighting) used with this random walk process. DAW and CAW random walk processes disallow termination of the photon due to absorption and instead deweight the photon along its trajectory and so "AW" in their name makes sense.

Thanks for your quick turn around!

[dcuccia]

Thanks.

Does everyone know that it is the ratio of the prior 2 plots?

Anyone reading the script will know. In the future I'd like to combine back in one, but suggest we first figure out how to resolve the colorbars issue.

fractional density matrix of incident light absorbed, A [1/cm3], in response to one unit of delivered power or energy.

So it's normalized for use with either power and energy, as we hypothesized, thanks!

Re: analog, I guess I was just saying it is binary absorption weighing: 0 if not terminated, 1 if so. But all fine w/ me re: naming!

[hayakawa16]

On Demo06FluenceOfRhoAndZTwoLayer if I set the number of rhos and zs to 50 instead of 100 the line is not there. It also isn't there if I set the "stop" for the rhos and zs to be 9.9 instead of 19.9. So it isn't a fineness of grid thing. Here is my guess, it is a floating point issue. If the rhos and zs are set to start:0.1, stop:19.1, number:100, then the rhos and zs values are: 0.1, 0.3, 0.5, 0.7, 0.899999999999991, 1.0999999999999999. This is right where problem is. If I set the rhos and zs to start:0.1, stop:20.0, number:100, then there is no line because no rhos and zs values have form x.x999999999999999. Do you buy this? I have pushed these fixes to Demo06 and Demo07 of ForwardSolvers so you can take a look.

[dcuccia]

@hayakawa16 thanks for looking into this. I appreciate the thought, but I'm not sure I buy it. This integer computation gets much closer to the double-precision limit of accuracy, but also produces the same line:

        var rhos = Enumerable.Range(1, 100).Select(i => ((double)(i * 2 - 1)) /10).ToArray(); // range of s-d separations in mm
        var zs = Enumerable.Range(1, 100).Select(i => ((double)(i * 2 - 1)) / 10).ToArray(); // range of depths in mm

(as does this, but it has the same double precision issue):

        var rhos = MathNet.Numerics.Generate.LinearSpaced(100, 0.1, 19.9); // range of s-d separations in mm
        var zs = MathNet.Numerics.Generate.LinearSpaced(100, 0.1, 19.9); // range of depths in mm

And if you zoom in, the data are "real" (and changing along rho!), but an order of magnitude larger than the rest:

[hayakawa16]

I have another idea. This is an SDA solution. One of the key assumptions about when you can apply an SDA solution is if mus'>>mua. The top layer OPs were set to mua=1, mus'=1. This violates this assumption. If you set the top layer OPs to mua=0.1, mus'=1, and bottom mua=0.01, mus'=1, then line is not there I think. What do you think about that?

[dcuccia]

Did a little more experimentation - it seems to only show up for me when the top layer optical properties are mua: 1, musp: 1. Change them slightly, and it will go away. Changing bottom layer properties doesn't matter, it seems.

[hayakawa16]

I reverted my rhos and zs modifications, and modified OPs so that line is not there. I pushed this fix.

[dcuccia]

Ok. Is it worth documenting the bug as a separate ticket? Some singularity with the TwoLayerSDAForwardSolver.StationaryFluence function?

[dcuccia]

And are we otherwise ready to approve/merge?

[dcuccia]

Missed the comment you made earlier, sorry. Yes I agree, if we see it pop up again for some other SDA-congruent optical properties, then maybe it's worth tracking as a bug.

[hayakawa16]

On the first comment, at first I thought we violated the assumption that the top layer needs to be greater than l because the embedded source at l deep needs to be in the top layer. This is well documented in the header of this class and there is code to check for that assumption. I'm not sure we can know a definite lower bound on mus'/mua when the solution becomes invalid, however I guess adding a check if mus'==mua (and I think this check should be done for both layers) could be done and adding documentation to the header in this regard. I can create an issue for that if that is what you had in mind.

[hayakawa16]

This check should be added to all SDA forward solvers if they don't already have it.

[hayakawa16]

Issue generated and assigned to me. I noticed that the check I made to ensure top layer thickness is thick enough is in each TwoLayerSDAForwardSolver method. Might be a better way to do this going forward. I might ask advice from those with stronger coding skills than mine as to best way to implement this check.

[dcuccia]

Thanks Carole! FWIW, I don't think it's related to musp==mua or limits. The original settings of (1,1) were the only thing I could get to generate the problem, and again only for the top layer.

[hayakawa16]

I think the check should be in the SDA forward solvers. Just because we don't see a visual aberration, doesn't mean the solution is valid.

[lmalenfant]

And are we otherwise ready to approve/merge?

Once you and Carole are good with the changes, I will take one last pass and we should be good to merge!

[dcuccia]

I think any other changes would be in a different branch. I think we're good here. Carole?

[hayakawa16]

My technical review is complete. I'll let Lisa approve after her final pass.

[lmalenfant]

I am unable to build the project, I think due to the static abstract in the IDemoScript interface.

This is the error: Using both 'static' and 'abstract' modifiers requires opting into preview features.

[lmalenfant]

@dcuccia can you look into this and then we should be ready to merge:

I am unable to build the project, I think due to the static abstract in the IDemoScript interface.

This is the error: Using both 'static' and 'abstract' modifiers requires opting into preview features.

[dcuccia]

Not sure what's going on here, as the preview features are enabled, per the .csproj file. I thought it might just be your machine, but right-clicking on the project (and running "dotnet build") and have the same problem. But, if I click the green "Play" button with the Vts.Scripting project as the active project, it works - and, I can see the Debug assets are built with a "just now" timestamp. Can you let me know if the same is true for you?

[dcuccia]

That said, I pushed a change that removes the contract for now, until we figure out what's going on. It's not essential to the use of this code, going forward, just a good guardrail for making new examples uniform (which will happen by copy-paste even without this).