Add function that uses automatic differentiation to calculate derivatives

[MichaelWedel]

This issue was originally TRAC 9782

As described in the following document:

https://github.com/mantidproject/documents/blob/master/Design/Autodiff.md

I would like to have a function interface that uses automatic differentiation instead of numerical differentiation, utilizing "adept":

http://www.met.rdg.ac.uk/clouds/adept/

Since this new function requires adept as a third party library (and two or three very small source code changes to one of the files in order to make it compile in windows), I would appreciate some advice on the best way to integrate this library.

I will use this branch to keep track of my local work, the code will not compile before the "library integration issue" is solved.

If you would like to try anyway, I can explain how to build and install adept (a target to make a shared object file has to be added to the makefile, it only builds a static library by default) on Linux and/or supply the changes required to build on windows.

[MichaelWedel]

@MichaelWedel (2014-07-01T15:23:28):
Re http://trac.mantidproject.org/mantid/ticket/9782. First rough draft

f3703ba

---

@MichaelWedel (2014-07-01T15:23:28):
Re http://trac.mantidproject.org/mantid/ticket/9782. Refined first draft, added tests

9ca66ea

---

@MichaelWedel (2014-07-03T06:47:00):
Example data with 20 Gaussian peaks, 5001 data points.

---

@MichaelWedel (2014-07-03T06:48:41):
Re http://trac.mantidproject.org/mantid/ticket/9782. Temporarily adding adept to MantidAPI CMakeLists.txt

05938a8

---

