diff --git a/doc/examples/srfit_parallel_test.py b/doc/examples/srfit_parallel_test.py new file mode 100644 index 00000000..32f35674 --- /dev/null +++ b/doc/examples/srfit_parallel_test.py @@ -0,0 +1,198 @@ +#!/usr/bin/env python +######################################################################## +# +# diffpy.srfit by DANSE Diffraction group +# Simon J. L. Billinge +# (c) 2009 The Trustees of Columbia University +# in the City of New York. All rights reserved. +# +# File coded by: Chris Farrow +# +# See AUTHORS.txt for a list of people who contributed. +# See LICENSE_DANSE.txt for license information. +# +######################################################################## + +"""Example of a PDF refinement using diffpy.structure and PDFGenerator. + +This is example of fitting the fcc nickel structure to measured PDF data. The +purpose of this example is to demonstrate and describe the classes in +configuraiton options involved with setting up a fit in this way. The main +benefit of using SrFit for PDF refinement is the flexibility of modifying the +PDF profile function for specific needs, adding restraints to a fit and the +ability to simultaneously refine a structure to PDF data and data from other +sources. This example demonstrates only the basic configuration. + +""" + +import numpy + +from diffpy.structure import Structure +from diffpy.srfit.pdf import PDFGenerator, PDFParser +from diffpy.srfit.fitbase import Profile +from diffpy.srfit.fitbase import FitContribution, FitRecipe +from diffpy.srfit.fitbase import FitResults + +from gaussianrecipe import scipyOptimize + +####### Example Code + +def makeRecipe(ciffile, datname): + """Create a fitting recipe for crystalline PDF data.""" + + ## The Profile + # This will be used to store the observed and calculated PDF profile. + profile = Profile() + + # Load data and add it to the Profile. Unlike in other examples, we use a + # class (PDFParser) to help us load the data. This class will read the data + # and relevant metadata from a two- to four-column data file generated + # with PDFGetX2 or PDFGetN. The metadata will be passed to the PDFGenerator + # when they are associated in the FitContribution, which saves some + # configuration steps. + parser = PDFParser() + parser.parseFile(datname) + profile.loadParsedData(parser) + profile.setCalculationRange(xmax = 20) + + ## The ProfileGenerator + # The PDFGenerator is for configuring and calculating a PDF profile. Here, + # we want to refine a Structure object from diffpy.structure. We tell the + # PDFGenerator that with the 'setStructure' method. All other configuration + # options will be inferred from the metadata that is read by the PDFParser. + # In particular, this will set the scattering type (x-ray or neutron), the + # Qmax value, as well as initial values for the non-structural Parameters. + generator = PDFGenerator("G") + stru = Structure() + stru.read(ciffile) + generator.setStructure(stru) + + ## The FitContribution + # Here we associate the Profile and ProfileGenerator, as has been done + # before. + contribution = FitContribution("nickel") + contribution.addProfileGenerator(generator) + contribution.setProfile(profile, xname = "r") + + ## Make the FitRecipe and add the FitContribution. + recipe = FitRecipe() + recipe.addContribution(contribution) + + ## Configure the fit variables + + # The PDFGenerator class holds the ParameterSet associated with the + # Structure passed above in a data member named "phase". (We could have + # given the ParameterSet a name other than "phase" when we added it to the + # PDFGenerator.) The ParameterSet in this case is a StructureParameterSet, + # the documentation for which is found in the + # diffpy.srfit.structure.diffpystructure module. + phase = generator.phase + + # We start by constraining the phase to the known space group. We could do + # this by hand, but there is a method in diffpy.srfit.structure named + # 'constrainAsSpaceGroup' for this purpose. The constraints will by default + # be applied to the sites, the lattice and to the ADPs. See the method + # documentation for more details. The 'constrainAsSpaceGroup' method may + # create new Parameters, which it returns in a SpaceGroupParameters object. + from diffpy.srfit.structure import constrainAsSpaceGroup + sgpars = constrainAsSpaceGroup(phase, "Fm-3m") + + # The SpaceGroupParameters object returned by 'constrainAsSpaceGroup' holds + # the free Parameters allowed by the space group constraints. Once a + # structure is constrained, we need (should) only use the Parameters + # provided in the SpaceGroupParameters, as the relevant structure + # Parameters are constrained to these. + # + # We know that the space group does not allow for any free sites because + # each atom is on a special position. There is one free (cubic) lattice + # parameter and one free (isotropic) ADP. We can access these Parameters in + # the xyzpars, latpars, and adppars members of the SpaceGroupParameters + # object. + for par in sgpars.latpars: + recipe.addVar(par) + for par in sgpars.adppars: + recipe.addVar(par, 0.005) + + # We now select non-structural parameters to refine. + # This controls the scaling of the PDF. + recipe.addVar(generator.scale, 1) + # This is a peak-damping resolution term. + recipe.addVar(generator.qdamp, 0.01) + # This is a vibrational correlation term that sharpens peaks at low-r. + recipe.addVar(generator.delta2, 5) + + # Give the recipe away so it can be used! + return recipe + +def plotResults(recipe): + """Plot the results contained within a refined FitRecipe.""" + + # All this should be pretty familiar by now. + r = recipe.nickel.profile.x + g = recipe.nickel.profile.y + gcalc = recipe.nickel.profile.ycalc + diffzero = -0.8 * max(g) * numpy.ones_like(g) + diff = g - gcalc + diffzero + + import pylab + pylab.plot(r,g,'bo',label="G(r) Data") + pylab.plot(r, gcalc,'r-',label="G(r) Fit") + pylab.plot(r,diff,'g-',label="G(r) diff") + pylab.plot(r,diffzero,'k-') + pylab.xlabel(r"$r (\AA)$") + pylab.ylabel(r"$G (\AA^{-2})$") + pylab.legend(loc=1) + + pylab.show() + return + + +def fit_wrapper(recipe): + """helper function""" + from scipy.optimize.minpack import leastsq + leastsq(recipe.residual, recipe.getValues()) + return recipe + +if __name__ == "__main__": + + # Make the data and the recipe + ciffile = "data/ni.cif" + data = "data/ni-q27r100-neutron.gr" + + # sanity check + print("==== Rw before refinements ====") + recipe = makeRecipe(ciffile, data) + recipe.clearFitHooks() + res = FitResults(recipe) + print(res.rw) + + # Make the recipe + recipe_list = [] + for i in range(5): + recipe = makeRecipe(ciffile, data) + recipe.clearFitHooks() + recipe_list.append(recipe) + + # Optimize + print("==== Rw: Sequencial refinements ====") + for recipe in recipe_list: + recipe = fit_wrapper(recipe) + res = FitResults(recipe) + print(res.rw) + + # Make the recipe + recipe_list = [] + for i in range(5): + recipe = makeRecipe(ciffile, data) + recipe.clearFitHooks() + recipe_list.append(recipe) + + # Optimize + print("==== Rw: Parallel refinements ====") + import multiprocessing as mp + n_process = 4 + with mp.Pool(n_process) as p: + rv = p.map(fit_wrapper, recipe_list) + for recipe in rv: + res = FitResults(recipe) + print(res.rw)