dkriegner / xrayutilities

xrayutilities - a package with useful scripts for X-ray diffraction physicists
http://xrayutilities.sourceforge.io
GNU General Public License v2.0
83 stars 29 forks source link

Convolutions for powder diffraction #32

Closed mendenmh closed 7 years ago

mendenmh commented 7 years ago

I have just recently discovered xrayutilities, and find it an exceedingly nice package for many purposes.

I do extensive powder diffraction work, and would like to integrate some extra features which would be broadly useful to the powder community. I have developed a pure-python code which implements the Fundamental Parameters Approach to compute peak shapes directly from basic machine and sample parameters. The paper on it is at: Journal of Research of NIST, and the code is provided as ancillary material, accessible at FPA Python Code. I think it would take very little effort for someone to wrap this code so that it fits in under the Powder experiment class as an alternative to the Gaussian convolver. This code includes essentially all currently documented aberrations to the peak shape from instrument geometry, crystal microstructure, and sample properties such as absorption. It is also extensible via simple plugins to handle extra aberrations.

I don't know what the development schedule for new bits of xrayutilities is. I probably do not have time to fully implement the needed wrapper for this, but would be glad to assist with the process.

dkriegner commented 7 years ago

generally every contribution to xrayutilities is welcome!

In particular the code you mention would fit very well to xrayutilities since the current Gaussian convolver is obviously totally oversimplified and most probably almost not used. I think your accurate models would be a valuable addition to this project.

I assume that writing some wrapper would not be too much effort, but will need some useful example scripts and a reasonable amount of documentation. The question is if you are willing to contribute your code in a way compatible with the GPL license? I think it would be most useful to include your code in the project and make it python 3 compatible so that it can be installed along with the current package.

Since we do not receive many "external" code contribution there is no strict procedure to follow in such a case. Feel free to contact me by email or any other means.

mendenmh commented 7 years ago

On Dec 15, 2016, at 3:47 PM, Dominik Kriegner notifications@github.com<mailto:notifications@github.com> wrote:

generally every contribution to xrayutilities is welcome!

In particular the code you mention would fit very well to xrayutilities since the current Gaussian convolver is obviously totally oversimplified and most probably almost not used. I think your accurate models would be a valuable addition to this project.

I assume that writing some wrapper would not be too much effort, but will need some useful example scripts and a reasonable amount of documentation. The question is if you are willing to contribute your code in a way compatible with the GPL license? I think it would be most useful to include your code in the project and make it python 3 compatible so that it can be installed

along with the current package.

I don't know if you use doxygen, but this code has intricate doxygen comments in it, so it generates a manual on itself that way. The function at the bottom, fourier_line_profile(), creates an instance of the class and populates it, so it shows moderately well how to use the class. Also, the if name=="main" section shows some usage. I will, of course, be glad to help with the wrapper, I just can't lead it. The class has features built into it to make it efficient, and I can help with how to take advantage of these, for example. I don't think the wrapper will be more than 100 lines of code.

The license is easy... this is written under U.S. Government support, and is in the public domain. I don't think there should be any big issues converting it to Python 3; it doesn't really do anything that uses deep python magic, so the 2to3 converter might just work. I haven't tried.

Since we do not receive many "external" code contribution there is no strict procedure to follow in such a case. Feel free to contact me by email or any other means.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHubhttps://github.com/dkriegner/xrayutilities/issues/32#issuecomment-267439575, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AXc4rCVEggxpfkmQMemQNhdGtFy9xzsMks5rIaddgaJpZM4LONwA.

dkriegner commented 7 years ago

This sounds good, I already had a look on some parts of the code yesterday. Indeed the conversion to python3 will not be complicated. For documentation we use sphinx but also here I think we can convert the existing documentation quite easy.

Once I have some free time I will look into this and ask you for your comments!

mendenmh commented 7 years ago

Thinking a bit in advance, I suspect we want to make a new Experiment subclass, for a diverging-beam powder diffractometer, into which all the machine parameters can be embedded, and maybe a powder subclass of Material into which the strain broadening and other things can be embedded. Then, these parameters would be trivially harvested by the convolver to make the peak shape.

On Dec 16, 2016, at 8:00 AM, Dominik Kriegner notifications@github.com<mailto:notifications@github.com> wrote:

This sounds good, I already had a look on some parts of the code yesterday. Indeed the conversion to python3 will not be complicated. For documentation we use sphinx but also here I think we can convert the existing documentation quite easy.

Once I have some free time I will look into this and ask you for your comments!

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHubhttps://github.com/dkriegner/xrayutilities/issues/32#issuecomment-267589245, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AXc4rCLEyoyOfsTVtL56yimJ2vz_ihHtks5rIotrgaJpZM4LONwA.

mendenmh commented 7 years ago

A question: are planning to move xrayutilities to python 3, or to maintain a full python 2 capability, also. I have not moved any code to python 3 yet, since a lot of the science packages are still not at 100%. I would like to start with the effort for this convolver in the python2 branch of xyayutilities, where I can give more help, and then handle the python 3 conversion as a second step.

On Dec 16, 2016, at 8:00 AM, Dominik Kriegner notifications@github.com<mailto:notifications@github.com> wrote:

