Compatibility with latest Mathematica version 13.3.1

[rudolfoldenbourg]

Hello Sjoerd Smit,

I am interested in Dual Numbers and Automatic differentiation for use in solving the inverse problem in non-linear imaging (polarized light microscopy). I am a member of the Wolfram Community, where you can learn more about my interests by exploring a notebook I created on polarized light ray tracing: https://community.wolfram.com/web/rudolfo/home

I came across your Wolfram Language paclet DualNumbers and installed the latest version (DualNumbers-1.2.0.1.paclet) for use in Mathematica's latest version 13.3.1 on my MacBokPro with M1 Max chip running the latest MacOS software version 13.6.

When using your example notebook Example_code.nb to test the installation, I am running into the following problems:

1. After installing the cloned repository, I load the paclet using <<DualNumbers`:

While the error message puzzles me, it seems the installation has not failed altogether. For example, it evaluates the following cell from the example section correctly:

However, in the earlier example of expanding a function f as a series in epsilon, DualSimplify[] does not evaluate correctly:

For testing if the current Mathematica version 13.3.1 might cause the problem, I installed Mathematica version 12.3.1, which is available for installation on a computer with the Apple M1 chip. Using this earlier version of Mathematica and your example notebook Example_code.nb, I reinstalled the repository and loaded the paclet using <<DualNumbers`. This time there was no error message about the SiegelTheta function and I was hopeful that the automatic differentiation of f[a+b epsilon] might work this time. However, DualSimplify gave the same Out[] line again, f[a+b epsilon].

I am very impressed by your programming skills and the amount of work you have contributed in creating the paclet, which by the way is a new programming paradigm for me. I hope your creation can be adapted to work with the current version of Mathematica on its various platforms.

Any information or hint you can give will be greatly appreciated!

Kind regards, Rudolf

[SjoerdSmitWolfram]

Hi Rudolf,

Thanks for checking out my code and reporting on your findings. A while ago I made a few changes that I never officially released and it seems like these address the problems you mentioned. I just released the 1.3.0 version that has all of these changes.

Please let me know if you find any other issues. I haven't really looked at this code in a while, but if new releases of Mathematica broke anything I'd happily take a look and fix them :)

[rudolfoldenbourg]

Hi Sjoerd, I was very happy to get your immediate response, with fixes to the code!, thank you so much! I upgraded the repository, downloaded the updated Example_code.nb and, indeed, find my issues resolved. All the Input cells in the Example section evaluate as expected, except that the Output of a Dual Number doesn’t look like Dual[0.900367, -0.321771], for example, as shown in the unevaluated notebook, but after my evaluation it looks like 0.900367 - 0.321771 epsilon.

I would prefer the representation of a Dual number in an Out[*] cell as Dual[a,b], instead of a+b epsilon. Is there a setting for this?

Apropos preference, I noticed that the In[] and Out[] numbers from the previous session are preserved, upon opening Example_code.nb. Generally, this is a bit confusing to me as I cannot tell if a cell has been evaluated or not in the current session. I would prefer the numbers to be erased upon saving a notebook, as, I think, is the default setting for new notebooks, but I don’t know how to change this setting.

Here are a few additional issues I find:

1. At the end of the Dual Arrays examples in Example_Code.nb, one can find the following entry:

The above entry is shown when opening the notebook. After loading the paclet and evaluating the same input cell, the Out[] cell looks different:

I don’t seem to be able to expand the short form of the DualArray. This might be related to my question of how to change the output form of Dual numbers.

2. When I evaluate all the input cells in the Symbols section, I’m getting a few error messages which I’m attaching here as screen shots. All other cells evaluate as expected:

3. In the section on Running tests, a test report is generated, together with some error messages: 

What is your take on these issues?

Thanks again for your quick response!

Best, Rudolf

[SjoerdSmitWolfram]

I would prefer the representation of a Dual number in an Out[*] cell as Dual[a,b], instead of a+b epsilon. Is there a setting for this?

Ah, I added that typesetting form at the request of another user, but never added a toggle for that. If it annoys you, you can evaluate

MakeBoxes[d : Dual[a_, b_], form : StandardForm] =.

to disable the custom typesetting. This will also disable the typesetting of packed dual arrays. Also, there's always InputForm and FullForm that you can use to inspect the internals of an expression with custom typesetting (e.g., InputForm[SparseArray[{0, 1, 3}]])

