jjhelmus / nmrglue

A module for working with NMR data in Python
BSD 3-Clause "New" or "Revised" License
211 stars 86 forks source link

Update proc_autophase.py #48

Open atomman opened 8 years ago

atomman commented 8 years ago

The Automatic phase correction did not work as intended. Especially the penalty function was calculated incorrectly. I think. Also there should be an option to set the derivative order and a penalty function scaling factor, gamma.

I would suggest that the function also returned the optimum phases. In the case where 100's of spectra should be processed it would be nice to be able to store the phases in case one wants to reprocess.

jjhelmus commented 8 years ago

This update seems to have broken the 'peak_minima' algorithm and neither algorithm is giving good results for some toy data I tried. @atomman can you provide an example of where this autops function works well?

atomman commented 8 years ago

OK. I see the problem. On the master of my branch(examples/processing/ACME test.ipynb) I have uploaded a notebook showing the new implementation compared to the current one. In order access the acme and peak_minima through autops one must pass the new gamma variable as a dummy variable. I dont like this hack and suggest that we create separate functions for acme and peak_minima auto phase functionality.

PS. Should we output a reference to the acme paper when the algorithm is used to remind the user of the citation? Or is it OK just to have it in the docs?

jjhelmus commented 8 years ago

PS. Should we output a reference to the acme paper when the algorithm is used to remind the user of the citation? Or is it OK just to have it in the docs?

I think keeping it it the docs only is acceptable.

atomman commented 8 years ago

Why did the CI test fail?

On aside note: What is the recommended strategy for testing a branch with modified code? If I download and install, it overwrites my current nmrglue. How to rename the package so it installs as say "tmpnmrglue"?

jjhelmus commented 8 years ago

@atomman I marked the lines which are causing the CI to fail.

There are a number of methods for setting up a development environment which allows the modified version of nmrglue to be used as the default version. The method I uses is as follows:

The source code in the new directory can be imported directly if the directory is added to the Python search path. I accomplish by adding a ".pth" file to the appropriate site-packages directory on my system. This file is a text file containing the path to the directory containing the nmrglue source . Another options is to add the source directory to the PYTHONPATH environmental variable.

chatcannon commented 8 years ago

To avoid overwriting your current install, you can install in a virtualenv:

virtualenv nmrglue-dev
source nmrglue-dev/bin/activate
pip install numpy scipy
pip install --editable [path-to-nmrglue-git-checkout]
coveralls commented 7 years ago

Coverage Status

Changes Unknown when pulling bf249861f9d694b3cf7f2d08b9b22437387da693 on atomman:patch-3 into \ on jjhelmus:master**.

jjhelmus commented 7 years ago

@mfitzp Can change you could take a look at these changes as you were to creator of the proc_autophase module. I do not know enough about the algorithms used to provide a good review of if these changes are reasonable.

mfitzp commented 7 years ago

@jjhelmus Sure, will take a look over the weekend.

@atomman the changes are difficult to assess as you've 'broken' the existing interface. The way the current algorithms were set up was to provide all the algorithms through a single function 'autops', with the algorithm selectable via fn. By removing this you remove this possibility.

Can you refactor based on the original interface?

If there are changes you would like to make to the API they would probably be better submitted as another PR (although I of course defer to @jjhelmus on this)

atomman commented 7 years ago

@mfitzp i will try and make the changes over the weekend. The reason I suggested the change in the first place was to be able to provide a user defined gamma. This variable is only relevant for acme and I could not come up with an good common interface for multiple autophase functions relying on different inputs. Any suggestions?

On Fri, 16 Sep 2016, 23:14 Martin Fitzpatrick, notifications@github.com wrote:

@jjhelmus https://github.com/jjhelmus Sure, will take a look over the weekend.

@atomman https://github.com/atomman the changes are difficult to assess as you've 'broken' the existing interface. The way the current algorithms were set up was to provide all the algorithm through a single function 'autops', with the algorithm selectable via fn. By removing this you remove this possibility.