This sounds good, I already had a look on some parts of the code yesterday. Indeed the conversion to python3 will not be complicated. For documentation we use sphinx but also here I think we can convert the existing documentation quite easy.

Once I have some free time I will look into this and ask you for your comments!

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHubhttps://github.com/dkriegner/xrayutilities/issues/32#issuecomment-267589245, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AXc4rCLEyoyOfsTVtL56yimJ2vz_ihHtks5rIotrgaJpZM4LONwA.

dkriegner commented 7 years ago

A question: are planning to move xrayutilities to python 3, or to maintain a full python 2 capability, also. I have not moved any code to python 3 yet, since a lot of the science packages are still not at 100%. I would like to start with the effort for this convolver in the python2 branch of xyayutilities, where I can give more help, and then handle the python 3 conversion as a second step.

So far all the code runs in Python-2.7.X and >3.2 simultaneously. I agree that its important to keep it like that for the near future!

Thinking a bit in advance, I suspect we want to make a new Experiment subclass, for a diverging-beam powder diffractometer, into which all the machine parameters can be embedded, and maybe a powder subclass of Material into which the strain broadening and other things can be embedded. Then, these parameters would be trivially harvested by the convolver to make the peak shape.

Yes I think a new Experiment subclass will be the best scenario. The machine parameters could be either embedded there or parsed from the config file. I would like the idea of the config file since it would allow to override the defaults via a local config file which one can place in the particular project directory. In addition to that of course changing parameters in the call to the Convolve function or Class constructor will be needed.

Regarding the Material class I will think a bit. But I imagine a dedicated subclass of SMaterial could do the job.

dkriegner commented 7 years ago

License issue seems indeed to be easy: see https://www.gnu.org/licenses/gpl-faq.en.html#GPLUSGovAdd

dkriegner commented 7 years ago

@mendenmh Please have a look on my recent commit 8dfba316e33d354bf58734e39a2e4305abff03fd

I transformed your code to comply with the rest of xrayutilities which tries to use pep8 source code formatting and wrote a basic wrapper which is now included in the Powder.Convolve() function

Basically the use of the code would be

import xrayutilities as xu
import numpy
p = xu.Powder(xu.materials.Si)
tt = numpy.arange(0,140,0.01)
inte = p.Convolve(tt)

The default parameters are obtained from the config file and can be also set by supplying a dictionary parameter of the class constructor. The used FP_Profile class can also be specified. Please let me know if you think this can be flexible enough to allow extensions.

There are few other things to be checked:

  1. Are you ok with the mentioning of the original license agreement and the copyright notice in general?

  2. I do not like the convolver name si_psd since not every linear detector is silicon based. Would you agree to rename this to only psd or linear_detector?

  3. The code now runs in both python2 and 3. In my tests I got the same results with both versions but possibly there could be integer division errors, since the default behavior for such divisions changed. Do you think this is an issue in your code?

  4. To obtain the full powder pattern of an material I used interpolation from the results obtained from the convolution. Do you think there is a better option for that?

  5. Is there any normalization of the output intensities necessary? In your example code you use

        #normalize both to integral=sum/dx
        dx1=(result.twotheta_deg[1]-result.twotheta_deg[0])
        peaknorm=40*result.peak/dx1
  6. Should this (and corresponding) debug code stay in the file?

    # used for debugging moments from FP_profile.axial_helper().
    moment_list = []
    # if this is *True*, compute and save moment errors
    collect_moment_errors = False
  7. in the axial convolver: slit_length_source can not be equal to slit_length_target. Is the following an expected behavior?

In [1]: import xrayutilities as xu; p = xu.Powder(xu.materials.Si)

In [2]: fpsettings = {'axial': {'slit_length_source': 10e-3, 'slit_length_target': 10e-3}}

In [3]: p = xu.Powder(xu.materials.Si, fpsettings=fpsettings)

In [4]: tt = arange(0,140,0.01)

In [5]: p.Convolve(tt)
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-5-d007978b440a> in <module>()
----> 1 p.Convolve(tt)

/home/dk/.local/lib64/python2.7/site-packages/xrayutilities-1.3.3-py2.7-linux-x86_64.egg/xrayutilities/powder.pyc in Convolve(self, twotheta, window_width)
   1909             # set parameters which are shared by many things
   1910             p.set_parameters(twotheta0_deg=twotheta_x)
-> 1911             result = p.compute_line_profile()
   1912 
   1913             outint += numpy.interp(tt, result.twotheta_deg,

/home/dk/.local/lib64/python2.7/site-packages/xrayutilities-1.3.3-py2.7-linux-x86_64.egg/xrayutilities/powder.pyc in compute_line_profile(self, convolver_names, compute_derivative)
   1468         # run through the name list, and call the convolver to harvest its
   1469         # result
-> 1470         conv_list = [self.convolver_funcs[x]() for x in convolver_names]
   1471 
   1472         # now, multiply everything together

/home/dk/.local/lib64/python2.7/site-packages/xrayutilities-1.3.3-py2.7-linux-x86_64.egg/xrayutilities/powder.pyc in conv_axial(self)
    985                 sollerDdeg=xx.angD_deg,
    986                 R=xx.diffractometer_radius,
--> 987                 twotheta=xx.twotheta0
    988             )
    989         axfn[:] = best_rfft(axbuf)