Cell labels

This one can be changed by evaluating SetOptions[InputNotebook[], CellLabelAutoDelete -> True] and saving the notebook. You can also set that option globally in the options inspector.

Symbols section

Hmmm, I guess that the usage message typesetting somehow clashes with the custom typesetting I added. Well, it's nothing serious. I might see if I can fix it sometime, but it's shouldn't affect anything.

Tests

Oh, that is weird. I thought I ran the tests successfully before releasing this version, but now something is failing. I'll investigate.

[SjoerdSmitWolfram]

Ok, I figured out what was happening with the tests. It was a global variable (g) that is set in the same notebook that causes this problem. I added a line that clears all global variables before running the tests to prevent this from happening. See Example_code.nb on main.

[rudolfoldenbourg]

Thank you, Sjoerd, for these fixes, which resolve all my issues, except for the error messages in the Symbols section, which are non-essential. Now comes the hard part for us in figuring out how the dual number frame work can be used to invert our forward model for polarized light ray tracing. If you have any comment or suggestion, we, i.e. a team of colleagues at the University of Chicago and the Marine Biological Laboratory in Woods Hole, MA, will be grateful to receive them at my address rudolfo@mbl.edu. Thanks again, Rudolf

[rudolfoldenbourg]

Hi Sjoerd, In my application for Dual numbers, I need to calculate the partial derivative of a multi-argument function, like you show in the DualNumbers_Example_code.nb:

My question is about how to get to the individual derivatives in front of b1 and b2, and, as in my case, of many more individual arguments (hundreds to thousands). I am looking for an efficient way of creating a list of the partial derivatives, as they are shown in the Out[4] cell.

In the Readme text you write: “So if you want to calculate the full gradient of g with dual numbers, you can either use symbolic nonstandard parts for the input arguments and then collect the coefficients afterwards or…” How do I efficiently collect the coefficients afterwards?

Thanks for any hint you can give.

Rudolf

[SjoerdSmitWolfram]

Ok, try this as an example:

bigFun[x_] := (SeedRandom[1]; x . RandomReal[1, Length[x]*{1, 1}] . Exp[x])
bigFun[RandomReal[1, 5]]
bigFun[Array[x, 5]] (* symbolic evaluation *)

Clear[d]
addSymbolicNonstandardPart[vec_List] := Dual[vec, Array[d, Length[vec]]]

input = addSymbolicNonstandardPart[Table[RandomReal[], 5]]
input // FullForm

poly = Expand @ NonStandard[bigFun[input]] (* get the symbolic gradient as a polynomial *)

(* Convert the polynomial to a vector *)
Last @ CoefficientArrays[
   poly,
   NonStandard[input]
] // Normal

[rudolfoldenbourg]

Thank you Sjoerd. It seems when passing a Dual number list into bigFun[x_List], the functions remains unevaluated and its NonStandard part is zero.

[ssmit1986]

Oh, my mistake. Shouldn't have over-edited that bit. Just replace that with bigFun[x_] := ...

[rudolfoldenbourg]

Thanks, it works now. Mathematica keeps surprising/impressing me with useful, high level functions like CoefficientArrays[]

[rudolfoldenbourg]

Hi Sjoerd, I prepared a simple notebook that illustrates a problem I have run into when using DualNumbers in matrix calculations of the Jones Calculus. As I mentioned before, we develop polarized light ray tracing methods and I hope to overcome this problem for implementing an optimization algorithm to reconstruct birefringence parameters in a 3D volume based on polarized light field images. I would be very grateful if you could have a look at this notebook and give me a hint that could help me overcome the problem. I am at a loss for why the code fails in producing the expected result when using DualNumbers, while standard numbers do fine. It would be great if you could allow me to upload this notebook which is less than 100kB. Thanks, Rudolf

[SjoerdSmitWolfram]

Feel free to upload the notebook here. I'll take a look when I have time.

[rudolfoldenbourg]

It seems you have to give me permission to attach a file

[SjoerdSmitWolfram]

Oh, I didn't know this. Can you share it with a Google drive or Dropbox link or something like that?

[rudolfoldenbourg]

Here it comes. Does it work?



On Nov 1, 2023, at 5:09 AM, Sjoerd Smit @.***> wrote:

Oh, I didn't know this. Can you share it with a Google drive or Dropbox link or something like that?