Can you refactor based on the original interface?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/jjhelmus/nmrglue/pull/48#issuecomment-247710667, or mute the thread https://github.com/notifications/unsubscribe-auth/AJwFa6PgXDkhBBAT3rMfm8CMhRoDkrGAks5qqwakgaJpZM4HpWkc .

Sent from my phone /Kresten

mfitzp commented 7 years ago

@atomman How about setting autops to take variable arguments and pass these into the acme function? e.g.

def autops(data, fn, p0=0.0, p1=0.0, *args):
    ...
    opt = scipy.optimize.fmin(fn, x0=opt, args=(data, )+args)

....
def _ps_acme_score(ph, data, gamma=5.0e-3):

This way any function can ave extra parameters sent to it (unfortunately fmin doesn't take keyword args which would be nicer).

What do you think?

The reason for the current API was to allow custom functions to be used easily. But there is probably an argument for method-specific functions (e.g. acmeps like you've done here) so they can be documented individually. But it would be nice to keep the general interface too.

coveralls commented 7 years ago

Coverage Status

Changes Unknown when pulling 2e2ae0e1c77456a54926a3d697310b7ac19fba65 on atomman:patch-3 into \ on jjhelmus:master**.

atomman commented 7 years ago

@mfitzp, I changed the functions back so they are compatible with theautops function. Now autops accepts keywords(gamma and order) via **kwargs which are passed to _ps_acme_score. They are also passed to _ps_peak_minima_score but are not used. This is the best solution I have found.

mfitzp commented 7 years ago

@atomman Thanks, but in the latest commit there is no autops function at all! :)

You shouldn't need to accept the arguments as kwargs in the _ps_acme_score function (they also aren't actually kwargs, but args! ...and to accept them into a dict you would need to use **kwargs). If the call from autops is structured as shown below:

opt = scipy.optimize.fmin(fn, x0=opt, args=(data, )+args)

...the individual arguments will be passed to _ps_peak_minima_score as normal function arguments. You should be able to pick them up as:

def _ps_acme_score(ph, data, gamma=5.0e-3, order=<something>):

The gamma, etc. above will accept the arguments as long as they are coming in the correct order. To accept and discard additional arguments you can add *args to the end, e.g. for _ps_peak_minima_score:

def _ps_peak_minima_score(ph, data, *args):

Sorry for the faff. If this isn't making any sense I can put together a complete example.

coveralls commented 7 years ago

Coverage Status

Changes Unknown when pulling e61f2ada5a286dd3f8f0d553fb07c62f9fbf3434 on atomman:patch-3 into \ on jjhelmus:master**.

atomman commented 7 years ago

@mfitzp Thank you for the tips. I think the current solution is quite nice. It would be nice if _ps_acme_score would default to an appropriate gamma and order parameter if they are not set by the user.

coveralls commented 7 years ago

Coverage Status

Changes Unknown when pulling 6f49c45fad4ac6fb452dca9b8d901815a471b019 on atomman:patch-3 into \ on jjhelmus:master**.

coveralls commented 7 years ago

Coverage Status

Changes Unknown when pulling ec4e92135b8d10ee8136b536caf7053f938fdd2d on atomman:patch-3 into \ on jjhelmus:master**.

jjhelmus commented 7 years ago

@atomman The changes made to the ACME algorithm seem to have broken the process. I'm using the following script as a test case:

import numpy as np
import matplotlib.pyplot as plt
import nmrglue as ng
from nmrglue.process.proc_autophase import autops

def complex_lorentzian(x, x_0, width, amp):
    return amp * (1. / (width + 1.j * (x-x_0)))

# simulate a spectrum with 2 lorentzians and noise
x = np.arange(-512, 512, dtype='complex64')
peak1 = complex_lorentzian(x, -200, 10, 1000)
peak2 = complex_lorentzian(x, 100, 15, 1500)
np.random.seed(0)
noise = np.random.uniform(size=(1024, ))
spectrum = peak1 + peak2 + noise

