Changeset 8a5f021 in sasmodels for example

Timestamp:

May 31, 2017 12:02:34 PM (8 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

8bc1a4b

Parents:

1f7792e

Message:

batch fit example: PEP8 cleanup; add .csv output; all pars on one plot; more flexible option handling

Location:

example

Files:

• batch_fit.py (modified) (1 diff, 1 prop)
• model_ellipsoid_hayter_msa.py (modified) (3 diffs)

example/batch_fit.py

@@ -1 @@
-'''
+#!/usr/bin/env python
+"""
 Script to run a batch fit in a series of files and plot the fitted parameters.
 
-Usage syntax:
+Usage syntax::
 
-    python batch_fit.py model.py "sample1.dat, sample2.dat, ..., other_sample.dat"
+    python batch_fit.py model.py sample1.dat sample2.dat ... other_sample.dat
     (files named sample1.dat, sample2.dat, ..., other_sample.dat)
-    
-    or if the file names are numbers (and the extension is .dat):
 
-                python batch_fit.py model.py 93190 93210 
-                (files named 093190.dat, 093191.dat, ..., 093210.dat)
-    
-    or for Grasp-like naming:
-    
-                python batch_fit.py model.py 93190 93210 200
-                (files named 093190_200.dat, 093191_201.dat, ..., 093210_202.dat)
-    
-    The script reads a series of files and fits the model defined by model.py.
-    E.g. python batch_fit.py  model_ellipsoid_hayter_msa.py fits a model 
-    consisting in an ellipsoid form factor multiplied by a Hayter MSA structure factor.  
-    
-    The file model.py must load the data using data = load_data('data.txt'), as the script
-    replaces 'data.txt' by the files given here.
-    
-    Modify the call to bumps and the options (minimizer, steps, etc.) as desired.
-    
-    For each file a directory named Fit_filename is created. There the file fit.par contains
-    the fitted parameters.
-    
-    Finally the fitted parameters are shown for the full series.
-    
-Example:     
+or if the file names are numbers (and the extension is .dat)::
+
+    python batch_fit.py model.py 93190 93210
+    (files named 093190.dat, 093191.dat, ..., 093210.dat)
+
+or for Grasp-like naming::
+
+    python batch_fit.py model.py 93190 93210 200
+    (files named 093190_200.dat, 093191_201.dat, ..., 093210_220.dat)
+
+The script reads a series of files and fits the model defined by model.py.
+For example model_ellipsoid_hayter_msa.py fits a model consisting in an
+ellipsoid form factor multiplied by a Hayter MSA structure factor.
+
+The file *model.py* must load the data using::
+
+    data = load_data(sys.argv[1])
+
+Include options to bumps (minimizer, steps, etc.) as desired.  For example::
+
+    python batch_fit.py model.py 93190 93210 200 --fit=lm --steps=200 --ftol=1.5e-8 --xtol=1.5e-8
+
+Fit options can come before or after the model and files.
+
+For each file a directory named Fit_filename is created. There the file
+model.par contains the fitted parameters.  These are gathered together into
+batch_fit.csv in the current directory.
+
+Finally the fitted parameters are plotted for the full series.
+
+Example::
 
     python batch_fit.py model_ellipsoid_hayter_msa.py 93191 93195 201
-    
-Note: 
 
-    If sasmodels is not in the path, edit the line sys.path.append to provide the 
-    right path to sasmodels.    
+Note:
 
-'''
+    If sasmodels, sasview or bumps are not in the path, use the PYTHONPATH
+    environment variable to set them.
+"""
+from __future__ import print_function
 
-from __future__ import print_function
 import sys
 import os
+
 import numpy as np
 import matplotlib.pyplot as plt
+from bumps.dream.views import tile_axes  # make a grid of plots
 
-''' GET INPUT AND ENSURE MODEL AND DATA FILES ARE DEFINED'''
+# GET INPUT AND ENSURE MODEL AND DATA FILES ARE DEFINED
 
-nargs = len(sys.argv) - 1
-if (nargs < 2) or (nargs > 4):
+fit_opts = [v for v in sys.argv[1:] if v.startswith('--')]
+args = [v for v in sys.argv[1:] if not v.startswith('--')]
+
+nargs = len(args)
+if nargs < 2:
     print ("Error in the list of arguments! \n")
     sys.exit()
 
-modelName = sys.argv[1]
-f = open(modelName, 'r')
-fileModel = f.read()
-f.close()
+model_file = args[0]
+if not model_file.endswith('.py'):
+    print("Expected model.py as the first argument")
+    sys.exit(1)
 
-if nargs == 2:
-    dataFiles = sys.argv[2].split(',')
-else: 
-    numorFirst = int(sys.argv[2])
-    numorLast = int(sys.argv[3])
-    dataFiles = []
+if '.' in args[1]:
+    data_files = args[1:]
+else:
+    first = int(args[1])
+    last = int(args[2])
+    count = last-first+1
+    data_files = []
     if nargs == 3:
-        for i in range(numorFirst, numorLast+1):
-            name = str(i).zfill(6) + '.dat'
-            dataFiles.append(name)
+        data_files = ['%06d.dat'%(first+i) for i in range(count)]
+    elif nargs == 4:
+        ext = int(args[3])
+        data_files = ['%06d_%d.dat'%(first+i, ext+i) for i in range(count)]
     else:
-        numorExt = int(sys.argv[4])    
-        for i in range(numorFirst, numorLast+1):
-            name = str(i).zfill(6) + '_' + str(numorExt) + '.dat'
-            numorExt += 1
-            dataFiles.append(name)
+        print("Unexpected arguments: " + " ".join(args[4:]))
+        sys.exit(1)
 
-for file in dataFiles:
-    if not os.path.isfile(file.strip()):
-        print ("File %s does not exist! \n" % file.strip())
-        sys.exit()
+# CHECK THAT DATA FILES EXIST
+missing = [filename for filename in data_files if not os.path.isfile(filename)]
+if missing:
+    print("Missing data files: %s" % ", ".join(missing))
+    sys.exit(1)
 
-        
-''' CALL TO BUMPS AND DEFINITION OF FITTING OPTIONS '''
+# STORE DIRECTORY FOR BUMPS FITS
+def fit_dir(filename):
+    "Return the store directory name for the given file"
+    return "Fit_" + os.path.splitext(filename)[0]
 
-msg0 = "python -m bumps.cli fit.py"
-options = " --fit=lm --steps=200 --ftol=1.5e-8 --xtol=1.5e-8 --batch"
+# LOOP OVER FILES AND CALL TO BUMPS FOR EACH OF THEM
+bumps_cmd = "python -m bumps.cli --batch"
+fit_opts = " ".join(fit_opts)
+for data_file in data_files:
+    store_opts = "--store=" + fit_dir(data_file)
+    cmd = " ".join((bumps_cmd, fit_opts, store_opts, model_file, data_file))
+    os.system(cmd)
 
+# GATHER RESULTS
+results = {}
+par_file = os.path.splitext(model_file)[0] + '.par'
+for data_file in data_files:
+    with open(os.path.join(fit_dir(data_file), par_file), 'r') as fid:
+        for line in fid:
+            parameter, value = line.split()
+            results.setdefault(parameter, []).append(float(value))
 
-''' LOOP OVER FILES AND CALL TO BUMPS FOR EACH OF THEM'''   
+# SAVE RESULTS INTO FILE
+with open('batch_fit.csv', 'w') as fid:
+    parameters = list(sorted(results.keys()))
+    values_by_file = zip(*(v for k, v in sorted(results.items())))
+    fid.write(','.join(['filename'] + parameters) + '\n')
+    for filename, values in zip(data_files, values_by_file):
+        fid.write(','.join([filename] + [str(v) for v in values]) + '\n')
 
-for file in dataFiles:
-    currentModel = fileModel.replace('data.txt', file.strip())
-    f = open('fit.py', 'w')
-    f.write(currentModel)
-    f.close()
-    store = " --store=Fit_" + file.strip()
-    msg = msg0 + options + store
-    os.system(msg)
-
-
-''' SHOW FITTED PARAMETERS '''    
-
-parDict = {}
-for file in dataFiles:
-    parFile = os.path.join('Fit_' + file.strip(), 'fit.par')
-    f = open(parFile, 'r')
-    lines = f.readlines()
-    for line in lines:
-        parName = line.split()[0]
-        if parName in parDict.keys():
-            parList = parDict[parName]
-            parList.append(line.split()[1])
-            parDict[parName] = parList
-        else:
-            parDict[parName] = [line.split()[1]]
-    
-for parName in parDict.keys():
-    values = np.array(map(float, parDict[parName]))   
-    plt.plot(values)
-    plt.title(parName)
-    plt.xlabel('Dataset #')
-    plt.ylabel('Value')
-    plt.show()
-
+# SHOW FITTED PARAMETERS
+nh, nw = tile_axes(len(results))
+ticks = np.arange(1, len(data_files)+1)
+labels = [os.path.splitext(filename)[0] for filename in data_files]
+for k, (parameter, values) in enumerate(sorted(results.items())):
+    plt.subplot(nh, nw, k+1)
+    plt.plot(ticks, values)
+    plt.xlim(ticks[0]-0.5, ticks[-1]+0.5)
+    if k%nh == nh-1:
+        #plt.xlabel('Dataset #')
+        plt.xticks(ticks, labels, rotation=30)
+    else:
+        plt.xticks(ticks, [' ']*len(labels))
+    plt.ylabel(parameter)
+plt.suptitle("Fit " + args[0])
+plt.show()

example/model_ellipsoid_hayter_msa.py

@@ -1 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-
-# To Sasview/documents/scripts
-
 import sys
 #sys.path.append('path_to_sasmodels')
+
+import numpy as np
 
 from bumps.names import *
@@ -11 +8 @@
 from sasmodels.bumps_model import Model, Experiment
 from sasmodels.data import load_data, plot_data
-import numpy as np
-import matplotlib.pyplot as plt
 
-
-""" IMPORT THE DATA USED """
-data = load_data('data.txt')
+# IMPORT THE DATA USED
+data = load_data(sys.argv[1])
 
 #setattr(data, 'qmin', 0.0)
 #setattr(data, 'qmax', 10.0)
 
-""" DEFINE THE MODEL """
-kernel = load_model('ellipsoid*hayter_msa', dtype="double")
+# DEFINE THE MODEL
+kernel = load_model('ellipsoid*hayter_msa')
 
-pars = dict(scale=6.4, background=0.06, sld=0.33, sld_solvent=2.15, radius_polar=14.0, 
-            radius_equatorial=24.0, volfraction=0.075, charge=66.373, temperature=298.0, 
+pars = dict(scale=6.4, background=0.06, sld=0.33, sld_solvent=2.15, radius_polar=14.0,
+            radius_equatorial=24.0, volfraction=0.075, charge=66.373, temperature=298.0,
             concentration_salt=0.001, dielectconst=71.0)
-            
+
 model = Model(kernel, **pars)
 
-""" PARAMETER RANGES (ONLY THOSE PARAMETERS ARE FITTED) """
+# PARAMETER RANGES (ONLY THOSE PARAMETERS ARE FITTED)
 model.scale.range(0, inf)
 model.background.range(-inf, inf)
@@ -43 +37 @@
 #model.dielectconst.range(0,inf)
 
-
 M = Experiment(data=data, model=model)
 
 problem = FitProblem(M)
-