@MichaelWedel (2014-07-03T06:52:21):
The testing algorithm can be used with the example data (it's just "proof of concept", so it only really does this one thing) in this way:

Load(Filename='WorkspaceData.nxs', OutputWorkspace='WorkspaceData')
AutoDiffTestAlg(InputWorkspace='WorkspaceData', OutputWorkspace='FitParameters', OutputWorkspacePlot='Plot')

Note that I will remove the testing algorithm and function before the ticket is finished.

---

@MichaelWedel (2014-07-03T11:24:51):
Re http://trac.mantidproject.org/mantid/ticket/9782. Cleaning up and documentation in IFunction1DAutoDiff

Added doxygen comments and some general documentation. Cleaned up the interface of AutoDiffJacobianAdapter and added error handling.

a94a28b

---

@MichaelWedel (2014-07-03T12:51:48):
Re http://trac.mantidproject.org/mantid/ticket/9782. Fixing incorrect operator in AutoDiffJacobianAdapter

9b6668c

---

@MichaelWedel (2014-07-03T14:08:48):
From my point of view this ticket is more or less done except for the unsolved question of how to integrate adept itself into Mantid.

---

@MichaelWedel (2014-07-11T09:17:02):
To sum up my packaging and cross-platform adventures with adept.

'''Linux'''

• Modified the makefiles of the original source archive to build an .so file instead of an .a file, so it can be linked dynamically. I thought this dynamic linking approach would be reasonable for Linux systems.
• Created .deb and .rpm packages. I'm not doing this too often, so I don't know whether the packages are okay the way I did it. One thing that could possibly be improved is that the .rpm package could be split up into adept and adept-dev, just as the debian package. On the other hand, there are actually only two file that need to be distributed (libadept.so and adept.h).
• Wrote a small test program and compiled it with -ladept after installing the created package on a clean system and make sure it produces the correct output.

'''Windows'''

• Found out that dynamic linking is not possible here, because __declspec(thread) and __declspec(dllexport) are not compatible, so on Windows there is (probably?) no alternative to static linking.
• Patched adept.h to use correct identifiers (__declspec(thread) instead of __thread, __restrict instead of __restrict__).
• Created a VS solution that builds adept as a static library (adept.lib)
• Created a VS solution that builds the above mentioned example program, linking adept.lib statically and checked the output.

'''Mac OS X'''

• Could not test on Mac OS X, since I have no machine with the proper environment around.

  ---

@MichaelWedel (2014-07-24T12:04:00):
Another update regarding the build system.

'''Windows'''

• Managed to build a DLL, sacrificing the possibility to make a thread safe build. This is not a problem either OpenMP is disabled or when the body of IFunction1DAutoDiff::functionWrapper is put into a PARALLEL_CRITICAL block (tested with a small example program).

'''General'''

• For easier cross-platform builds I decided to create a CMake based build system for the library and the included tests (only those that don't require other libraries, e.g. benchmarks).
• The CMake system works in Linux and Windows with Visual Studio.
• I created a Findadept.cmake file for easy inclusion into the Mantid build system and tested that it actually finds the library (at least on Linux, it does).
• Since these changes added up, I decided for now to have this as a fork of the library. The library itself is only two files, and by putting this into a github repository, the original author can maybe profit from the changes (already contacted him about this).

'''What's left to do'''

• Verify that it builds on Mac OS X. I tried building with old gcc versions agains old standard libraries, this worked okay, so I hope it should be okay on Mac OS X.
• Put readme-files, docs and other additional stuff into the fork repo.
• Adjust packaging scripts (should be much easier now, as there are no patches required in the packaging process itself).

  ---

@MichaelWedel (2014-07-25T08:02:08):
Refs http://trac.mantidproject.org/mantid/ticket/9782. Added CMake find module for adept

244715a

---

@MichaelWedel (2014-07-25T08:02:08):
Refs http://trac.mantidproject.org/mantid/ticket/9782. Added PARALLEL_CRITICAL section.

62cc1e6

---

@MichaelWedel (2014-07-25T08:20:12):
Packaging scripts are adjusted as well. I built packages for Ubuntu 13.10 and Fedora 13 (Don't have RHEL6 available, so I chose an old version - if it builds there, it should work on newer versions as well). CMake build system works on Windows as well.

Only Mac OS X left now.

---

@MichaelWedel (2014-07-31T08:43:24):
Refs http://trac.mantidproject.org/mantid/ticket/9782. Added Gaussian function with numerical derivatives

The same as GaussianAutoDiff, for comparison.

1fbcda1

---

@MichaelWedel (2014-07-31T08:43:24):
Refs http://trac.mantidproject.org/mantid/ticket/9782. Comparison capability on AutoDiffTestAlg

Roman's changes to make a comparison possible between numerical and auto diff.

8054491

---

@MichaelWedel (2014-07-31T08:43:24):
Refs http://trac.mantidproject.org/mantid/ticket/9782. Modifying fit to output number of iterations.

d15ce05

---

@MichaelWedel (2014-08-13T07:09:32):
Refs http://trac.mantidproject.org/mantid/ticket/9782. Added gaussian with hand-coded derivatives

1538ab6

---

@MichaelWedel (2014-08-13T07:17:13):
Refs http://trac.mantidproject.org/mantid/ticket/9782. Removed hard-coded number from AutoDiffTestAlg

0acd294

---

@MichaelWedel (2014-08-13T08:03:51):
'''Update'''

Performance tests conducted by Roman and me showed that (in contrast to my expectation), automatic differentiation is slightly slower for the cases we tested, even with a fair amount of parameters and spectra to fit. I also tested different approaches regarding the integration of adept into the function fitting framework, but none of that lead to significant speedups. Compiling adept with a larger initial stack size helped a bit, since it does not need to grow the stack so often anymore (it involves re-allocation of memory and is rather expensive).

Another argument generally made in favor of automatic differentiation is the better precision of the resulting derivatives, so I wanted to test this as well.

I modified AutoDiffTestAlg to also include a function with hand-coded derivatives and make it more suitable for systematic testing.

The algorithm expects a MatrixWorkspace with 4 spectra, which contain the following data:

1. f(x)
2. df/d(centre)
3. df/d(height)
4. df/d(sigma)

It produces an equivalent workspace with the differences (reference - calculated) from calls to IFunction::function and IFunction::functionDeriv, using the specified method (num, adept, handcoded).

In attachment:Gaussian_y_J_053.nxs there are data for a Gaussian (Centre: 0.0, Height: 4.0, Sigma: 3.5) on the range [-5,5](Produced with Matlab symbolic toolbox). Assuming the file is located in the data search path, the following script generates workspaces for all three methods:

Load(Filename='Gaussian_y_J_053.nxs', OutputWorkspace='WorkspaceData', LoaderName='LoadNexusProcessed', LoaderVersion=1)
AutoDiffTestAlg(InputWorkspace='WorkspaceData', OutputWorkspace="HandCoded",DerivativeType="handcoded")
AutoDiffTestAlg(InputWorkspace='WorkspaceData', OutputWorkspace="Numerical",DerivativeType="num")
AutoDiffTestAlg(InputWorkspace='WorkspaceData', OutputWorkspace="AutoDiff",DerivativeType="adept")

For hand-coded and automatic derivatives the differences are all in the range [-5e-16, 5e-16], but not identical. I suspect differences on that level are not significant taking into account floating point precision (the values are on the order of 1e0 or 1e-1, with 15-17 significant digits I would classify these differences as "tolerable").

In the numerical case all three derivatives behave differently. For df/d(centre) the differences to the correct value are largest with values on the interval [-0.03, 0.03], where the derivatives themselves are on a scale of ca. [-0.8, 0.8]. df/d(height) looks reasonable, with differences of [-5e-14, 5e-14] with a distribution that seems more or less random. The derivations of df/ds are on the interval of [0, 7e-4], with values [0, 0.9]. Here the differences look very systematic, the difference curve looks very similar to the derivative, just 3 orders of magnitudes smaller.

With the script above and the attached data these results should be reproducible.

Overall, the hand-coded and automatic derivatives are comparable and do not show significant derivations from the reference results. Numeric derivatives show differences which can become large in some cases. I know that this is expected, because of truncation and rounding of floating point numbers and depends on many factors such as parameter value range, function value range and possibly other factors that I am forgetting now. I tried to choose "reasonable" example parameter ranges, but I am willing to try some more functions/parameter combinations as well.

---

@MichaelWedel (2014-08-13T12:55:29):
Refs http://trac.mantidproject.org/mantid/ticket/9782. Changed AutoDiffTestAlg

The algorithm now takes function parameters (for calculating values) and measures performance of derivative calculation.

e1440db

---

@MichaelWedel (2014-08-13T12:55:29):
Refs http://trac.mantidproject.org/mantid/ticket/9782. Changed implementation of Gaussians to not use pow

I was not aware how incredibly slow pow() is for squaring a value.

5e7de3e

---

@MichaelWedel (2014-08-13T13:07:41):
I decided to spend some more time on profiling, because it seemed odd that hand-coded derivatives would sometimes be slower than numerical and discovered that using "pow" for squaring is not a good idea, especially with automatic differentiation. This is why numerical derivatives were getting away quite cheaply - the derivatives only called "pow" once, and only the double version, as opposed to the adept-version. Now the relative performance of calculating derivatives directly through IFunction::functionDeriv() is ca. 1 : 2.5 : 2.8 (hand-coded, adept, numerical).

---

@MichaelWedel (2014-08-13T14:37:13):
Refs http://trac.mantidproject.org/mantid/ticket/9782. Added Lorentz-profiles

They all live in a namespace called "Lorentzians".

a20c5fc

---

@MichaelWedel (2014-08-13T14:37:13):
Refs http://trac.mantidproject.org/mantid/ticket/9782. AutoDiffTestAlg uses Lorentzians instead of Gaussians

8b0f0e5

---

@MichaelWedel (2014-08-13T14:44:54):
The Lorentzian is an interesting test-case, because there isn't even a call to exp, it's just basic arithmetic operations. In this case the relative performance is 1 : 1.94 : 6.97, so roughly 1 : 2 : 7 (hand-coded, numerical, adept). It's not what I expected though (what I had written previously was caused by wrong ordering of my output...), so I will investigate further

---

@MichaelWedel (2014-08-18T08:52:03):
Refs http://trac.mantidproject.org/mantid/ticket/9782. AutoDiffTestAlg improved

f620746

---

@MichaelWedel (2014-08-18T08:52:03):
Refs http://trac.mantidproject.org/mantid/ticket/9782. Adding necessary compile time defines

bb5616b

---

@MichaelWedel (2014-08-18T08:52:03):
Refs http://trac.mantidproject.org/mantid/ticket/9782. Const string reference

d96c282

---

@MichaelWedel (2014-08-18T08:52:03):
Refs http://trac.mantidproject.org/mantid/ticket/9782. Added Pearson-VII, changed AutoDiffTestAlg

3b6f16b

---

@MichaelWedel (2014-08-18T09:11:19):
In the last commit I added the Pearson-VII function (as defined here http://pd.chem.ucl.ac.uk/pdnn/peaks/pvii.htm), which has 4 parameters, definitely requires the use of pow and has fairly complex derivatives.

With the attachment data, you can do:

Load(Filename='PearsonVII_y_J_053.nxs', OutputWorkspace='WorkspaceData', LoaderName='LoadNexusProcessed', LoaderVersion=1)
AutoDiffTestAlg(InputWorkspace='WorkspaceData', OutputWorkspace="HandCoded",DerivativeType="handcoded", GaussianParameters="0.665 200.0 0.0007 2.0")
AutoDiffTestAlg(InputWorkspace='WorkspaceData', OutputWorkspace="Numerical",DerivativeType="num", GaussianParameters="0.665 200.0 0.0007 2.0")
AutoDiffTestAlg(InputWorkspace='WorkspaceData', OutputWorkspace="AutoDiff",DerivativeType="adept", GaussianParameters="0.665 200.0 0.0007 2.0")

The derivatives are not as precise as the ones for Gaussian and Lorentzian, not even the hand-coded ones. I guess it has something to do with the implementation of the mathematical functions that are used and in the end, also floating point accuracy.

This is a summary of the performance (on my Ubuntu machine). It's always Hand-coded : Numerical : Adept.

'''Gaussian'''

• Function: 1 : 1 : 1.3446
• Derivatives: 1 : 2.8005 : 1.9767
• Fit (iterations): 6 : 4 : 6

'''Lorentzian'''

• Function: 1 : 1 : 3.3559
• Derivatives: 1 : 1.9318 : 5.2510
• Fit (iterations): 7 : 6 : 7

'''Pearson-VII'''

• Function: 1 : 1 : 1.0908
• Derivatives: 1 : 2.2272 : 1.0737
• Fit (iterations): 6 : 7 : 6

Based on these findings I think it's really worthwhile to include this feature into Mantid.

---

@MichaelWedel (2014-10-31T15:49:40):
Local checkpoint Refs http://trac.mantidproject.org/mantid/ticket/9782

ad174be

---

@MichaelWedel (2014-10-31T15:49:40):
Refs http://trac.mantidproject.org/mantid/ticket/9782. Some cleanup and comments in AutoDiffTestAlg

7f73184

---

@MichaelWedel (2014-10-31T16:26:28):
The main adept integration is in Mantid::API::IFunctionAutoDiff. Because of the way parameters for automatic differentiation work they can not be used as usual. There are examples how to implement concrete functions in Mantid::CurveFitting (GaussianAutoDiff, Lorentzians, Pearsons).

In Mantid::CurveFitting there is a test algorithm, AutoDiffTestAlg, which can be invoked by this script:

Load(Filename='/home/mwedel/Data/TestData/AutoDiff/Gaussian_y_J_d.nxs', OutputWorkspace='WorkspaceData', LoaderName='LoadNexusProcessed', LoaderVersion=1)
AutoDiffTestAlg(InputWorkspace='WorkspaceData', OutputWorkspace="G_HandCoded", FunctionName="Gaussian",DerivativeType="analytical", FunctionParameters="0.665 200.0 0.0007")
AutoDiffTestAlg(InputWorkspace='WorkspaceData', OutputWorkspace="G_Numerical", FunctionName="Gaussian",DerivativeType="num", FunctionParameters="0.665 200.0 0.0007")
AutoDiffTestAlg(InputWorkspace='WorkspaceData', OutputWorkspace="G_AutoDiff", FunctionName="Gaussian",DerivativeType="adept", FunctionParameters="0.665 200.0 0.0007")

Load(Filename='/home/mwedel/Data/TestData/AutoDiff/Lorentzian_y_J_d.nxs', OutputWorkspace='WorkspaceData', LoaderName='LoadNexusProcessed', LoaderVersion=1)
AutoDiffTestAlg(InputWorkspace='WorkspaceData', OutputWorkspace="L_HandCoded", FunctionName="Lorentzian",DerivativeType="analytical", FunctionParameters="0.665 200.0 0.0007")
AutoDiffTestAlg(InputWorkspace='WorkspaceData', OutputWorkspace="L_Numerical", FunctionName="Lorentzian",DerivativeType="num", FunctionParameters="0.665 200.0 0.0007")
AutoDiffTestAlg(InputWorkspace='WorkspaceData', OutputWorkspace="L_AutoDiff", FunctionName="Lorentzian",DerivativeType="adept", FunctionParameters="0.665 200.0 0.0007")

Load(Filename='/home/mwedel/Data/TestData/AutoDiff/PearsonVII_y_J_d.nxs', OutputWorkspace='WorkspaceData', LoaderName='LoadNexusProcessed', LoaderVersion=1)
AutoDiffTestAlg(InputWorkspace='WorkspaceData', OutputWorkspace="P_HandCoded", FunctionName="Pearson",DerivativeType="analytical", FunctionParameters="0.665 200.0 0.0007 2.0")
AutoDiffTestAlg(InputWorkspace='WorkspaceData', OutputWorkspace="P_Numerical", FunctionName="Pearson",DerivativeType="num", FunctionParameters="0.665 200.0 0.0007 2.0")
AutoDiffTestAlg(InputWorkspace='WorkspaceData', OutputWorkspace="P_AutoDiff", FunctionName="Pearson",DerivativeType="adept", FunctionParameters="0.665 200.0 0.0007 2.0")

It records timing information and generates output with accuracy information for the partial derivatives. To use it you need to download the attached data files.

---

@MichaelWedel (2014-12-05T12:14:43):
Unfortunately I have to postpone this one to Release 3.4

[MichaelWedel]

http://trac.mantidproject.org/mantid/raw-attachment/ticket/9782/WorkspaceData.nxs
(uploaded by @MichaelWedel on 2014-07-03T06:47:00)

---

http://trac.mantidproject.org/mantid/raw-attachment/ticket/9782/6904_peaks_table_fortran.nxs
(uploaded by @MichaelWedel on 2014-07-31T11:28:18)

---

http://trac.mantidproject.org/mantid/raw-attachment/ticket/9782/6904_calc_peaks1-4_fortran.nxs
(uploaded by @MichaelWedel on 2014-07-31T11:28:32)

---

http://trac.mantidproject.org/mantid/raw-attachment/ticket/9782/Gaussian_y_J_053.nxs
(uploaded by @MichaelWedel on 2014-08-13T07:10:01)

---

http://trac.mantidproject.org/mantid/raw-attachment/ticket/9782/Lorentzian_y_J_d.nxs
(uploaded by @MichaelWedel on 2014-08-13T14:35:16)

---

http://trac.mantidproject.org/mantid/raw-attachment/ticket/9782/PearsonVII_y_J_d.nxs
(uploaded by @MichaelWedel on 2014-08-18T08:52:38)

---

http://trac.mantidproject.org/mantid/raw-attachment/ticket/9782/Gaussian_y_J_d.nxs
(uploaded by @MichaelWedel on 2014-10-31T15:43:32)

---

[mantid-builder]

http://trac.mantidproject.org/mantid/raw-attachment/ticket/9782/WorkspaceData.nxs
(uploaded by @MichaelWedel on 2014-07-03T06:47:00)

---

[mantid-builder]

http://trac.mantidproject.org/mantid/raw-attachment/ticket/9782/6904_peaks_table_fortran.nxs
(uploaded by @MichaelWedel on 2014-07-31T11:28:18)

---

[mantid-builder]

http://trac.mantidproject.org/mantid/raw-attachment/ticket/9782/6904_calc_peaks1-4_fortran.nxs
(uploaded by @MichaelWedel on 2014-07-31T11:28:32)

---

[mantid-builder]

http://trac.mantidproject.org/mantid/raw-attachment/ticket/9782/Gaussian_y_J_053.nxs
(uploaded by @MichaelWedel on 2014-08-13T07:10:01)

---

[mantid-builder]

http://trac.mantidproject.org/mantid/raw-attachment/ticket/9782/Lorentzian_y_J_d.nxs
(uploaded by @MichaelWedel on 2014-08-13T14:35:16)

---

[mantid-builder]

http://trac.mantidproject.org/mantid/raw-attachment/ticket/9782/PearsonVII_y_J_d.nxs
(uploaded by @MichaelWedel on 2014-08-18T08:52:38)

---

[mantid-builder]

http://trac.mantidproject.org/mantid/raw-attachment/ticket/9782/Gaussian_y_J_d.nxs
(uploaded by @MichaelWedel on 2014-10-31T15:43:32)

---

[stale]

This issue has been automatically marked as stale because it has not had recent activity. It will be closed in 7 days if no further activity occurs. If you feel this is incorrect please comment to keep it alive, with a reason why.

To prevent closure, e.g. for long-term planning issues, add the "Never Stale" label.

[stale]

This issue has been closed automatically. If this still affects you please re-open this issue with a comment so we can look into resolving it.