# apply linear phase to the spectrum to create 'raw' un-phased spectrum
p0 = 48.
p1 = 3.
raw = ng.proc_base.ps(spectrum, p0=p0, p1=p1, inv=True)

# phase raw spectrum using various techniques
manual = ng.proc_base.ps(raw, p0=p0, p1=p1)
_, acme = autops(raw, 'acme', order=1, gamma=0)
_, peak_min = autops(raw, 'peak_minima', order=2, gamma=2)

# plot spectrums
p1 = plt.plot(raw.real, 'r-')
p2 = plt.plot(manual.real, 'k-')
p3 = plt.plot(acme.real, 'b-')
p4 = plt.plot(peak_min.real, 'g-')
plt.legend(['Raw', 'Manual', 'ACME', 'Peak Minima'])
plt.show()

With these changes I am get the following results. Notice that that the blue line phased using the ACME method is not correct.

figure_1-1

Without these changes I get a reasonable phasing of the spectrum using the ACME method:

figure_1-2

Also, please correct the PEP8 style issues in the new code and add reasonable defaults for the order and gamma parameters.

atomman commented 7 years ago

Dear @jjhelmus. I will make the suggested corrections. The reason it does not work in your example is that gamma is set to zero. Try setting gamma to a value different from zero, ex. 5e-3, otherwise negative amplitudes in the spectrum will not be penalized and the correct phase will not be found.

On Wed, 5 Oct 2016, 18:37 Jonathan J. Helmus, notifications@github.com wrote:

@atomman https://github.com/atomman The changes made to the ACME algorithm seem to have broken the process. I'm using the following script as a test case:

import numpy as npimport matplotlib.pyplot as pltimport nmrglue as ngfrom nmrglue.process.proc_autophase import autops def complex_lorentzian(x, x_0, width, amp): return amp * (1. / (width + 1.j * (x-x_0)))

simulate a spectrum with 2 lorentzians and noise

x = np.arange(-512, 512, dtype='complex64') peak1 = complex_lorentzian(x, -200, 10, 1000) peak2 = complex_lorentzian(x, 100, 15, 1500) np.random.seed(0) noise = np.random.uniform(size=(1024, )) spectrum = peak1 + peak2 + noise

apply linear phase to the spectrum to create 'raw' un-phased spectrum

p0 = 48. p1 = 3. raw = ng.proc_base.ps(spectrum, p0=p0, p1=p1, inv=True)

phase raw spectrum using various techniques

manual = ng.procbase.ps(raw, p0=p0, p1=p1) , acme = autops(raw, 'acme', order=1, gamma=0) _, peak_min = autops(raw, 'peak_minima', order=2, gamma=2)

plot spectrums

p1 = plt.plot(raw.real, 'r-') p2 = plt.plot(manual.real, 'k-') p3 = plt.plot(acme.real, 'b-') p4 = plt.plot(peak_min.real, 'g-') plt.legend(['Raw', 'Manual', 'ACME', 'Peak Minima']) plt.show()

With these changes I am get the following results. Notice that that the blue line phased using the ACME method is not correct.

[image: figure_1-1] https://cloud.githubusercontent.com/assets/1050278/19122274/56a10e74-8aef-11e6-86ee-bbc350d1650b.png

Without these changes I get a reasonable phasing of the spectrum using the ACME method:

[image: figure_1-2] https://cloud.githubusercontent.com/assets/1050278/19122352/b335d66a-8aef-11e6-8246-7cf3e4436af8.png

Also, please correct the PEP8 https://www.python.org/dev/peps/pep-0008/ style issues in the new code and add reasonable defaults for the order and gamma parameters.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/jjhelmus/nmrglue/pull/48#issuecomment-251728542, or mute the thread https://github.com/notifications/unsubscribe-auth/AJwFa952kCtvRkC6QsiUajddsy-8pGUuks5qw9IqgaJpZM4HpWkc .

Sent from my phone /Kresten