Changeset 8a5f021 in sasmodels


Ignore:
Timestamp:
May 31, 2017 10:02:34 AM (7 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

Files:
3 edited

Legend:

Unmodified
Added
Removed
  • .gitignore

    r2a55a6f r8a5f021  
    2020/.idea 
    2121/sasmodels.egg-info/ 
     22/example/Fit_*/ 
     23/example/batch_fit.csv 
  • example/batch_fit.py

    • Property mode changed from 100644 to 100755
    r93d0ea7 r8a5f021  
    1 ''' 
     1#!/usr/bin/env python 
     2""" 
    23Script to run a batch fit in a series of files and plot the fitted parameters. 
    34 
    4 Usage syntax: 
     5Usage syntax:: 
    56 
    6     python batch_fit.py model.py "sample1.dat, sample2.dat, ..., other_sample.dat" 
     7    python batch_fit.py model.py sample1.dat sample2.dat ... other_sample.dat 
    78    (files named sample1.dat, sample2.dat, ..., other_sample.dat) 
    8      
    9     or if the file names are numbers (and the extension is .dat): 
    109 
    11                 python batch_fit.py model.py 93190 93210  
    12                 (files named 093190.dat, 093191.dat, ..., 093210.dat) 
    13      
    14     or for Grasp-like naming: 
    15      
    16                 python batch_fit.py model.py 93190 93210 200 
    17                 (files named 093190_200.dat, 093191_201.dat, ..., 093210_202.dat) 
    18      
    19     The script reads a series of files and fits the model defined by model.py. 
    20     E.g. python batch_fit.py  model_ellipsoid_hayter_msa.py fits a model  
    21     consisting in an ellipsoid form factor multiplied by a Hayter MSA structure factor.   
    22      
    23     The file model.py must load the data using data = load_data('data.txt'), as the script 
    24     replaces 'data.txt' by the files given here. 
    25      
    26     Modify the call to bumps and the options (minimizer, steps, etc.) as desired. 
    27      
    28     For each file a directory named Fit_filename is created. There the file fit.par contains 
    29     the fitted parameters. 
    30      
    31     Finally the fitted parameters are shown for the full series. 
    32      
    33 Example:      
     10or if the file names are numbers (and the extension is .dat):: 
     11 
     12    python batch_fit.py model.py 93190 93210 
     13    (files named 093190.dat, 093191.dat, ..., 093210.dat) 
     14 
     15or for Grasp-like naming:: 
     16 
     17    python batch_fit.py model.py 93190 93210 200 
     18    (files named 093190_200.dat, 093191_201.dat, ..., 093210_220.dat) 
     19 
     20The script reads a series of files and fits the model defined by model.py. 
     21For example model_ellipsoid_hayter_msa.py fits a model consisting in an 
     22ellipsoid form factor multiplied by a Hayter MSA structure factor. 
     23 
     24The file *model.py* must load the data using:: 
     25 
     26    data = load_data(sys.argv[1]) 
     27 
     28Include options to bumps (minimizer, steps, etc.) as desired.  For example:: 
     29 
     30    python batch_fit.py model.py 93190 93210 200 --fit=lm --steps=200 --ftol=1.5e-8 --xtol=1.5e-8 
     31 
     32Fit options can come before or after the model and files. 
     33 
     34For each file a directory named Fit_filename is created. There the file 
     35model.par contains the fitted parameters.  These are gathered together into 
     36batch_fit.csv in the current directory. 
     37 
     38Finally the fitted parameters are plotted for the full series. 
     39 
     40Example:: 
    3441 
    3542    python batch_fit.py model_ellipsoid_hayter_msa.py 93191 93195 201 
    36      
    37 Note:  
    3843 
    39     If sasmodels is not in the path, edit the line sys.path.append to provide the  
    40     right path to sasmodels.     
     44Note: 
    4145 
    42 ''' 
     46    If sasmodels, sasview or bumps are not in the path, use the PYTHONPATH 
     47    environment variable to set them. 
     48""" 
     49from __future__ import print_function 
    4350 
    44 from __future__ import print_function 
    4551import sys 
    4652import os 
     53 
    4754import numpy as np 
    4855import matplotlib.pyplot as plt 
     56from bumps.dream.views import tile_axes  # make a grid of plots 
    4957 
    50 ''' GET INPUT AND ENSURE MODEL AND DATA FILES ARE DEFINED''' 
     58# GET INPUT AND ENSURE MODEL AND DATA FILES ARE DEFINED 
    5159 
    52 nargs = len(sys.argv) - 1 
    53 if (nargs < 2) or (nargs > 4): 
     60fit_opts = [v for v in sys.argv[1:] if v.startswith('--')] 
     61args = [v for v in sys.argv[1:] if not v.startswith('--')] 
     62 
     63nargs = len(args) 
     64if nargs < 2: 
    5465    print ("Error in the list of arguments! \n") 
    5566    sys.exit() 
    5667 
    57 modelName = sys.argv[1] 
    58 f = open(modelName, 'r') 
    59 fileModel = f.read() 
    60 f.close() 
     68model_file = args[0] 
     69if not model_file.endswith('.py'): 
     70    print("Expected model.py as the first argument") 
     71    sys.exit(1) 
    6172 
    62 if nargs == 2: 
    63     dataFiles = sys.argv[2].split(',') 
    64 else:  
    65     numorFirst = int(sys.argv[2]) 
    66     numorLast = int(sys.argv[3]) 
    67     dataFiles = [] 
     73if '.' in args[1]: 
     74    data_files = args[1:] 
     75else: 
     76    first = int(args[1]) 
     77    last = int(args[2]) 
     78    count = last-first+1 
     79    data_files = [] 
    6880    if nargs == 3: 
    69         for i in range(numorFirst, numorLast+1): 
    70             name = str(i).zfill(6) + '.dat' 
    71             dataFiles.append(name) 
     81        data_files = ['%06d.dat'%(first+i) for i in range(count)] 
     82    elif nargs == 4: 
     83        ext = int(args[3]) 
     84        data_files = ['%06d_%d.dat'%(first+i, ext+i) for i in range(count)] 
    7285    else: 
    73         numorExt = int(sys.argv[4])     
    74         for i in range(numorFirst, numorLast+1): 
    75             name = str(i).zfill(6) + '_' + str(numorExt) + '.dat' 
    76             numorExt += 1 
    77             dataFiles.append(name) 
     86        print("Unexpected arguments: " + " ".join(args[4:])) 
     87        sys.exit(1) 
    7888 
    79 for file in dataFiles: 
    80     if not os.path.isfile(file.strip()): 
    81         print ("File %s does not exist! \n" % file.strip()) 
    82         sys.exit() 
     89# CHECK THAT DATA FILES EXIST 
     90missing = [filename for filename in data_files if not os.path.isfile(filename)] 
     91if missing: 
     92    print("Missing data files: %s" % ", ".join(missing)) 
     93    sys.exit(1) 
    8394 
    84          
    85 ''' CALL TO BUMPS AND DEFINITION OF FITTING OPTIONS ''' 
     95# STORE DIRECTORY FOR BUMPS FITS 
     96def fit_dir(filename): 
     97    "Return the store directory name for the given file" 
     98    return "Fit_" + os.path.splitext(filename)[0] 
    8699 
    87 msg0 = "python -m bumps.cli fit.py" 
    88 options = " --fit=lm --steps=200 --ftol=1.5e-8 --xtol=1.5e-8 --batch" 
     100# LOOP OVER FILES AND CALL TO BUMPS FOR EACH OF THEM 
     101bumps_cmd = "python -m bumps.cli --batch" 
     102fit_opts = " ".join(fit_opts) 
     103for data_file in data_files: 
     104    store_opts = "--store=" + fit_dir(data_file) 
     105    cmd = " ".join((bumps_cmd, fit_opts, store_opts, model_file, data_file)) 
     106    os.system(cmd) 
    89107 
     108# GATHER RESULTS 
     109results = {} 
     110par_file = os.path.splitext(model_file)[0] + '.par' 
     111for data_file in data_files: 
     112    with open(os.path.join(fit_dir(data_file), par_file), 'r') as fid: 
     113        for line in fid: 
     114            parameter, value = line.split() 
     115            results.setdefault(parameter, []).append(float(value)) 
    90116 
    91 ''' LOOP OVER FILES AND CALL TO BUMPS FOR EACH OF THEM'''    
     117# SAVE RESULTS INTO FILE 
     118with open('batch_fit.csv', 'w') as fid: 
     119    parameters = list(sorted(results.keys())) 
     120    values_by_file = zip(*(v for k, v in sorted(results.items()))) 
     121    fid.write(','.join(['filename'] + parameters) + '\n') 
     122    for filename, values in zip(data_files, values_by_file): 
     123        fid.write(','.join([filename] + [str(v) for v in values]) + '\n') 
    92124 
    93 for file in dataFiles: 
    94     currentModel = fileModel.replace('data.txt', file.strip()) 
    95     f = open('fit.py', 'w') 
    96     f.write(currentModel) 
    97     f.close() 
    98     store = " --store=Fit_" + file.strip() 
    99     msg = msg0 + options + store 
    100     os.system(msg) 
    101  
    102  
    103 ''' SHOW FITTED PARAMETERS '''     
    104  
    105 parDict = {} 
    106 for file in dataFiles: 
    107     parFile = os.path.join('Fit_' + file.strip(), 'fit.par') 
    108     f = open(parFile, 'r') 
    109     lines = f.readlines() 
    110     for line in lines: 
    111         parName = line.split()[0] 
    112         if parName in parDict.keys(): 
    113             parList = parDict[parName] 
    114             parList.append(line.split()[1]) 
    115             parDict[parName] = parList 
    116         else: 
    117             parDict[parName] = [line.split()[1]] 
    118      
    119 for parName in parDict.keys(): 
    120     values = np.array(map(float, parDict[parName]))    
    121     plt.plot(values) 
    122     plt.title(parName) 
    123     plt.xlabel('Dataset #') 
    124     plt.ylabel('Value') 
    125     plt.show() 
    126  
     125# SHOW FITTED PARAMETERS 
     126nh, nw = tile_axes(len(results)) 
     127ticks = np.arange(1, len(data_files)+1) 
     128labels = [os.path.splitext(filename)[0] for filename in data_files] 
     129for k, (parameter, values) in enumerate(sorted(results.items())): 
     130    plt.subplot(nh, nw, k+1) 
     131    plt.plot(ticks, values) 
     132    plt.xlim(ticks[0]-0.5, ticks[-1]+0.5) 
     133    if k%nh == nh-1: 
     134        #plt.xlabel('Dataset #') 
     135        plt.xticks(ticks, labels, rotation=30) 
     136    else: 
     137        plt.xticks(ticks, [' ']*len(labels)) 
     138    plt.ylabel(parameter) 
     139plt.suptitle("Fit " + args[0]) 
     140plt.show() 
  • example/model_ellipsoid_hayter_msa.py

    r1f7792e r8a5f021  
    1 #!/usr/bin/env python 
    2 # -*- coding: utf-8 -*- 
    3  
    4 # To Sasview/documents/scripts 
    5  
    61import sys 
    72#sys.path.append('path_to_sasmodels') 
     3 
     4import numpy as np 
    85 
    96from bumps.names import * 
     
    118from sasmodels.bumps_model import Model, Experiment 
    129from sasmodels.data import load_data, plot_data 
    13 import numpy as np 
    14 import matplotlib.pyplot as plt 
    1510 
    16  
    17 """ IMPORT THE DATA USED """ 
    18 data = load_data('data.txt') 
     11# IMPORT THE DATA USED 
     12data = load_data(sys.argv[1]) 
    1913 
    2014#setattr(data, 'qmin', 0.0) 
    2115#setattr(data, 'qmax', 10.0) 
    2216 
    23 """ DEFINE THE MODEL """ 
    24 kernel = load_model('ellipsoid*hayter_msa', dtype="double") 
     17# DEFINE THE MODEL 
     18kernel = load_model('ellipsoid*hayter_msa') 
    2519 
    26 pars = dict(scale=6.4, background=0.06, sld=0.33, sld_solvent=2.15, radius_polar=14.0,  
    27             radius_equatorial=24.0, volfraction=0.075, charge=66.373, temperature=298.0,  
     20pars = dict(scale=6.4, background=0.06, sld=0.33, sld_solvent=2.15, radius_polar=14.0, 
     21            radius_equatorial=24.0, volfraction=0.075, charge=66.373, temperature=298.0, 
    2822            concentration_salt=0.001, dielectconst=71.0) 
    29              
     23 
    3024model = Model(kernel, **pars) 
    3125 
    32 """ PARAMETER RANGES (ONLY THOSE PARAMETERS ARE FITTED) """ 
     26# PARAMETER RANGES (ONLY THOSE PARAMETERS ARE FITTED) 
    3327model.scale.range(0, inf) 
    3428model.background.range(-inf, inf) 
     
    4337#model.dielectconst.range(0,inf) 
    4438 
    45  
    4639M = Experiment(data=data, model=model) 
    4740 
    4841problem = FitProblem(M) 
    49  
Note: See TracChangeset for help on using the changeset viewer.