jjhelmus / nmrglue

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

Matrix discrepancies between nmrpipe and nmrglue. #143

Closed andrealorenzon closed 3 years ago

andrealorenzon commented 3 years ago

I was trying to replicate a pipeline with the 2 tools and, step by step, compared results.

I've noticed a large discrepancy in matrix values for most functions, that I cannot explain with numeric issues. I am not an NMR expert or a physicist, so there could be some things I forgot to consider.

I am opening this issue to ask for any tip or assistance in sorting out this problem, or any information that can help.

The nmrpipe pipeline I used is:

nmrPipe -in test.fid \
| nmrPipe -fn SP -off 0.5 -end 0.95 -pow 2 -elb 0.0 -glb 0.0 -c 0.5 \
| nmrPipe -fn ZF -zf 1 -auto \
| nmrPipe -fn FT -verb \
| nmrPipe -fn PS -p0 -41.5 -p1 -8.0 -di \
  -out test.ft2 -ov

The nmrglue pipeline I used is:

dicFT1,  dataFT1  = ng.pipe.read("./transpose/test.fid")  # same origin fid file
dicPIPE, dataPIPE = ng.pipe.read("./transpose/testPS.ft2") # nmrpipe data, to be compared, see below for results
dicGLUE, dataGLUE = ng.pipe_proc.sp(dic, data,
                                off=0.5,
                                end=0.95,
                                pow=2,
                                c=0.5)
dicGLUE, dataGLUE = ng.pipe_proc.zf(dicGLUE, dataGLUE, 
                                zf=1, 
                                auto=True)

dicGLUE, dataGLUE = ng.pipe_proc.ft(dicGLUE, dataGLUE) 

dicGLUE, dataGLUE = ng.pipe_proc.ps(dicGLUE, dataGLUE, 
                                p0=-41.5, 
                                p1=-8.,)

Results follow: Note: the first column of the matrix is more or less similar, while the others are completely different.

image

kaustubhmote commented 3 years ago

@andrealorenzon Is this a dataset recorded on a Bruker machine and then converted to pipe format using NMRpipe's standard interface? If yes, then this issue is related to the process with which the digital filter is removed. At some point in the past (not sure exactly when), nmrpipe switched to removing the digital filter for Bruker datasets after processing by default, so your nmrpipe FID (test.fid) most likely still has the digital filter. When you process using NMRPipe, this is automatically detected and taken care of in the FT function, but in nmrglue, you will have to do this manually.

FYI, I can reproduce this issue with practically any Bruker dataset. If you click the button saying During conversion (Normal FID) in the nmrpipe conversion window, then this issue goes away. However, that is not an ideal solution.

Please drop a line if this explains the discrepancies for your dataset, and I can try and come up with a solution using the available information in dicGLUE. If this is not the case, it would be great if you can share how your first Fourier transformed FID looks like. It will help in debugging this.

andrealorenzon commented 3 years ago

I can confirm, these are Bruker data. Yes, that could explain the behavior.. at the first line of data import (nmrpipe vs nmrglue) the matrices are the same.

nmrPipe -in test.fid -out test2.fid == ng.pipe.read("test.fid")

I'll discuss your reply with the scientist using nmrpipe and follow up with any information. If it's ok for you, I can prepare and upload somewhere both sequences of transformed matrices, to be compared. May this help?

Meanwhile: to replicate, in nmrglue, the same pipeline, how can I manually take care of the digital filter? Thank you for your support!

EDIT: I tried my nmrglue pipeline both with filtered data as before, and without digital filter. Using np.allclose, I compared the outputs with or without filters in the 5 steps (reading, SP, ZF, FT, PS), with the following results (True if the values are comparable with the same nmrpipe process):

