Skip to content

Feature/197 powerlawscatterer option to define lambda0 - #200

Merged
hayakawa16 merged 14 commits into
masterfrom
feature/197-powerlawscatterer-option-to-define-lambda0
Aug 11, 2025
Merged

Feature/197 powerlawscatterer option to define lambda0#200
hayakawa16 merged 14 commits into
masterfrom
feature/197-powerlawscatterer-option-to-define-lambda0

Conversation

@hayakawa16

Copy link
Copy Markdown
Member

I made an attempt at solving this by adding Property to PowerLawScatterer and overload to constructor. In Tissue I added overload to GetMusp with parameter lambda0 which checks if Scatterer is PowerLawScatterer and if so call the overload with parameter lambda0. Also needed to add 2 overloads to GetOpticalProperties, one for single wavelength and one for array of wavelengths.

I have not added unit tests. Will do this once the code looks okay.

…ly added lamba0 parameter only needed by PowerLawScatterer (other IScatterers don't need). If this looks good, I will add unit tests.
@hayakawa16
hayakawa16 requested a review from lmalenfant July 1, 2025 23:07
@hayakawa16 hayakawa16 self-assigned this Jul 1, 2025
… modified subsequent constructors to call top constructor with fixed values. Same for GetMusp method. In Tissue cleaned up code.

@lmalenfant lmalenfant left a comment

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I didn't notice these methods before, hopefully I'm correct about being able to have the optional parameter.

Comment thread src/Vts/Modeling/Spectroscopy/Tissue.cs Outdated
/// <param name="wavelength">Wavelength</param>
/// <param name="lambda0">Wavelength normalization factor</param>
/// <returns>The optical properties</returns>
public OpticalProperties GetOpticalProperties(double wavelength, double lambda0)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since GetOpticalProperties is not in the interface, I think you could combine these 2 methods by passing the lambda0 as an optional parameter like we did for GetMusp initially. Then if lambda0 is NaN we call GetMusp without that parameter, otherwise we call it with the parameter.

Comment thread src/Vts/Modeling/Spectroscopy/Tissue.cs Outdated
/// <param name="wavelengths">Wavelength</param>
/// <param name="lambda0">Wavelength normalization factor</param>
/// <returns>The optical properties</returns>
public OpticalProperties[] GetOpticalProperties(double[] wavelengths, double lambda0)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same comment as above, since GetOpticalProperties is not in the interface, I think you could combine these 2 methods by passing the lambda0 as an optional parameter like we did for GetMusp initially. Then if lambda0 is NaN we call GetMusp without that parameter, otherwise we call it with the parameter.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thank you @lmalenfant. I think I understand what you mean. I will update.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have pushed my updates.

Comment thread src/Vts/Modeling/Spectroscopy/Tissue.cs Outdated
{
var opArray = new OpticalProperties[wavelengths.Length];
for (int i = 0; i < wavelengths.Length; i++)
if (double.IsNaN(lambda0))

@lmalenfant lmalenfant Jul 9, 2025

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We could move this if statement inside the loop, it would be even less duplicated code.

I asked CoPilot if it thought that was a good solution but it said not, that it was good the way it was written, so I asked about duplication and it came up with this solution:

for (var i = 0; i < wavelengths.Length; i++)
{
    opArray[i] = double.IsNaN(lambda0)
        ? GetOpticalProperties(wavelengths[i])
        : GetOpticalProperties(wavelengths[i], lambda0);
}

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

After giving CoPilot both options, it finally agreed with me that the code above is the best option:
image

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @lmalenfant, I see that your proposed change is better for code duplication, but it is better for performance? I considered performance when I put the "if" outside the "for loop". Going through an "if" within a for loop when the outcome of the "if" is the same for each instance of the "for" seems like additional processing. Does CoPilot consider performance? I'd be interested what it says.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It does, if you look at the Recommendation, it is taking that into account and it said the difference is insignificant. However, if you think this could affect performance in a negative way due to the way we utilize the code, I am good leaving it as is. We could do some performance checks of our own later when we have more time.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry missed that screen shot before my last comment. Sounds good. I might run some benchmarks just for fun.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Tissue is not referenced by any class in the Vts except the unit tests. It is referenced by the GUI though. I think the timing difference will show when numerous wavelengths on Spectral Panel specified. I'll give this a try.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I found no difference in performance with the GUI. In addition, I found documentation that branch prediction will usually make sure that there is no difference at all in these cases, and the code is more readable.
https://stackoverflow.com/questions/38741599/using-an-if-statement-inside-loop-vs-outside-the-loop
So I have pushed your suggested code modification. Thank you @lmalenfant !

…ith code that has a single "for" and the "if" statement inside. This:

https://stackoverflow.com/questions/38741599/using-an-if-statement-inside-loop-vs-outside-the-loop
states that branch prediction will usually make sure that there is no difference at all in these caess.  And the code is more readable.
… loops with code that has a single "for" and the "if" statement inside. This:"

This reverts commit d754ed2.
…ith code that has a single "for" and the "if" statement inside. This:

https://stackoverflow.com/questions/38741599/using-an-if-statement-inside-loop-vs-outside-the-loop
states that branch prediction will usually make sure that there is no difference at all in these caess.  And the code is more readable.
@hayakawa16

Copy link
Copy Markdown
Member Author

Arg! I inadvertently committed from Vts that was integrated into the Vts.Gui.Wpf. So I had two sources committing to same branch. I reverted commit on Wpf vts, cloned a clean clone of Vts, updated to branch, reapplied commit, and pushed. It worked!

@dcuccia

dcuccia commented Jul 15, 2025

Copy link
Copy Markdown
Contributor

Late to the game here, and sorry if I'm wrong about this, but power laws are scale invariant - no need to pollute the per-call API here at all, just take in the constructor of PowerLawScatterer a lambda0 once (use the current assumption as the default constructor parameter), and then convert its A to the "desired" reference lambda based on the similarity relationship:

A'(λ/λ1)^-b = A(λ/λ0)^-b where A' = A*(λ0​/λ1​)b

Carole, does this track?

@hayakawa16

Copy link
Copy Markdown
Member Author

Hi @dcuccia, thank you for chiming in. If I understand your comment correctly, you feel that there is no need for a lambda0 additional parameter because it can be taken into account by modifying the up-front A coefficient? If so, I see your point.
I can think of two reasons to allow users to specify this lambda0:

  1. The formulation of the Power Law scatterer is often defined with this wavelength normalization factor. Our code in PowerLawScatterer hardcodes lambda0=1000nm and we cite OMLC because OMLC original web pages provided this value? However this value is not listed on the current page. This publication: https://omlc.org/news/dec14/Jacques_PMB2013/Jacques_PMB2013.pdf
    uses 500nm as the normalization factor in Equations (1) and (2). So this form:
    mus'(lambda)=A(lambda/lambda0)^(-b)
    of the equation is standard, but the value of lambda0 is not. In addition, if a user wanted to use Jacques' equations and specify lambda0=500nm, they could not directly, they'd have to use your scaling.
  2. This modification arose because both @janakarana and I have been solving the inverse problem which tries to determine the coefficients A and b or A1, b1, A2, b2 in the extended Power Law:
    A1(lambda/lambda0)^(b1)+A2(lambda/lambda0)^b2.
    What we have found is that if you select lambda0 to be around the midpoint of the wavelengths you are using in the inversion, the solution has better success. Often since the measured data is not perfect (it does not obey the Power Law perfectly), altering lambda0 has influence over the recovered A and b and so playing around with various lambda0 has merit. Now we could determine the appropriate A1' and A2', however it is more straight-forward to use another lambda0 value and solve for A and b directly because initial guesses and final converged values make more intuitive sense since they are consistent with the standard form of the Power Law scatterer equation.

@hayakawa16

Copy link
Copy Markdown
Member Author

Thanks @lmalenfant for your approval! @dcuccia, any responses to my comments above?

@dcuccia

dcuccia commented Aug 6, 2025

Copy link
Copy Markdown
Contributor

Hey Carole, I'm sorry I hadn't responded. I was hoping to provide an example to make it easy, but I haven't had time to dedicate to this, as I'd been prepping for customer deliveries and conference presentations. My hypothesis was that we could improve the constructors of the solvers, such that the interface could remain unchanged - the necessary configuration/normalization parameters could be stored as values within the class, and used within the methods, without changing/polluting the per-call function signatures.

@lmalenfant

lmalenfant commented Aug 6, 2025

Copy link
Copy Markdown
Member

@dcuccia we didn't change the interface, we added a constructor for PowerLaw only.

@dcuccia

dcuccia commented Aug 6, 2025

Copy link
Copy Markdown
Contributor

@lmalenfant that's great! And, is there a uniform implementation for all uses of PowerLawScatterer functions, with the class-wide normalization parameters used in the calculations (default center wavelength, or custom)?

@hayakawa16

Copy link
Copy Markdown
Member Author

Hi @dcuccia, in PowerLawScatterer class we have added overloads to the constructor and GetMusp to enable default center wavelength or custom. Within the Vts the only class calling PowerLawScatterer is the Tissue class (within the SpectralMapping namespace). We have added overloads to the methods in Tissue class for instance, GetOpticalProperties and GetMusp, that enable default center wavelength or custom.
Is that what you were interested in?

@dcuccia

dcuccia commented Aug 7, 2025

Copy link
Copy Markdown
Contributor

My hypothesis is that, assuming private class fields of wavelength normalization points and prefactors are appropriately captured (via the single constructor, with default parameters being set to the current values and being customizable if desired), we should not need separate overloads of any instance methods. We should be able to calculate mus' based on the A and b and Lambda0 provided - no separate code paths, no separate overloads. As long as we have these three inputs (A, b, lambda0) we should be able to GetMusp regardless of the chosen values.

@hayakawa16

Copy link
Copy Markdown
Member Author

@dcuccia, this code is on a branch and you can see from this PR what code changes we've made. If you feel you have a better solution, please edit code on this branch. I'd be interested in seeing what you have in mind.

@hayakawa16

Copy link
Copy Markdown
Member Author

This is example python code I wrote that drove some of the changes:
# Define parameters to fit [Hb,HbO2,A,b]
measuredData = [28.4, 22.4, 1.2, 1.42]
# Construct absorber
chromophoresMeasuredData = Array.CreateInstance(IChromophoreAbsorber, 2)
chromophoresMeasuredData[0] = ChromophoreAbsorber(ChromophoreType.HbO2, measuredData[0])
chromophoresMeasuredData[1] = ChromophoreAbsorber(ChromophoreType.Hb, measuredData[1])
# Construct scatterer
scattererMeasuredData = PowerLawScatterer(measuredData[2], measuredData[3])
# len(opsMeasured)=4 per number of wavelengths
lambda0=750.0
opsMeasured = Tissue(chromophoresMeasuredData, scattererMeasuredData, "", n=1.4).GetOpticalProperties(wavelengths, lambda0)

I needed to be able to specify lambda0 and also obtain the optical properties via Tissue. GetOpticalProperties calls GetMusp.

@dcuccia

dcuccia commented Aug 7, 2025

Copy link
Copy Markdown
Contributor

I might have time this weekend to chime in further, but my main point is to take lambda0 in the constructor, and leave it out of GetOP

@hayakawa16

Copy link
Copy Markdown
Member Author

I feel your suggestion. I created fix from bottom up. I think you are suggesting top down addition of lambda0. By taking lambda0 in constructor, did you mean PowerLawScatterer constructor and Tissue constructor? I thought taking lambda0 in Tissue constructor didn't make sense because IScatterer could be other than PowerLawScatterer. This is where I get stuck.

@dcuccia

dcuccia commented Aug 7, 2025

Copy link
Copy Markdown
Contributor

From the C# API perspective:

IChromophoreAbsorber[] abs = [
    new ChromophoreAbsorber(ChromophoreType.HbO2, measuredData[0]),
    new ChromophoreAbsorber(ChromophoreType.Hb, measuredData[1]),
];
var (A, b, lambda0) = (measuredData[2], measuredData[3], 750f);
var scat = new PowerLawScatterer(A, b, lambda0); // lambda0 default for constructor can be whatever you want (sometimes best at 1nm for historical compatibility...)
var tiss = new Tissue(abs, scat, "", 1.4);
tiss.GetOpticalProperties(wavelengths); // calls scat.GetOpticalProperties(wavelengths) under the hood, "doesn't care" how scat does it's job...

@hayakawa16

Copy link
Copy Markdown
Member Author

@dcuccia, thanks for your help! I understand concept. I'll try again.

…LawScatterer. I'm not sure why I didn't see this solution earlier. Updated unit tests that test using default lambda0 value (1000nm) and user specified value.
…bda/lambda0)^-b, i.e. simple not extended power law definition,
@hayakawa16

Copy link
Copy Markdown
Member Author

@dcuccia, this is a much improved solution. Now the dependency of PowerLawScatterer on lambda0 is injected upon instantiation of this class as it should be. Overloads in Tissue removed. Thanks for your help!

@hayakawa16

Copy link
Copy Markdown
Member Author

Thank you @lmalenfant!

@dcuccia

dcuccia commented Aug 8, 2025

Copy link
Copy Markdown
Contributor

No problem. :) Nitpick suggestion would be to consolidate all the PowerLawScatterer overloads and instead use default parameters.