/home/dk/.local/lib64/python2.7/site-packages/xrayutilities-1.3.3-py2.7-linux-x86_64.egg/xrayutilities/powder.pyc in full_axdiv_I3(self, Lx, Ls, Lr, R, twotheta, epsvals, sollerIdeg, sollerDdeg, nsteps, axDiv)
    907             eps, idxmin, idxmax, I2p, I2m = self.full_axdiv_I2(
    908                 Lx=Lx, Lr=Lr, Ls=Ls, beta=beta, R=R,
--> 909                 twotheta=twotheta, epsvals=epsvals)
    910 
    911             # after eq. 26 in Ch&Co

/home/dk/.local/lib64/python2.7/site-packages/xrayutilities-1.3.3-py2.7-linux-x86_64.egg/xrayutilities/powder.pyc in full_axdiv_I2(self, Lx, Ls, Lr, R, twotheta, beta, epsvals)
    840             indices += F2(dst=I2m, lower=ed, upper=ec, eea=ed)
    841         elif rng == 3:
--> 842             indices += F4(dst=I2m, lower=eb, upper=ea, eea=ea)
    843             indices += F1(dst=I2m, lower=ec, upper=eb, eea=ec, eeb=ed)
    844             indices += F2(dst=I2m, lower=ed, upper=ec, eea=ed)

/home/dk/.local/lib64/python2.7/site-packages/xrayutilities-1.3.3-py2.7-linux-x86_64.egg/xrayutilities/powder.pyc in F4(dst, lower, upper, eea)
    821                                      innerbound=upper, outerbound=lower,
    822                                      epsvals=epsvals, peakpos=eps0,
--> 823                                      k=-sqrt(abs(eps0 - eea)), y0=+1)
    824 
    825         I2p = self._I2p

/home/dk/.local/lib64/python2.7/site-packages/xrayutilities-1.3.3-py2.7-linux-x86_64.egg/xrayutilities/powder.pyc in axial_helper(self, outerbound, innerbound, epsvals, destination, peakpos, y0, k)
    644             deps[:] = peakpos
    645             deps -= epsvals[idx0:idx1]
--> 646             deps[-1] = peakpos - min(innerbound, peakpos)
    647             deps[0] = peakpos - outerbound
    648         else:

IndexError: index -1 is out of bounds for axis 0 with size 0
  1. most important: I think that there should be some documentation about the possible parameters of each convolver in the code. This seems to be missing. Do you think one line of description for every parameter will help or is the user expected to read your paper?

I am aware that few things (especially the sample related parameters) can so far not be set in a useful way. As I said above I think this should be solved by another dedicated class which works with SMaterial or lists of such materials which can also contain the sample parameters.

mendenmh commented 7 years ago

Wow, this is fantastic that you got this going so quickly!

I will look over it in detail.

The documentation really should be the paper. There is just no way to do Fundamental Parameters without having read the papers on it, since it is fairly intricate. Since the parameters are dependent on the actual physical setup and material the user has, there is no sensible concept of 'recommended' values.

answers to specific numbered questions:

