Skip to content

Add convoluted Gaussian and Exponential Component#228

Open
jlaehne wants to merge 5 commits intoLumiSpy:mainfrom
jlaehne:tr-expression
Open

Add convoluted Gaussian and Exponential Component#228
jlaehne wants to merge 5 commits intoLumiSpy:mainfrom
jlaehne:tr-expression

Conversation

@jlaehne
Copy link
Contributor

@jlaehne jlaehne commented Apr 1, 2025

Description of the change

Add a component that contains the convolution of an exponential decay with a Gaussian instrument response function.

Progress of the PR

  • Change implemented (can be split into several points),
  • docstring updated (if appropriate),
  • update user guide (if appropriate),
  • added tests,
  • add a changelog entry in the upcoming_changes folder (see upcoming_changes/README.rst),
  • Check formatting of the changelog entry (and eventual user guide changes) in the docs/readthedocs.org:lumispy build of this PR (link in github checks),
  • ready for review.

Minimal example of the bug fix or the new feature

import lumispy as lum
m = s.create_model()
g1 = lum.components.ConvGaussianExponential()
m.append(g1)
m.fit()

@jlaehne
Copy link
Contributor Author

jlaehne commented Apr 1, 2025

@aidanc151 example of adding a component to lumispy :-)

@jlaehne jlaehne force-pushed the tr-expression branch 2 times, most recently from c1675b7 to 3237cfb Compare April 2, 2025 21:17
Copy link
Contributor

@ericpre ericpre left a comment

Choose a reason for hiding this comment

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

It should be good to add this component to the document - I was looking for it to check for the rendering of the analytical expression!

np.testing.assert_allclose(g.function(35), 1.3874364)


class TestConvGaussExp:
Copy link
Contributor

Choose a reason for hiding this comment

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

Add a test that check that there is no normalisation issue: for example by making a dataset where the Gaussian convolution is negligible and checking that the fitting parameter are very close for the ConvGaussExp and Exponential components.

):
super().__init__(
expression="1/2*height*exp(-((x-t0) - sigma**2/(2*tau))/tau)*\
(1 + erf(((x-t0) - sigma**2/tau)/(sqrt(2)*sigma)))",
Copy link
Contributor

Choose a reason for hiding this comment

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

I didn't check the math but from a quick look (based on https://math.stackexchange.com/questions/685041/convolute-exponential-with-a-gaussian/872561#872561), it seems that it should be a "-"?

Suggested change
(1 + erf(((x-t0) - sigma**2/tau)/(sqrt(2)*sigma)))",
(1 - erf(((x-t0) - sigma**2/tau)/(sqrt(2)*sigma)))",

Could it be replaced by the erfc function?


.. math::

f(x) = \frac{h}{2} \cdot \exp\left[(x - t_0 - \frac{\sigma^2}{2\tau})\tau^{-1}\right]
Copy link
Contributor

Choose a reason for hiding this comment

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

Add a reference where the math comes from?
For example, https://math.stackexchange.com/questions/685041/convolute-exponential-with-a-gaussian/872561#872561 has good pointers?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

So far, it is based on https://zenodo.org/records/6384277 - but we are checking the math again with mathematica. As perspective we will at least add a component with double exponential as well.

Copy link
Contributor

Choose a reason for hiding this comment

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

Out of curiosity, are you using mathematica because sympy can't do it or just haven't tried with sympy?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

A student already did the calculations in mathematica and I did not have the time yet to copy that in sympy :-)

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.

2 participants