Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update proc_autophase.py #48

Open
wants to merge 12 commits into
base: master
Choose a base branch
from
95 changes: 39 additions & 56 deletions nmrglue/process/proc_autophase.py
Original file line number Diff line number Diff line change
Expand Up @@ -14,46 +14,46 @@
from .proc_base import ps


def autops(data, fn, p0=0.0, p1=0.0):
def psacme(data, p0=0.0, p1=0.0, gamma=5.0e-3):
Copy link
Owner

Choose a reason for hiding this comment

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

I do not think the name of this function should be changed. If it is the name in the documentation must also be changed.

"""
Automatic linear phase correction
Phase correction using ACME algorithm by Chen Li et al.
Journal of Magnetic Resonance 158 (2002) 164-168

Parameters
----------
ph : tuple
Current p0 and p1 values
data : ndarray
Array of NMR data.
fn : str or function
Algorithm to use for phase scoring. Built in functions can be
specified by one of the following strings: "acme", "peak_minima"
p0 : float
Initial zero order phase in degrees.
p1 : float
Initial first order phase in degrees.
gamma : float
scale factor set to balance the penalty and entropy
optreturn : boolean
option to return the optimum p0 and p1

Returns
-------
ndata : ndarray
Phased NMR data.
score : float
Value of the objective function (phase score)
opt : tuple
phases returned by algoritm in a tuple

"""
if not callable(fn):
fn = {
'peak_minima': _ps_peak_minima_score,
'acme': _ps_acme_score,
}[fn]

opt = [p0, p1]
opt = scipy.optimize.fmin(fn, x0=opt, args=(data, ))
opt = scipy.optimize.fmin(_ps_acme_score, x0=opt, args=(data, gamma))

phasedspc = ps(data, p0=opt[0], p1=opt[1])

return phasedspc
return tuple(opt), phasedspc


def _ps_acme_score(ph, data):
def pspmin(data, p0=0.0, p1=0.0):
"""
Phase correction using ACME algorithm by Chen Li et al.
Journal of Magnetic Resonance 158 (2002) 164-168
Phase correction using simple minima-minimisation around highest peak

This is a naive approach but is quick and often achieves reasonable
results. The optimisation is performed by finding the highest peak in the
spectra (e.g. TMSP) and then attempting to reduce minima surrounding it.

Parameters
----------
Expand All @@ -66,60 +66,43 @@ def _ps_acme_score(ph, data):
-------
score : float
Value of the objective function (phase score)
opt : tuple
phases returned by algoritm in a tuple

"""
stepsize = 1

opt = [p0, p1]
opt = scipy.optimize.fmin(_ps_acme_score, x0=opt, args=(data,))

phasedspc = ps(data, p0=opt[0], p1=opt[1])

return tuple(opt), phasedspc


def _ps_acme_score(ph, data, gamma):
phc0, phc1 = ph

s0 = ps(data, p0=phc0, p1=phc1)
data = np.real(s0)

# Calculation of first derivatives
ds1 = np.abs((data[1:]-data[:-1]) / (stepsize*2))
p1 = ds1 / np.sum(ds1)
ds1 = np.diff(data, 1)
p1 = np.abs(ds1) / np.abs(np.sum(ds1))

# Calculation of entropy
p1[p1 == 0] = 1

h1 = -p1 * np.log(p1)
h1s = np.sum(h1)

# Calculation of penalty
pfun = 0.0
as_ = data - np.abs(data)
sumas = np.sum(as_)

if sumas < 0:
pfun = pfun + np.sum((as_/2) ** 2)
fr = data
Copy link
Owner

Choose a reason for hiding this comment

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

This needs to be fr = data.copy() or else the next few lines changes the values in the original data array.

Copy link
Owner

Choose a reason for hiding this comment

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

Looking into this a bit more a boolean comparison can be used here. I believe the following should suffice.

    # Calculation of penalty
    fr = data < 0   # 0 if y >= 0, 1 if y < 0
    pr = kwargs['gamma'] * np.nansum(fr * data**2)

fr[fr >= 0] = 0
fr[fr < 0] = 1
pr = gamma * np.sum(fr * data**2)

p = 1000 * pfun

return h1s + p
return h1s + pr


def _ps_peak_minima_score(ph, data):
"""
Phase correction using simple minima-minimisation around highest peak

This is a naive approach but is quick and often achieves reasonable
results. The optimisation is performed by finding the highest peak in the
spectra (e.g. TMSP) and then attempting to reduce minima surrounding it.

Parameters
----------
pd : tuple
Current p0 and p1 values
data : ndarray
Array of NMR data.

Returns
-------
score : float
Value of the objective function (phase score)

"""

phc0, phc1 = ph

s0 = ps(data, p0=phc0, p1=phc1)
Expand Down