— Reply to this email directly, view it on GitHub https://github.com/ssmit1986/DualNumbers/issues/6#issuecomment-1788635281, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAWCVQTRR5Q7PZQMREVO6OLYCIGTTAVCNFSM6AAAAAA5G2TB7WVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTOOBYGYZTKMRYGE. You are receiving this because you authored the thread.

[SjoerdSmitWolfram]

That's not a link? I'm not sure what that is.

[rudolfoldenbourg]

Try this:

https://www.dropbox.com/scl/fi/myr0rnenjcrwez55cbcn2/Note-to-Sjoerd-Smit.nb?rlkey=i890w9eepopk3z80rnqraynkx&dl=0 Note to Sjoerd Smit.nb dropbox.com

On Nov 1, 2023, at 5:18 AM, Sjoerd Smit @.***> wrote:

That's not a link? I'm not sure what that is.

— Reply to this email directly, view it on GitHub https://github.com/ssmit1986/DualNumbers/issues/6#issuecomment-1788646177, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAWCVQS4AF7JGHZO2NRNEE3YCIHWLAVCNFSM6AAAAAA5G2TB7WVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTOOBYGY2DMMJXG4. You are receiving this because you authored the thread.

[SjoerdSmitWolfram]

A couple of observations:

1. The function Standard is only for applying directly to dual numbers. You can use StandardAll to replace all instances of dual numbers with their standard part in an expression. However, the fact that the returned expression is not manifestly a Dual number indicates some kind of problem.
2. The Re function is not differentiable! I don't know what exactly you need to do here to make things work the way they should, to be honest because this is a problem of setting up the mathematics correctly. The problem is illustrated by this simple example:

D[Re[x], x]
(* Re'[x] *) 

Re[Dual[x, 1]]
(* returns unevaluated *)

As you can see, Re refuses to propagate dual numbers and I honestly don't know what would be the correct way to do it. I recommend you try to re-formulate the formula for jMRet in a way that avoids non-differentiable functions like Re, Im, Conjugate and Arg. If the problem is well-defined, you should be able to do this by placing appropriate conditions on the input arguments.

[SjoerdSmitWolfram]

One more note: for convenience I did implement the relationship Abs[Dual[x, y]] -> Dual[Abs[x], y Sign[x]], but it should be noted that this is only legal for real numbers. So I highly recommend staying within the realm of real numbers and be very very careful if you need to move to complex numbers.

[rudolfoldenbourg]

Thank you, Sjoerd, these are helpful observations and comments, which lead me to explore another approach for polarized light ray tracing. We have mostly used what’s called the Jones Calculus, but there is an alternative, the Mueller Calculus. While the Mueller Calculus uses 4 by 4 matrices for representing the anisotropic optical properties of a slab or volume element, for example, instead of the 2 by 2 matrices of the Jones Calculus, the Mueller matrix elements are all real valued, no imaginary or complex numbers, and so are the Stokes vectors which represent the polarization of light rays in the Mueller Calculus. So I think this will be easier to adapt to the DualNumber concept for automatic differentiation.

By the way, we successfully used DualNumbers for significantly speeding up the optimization routine for simulating fluorescence light field imaging, which only requires one real valued parameter per ray, namely its intensity.

Thanks again, Rudolf

On Nov 1, 2023, at 6:04 AM, Sjoerd Smit @.***> wrote:

One more note: for convenience I did implement the relationship Abs[Dual[x, y]] -> Dual[Abs[x], y Sign[x]], but it should be noted that this is only legal for real numbers. So I highly recommend staying within the realm of real numbers and be very very careful if you need to move to complex numbers.

— Reply to this email directly, view it on GitHub https://github.com/ssmit1986/DualNumbers/issues/6#issuecomment-1788702011, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAWCVQXC6I3P75N4YHFW2LDYCINDPAVCNFSM6AAAAAA5G2TB7WVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTOOBYG4YDEMBRGE. You are receiving this because you authored the thread.

[rebcabin]

Just a friendly confirmation that DualNumbers-1.3.0 loads and passes ad-hoc tests on Mathematica 14.0.

[rebcabin]

However, DualNumbers-1.3.0 did not pass the tests included in the paclet. Here is the generated Notebook and the Test Report Archive.zip

[SjoerdSmitWolfram]

Hi all. Thanks for all the input. Unfortunately I currently do not have much time to work on this...