JeschkeLab / DeerLab

Python package for data analysis for dipolar EPR spectroscopy
https://jeschkelab.github.io/DeerLab
MIT License
14 stars 10 forks source link

Normalization in dd_models._multigaussfun #481

Closed sergeqzin closed 1 day ago

sergeqzin commented 4 months ago

Dear all,

The function dd_models._multigaussfun has a misleading normalization step, in my opinion

if not np.all(P==0):
        # Normalization
        P = np.squeeze(P)/np.sum([np.trapz(c,r) for c in P.T])

In K-Gaussian models, P is a matrix N \times K, where N is the size of the input r-axis. The columns are individual Gaussians normalized to the unit area. After the current normalization step, the integral of each Gaussian is 1/K. This is at least inconsistent with the model description in the documentation. E.g. from https://jeschkelab.github.io/DeerLab/_autosummary/deerlab.dd_gauss2.html#deerlab.dd_gauss2 I conclude that the condition (a1 + a2 = 1) yields a normalized distribution which is not true in reality.

The code below tests the normalization

import deerlab as dl
import numpy as np

r = np.linspace(1, 5, 200)

mean = 3
std = 0.5

P1 = dl.dd_gauss(r, 
                 mean, std
                 ) # expected: area = 1
P2 = dl.dd_gauss2(r, 
                  mean, std, 
                  mean, std, 
                  1/2, 1/2) # expected: area = 1
P3 = dl.dd_gauss3(r, 
                  mean, std, 
                  mean, std, 
                  mean, std, 
                  1/3, 1/3, 1/3) # expected: area = 1

print(np.trapz(P1, r)) # output: 1
print(np.trapz(P2, r)) # output: 0.5
print(np.trapz(P3, r)) # output: 0.3333

Thanks!

sergeqzin commented 4 months ago

Issue 2:

in dd_model.py, in the Section dd_gauss3, the variable notes is misspelled as ntoes. Because of this, the documentation on dd_gauss3 displays a LaTeX formula of dd_gauss2 (as it precedes dd_gauss3 in dd_models.py).

Please, fix this as well. Thank you!

HKaras commented 1 day ago

Upon further investigation the normalisation criteria is working as expected. The amplitude for the all the individual gaussian need to be set to 1. These are amplitudes not weights.

Since changing this would implement a breaking change, I do not propose making the change.

sergeqzin commented 1 day ago

Dear @HKaras

I understand what you are trying to say. Nevertheless, I must point out again that the documentation on these functions disagrees with their implementation. I kindly suggest to fix either of them.