step    AMX DMX      tolerance =1e-10
0   True    True    OK
1   False   False   !!
2   False   False   !!
3   False   False   !!
4   False   False   !!
step    AMX DMX      tolerance =1e-08
0   True    True    OK
1   False   False   !!
2   False   False   !!
3   False   False   !!
4   False   False   !!
step    AMX DMX      tolerance =1e-06
0   True    True    OK
1   True    True    OK
2   True    True    OK
3   False   False   !!
4   False   False   !!
step    AMX DMX      tolerance =0.0001
0   True    True    OK
1   True    True    OK
2   True    True    OK
3   False   False   !!
4   False   False   !!
step    AMX DMX      tolerance =0.01
0   True    True    OK
1   True    True    OK
2   True    True    OK
3   False   False   !!
4   False   False   !!
step    AMX DMX      tolerance =1
0   True    True    OK
1   True    True    OK
2   True    True    OK
3   False   True    !!
4   False   False   !!
step    AMX DMX      tolerance =100
0   True    True    OK
1   True    True    OK
2   True    True    OK
3   True    True    OK
4   False   False   !!
step    AMX DMX      tolerance =10000
0   True    True    OK
1   True    True    OK
2   True    True    OK
3   True    True    OK
4   False   False   !!
kaustubhmote commented 3 years ago

Great. At this point, I can't seem to find how/where NMRPipe stores the value of GRPDLY. Even if I can find where nmrpipe stores this, the current nmrglue version does not seem to be able to read it.

There are 2 ways you can do this:

  1. Read in the raw data from Bruker into nmrglue and process it. However, you will have to use proc_base for all processing functions, and some more customization of each function will be required to exactly match the nmrpipe functions.

  2. Read in raw data from Bruker into nmrglue. Then read in the test.fid separately. Then you can use the read the value of GPDLY and use that to correct your data. (This should have worked with the ng.bruker.remove_digital_filter() function, but I am seeing some inconsistencies with this function and how pipe stores data).

I would recommend (1) if you want to just process data. Examples for this are given here. If you want to compare with nmrpipe, (2) is currently the only way I can think of doing this.

For (2), this works for me on a Bruker 1D dataset. By "works", I mean the arrays are virtually identical.

dicBRUK, dataBRUK = ng.bruker.read("../../") # raw bruker data
dicPIPE, dataPIPE = ng.pipe.read("./test.ft")  # nmrpipe data, to be compared, see below for results
dicFT1, dataFT1 = ng.pipe.read("../../test.fid")  # same origin fid file

dicGLUE, dataGLUE = ng.pipe_proc.sp(dicFT1, dataFT1, off=0.5, end=0.95, pow=2, c=0.5)
dicGLUE, dataGLUE = ng.pipe_proc.zf(dicGLUE, dataGLUE, zf=1, auto=True)
dicGLUE, dataGLUE = ng.pipe_proc.ft(dicGLUE, dataGLUE)
dicGLUE, dataGLUE = ng.pipe_proc.ps(dicGLUE, dataGLUE, p0=0, p1=-360*dicBRUK["acqus"]["GRPDLY"],)
dicGLUE, dataGLUE = ng.pipe_proc.ps(dicGLUE, dataGLUE, p0=-60.5, p1=0)

Note that removing the digital filter is equivalent to applying a first-order phase correction that depends on the GRPDLY parameter. This is discussed in #67.

kaustubhmote commented 3 years ago

@andrealorenzon If I am reading the AMX and DMX numbers correctly, it seems that there is no difference between these two modes? I am not quite sure why this is. At the very least, in the DMX mode, the arrays should be quite different if the filter is not removed/corrected. Also not quite sure why phase correction should make the arrays different. Most likely I am missing something in your post. If you still have trouble with this, you can mail me the matrices as NumPy arrays, and I will be happy to look at them.

andrealorenzon commented 3 years ago

Yes, I assumed them to be very different, too... something's wrong somewhere.. Probably I lack some topic knowledge to explain properly, and I thank you for your extreme patience and support, since I'm a newbie in NMR processing.

I've created a repository with a deeper explanation of my actions, code, and folders with raw and processed data for each step. I pushed everything at https://github.com/andrealorenzon/nmr-pipelines.

kaustubhmote commented 3 years ago

@andrealorenzon No problem at all. I just opened an issue here. As far as I can see, the nmrglue and nmrpipe seem quite equivalent. Let me know if you see some other issues.

andrealorenzon commented 3 years ago

Great, thank you. Now I just need to know how to find that "76" procedurally and I should be done.

EDIT: done, thanks to @jjhelmus . I'll close this issue, thank you everyone!