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
Open

Update proc_autophase.py #48

wants to merge 12 commits into from

Conversation

atomman
Copy link
Contributor

@atomman atomman commented Mar 4, 2016

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.

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.
Needed to take the absolute of the n-derivative of the spectrum
@@ -68,35 +68,42 @@ def _ps_acme_score(ph, data):
Value of the objective function (phase score)

"""
stepsize = 1

"""
Copy link
Contributor

Choose a reason for hiding this comment

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

I think this is why the Travis CI tests are failing.

Edited: there are two many parentheses for the doctoring. One of the parenthesis is not indented correctly.

I have made changes to the autops function to accept the argument for the gamma constant.
Also I have set the options to return the phases in a tuple.

I'm unsure about how to describe these changes in the doc?
@jjhelmus
Copy link
Owner

jjhelmus commented Mar 7, 2016

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
Copy link
Contributor Author

atomman commented Mar 9, 2016

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
Copy link
Owner

jjhelmus commented Mar 9, 2016

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.

Replaced autops with functions psacme and pspmin.
As the two methods relies on different input the should have their own functions I think.
@atomman
Copy link
Contributor Author

atomman commented Mar 12, 2016

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"?

@atomman atomman closed this Mar 12, 2016
@atomman atomman reopened this Mar 12, 2016
@@ -14,46 +14,49 @@
from .proc_base import ps


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

Choose a reason for hiding this comment

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

This syntax is not allowed.

@jjhelmus
Copy link
Owner

@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:

  • Remove any previous installations of nmrglue.
  • Fork the nmrglue repository on GitHub using the top right "Fork" button".
  • In a development directory use git clone address_of_your_fork to clone your fork of nmrglue.

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
Copy link
Contributor

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
Copy link

coveralls commented Sep 16, 2016

Coverage Status

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

@@ -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.

@jjhelmus
Copy link
Owner

@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
Copy link
Contributor

mfitzp commented Sep 16, 2016

@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
Copy link
Contributor Author

atomman commented Sep 17, 2016

@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, [email protected]
wrote:

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

@atomman https:/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
#48 (comment), or mute
the thread
https:/notifications/unsubscribe-auth/AJwFa6PgXDkhBBAT3rMfm8CMhRoDkrGAks5qqwakgaJpZM4HpWkc
.

Sent from my phone
/Kresten

@mfitzp
Copy link
Contributor

mfitzp commented Sep 17, 2016

@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.

@mfitzp, I changed the functions back so they are compatible with the 'autops' 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.
@coveralls
Copy link

coveralls commented Sep 20, 2016

Coverage Status

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

@atomman
Copy link
Contributor Author

atomman commented Sep 20, 2016

@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
Copy link
Contributor

mfitzp commented Sep 20, 2016

@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.

forgot autops function in last commit
@coveralls
Copy link

coveralls commented Sep 20, 2016

Coverage Status

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

@atomman
Copy link
Contributor Author

atomman commented Sep 20, 2016

@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.

Added a few white spaces and changed a np.sum to np.nansum
@coveralls
Copy link

coveralls commented Sep 20, 2016

Coverage Status

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


return phasedspc
phasedspc = ng.proc_base.ps(data, p0=opt[0], p1=opt[1])
Copy link
Owner

Choose a reason for hiding this comment

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

ng.proc_base.ps is not valid here. This should be ps


def _ps_acme_score(ph, data):
def _ps_acme_score(ph, data, kwargs):
Copy link
Owner

Choose a reason for hiding this comment

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

**kwargs?

phc0, phc1 = ph
s = ng.proc_base.ps(data, p0=phc0, p1=phc1)
Copy link
Owner

Choose a reason for hiding this comment

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

ps

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I will fix this



def _ps_peak_minima_score(ph, data):
def _ps_peak_minima_score(ph, data, kwargs):
Copy link
Owner

Choose a reason for hiding this comment

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

**kwargs

Copy link
Contributor Author

Choose a reason for hiding this comment

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

@jjhelmus I'm a little unsure how to do this. The main functions 'autops' accepts **kwargs. Then these are passed to _ps_acme_score as kwargs. The _ps_acme_score needs additional parameters; here implemented as kwargs['gamma'] and kwargs['order']

@@ -122,7 +123,7 @@ def _ps_peak_minima_score(ph, data):

phc0, phc1 = ph

s0 = ps(data, p0=phc0, p1=phc1)
s0 = ng.proc_base.ps(data, p0=phc0, p1=phc1)
Copy link
Owner

Choose a reason for hiding this comment

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

Just ps

@coveralls
Copy link

coveralls commented Sep 23, 2016

Coverage Status

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

@jjhelmus
Copy link
Owner

jjhelmus commented Oct 5, 2016

@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
Copy link
Contributor Author

atomman commented Oct 5, 2016

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, [email protected]
wrote:

@atomman https:/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.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.

[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
#48 (comment), or mute
the thread
https:/notifications/unsubscribe-auth/AJwFa952kCtvRkC6QsiUajddsy-8pGUuks5qw9IqgaJpZM4HpWkc
.

Sent from my phone
/Kresten

p = 1000 * pfun

return h1s + p
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)

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.

6 participants