1: any 'open' license wording is fine, but the original attribution should be maintained.

  1. the name si_psd is quite specific; the other kind of psd (gas) is not handled by this. I do not know of any solid-state psd devices that are not silicon, although they could exist. We could name it solid_state_psd or strip_detector_psd or something like that.
  2. I think I was very careful to use // for integer division, and to always cast floats where needed in single-/ division. I've been watching this since the first creation of the // operator. I don't think there should be any problems on this front.
  3. Interpolation is probably the best bet. I would also suggest a significant oversampling factor (note that in the fourier_line_profile function, the value is set to 10; this is bigger than usually necessary, but 4 isn't crazy). This results in line profiles which are more accurate at the sub-pixel level. For very high precision work, (lattice constant measurement, typically), this makes a difference.
  4. The intensities should be normalized in a way that is consistent as angles are varied. I think that multiplying these lines by the appropriate scattering strength, Lorentz and polarization factor, should result in peak intensities that can be properly compared. Note that I have never tested this in detail, since I have only used this for Pawley-type fits, in which the intensity of each line is independent. Now that it is being integrated into xrayutilities, it will be easy for me to test it against full Rietveld-type fits, since you have the scattering factors.
  5. The debug code doesn't hurt; I think . would leave it in.
  6. The problem with the two lengths being identical is a known bug; I just set them a tiny bit apart if they really happen to be the same. Someday I will get around to making the fix for it.

On Jan 2, 2017, at 8:42 AM, Dominik Kriegner notifications@github.com<mailto:notifications@github.com> wrote:

@mendenmhhttps://github.com/mendenmh Please have a look on my recent commit 8dfba31https://github.com/dkriegner/xrayutilities/commit/8dfba316e33d354bf58734e39a2e4305abff03fd

I transformed your code to comply with the rest of xrayutilities which tries to use pep8 source code formatting and wrote a basic wrapper which is now included in the Powder.Convolve() functionhttps://github.com/dkriegner/xrayutilities/blob/fpamodel/xrayutilities/powder.py#L1870

Basically the use of the code would be

import xrayutilities as xu import numpy p = xu.Powder(xu.materials.Si) tt = numpy.arange(0,140,0.01) inte = p.Convolve(tt)

The default parameters are obtained from the config file and can be also set by supplying a dictionary parameter of the class constructor. The used FP_Profile class can also be specified. Please let me know if you think this can be flexible enough to allow extensions.

There are few other things to be checked:

  1. Are you ok with the mentioning of the original license agreement and the copyright notice in general?

  2. I do not like the convolver name si_psd since not every linear detector is silicon based. Would you agree to rename this to only psd or linear_detector?

  3. The code now runs in both python2 and 3. In my tests I got the same results with both versions but possibly there could be integer division errors, since the default behavior for such divisions changedhttps://www.python.org/dev/peps/pep-0238/. Do you think this is an issue in your code?

  4. To obtain the full powder pattern of an material I used interpolation from the results obtained from the convolution. Do you think there is a better option for that?

  5. Is there any normalization of the output intensities necessary? In your example code you use

    #normalize both to integral=sum/dx
    dx1=(result.twotheta_deg[1]-result.twotheta_deg[0])
    peaknorm=40*result.peak/dx1
  6. Should this (and corresponding) debug code stay in the file?

used for debugging moments from FP_profile.axial_helper().

moment_list = []

if this is True, compute and save moment errors

collect_moment_errors = False

  1. in the axial convolver: slit_length_source can not be equal to slit_length_target. Is the following an expected behavior?

In [1]: import xrayutilities as xu; p = xu.Powder(xu.materials.Si)

In [2]: fpsettings = {'axial': {'slit_length_source': 10e-3, 'slit_length_target': 10e-3}}

In [3]: p = xu.Powder(xu.materials.Si, fpsettings=fpsettings)

In [4]: tt = arange(0,140,0.01)

In [5]: p.Convolve(tt)

IndexError Traceback (most recent call last)

in () ----> 1 p.Convolve(tt) /home/dk/.local/lib64/python2.7/site-packages/xrayutilities-1.3.3-py2.7-linux-x86_64.egg/xrayutilities/powder.pyc in Convolve(self, twotheta, window_width) 1909 # set parameters which are shared by many things 1910 p.set_parameters(twotheta0_deg=twotheta_x) -> 1911 result = p.compute_line_profile() 1912 1913 outint += numpy.interp(tt, result.twotheta_deg, /home/dk/.local/lib64/python2.7/site-packages/xrayutilities-1.3.3-py2.7-linux-x86_64.egg/xrayutilities/powder.pyc in compute_line_profile(self, convolver_names, compute_derivative) 1468 # run through the name list, and call the convolver to harvest its 1469 # result -> 1470 conv_list = [self.convolver_funcs[x]() for x in convolver_names] 1471 1472 # now, multiply everything together /home/dk/.local/lib64/python2.7/site-packages/xrayutilities-1.3.3-py2.7-linux-x86_64.egg/xrayutilities/powder.pyc in conv_axial(self) 985 sollerDdeg=xx.angD_deg, 986 R=xx.diffractometer_radius, --> 987 twotheta=xx.twotheta0 988 ) 989 axfn[:] = best_rfft(axbuf) /home/dk/.local/lib64/python2.7/site-packages/xrayutilities-1.3.3-py2.7-linux-x86_64.egg/xrayutilities/powder.pyc in full_axdiv_I3(self, Lx, Ls, Lr, R, twotheta, epsvals, sollerIdeg, sollerDdeg, nsteps, axDiv) 907 eps, idxmin, idxmax, I2p, I2m = self.full_axdiv_I2( 908 Lx=Lx, Lr=Lr, Ls=Ls, beta=beta, R=R, --> 909 twotheta=twotheta, epsvals=epsvals) 910 911 # after eq. 26 in Ch&Co /home/dk/.local/lib64/python2.7/site-packages/xrayutilities-1.3.3-py2.7-linux-x86_64.egg/xrayutilities/powder.pyc in full_axdiv_I2(self, Lx, Ls, Lr, R, twotheta, beta, epsvals) 840 indices += F2(dst=I2m, lower=ed, upper=ec, eea=ed) 841 elif rng == 3: --> 842 indices += F4(dst=I2m, lower=eb, upper=ea, eea=ea) 843 indices += F1(dst=I2m, lower=ec, upper=eb, eea=ec, eeb=ed) 844 indices += F2(dst=I2m, lower=ed, upper=ec, eea=ed) /home/dk/.local/lib64/python2.7/site-packages/xrayutilities-1.3.3-py2.7-linux-x86_64.egg/xrayutilities/powder.pyc in F4(dst, lower, upper, eea) 821 innerbound=upper, outerbound=lower, 822 epsvals=epsvals, peakpos=eps0, --> 823 k=-sqrt(abs(eps0 - eea)), y0=+1) 824 825 I2p = self._I2p /home/dk/.local/lib64/python2.7/site-packages/xrayutilities-1.3.3-py2.7-linux-x86_64.egg/xrayutilities/powder.pyc in axial_helper(self, outerbound, innerbound, epsvals, destination, peakpos, y0, k) 644 deps[:] = peakpos 645 deps -= epsvals[idx0:idx1] --> 646 deps[-1] = peakpos - min(innerbound, peakpos) 647 deps[0] = peakpos - outerbound 648 else: IndexError: index -1 is out of bounds for axis 0 with size 0 1. most important: I think that there should be some documentation about the possible parameters of each convolver in the code. This seems to be missing. Do you think one line of description for every parameter will help or is the user expected to read your paper? I am aware that few things (especially the sample related parameters) can so far not be set in a useful way. As I said above I think this should be solved by another dedicated class which works with SMaterial or lists of such materials which can also contain the sample parameters. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.
dkriegner commented 7 years ago

Thanks for your reply, most points seem solved by your answers. According to your explaination I kept "si_psd". Regarding point 5 I will do some tests myself later. I think this will be interesting but will also reveal some other weaknesses of my code (like issue #7 ) and that for a more speedy calculation the database needs to be speed up (issue #33 )

I am now working on another bigger restructuring of the code to enable a sensible way to set the sample parameters. I hope to commit/push later today.

Meanwhile I discovered the there are some unused parameters in fourier_line_profile: these are:

mendenmh commented 7 years ago

On Jan 4, 2017, at 8:31 AM, Dominik Kriegner notifications@github.com<mailto:notifications@github.com> wrote:

Thanks for your reply, most points seem solved by your answers. According to your explaination I kept "si_psd". Regarding point 5 I will do some tests myself later. I think this will be interesting but will also reveal some other weaknesses of my code (like issue #7https://github.com/dkriegner/xrayutilities/issues/7 ) and that for a more speedy calculation the database needs to be speed up (issue #33https://github.com/dkriegner/xrayutilities/issues/33 )

I am now working on another bigger restructuring of the code to enable a sensible way to set the sample parameters. I hope to commit/push later today.

Meanwhile I discovered the there are some unused parameters in fourier_line_profile: these are:

the function fourier_line_profile is a convenience function, and originally was written to have some aberrations that were never implemented. I really only use it in testing code, not in production. I would ignore it, or delete it from your production code. Using the full class enables easy parameter setting/getting, and also causes a lot of cacheing, so that if one parameter is changed, very little gets recalculated. This results in its usability in fitting procedures, in which parameters get shifted to compute numerical derivatives.

The error when the two lengths are identical is just a lack of an if statement somewhere that the Fundamental Parameters methods end up cutting off a boundary wrong. Maybe I should fix that today.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://github.com/dkriegner/xrayutilities/issues/32#issuecomment-270370505, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AXc4rC5BQF6oeXmpR6jnlJY9M3EKCgU_ks5rO58hgaJpZM4LONwA.

mendenmh commented 7 years ago

I have changed two places in the code here, to avoid singular conditions, I think. One was the problems with the lengths being the same. The other was that attempting to compute a peak at exactly 90 degrees (two theta) would crash, since there is a tangent in the axial divergence model. I now have (I think) put in conditions for both of them.\

You can diff this to the version I sent before to see the lines changed, and adapt your copy appropriately. Let me know if this resolves that issue.

Thanks.

On Jan 4, 2017, at 8:37 AM, Mendenhall, Marcus marcus.mendenhall@nist.gov<mailto:marcus.mendenhall@nist.gov> wrote:

On Jan 4, 2017, at 8:31 AM, Dominik Kriegner notifications@github.com<mailto:notifications@github.com> wrote:

Thanks for your reply, most points seem solved by your answers. According to your explaination I kept "si_psd". Regarding point 5 I will do some tests myself later. I think this will be interesting but will also reveal some other weaknesses of my code (like issue #7https://github.com/dkriegner/xrayutilities/issues/7 ) and that for a more speedy calculation the database needs to be speed up (issue #33https://github.com/dkriegner/xrayutilities/issues/33 )

I am now working on another bigger restructuring of the code to enable a sensible way to set the sample parameters. I hope to commit/push later today.

Meanwhile I discovered the there are some unused parameters in fourier_line_profile: these are:

the function fourier_line_profile is a convenience function, and originally was written to have some aberrations that were never implemented. I really only use it in testing code, not in production. I would ignore it, or delete it from your production code. Using the full class enables easy parameter setting/getting, and also causes a lot of cacheing, so that if one parameter is changed, very little gets recalculated. This results in its usability in fitting procedures, in which parameters get shifted to compute numerical derivatives.

The error when the two lengths are identical is just a lack of an if statement somewhere that the Fundamental Parameters methods end up cutting off a boundary wrong. Maybe I should fix that today.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://github.com/dkriegner/xrayutilities/issues/32#issuecomment-270370505, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AXc4rC5BQF6oeXmpR6jnlJY9M3EKCgU_ks5rO58hgaJpZM4LONwA.

dkriegner commented 7 years ago

I have changed two places in the code here, to avoid singular conditions, I think. One was the problems with the lengths being the same. The other was that attempting to compute a peak at exactly 90 degrees (two theta) would crash, since there is a tangent in the axial divergence model. I now have (I think) put in conditions for both of them.\

You can diff this to the version I sent before to see the lines changed, and adapt your copy appropriately. Let me know if this resolves that issue.

I am of course happy to include fixes in the code, I could however not find the new version of the code you mentioned.

In commit 8bf6898deff5902cec55a8dbc85ba58a492bd861 I have finished some bigger movement and began some documentation/example code. As you will see I have recreated the original PowderExperiment class (without the primitive convolver) and moved the PowderDiffraction class to the simpack subpackage where it fits better. PowderDiffraction (similar to the old Powder calculating the convolution for one material) and PowderModel (calculating the convolution for N materials) which are in simpack should be used when one attempts to use the fundamental parameter convolution. I also changed a bit the default values in the configuration and removed sample specific values from there. these are now supplied to the xrayutilities.simpack.Powder class which should be used to define materials. See here for a working example script.

A question I have at this point is: I attempted to convolute the signal from 0 to 180 degree as in

pd = xu.simpack.PowderDiffraction(xu.materials.LaB6)
tt = arange(0,180,0.01)
inte = pd.Convolve(tt)

When I do this the last peak at highest angle is missing after the convolution and I get some nan values around its position. Do you know what can cause this? Is this an edgecase at such high angles (~172deg)?

From my side there are only few more smaller problems which I will work on (performance and setting absorption from the material properties) but otherwise it seems to me its getting complete.

dkriegner commented 7 years ago

Now I just found out that with my implementation the convolution history is effectively not working because I call set_window before every calculation. Is it true that one should create an instance of FP_profile for every peak in the powder pattern? Like this one can avoid to call set_window unless the lattice parameters/peak position is changed by the user/fit function.

mendenmh commented 7 years ago

On Jan 4, 2017, at 6:10 PM, Dominik Kriegner notifications@github.com<mailto:notifications@github.com> wrote:

Now I just found out that with my implementation the convolution history is effectively not working because I call set_window before every calculation. Is it true that one should create an instance of FP_profile for every peak in the powder pattern? Like this one can avoid to call set_window unless the lattice parameters/peak position is changed by the user/fit function.

Yes, that is how to use it if you want optimum cache usage.

Once you have this working generally, I will send along a little module which uses the python multiprocessing module to manage peak computation in parallel. I think it is best, though, to first get it working in single-threaded mode.

I have re-attached the python code. If it disappears again, I will find a different way to send it. It is possible one of the mailers doesn't like code attachments.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://github.com/dkriegner/xrayutilities/issues/32#issuecomment-270514264, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AXc4rEicCV32Och_ToDLNLy-50I4zEu1ks5rPCbNgaJpZM4LONwA.

dkriegner commented 7 years ago

Your attachment did again not make it! I guess it does not work to attach files by email. you either have to do it via the webinterface (where it should work) or drop me a classic email.

Good news from my side is that as of commit d2d5f01e2e17a8963807b6ddaa0fa644b2007608 the buffering works in the Convolve function. The Calculate function on the other hand avoids the cache and recalculates everything. A more clever mechanism to decide if this is really necessary needs to be thought of but this needs to wait until more clever ways to change/update material parameters are developed.

Btw.: I know it will not be a sufficient description but could you suggest a line of text for the parameters in these lines L321 L322

dkriegner commented 7 years ago

@mendenmh I have now included an example script for very basic refinement. Unfortunately some important parameters are currently missing (lattice parameters and other material parameters like atom positions, occupation, Debye Waller factors). Also it runs very slow since it currently takes no benefit from caching. I need to redesign the Crystal class and include some sort of checksum to allow to detect when the peak positions changed and only then use FP_profile.set_window(). For the fitting I use lmfit since I think it offers lots of features which will be useful for this sort of refinement.

But I would appreciate your feedback on the implementation in its current state and if possible suggest a phrasing for these lines L321 L322

mendenmh commented 7 years ago

Generally, the window should start out wide enough to cover routine, small shifts in peak positions which result from refinement. There is no need to have the window exactly track the peak. I will look at the example later today.

Marcus

On Jan 13, 2017, at 3:20 AM, Dominik Kriegner notifications@github.com<mailto:notifications@github.com> wrote:

@mendenmhhttps://github.com/mendenmh I have now included an example scripthttps://github.com/dkriegner/xrayutilities/blob/fpamodel/examples/simpack_powdermodel.py for very basic refinement. Unfortunately some important parameters are currently missing (lattice parameters and other material parameters like atom positions, occupation, Debye Waller factors). Also it runs very slow since it currently takes no benefit from caching. I need to redesign the Crystalhttps://github.com/dkriegner/xrayutilities/blob/master/xrayutilities/materials/material.py#L463 class and include some sort of checksum to allow to detect when the peak positions changed and only then use FP_profile.set_window(). For the fitting I use lmfithttps://lmfit.github.io/lmfit-py/ since I think it offers lots of features which will be useful for this sort of refinement.

But I would appreciate your feedback on the implementation in its current state and if possible suggest a phrasing for these lines L321https://github.com/dkriegner/xrayutilities/blob/fpamodel/xrayutilities/simpack/smaterials.py#L321 L322https://github.com/dkriegner/xrayutilities/blob/fpamodel/xrayutilities/simpack/smaterials.py#L322

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://github.com/dkriegner/xrayutilities/issues/32#issuecomment-272386465, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AXc4rHrc6syroEtFH9NRcekNtJr1LslZks5rRzPLgaJpZM4LONwA.

mendenmh commented 7 years ago

On Jan 13, 2017, at 3:20 AM, Dominik Kriegner notifications@github.com<mailto:notifications@github.com> wrote:

@mendenmhhttps://github.com/mendenmh I have now included an example scripthttps://github.com/dkriegner/xrayutilities/blob/fpamodel/examples/simpack_powdermodel.py for very basic refinement. Unfortunately some important parameters are currently missing (lattice parameters and other material parameters like atom positions, occupation, Debye Waller factors). Also it runs very slow since it currently takes no benefit from caching. I need to redesign the Crystalhttps://github.com/dkriegner/xrayutilities/blob/master/xrayutilities/materials/material.py#L463 class and include some sort of checksum to allow to detect when the peak positions changed and only then use FP_profile.set_window(). For the fitting I use lmfithttps://lmfit.github.io/lmfit-py/ since I think it offers lots of features which will be useful for this sort of refinement.

But I would appreciate your feedback on the implementation in its current state and if possible suggest a phrasing for these lines L321https://github.com/dkriegner/xrayutilities/blob/fpamodel/xrayutilities/simpack/smaterials.py#L321 L322https://github.com/dkriegner/xrayutilities/blob/fpamodel/xrayutilities/simpack/smaterials.py#L322

OK, the comments that go with the strain_lor and strain_gauss should say something like:

Extra peak width proportional to tan(theta), typically interpreted as microstrain broadening.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://github.com/dkriegner/xrayutilities/issues/32#issuecomment-272386465, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AXc4rHrc6syroEtFH9NRcekNtJr1LslZks5rRzPLgaJpZM4LONwA.

mendenmh commented 7 years ago

A few comments about the optimizer:

For standard LaB6, the crystal size and strain are very hard to refine, because the crystals are large and they have basically no strain. This is why this material is well-liked for calibration. Also, the sample attenuation is very high, making it hard to refine, too.

For this material, the sample displacement and zero offset can be refined. I would at least start the crystallite size frozen at something like 1e-6 (meters), which is physically plausible. Also, an attenuation of about 3e4 m^-1 (300 cm^-1) is physically realistic for the typical SRM660x powder that is used. Freezing it here (at least for a few passes of the fit) is a good starting example. Note that simultaneously refining the zero offset and attenuation of a material is rarely stable for materials with a high attenuation, since the peak shifts are very small and the correlation high. For a high-attenuation sample, there is no chance of getting convergence on the sample thickness. Except for very-low-absorption materials, or very thin specimens, this is usually left out, since the sample is very many absorption thicknesses deep, and the back side is invisible. It is certainly the case for a high-Z, high-attenuation material such as LaB6.

It sort of depends whether we want this example to be just a let-it-run example to show how all parameters can vary, or one with some physical realism and good practice in it. I'm not sure which is more appropriate here.

On Jan 13, 2017, at 3:20 AM, Dominik Kriegner notifications@github.com<mailto:notifications@github.com> wrote:

@mendenmhhttps://github.com/mendenmh I have now included an example scripthttps://github.com/dkriegner/xrayutilities/blob/fpamodel/examples/simpack_powdermodel.py for very basic refinement. Unfortunately some important parameters are currently missing (lattice parameters and other material parameters like atom positions, occupation, Debye Waller factors). Also it runs very slow since it currently takes no benefit from caching. I need to redesign the Crystalhttps://github.com/dkriegner/xrayutilities/blob/master/xrayutilities/materials/material.py#L463 class and include some sort of checksum to allow to detect when the peak positions changed and only then use FP_profile.set_window(). For the fitting I use lmfithttps://lmfit.github.io/lmfit-py/ since I think it offers lots of features which will be useful for this sort of refinement.

But I would appreciate your feedback on the implementation in its current state and if possible suggest a phrasing for these lines L321https://github.com/dkriegner/xrayutilities/blob/fpamodel/xrayutilities/simpack/smaterials.py#L321 L322https://github.com/dkriegner/xrayutilities/blob/fpamodel/xrayutilities/simpack/smaterials.py#L322

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://github.com/dkriegner/xrayutilities/issues/32#issuecomment-272386465, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AXc4rHrc6syroEtFH9NRcekNtJr1LslZks5rRzPLgaJpZM4LONwA.

dkriegner commented 7 years ago

Generally, the window should start out wide enough to cover routine, small shifts in peak positions which result from refinement. There is no need to have the window exactly track the peak. I will look at the example later today.

Ah then I understood something wrong. This means I can call FP_profile.set_parameters(twotheta0_deg=ttpeak) without resetting the window?! I will consider that to improve the speed of the calculation.

Still I need to make changes to the Crystal objects to enable semi-automatic detection/declaration of fitting parameters. I will look into that, hopefully next week but from my side most of the remaining wishes go beyond the original aim of this issue and when you do not see major problems in the current design of PowderDiffraction and PowderModel I suggest to make separate issues for the remaining glitches.

A few comments about the optimizer: ...

I agree with your comments. I think the example should be as realistic as possible. I believe that all parameters which you propose to fix do not make a big impact on the quality of the refinement.

mendenmh commented 7 years ago

On Jan 13, 2017, at 8:10 AM, Dominik Kriegner notifications@github.com<mailto:notifications@github.com> wrote:

Generally, the window should start out wide enough to cover routine, small shifts in peak positions which result from refinement. There is no need to have the window exactly track the peak. I will look at the example later today.

Ah then I understood something wrong. This means I can call FP_profile.set_parameters(twotheta0_deg=ttpeak) without resetting the window?!

Yes, absolutely. As long as the peak stays close enough to the center of the window that its tails are not hanging off the end, this is the appropriate usage.

I will consider that to improve the speed of the calculation.

Still I need to make changes to the Crystal objects to enable semi-automatic detection/declaration of fitting parameters. I will look into that, hopefully next week but from my side most of the remaining wishes go beyond the original aim of this issue and when you do not see major problems in the current design of PowderDiffractionhttps://github.com/dkriegner/xrayutilities/blob/fpamodel/xrayutilities/simpack/powder.py#L1601 and PowderModelhttps://github.com/dkriegner/xrayutilities/blob/fpamodel/xrayutilities/simpack/powdermodel.py#L38 I suggest to make separate issues for the remaining glitches.

A few comments about the optimizer: ...

I agree with your comments. I think the example should be as realistic as possible. I believe that all parameters which you propose to fix do not make a big impact on the quality of the refinement.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://github.com/dkriegner/xrayutilities/issues/32#issuecomment-272439566, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AXc4rBF17nkC77pEDyA2A6EWWAjeEeRYks5rR3fGgaJpZM4LONwA.

dkriegner commented 7 years ago

@mendenmh with the latest commit I have fixed some errors in the fitting and updated the example script to first only fit the displacement and the zero error. in another optimization step I have then optimized the displacement and absorption. I now also avoid calling set_window when its not necessary. Now it would be quite interesting to use multiprocessing to speed up the calculation of one powder pattern in PowderDiffraction.Convolve() since most of the heavy lifting can be done for every Bragg peak independently.

I so far did not use multiprocessing, so i am not experienced in this! But I do not understand how this will work efficiently since upon unpickling the set_window function is called and will therefore destroy any caching?

mendenmh commented 7 years ago

I am attaching the code I use to do multiprocessing on this. It has never really been prepared for release, but it is straightforward to use. It basically makes an object that looks a lot like a lit of peak generators, and distributes this list over the multiprocessing back end transparently.

On Jan 16, 2017, at 11:09 AM, Dominik Kriegner notifications@github.com<mailto:notifications@github.com> wrote:

@mendenmhhttps://github.com/mendenmh with the latest commithttps://github.com/dkriegner/xrayutilities/commit/e44f0c751a27fa510161e91f50bebd81756fd2c1 I have fixed some errors in the fitting and updated the example script to first only fit the displacement and the zero error. in another optimization step I have then optimized the displacement and absorption. I now also avoid calling set_window when its not necessary. Now it would be quite interesting to use multiprocessing to speed up the calculation of one powder pattern in PowderDiffraction.Convolve() since most of the heavy lifting can be done for every Bragg peak independently.

I so far did not use multiprocessing, so i am not experienced in this! But I do not understand how this will work efficiently since upon unpickling the set_window function is called and will therefore destroy any caching?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://github.com/dkriegner/xrayutilities/issues/32#issuecomment-272902371, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AXc4rFhTj0AIc8m3_V0jeomvVzfDWwJNks5rS5ZKgaJpZM4LONwA.

dkriegner commented 7 years ago

@mendenmh If you do not object I will close this issue and merge the fpamodel branch back into master. For the remaining problems I have identified so far (all of them in principle limitations of my implementation and not of FP_profile) I have opened new issues:

Please see #34, #35, #36 and add comments or issues as you see the need for them.

mendenmh commented 7 years ago

This seems a reasonable plan.

Once everything looks as if it is working, I can also contribute some other bits I have to do new FPA convolutions. I have code to do Warren model stacking and deformation faulting in hexagonal crystals, and code to do log-normal crystallite size distributions for a few important geometries. These are all trivial plug-in extensions to the basic FPA class (i.e. derive the new class from the one you have with these inherited, too.)

I haven't sent these all at once, since I thought it would get confusing. Let me know when you are ready for the extra bits. They should take minutes to integrate, but a while to document, since they have not been prepared for release at all.

Marcus Mendenhall

Materials Measurement Science Division National Institute of Standards and Technology 100 Bureau Dr. stop 8370 (217/B115) Gaithersburg MD 20899 USA phone: +1-301-975-8631

On Jan 18, 2017, at 9:27 AM, Dominik Kriegner notifications@github.com<mailto:notifications@github.com> wrote:

@mendenmhhttps://github.com/mendenmh If you do not object I will close this issue and merge the fpamodelhttps://github.com/dkriegner/xrayutilities/tree/fpamodel branch back into master. For the remaining problems I have identified so far (all of them in principle limitations of my implementation and not of FP_profile) I have opened new issues:

Please see #34https://github.com/dkriegner/xrayutilities/issues/34, #35https://github.com/dkriegner/xrayutilities/issues/35, #36https://github.com/dkriegner/xrayutilities/issues/36 and add comments or issues as you see the need for them.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://github.com/dkriegner/xrayutilities/issues/32#issuecomment-273488603, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AXc4rPNuahawdRhl42YJD6DDJuuPJ8f9ks5rTiFTgaJpZM4LONwA.

dkriegner commented 7 years ago

I am looking forward for these extra bits once the current code is ready for them. Thanks for the input so far!