@hayakawa16

Copy link
Copy Markdown
Member Author

Like the suggestion. Will do.

…eeded to cover: 1. mus' = A1*(lamda/lambda0)(-b1) + A2*(lambda/lambda0)(-b2) with an optional lambda0, i.e. public PowerLawScatterer(double a, double b, double c, double d, double lambda0 = 1000.0). 2. mus' = a*(lamda/lambda0)^-b with an optional lambda0, i.e.

 public PowerLawScatterer(double a, double b, double lambda0 = 1000.0). I couldn't combine without parameter ambiguity.
@hayakawa16

Copy link
Copy Markdown
Member Author

I needed one overload. There are two use cases I needed to cover:

  1. mus' = A1*(lamda/lambda0)(-b1) + A2*(lambda/lambda0)(-b2) with an optional lambda0, i.e.
    public PowerLawScatterer(double a, double b, double c, double d, double lambda0 = 1000.0)
  2. mus' = a*(lamda/lambda0)^-b with an optional lambda0, i.e.
    public PowerLawScatterer(double a, double b, double lambda0 = 1000.0)
    I couldn't combine without parameter ambiguity.

@sonarqubecloud

Copy link
Copy Markdown

@dcuccia

dcuccia commented Aug 10, 2025

Copy link
Copy Markdown
Contributor

Sounds great

@hayakawa16
hayakawa16 merged commit daa85fa into master Aug 11, 2025
3 checks passed
@hayakawa16
hayakawa16 deleted the feature/197-powerlawscatterer-option-to-define-lambda0 branch August 11, 2025 18:27
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants