Changeset ec7e360 in sasmodels


Ignore:
Timestamp:
Dec 23, 2015 12:17:49 PM (8 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
e21cc31
Parents:
ce166d3
Message:

refactor option processing for compare.py, allowing more flexible selection of calculation engines

Files:
4 edited

Legend:

Unmodified
Added
Removed
  • compare.sh

    r5753e4e rec7e360  
    1212#SAS_OPENMP=1; export SAS_OPENMP 
    1313 
    14 python -m sasmodels.compare $* 
     14${PYTHON:-python} -m sasmodels.compare $* 
  • sasmodels/bumps_model.py

    r9404dd3 rec7e360  
    1212""" 
    1313 
    14 import datetime 
    1514import warnings 
    1615 
    1716import numpy as np 
    1817 
    19 from . import sesans 
    20 from . import weights 
    2118from .data import plot_theory 
    2219from .direct_model import DataMixin 
     
    3128    return experiment 
    3229 
     30def create_parameters(model_info, **kwargs): 
     31    # lazy import; this allows the doc builder and nosetests to run even 
     32    # when bumps is not on the path. 
     33    from bumps.names import Parameter 
     34 
     35    partype = model_info['partype'] 
     36 
     37    pars = {} 
     38    for p in model_info['parameters']: 
     39        name, default, limits = p[0], p[2], p[3] 
     40        value = kwargs.pop(name, default) 
     41        pars[name] = Parameter.default(value, name=name, limits=limits) 
     42    for name in partype['pd-2d']: 
     43        for xpart, xdefault, xlimits in [ 
     44            ('_pd', 0., limits), 
     45            ('_pd_n', 35., (0, 1000)), 
     46            ('_pd_nsigma', 3., (0, 10)), 
     47        ]: 
     48            xname = name + xpart 
     49            xvalue = kwargs.pop(xname, xdefault) 
     50            pars[xname] = Parameter.default(xvalue, name=xname, limits=xlimits) 
     51 
     52    pd_types = {} 
     53    for name in partype['pd-2d']: 
     54        xname = name + '_pd_type' 
     55        xvalue = kwargs.pop(xname, 'gaussian') 
     56        pd_types[xname] = xvalue 
     57 
     58    if kwargs:  # args not corresponding to parameters 
     59        raise TypeError("unexpected parameters: %s" 
     60                        % (", ".join(sorted(kwargs.keys())))) 
     61 
     62    return pars, pd_types 
    3363 
    3464class Model(object): 
    35     def __init__(self, model, **kw): 
    36         # lazy import; this allows the doc builder and nosetests to run even 
    37         # when bumps is not on the path. 
    38         from bumps.names import Parameter 
    39  
     65    def __init__(self, model, **kwargs): 
    4066        self._sasmodel = model 
    41         partype = model.info['partype'] 
    42  
    43         pars = [] 
    44         for p in model.info['parameters']: 
    45             name, default, limits = p[0], p[2], p[3] 
    46             value = kw.pop(name, default) 
    47             setattr(self, name, Parameter.default(value, name=name, limits=limits)) 
    48             pars.append(name) 
    49         for name in partype['pd-2d']: 
    50             for xpart, xdefault, xlimits in [ 
    51                 ('_pd', 0, limits), 
    52                 ('_pd_n', 35, (0, 1000)), 
    53                 ('_pd_nsigma', 3, (0, 10)), 
    54                 ('_pd_type', 'gaussian', None), 
    55                 ]: 
    56                 xname = name + xpart 
    57                 xvalue = kw.pop(xname, xdefault) 
    58                 if xlimits is not None: 
    59                     xvalue = Parameter.default(xvalue, name=xname, limits=xlimits) 
    60                     pars.append(xname) 
    61                 setattr(self, xname, xvalue) 
    62         self._parameter_names = pars 
    63         if kw: 
    64             raise TypeError("unexpected parameters: %s" 
    65                             % (", ".join(sorted(kw.keys())))) 
     67        pars, pd_types = create_parameters(model.info, **kwargs) 
     68        for k,v in pars.items(): 
     69            setattr(self, k, v) 
     70        for k,v in pd_types.items(): 
     71            setattr(self, k, v) 
     72        self._parameter_names = list(pars.keys()) 
     73        self._pd_type_names = list(pd_types.keys()) 
    6674 
    6775    def parameters(self): 
     
    7179        return dict((k, getattr(self, k)) for k in self._parameter_names) 
    7280 
     81    def state(self): 
     82        pars = dict((k, getattr(self, k).value) for k in self._parameter_names) 
     83        pars.update((k, getattr(self, k)) for k in self._pd_type_names) 
     84        return pars 
    7385 
    7486class Experiment(DataMixin): 
     
    113125    def theory(self): 
    114126        if 'theory' not in self._cache: 
    115             pars = dict((k, v.value) for k,v in self.model.parameters().items()) 
     127            pars = self.model.state() 
    116128            self._cache['theory'] = self._calc_theory(pars, cutoff=self.cutoff) 
    117             """ 
    118             if self._fn is None: 
    119                 q_input = self.model.kernel.make_input(self._kernel_inputs) 
    120                 self._fn = self.model.kernel(q_input) 
    121  
    122             fixed_pars = [getattr(self.model, p).value for p in self._fn.fixed_pars] 
    123             pd_pars = [self._get_weights(p) for p in self._fn.pd_pars] 
    124             #print(fixed_pars,pd_pars) 
    125             Iq_calc = self._fn(fixed_pars, pd_pars, self.cutoff) 
    126             #self._theory[:] = self._fn.eval(pars, pd_pars) 
    127             if self.data_type == 'sesans': 
    128                 result = sesans.hankel(self.data.x, self.data.lam * 1e-9, 
    129                                        self.data.sample.thickness / 10, 
    130                                        self._kernel_inputs[0], Iq_calc) 
    131                 self._cache['theory'] = result 
    132             else: 
    133                 Iq = self.resolution.apply(Iq_calc) 
    134                 self._cache['theory'] = Iq 
    135             """ 
    136129        return self._cache['theory'] 
    137130 
     
    162155        pass 
    163156 
    164     def remove_get_weights(self, name): 
    165         """ 
    166         Get parameter dispersion weights 
    167         """ 
    168         info = self.model.kernel.info 
    169         relative = name in info['partype']['pd-rel'] 
    170         limits = info['limits'][name] 
    171         disperser, value, npts, width, nsigma = [ 
    172             getattr(self.model, name + ext) 
    173             for ext in ('_pd_type', '', '_pd_n', '_pd', '_pd_nsigma')] 
    174         value, weight = weights.get_weights( 
    175             disperser, int(npts.value), width.value, nsigma.value, 
    176             value.value, limits, relative) 
    177         return value, weight / np.sum(weight) 
    178  
    179157    def __getstate__(self): 
    180158        # Can't pickle gpu functions, so instead make them lazy 
  • sasmodels/compare.py

    rf4f3919 rec7e360  
    6666 
    6767 
    68 def sasview_model(model_definition, **pars): 
    69     """ 
    70     Load a sasview model given the model name. 
    71     """ 
     68def parameter_range(p, v): 
     69    """ 
     70    Choose a parameter range based on parameter name and initial value. 
     71    """ 
     72    if p.endswith('_pd_n'): 
     73        return [0, 100] 
     74    elif p.endswith('_pd_nsigma'): 
     75        return [0, 5] 
     76    elif p.endswith('_pd_type'): 
     77        return v 
     78    elif any(s in p for s in ('theta','phi','psi')): 
     79        # orientation in [-180,180], orientation pd in [0,45] 
     80        if p.endswith('_pd'): 
     81            return [0,45] 
     82        else: 
     83            return [-180, 180] 
     84    elif 'sld' in p: 
     85        return [-0.5, 10] 
     86    elif p.endswith('_pd'): 
     87        return [0, 1] 
     88    elif p in ['background', 'scale']: 
     89        return [0, 1e3] 
     90    else: 
     91        return [0, (2*v if v>0 else 1)] 
     92 
     93def _randomize_one(p, v): 
     94    """ 
     95    Randomizing parameter. 
     96    """ 
     97    if any(p.endswith(s) for s in ('_pd_n','_pd_nsigma','_pd_type')): 
     98        return v 
     99    else: 
     100        return np.random.uniform(*parameter_range(p, v)) 
     101 
     102def randomize_pars(pars, seed=None): 
     103    np.random.seed(seed) 
     104    # Note: the sort guarantees order `of calls to random number generator 
     105    pars = dict((p,_randomize_one(p,v)) 
     106                for p,v in sorted(pars.items())) 
     107    return pars 
     108 
     109def constrain_pars(model_definition, pars): 
     110    """ 
     111    Restrict parameters to valid values. 
     112    """ 
     113    name = model_definition.name 
     114    if name == 'capped_cylinder' and pars['cap_radius'] < pars['radius']: 
     115        pars['radius'],pars['cap_radius'] = pars['cap_radius'],pars['radius'] 
     116    if name == 'barbell' and pars['bell_radius'] < pars['radius']: 
     117        pars['radius'],pars['bell_radius'] = pars['bell_radius'],pars['radius'] 
     118 
     119    # Limit guinier to an Rg such that Iq > 1e-30 (single precision cutoff) 
     120    if name == 'guinier': 
     121        #q_max = 0.2  # mid q maximum 
     122        q_max = 1.0  # high q maximum 
     123        rg_max = np.sqrt(90*np.log(10) + 3*np.log(pars['scale']))/q_max 
     124        pars['rg'] = min(pars['rg'],rg_max) 
     125 
     126def parlist(pars): 
     127    return "\n".join("%s: %s"%(p,v) for p,v in sorted(pars.items())) 
     128 
     129def suppress_pd(pars): 
     130    """ 
     131    Suppress theta_pd for now until the normalization is resolved. 
     132 
     133    May also suppress complete polydispersity of the model to test 
     134    models more quickly. 
     135    """ 
     136    pars = pars.copy() 
     137    for p in pars: 
     138        if p.endswith("_pd"): pars[p] = 0 
     139    return pars 
     140 
     141def eval_sasview(model_definition, data): 
     142    # importing sas here so that the error message will be that sas failed to 
     143    # import rather than the more obscure smear_selection not imported error 
     144    import sas 
     145    from sas.models.qsmearing import smear_selection 
     146 
    72147    # convert model parameters from sasmodel form to sasview form 
    73148    #print("old",sorted(pars.items())) 
    74     modelname, pars = revert_model(model_definition, pars) 
     149    modelname, pars = revert_model(model_definition, {}) 
    75150    #print("new",sorted(pars.items())) 
    76151    sas = __import__('sas.models.'+modelname) 
     
    79154        raise ValueError("could not find model %r in sas.models"%modelname) 
    80155    model = ModelClass() 
    81  
    82     for k,v in pars.items(): 
    83         parts = k.split('.')  # polydispersity components 
    84         if len(parts) == 2: 
    85             model.dispersion[parts[0]][parts[1]] = v 
     156    smearer = smear_selection(data, model=model) 
     157 
     158    if hasattr(data, 'qx_data'): 
     159        q = np.sqrt(data.qx_data**2 + data.qy_data**2) 
     160        index = ((~data.mask) & (~np.isnan(data.data)) 
     161                 & (q >= data.qmin) & (q <= data.qmax)) 
     162        if smearer is not None: 
     163            smearer.model = model  # because smear_selection has a bug 
     164            smearer.accuracy = data.accuracy 
     165            smearer.set_index(index) 
     166            theory = lambda: smearer.get_value() 
    86167        else: 
    87             model.setParam(k, v) 
    88     return model 
    89  
    90 def randomize(p, v): 
    91     """ 
    92     Randomizing parameter. 
    93  
    94     Guess the parameter type from name. 
    95     """ 
    96     if any(p.endswith(s) for s in ('_pd_n','_pd_nsigma','_pd_type')): 
    97         return v 
    98     elif any(s in p for s in ('theta','phi','psi')): 
    99         # orientation in [-180,180], orientation pd in [0,45] 
    100         if p.endswith('_pd'): 
    101             return 45*np.random.rand() 
    102         else: 
    103             return 360*np.random.rand() - 180 
    104     elif 'sld' in p: 
    105         # sld in in [-0.5,10] 
    106         return 10.5*np.random.rand() - 0.5 
    107     elif p.endswith('_pd'): 
    108         # length pd in [0,1] 
    109         return np.random.rand() 
    110     else: 
    111         # values from 0 to 2*x for all other parameters 
    112         return 2*np.random.rand()*(v if v != 0 else 1) 
    113  
    114 def randomize_model(pars, seed=None): 
    115     if seed is None: 
    116         seed = np.random.randint(1e9) 
    117     np.random.seed(seed) 
    118     # Note: the sort guarantees order of calls to random number generator 
    119     pars = dict((p,randomize(p,v)) for p,v in sorted(pars.items())) 
    120  
    121     return pars, seed 
    122  
    123 def constrain_pars(model_definition, pars): 
    124     """ 
    125     Restrict parameters to valid values. 
    126     """ 
    127     name = model_definition.name 
    128     if name == 'capped_cylinder' and pars['cap_radius'] < pars['radius']: 
    129         pars['radius'],pars['cap_radius'] = pars['cap_radius'],pars['radius'] 
    130     if name == 'barbell' and pars['bell_radius'] < pars['radius']: 
    131         pars['radius'],pars['bell_radius'] = pars['bell_radius'],pars['radius'] 
    132  
    133     # Limit guinier to an Rg such that Iq > 1e-30 (single precision cutoff) 
    134     if name == 'guinier': 
    135         #q_max = 0.2  # mid q maximum 
    136         q_max = 1.0  # high q maximum 
    137         rg_max = np.sqrt(90*np.log(10) + 3*np.log(pars['scale']))/q_max 
    138         pars['rg'] = min(pars['rg'],rg_max) 
    139  
    140 def parlist(pars): 
    141     return "\n".join("%s: %s"%(p,v) for p,v in sorted(pars.items())) 
    142  
    143 def suppress_pd(pars): 
    144     """ 
    145     Suppress theta_pd for now until the normalization is resolved. 
    146  
    147     May also suppress complete polydispersity of the model to test 
    148     models more quickly. 
    149     """ 
    150     pars = pars.copy() 
    151     for p in pars: 
    152         if p.endswith("_pd"): pars[p] = 0 
    153     return pars 
    154  
    155 def eval_sasview(model_definition, pars, data, Nevals=1): 
    156     # importing sas here so that the error message will be that sas failed to 
    157     # import rather than the more obscure smear_selection not imported error 
    158     import sas 
    159     from sas.models.qsmearing import smear_selection 
    160     model = sasview_model(model_definition, **pars) 
    161     smearer = smear_selection(data, model=model) 
    162     value = None  # silence the linter 
    163     toc = tic() 
    164     for _ in range(max(Nevals, 1)):  # make sure there is at least one eval 
    165         if hasattr(data, 'qx_data'): 
    166             q = np.sqrt(data.qx_data**2 + data.qy_data**2) 
    167             index = ((~data.mask) & (~np.isnan(data.data)) 
    168                      & (q >= data.qmin) & (q <= data.qmax)) 
    169             if smearer is not None: 
    170                 smearer.model = model  # because smear_selection has a bug 
    171                 smearer.accuracy = data.accuracy 
    172                 smearer.set_index(index) 
    173                 value = smearer.get_value() 
     168            theory = lambda: model.evalDistribution([data.qx_data[index], data.qy_data[index]]) 
     169    elif smearer is not None: 
     170        theory = lambda: smearer(model.evalDistribution(data.x)) 
     171    else: 
     172        theory = lambda: model.evalDistribution(data.x) 
     173 
     174    def calculator(**pars): 
     175        # paying for parameter conversion each time to keep life simple, if not fast 
     176        _, pars = revert_model(model_definition, pars) 
     177        for k,v in pars.items(): 
     178            parts = k.split('.')  # polydispersity components 
     179            if len(parts) == 2: 
     180                model.dispersion[parts[0]][parts[1]] = v 
    174181            else: 
    175                 value = model.evalDistribution([data.qx_data[index], data.qy_data[index]]) 
    176         else: 
    177             value = model.evalDistribution(data.x) 
    178             if smearer is not None: 
    179                 value = smearer(value) 
    180     average_time = toc()*1000./Nevals 
    181     return value, average_time 
    182  
    183 def eval_opencl(model_definition, pars, data, dtype='single', Nevals=1, 
    184                 cutoff=0., fast=False): 
     182                model.setParam(k, v) 
     183        return theory() 
     184 
     185    calculator.engine = "sasview" 
     186    return calculator 
     187 
     188DTYPE_MAP = { 
     189    'half': '16', 
     190    'fast': 'fast', 
     191    'single': '32', 
     192    'double': '64', 
     193    'quad': '128', 
     194    'f16': '16', 
     195    'f32': '32', 
     196    'f64': '64', 
     197    'longdouble': '128', 
     198} 
     199def eval_opencl(model_definition, data, dtype='single', cutoff=0.): 
    185200    try: 
    186         model = core.load_model(model_definition, dtype=dtype, 
    187                                 platform="ocl", fast=fast) 
     201        model = core.load_model(model_definition, dtype=dtype, platform="ocl") 
    188202    except Exception as exc: 
    189203        print(exc) 
    190204        print("... trying again with single precision") 
    191         model = core.load_model(model_definition, dtype='single', 
    192                                 platform="ocl", fast=fast) 
     205        dtype = 'single' 
     206        model = core.load_model(model_definition, dtype=dtype, platform="ocl") 
    193207    calculator = DirectModel(data, model, cutoff=cutoff) 
    194     # Run the calculator once before starting the timing loop 
     208    calculator.engine = "OCL%s"%DTYPE_MAP[dtype] 
     209    return calculator 
     210 
     211def eval_ctypes(model_definition, data, dtype='double', cutoff=0.): 
     212    if dtype=='quad': 
     213        dtype = 'longdouble' 
     214    model = core.load_model(model_definition, dtype=dtype, platform="dll") 
     215    calculator = DirectModel(data, model, cutoff=cutoff) 
     216    calculator.engine = "OMP%s"%DTYPE_MAP[dtype] 
     217    return calculator 
     218 
     219def time_calculation(calculator, pars, Nevals=1): 
     220    # initialize the code so time is more accurate 
    195221    value = calculator(**suppress_pd(pars)) 
    196     # Now run the timing loop 
    197222    toc = tic() 
    198     for _ in range(max(Nevals, 1)):  # force at least one eval 
     223    for _ in range(max(Nevals, 1)):  # make sure there is at least one eval 
    199224        value = calculator(**pars) 
    200225    average_time = toc()*1000./Nevals 
    201226    return value, average_time 
    202227 
    203  
    204 def eval_ctypes(model_definition, pars, data, dtype='double', Nevals=1, cutoff=0.): 
    205     model = core.load_model(model_definition, dtype=dtype, platform="dll") 
    206     calculator = DirectModel(data, model, cutoff=cutoff) 
    207     # Run the calculator once before starting the timing loop 
    208     value = calculator(**suppress_pd(pars)) 
    209     # Now run the timing loop 
    210     toc = tic() 
    211     for _ in range(max(Nevals, 1)):  # force at least one eval 
    212         value = calculator(**pars) 
    213     average_time = toc()*1000./Nevals 
    214     return value, average_time 
    215  
    216  
    217 def make_data(qmax, is2D, Nq=128, resolution=0.0, accuracy='Low', view='log'): 
    218     if is2D: 
    219         data = empty_data2D(np.linspace(-qmax, qmax, Nq), resolution=resolution) 
    220         data.accuracy = accuracy 
     228def make_data(opts): 
     229    qmax, nq, res = opts['qmax'], opts['nq'], opts['res'] 
     230    if opts['is2d']: 
     231        data = empty_data2D(np.linspace(-qmax, qmax, nq), resolution=res) 
     232        data.accuracy = opts['accuracy'] 
    221233        set_beam_stop(data, 0.004) 
    222234        index = ~data.mask 
    223235    else: 
    224         if view == 'log': 
     236        if opts['view'] == 'log': 
    225237            qmax = math.log10(qmax) 
    226             q = np.logspace(qmax-3, qmax, Nq) 
     238            q = np.logspace(qmax-3, qmax, nq) 
    227239        else: 
    228             q = np.linspace(0.001*qmax, qmax, Nq) 
    229         data = empty_data1D(q, resolution=resolution) 
     240            q = np.linspace(0.001*qmax, qmax, nq) 
     241        data = empty_data1D(q, resolution=res) 
    230242        index = slice(None, None) 
    231243    return data, index 
    232244 
    233 def compare(name, pars, Ncomp, Nbase, opts, set_pars): 
    234     model_definition = core.load_model_definition(name) 
    235  
    236     view = ('linear' if '-linear' in opts 
    237             else 'log' if '-log' in opts 
    238             else 'q4' if '-q4' in opts 
    239             else 'log') 
    240  
    241     opt_values = dict(split 
    242                       for s in opts for split in ((s.split('='),)) 
    243                       if len(split) == 2) 
    244     # Sort out data 
    245     qmax = (10.0 if '-exq' in opts 
    246             else 1.0 if '-highq' in opts 
    247             else 0.2 if '-midq' in opts 
    248             else 0.05) 
    249     Nq = int(opt_values.get('-Nq', '128')) 
    250     res = float(opt_values.get('-res', '0')) 
    251     accuracy = opt_values.get('-accuracy', 'Low') 
    252     is2D = "-2d" in opts 
    253     data, index = make_data(qmax, is2D, Nq, res, accuracy, view=view) 
    254  
    255  
    256     # modelling accuracy is determined by dtype and cutoff 
    257     dtype = ('longdouble' if '-quad' in opts 
    258              else 'double' if '-double' in opts 
    259              else 'half' if '-half' in opts 
    260              else 'single') 
    261     cutoff = float(opt_values.get('-cutoff','1e-5')) 
    262     fast = "-fast" in opts and dtype is 'single' 
    263  
    264     # randomize parameters 
    265     #pars.update(set_pars)  # set value before random to control range 
    266     if '-random' in opts or '-random' in opt_values: 
    267         seed = int(opt_values['-random']) if '-random' in opt_values else None 
    268         pars, seed = randomize_model(pars, seed=seed) 
    269         print("Randomize using -random=%i"%seed) 
    270     pars.update(set_pars)  # set value after random to control value 
    271     constrain_pars(model_definition, pars) 
    272     constrain_new_to_old(model_definition, pars) 
    273  
    274     # parameter selection 
    275     if '-mono' in opts: 
    276         pars = suppress_pd(pars) 
    277     if '-pars' in opts: 
    278         print("pars "+str(parlist(pars))) 
     245def make_engine(model_definition, data, dtype, cutoff): 
     246    if dtype == 'sasview': 
     247        return eval_sasview(model_definition, data) 
     248    elif dtype.endswith('!'): 
     249        return eval_ctypes(model_definition, data, dtype=dtype[:-1], 
     250                           cutoff=cutoff) 
     251    else: 
     252        return eval_opencl(model_definition, data, dtype=dtype, 
     253                           cutoff=cutoff) 
     254 
     255def compare(opts): 
     256    Nbase, Ncomp = opts['N1'], opts['N2'] 
     257    pars = opts['pars'] 
     258    data = opts['data'] 
    279259 
    280260    # Base calculation 
    281     if 0: 
    282         from sasmodels.models import sphere as target 
    283         base_name = target.name 
    284         base, base_time = eval_ctypes(target, pars, data, 
    285                 dtype='longdouble', cutoff=0., Nevals=Ncomp) 
    286     elif Nbase > 0 and "-ctypes" in opts and "-sasview" in opts: 
     261    if Nbase > 0: 
     262        base = opts['engines'][0] 
    287263        try: 
    288             base, base_time = eval_sasview(model_definition, pars, data, Ncomp) 
    289             base_name = "sasview" 
    290             #print("base/sasview", (base-pars['background'])/(comp-pars['background'])) 
    291             print("sasview t=%.1f ms, intensity=%.0f"%(base_time, sum(base))) 
    292             #print("sasview",comp) 
     264            base_value, base_time = time_calculation(base, pars, Nbase) 
     265            print("%s t=%.1f ms, intensity=%.0f"%(base.engine, base_time, sum(base_value))) 
    293266        except ImportError: 
    294267            traceback.print_exc() 
    295268            Nbase = 0 
    296     elif Nbase > 0: 
    297         base, base_time = eval_opencl(model_definition, pars, data, 
    298                 dtype=dtype, cutoff=cutoff, Nevals=Nbase, fast=fast) 
    299         base_name = "ocl" 
    300         print("opencl t=%.1f ms, intensity=%.0f"%(base_time, sum(base))) 
    301         #print("base " + base) 
    302         #print(max(base), min(base)) 
    303269 
    304270    # Comparison calculation 
    305     if Ncomp > 0 and "-ctypes" in opts: 
    306         comp, comp_time = eval_ctypes(model_definition, pars, data, 
    307                 dtype=dtype, cutoff=cutoff, Nevals=Ncomp) 
    308         comp_name = "ctypes" 
    309         print("ctypes t=%.1f ms, intensity=%.0f"%(comp_time, sum(comp))) 
    310     elif Ncomp > 0: 
     271    if Ncomp > 0: 
     272        comp = opts['engines'][1] 
    311273        try: 
    312             comp, comp_time = eval_sasview(model_definition, pars, data, Ncomp) 
    313             comp_name = "sasview" 
    314             #print("base/sasview", (base-pars['background'])/(comp-pars['background'])) 
    315             print("sasview t=%.1f ms, intensity=%.0f"%(comp_time, sum(comp))) 
    316             #print("sasview",comp) 
     274            comp_value, comp_time = time_calculation(comp, pars, Ncomp) 
     275            print("%s t=%.1f ms, intensity=%.0f"%(comp.engine, comp_time, sum(comp_value))) 
    317276        except ImportError: 
    318277            traceback.print_exc() 
     
    322281    if Nbase > 0 and Ncomp > 0: 
    323282        #print("speedup %.2g"%(comp_time/base_time)) 
    324         #print("max |base/comp|", max(abs(base/comp)), "%.15g"%max(abs(base)), "%.15g"%max(abs(comp))) 
    325         #comp *= max(base/comp) 
    326         resid = (base - comp) 
    327         relerr = resid/comp 
    328         #bad = (relerr>1e-4) 
    329         #print(relerr[bad],comp[bad],base[bad],data.qx_data[bad],data.qy_data[bad]) 
    330         _print_stats("|%s-%s|"%(base_name,comp_name)+(" "*(3+len(comp_name))), resid) 
    331         _print_stats("|(%s-%s)/%s|"%(base_name,comp_name,comp_name), relerr) 
     283        #print("max |base/comp|", max(abs(base_value/comp_value)), "%.15g"%max(abs(base_value)), "%.15g"%max(abs(comp_value))) 
     284        #comp *= max(base_value/comp_value) 
     285        resid = (base_value - comp_value) 
     286        relerr = resid/comp_value 
     287        _print_stats("|%s - %s|"%(base.engine,comp.engine)+(" "*(3+len(comp.engine))), resid) 
     288        _print_stats("|(%s - %s) / %s|"%(base.engine,comp.engine,comp.engine), relerr) 
    332289 
    333290    # Plot if requested 
    334     if '-noplot' in opts: return 
     291    if not opts['plot'] and not opts['explore']: return 
     292    view = opts['view'] 
    335293    import matplotlib.pyplot as plt 
     294    if Nbase > 0: 
     295        if Ncomp > 0: plt.subplot(131) 
     296        plot_theory(data, base_value, view=view, plot_data=False) 
     297        plt.title("%s t=%.1f ms"%(base.engine, base_time)) 
     298        #cbar_title = "log I" 
    336299    if Ncomp > 0: 
    337         if Nbase > 0: plt.subplot(131) 
    338         plot_theory(data, comp, view=view, plot_data=False) 
    339         plt.title("%s t=%.1f ms"%(comp_name,comp_time)) 
    340         #cbar_title = "log I" 
    341     if Nbase > 0: 
    342         if Ncomp > 0: plt.subplot(132) 
    343         plot_theory(data, base, view=view, plot_data=False) 
    344         plt.title("%s t=%.1f ms"%(base_name,base_time)) 
     300        if Nbase > 0: plt.subplot(132) 
     301        plot_theory(data, comp_value, view=view, plot_data=False) 
     302        plt.title("%s t=%.1f ms"%(comp.engine,comp_time)) 
    345303        #cbar_title = "log I" 
    346304    if Ncomp > 0 and Nbase > 0: 
     
    363321        v[v==0] = 0.5*np.min(np.abs(v[v!=0])) 
    364322        plt.hist(np.log10(np.abs(v)), normed=1, bins=50); 
    365         plt.xlabel('log10(err), err = | F(q) single - F(q) double| / | F(q) double |'); 
     323        plt.xlabel('log10(err), err = |(%s - %s) / %s|'%(base.engine, comp.engine, comp.engine)); 
    366324        plt.ylabel('P(err)') 
    367         plt.title('Comparison of single and double precision models for %s'%name) 
    368  
    369     plt.show() 
     325        plt.title('Distribution of relative error between calculation engines') 
     326 
     327    if not opts['explore']: 
     328        plt.show() 
    370329 
    371330def _print_stats(label, err): 
     
    387346# 
    388347USAGE=""" 
    389 usage: compare.py model [Nopencl] [Nsasview] [options...] [key=val] 
     348usage: compare.py model N1 N2 [options...] [key=val] 
    390349 
    391350Compare the speed and value for a model between the SasView original and the 
    392 OpenCL rewrite. 
     351sasmodels rewrite. 
    393352 
    394353model is the name of the model to compare (see below). 
    395 Nopencl is the number of times to run the OpenCL model (default=5) 
    396 Nsasview is the number of times to run the Sasview model (default=1) 
     354N1 is the number of times to run sasmodels (default=1). 
     355N2 is the number times to run sasview (default=1). 
    397356 
    398357Options (* for default): 
    399358 
    400359    -plot*/-noplot plots or suppress the plot of the model 
    401     -half/-single*/-double/-quad/-fast sets the calculation precision 
    402360    -lowq*/-midq/-highq/-exq use q values up to 0.05, 0.2, 1.0, 10.0 
    403     -Nq=128 sets the number of Q points in the data set 
     361    -nq=128 sets the number of Q points in the data set 
    404362    -1d*/-2d computes 1d or 2d data 
    405363    -preset*/-random[=seed] preset or random parameters 
    406364    -mono/-poly* force monodisperse/polydisperse 
    407     -ctypes/-sasview* selects gpu:cpu, gpu:sasview, or sasview:cpu if both 
    408365    -cutoff=1e-5* cutoff value for including a point in polydispersity 
    409366    -pars/-nopars* prints the parameter set or not 
    410367    -abs/-rel* plot relative or absolute error 
    411     -linear/-log/-q4 intensity scaling 
     368    -linear/-log*/-q4 intensity scaling 
    412369    -hist/-nohist* plot histogram of relative error 
    413370    -res=0 sets the resolution width dQ/Q if calculating with resolution 
    414371    -accuracy=Low accuracy of the resolution calculation Low, Mid, High, Xhigh 
    415  
    416 Key=value pairs allow you to set specific values to any of the model 
    417 parameters. 
     372    -edit starts the parameter explorer 
     373 
     374Any two calculation engines can be selected for comparison: 
     375 
     376    -single/-double/-half/-fast sets an OpenCL calculation engine 
     377    -single!/-double!/-quad! sets an OpenMP calculation engine 
     378    -sasview sets the sasview calculation engine 
     379 
     380The default is -single -sasview 
     381 
     382Key=value pairs allow you to set specific values for the model parameters. 
    418383 
    419384Available models: 
     
    423388NAME_OPTIONS = set([ 
    424389    'plot', 'noplot', 
    425     'half', 'single', 'double', 'quad', 'fast', 
     390    'half', 'fast', 'single', 'double', 
     391    'single!', 'double!', 'quad!', 'sasview', 
    426392    'lowq', 'midq', 'highq', 'exq', 
    427393    '2d', '1d', 
    428394    'preset', 'random', 
    429395    'poly', 'mono', 
    430     'sasview', 'ctypes', 
    431396    'nopars', 'pars', 
    432397    'rel', 'abs', 
    433398    'linear', 'log', 'q4', 
    434399    'hist', 'nohist', 
     400    'edit', 
    435401    ]) 
    436402VALUE_OPTIONS = [ 
    437403    # Note: random is both a name option and a value option 
    438     'cutoff', 'random', 'Nq', 'res', 'accuracy', 
     404    'cutoff', 'random', 'nq', 'res', 'accuracy', 
    439405    ] 
    440406 
     
    453419def get_demo_pars(model_definition): 
    454420    info = generate.make_info(model_definition) 
     421    # Get the default values for the parameters 
    455422    pars = dict((p[0],p[2]) for p in info['parameters']) 
     423 
     424    # Fill in default values for the polydispersity parameters 
     425    for p in info['parameters']: 
     426        if p[4] in ('volume', 'orientation'): 
     427            pars[p[0]+'_pd'] = 0.0 
     428            pars[p[0]+'_pd_n'] = 0 
     429            pars[p[0]+'_pd_nsigma'] = 3.0 
     430            pars[p[0]+'_pd_type'] = "gaussian" 
     431 
     432    # Plug in values given in demo 
    456433    pars.update(info['demo']) 
    457434    return pars 
    458435 
    459 def main(): 
    460     opts = [arg for arg in sys.argv[1:] if arg.startswith('-')] 
    461     popts = [arg for arg in sys.argv[1:] if not arg.startswith('-') and '=' in arg] 
     436def parse_opts(): 
     437    flags = [arg for arg in sys.argv[1:] if arg.startswith('-')] 
     438    values = [arg for arg in sys.argv[1:] if not arg.startswith('-') and '=' in arg] 
    462439    args = [arg for arg in sys.argv[1:] if not arg.startswith('-') and '=' not in arg] 
    463440    models = "\n    ".join("%-15s"%v for v in MODELS) 
     
    472449        print("expected parameters: model Nopencl Nsasview") 
    473450 
    474     invalid = [o[1:] for o in opts 
     451    invalid = [o[1:] for o in flags 
    475452               if o[1:] not in NAME_OPTIONS 
    476453                  and not any(o.startswith('-%s='%t) for t in VALUE_OPTIONS)] 
     
    479456        sys.exit(1) 
    480457 
     458 
     459    # Interpret the flags 
     460    opts = { 
     461        'plot'      : True, 
     462        'view'      : 'log', 
     463        'is2d'      : False, 
     464        'qmax'      : 0.05, 
     465        'nq'        : 128, 
     466        'res'       : 0.0, 
     467        'accuracy'  : 'Low', 
     468        'cutoff'    : 1e-5, 
     469        'seed'      : -1,  # default to preset 
     470        'mono'      : False, 
     471        'show_pars' : False, 
     472        'show_hist' : False, 
     473        'rel_err'   : True, 
     474        'explore'   : False, 
     475    } 
     476    engines = [] 
     477    for arg in flags: 
     478        if arg == '-noplot':    opts['plot'] = False 
     479        elif arg == '-plot':    opts['plot'] = True 
     480        elif arg == '-linear':  opts['view'] = 'linear' 
     481        elif arg == '-log':     opts['view'] = 'log' 
     482        elif arg == '-q4':      opts['view'] = 'q4' 
     483        elif arg == '-1d':      opts['is2d'] = False 
     484        elif arg == '-2d':      opts['is2d'] = True 
     485        elif arg == '-exq':     opts['qmax'] = 10.0 
     486        elif arg == '-highq':   opts['qmax'] = 1.0 
     487        elif arg == '-midq':    opts['qmax'] = 0.2 
     488        elif arg == '-loq':     opts['qmax'] = 0.05 
     489        elif arg.startswith('-nq='):       opts['nq'] = int(arg[4:]) 
     490        elif arg.startswith('-res='):      opts['res'] = float(arg[5:]) 
     491        elif arg.startswith('-accuracy='): opts['accuracy'] = arg[10:] 
     492        elif arg.startswith('-cutoff='):   opts['cutoff'] = float(arg[8:]) 
     493        elif arg.startswith('-random='):   opts['seed'] = int(arg[8:]) 
     494        elif arg == '-random':  opts['seed'] = np.random.randint(1e6) 
     495        elif arg == '-preset':  opts['seed'] = -1 
     496        elif arg == '-mono':    opts['mono'] = True 
     497        elif arg == '-poly':    opts['mono'] = False 
     498        elif arg == '-pars':    opts['show_pars'] = True 
     499        elif arg == '-nopars':  opts['show_pars'] = False 
     500        elif arg == '-hist':    opts['show_hist'] = True 
     501        elif arg == '-nohist':  opts['show_hist'] = False 
     502        elif arg == '-rel':     opts['rel_err'] = True 
     503        elif arg == '-abs':     opts['rel_err'] = False 
     504        elif arg == '-half':    engines.append(arg[1:]) 
     505        elif arg == '-fast':    engines.append(arg[1:]) 
     506        elif arg == '-single':  engines.append(arg[1:]) 
     507        elif arg == '-double':  engines.append(arg[1:]) 
     508        elif arg == '-single!': engines.append(arg[1:]) 
     509        elif arg == '-double!': engines.append(arg[1:]) 
     510        elif arg == '-quad!':   engines.append(arg[1:]) 
     511        elif arg == '-sasview': engines.append(arg[1:]) 
     512        elif arg == '-edit':    opts['explore'] = True 
     513 
     514    if len(engines) == 0: 
     515        engines.extend(['single','sasview']) 
     516    elif len(engines) == 1: 
     517        if engines[0][0] != 'sasview': 
     518            engines.append('sasview') 
     519        else: 
     520            engines.append('single') 
     521    elif len(engines) > 2: 
     522        del engines[2:] 
     523 
     524    name = args[0] 
     525    model_definition = core.load_model_definition(name) 
     526 
     527    N1 = int(args[1]) if len(args) > 1 else 1 
     528    N2 = int(args[2]) if len(args) > 2 else 1 
     529 
    481530    # Get demo parameters from model definition, or use default parameters 
    482531    # if model does not define demo parameters 
    483     name = args[0] 
    484     model_definition = core.load_model_definition(name) 
    485532    pars = get_demo_pars(model_definition) 
    486533 
    487     Ncomp = int(args[1]) if len(args) > 1 else 5 
    488     Nbase = int(args[2]) if len(args) > 2 else 1 
    489  
    490     # Fill in default polydispersity parameters 
    491     pds = set(p.split('_pd')[0] for p in pars if p.endswith('_pd')) 
    492     for p in pds: 
    493         if p+"_pd_nsigma" not in pars: pars[p+"_pd_nsigma"] = 3 
    494         if p+"_pd_type" not in pars: pars[p+"_pd_type"] = "gaussian" 
    495  
    496534    # Fill in parameters given on the command line 
    497     set_pars = {} 
    498     for arg in popts: 
     535    presets = {} 
     536    for arg in values: 
    499537        k,v = arg.split('=',1) 
    500538        if k not in pars: 
    501             # extract base name without distribution 
     539            # extract base name without polydispersity info 
    502540            s = set(p.split('_pd')[0] for p in pars) 
    503541            print("%r invalid; parameters are: %s"%(k,", ".join(sorted(s)))) 
    504542            sys.exit(1) 
    505         set_pars[k] = float(v) if not v.endswith('type') else v 
    506  
    507     compare(name, pars, Ncomp, Nbase, opts, set_pars) 
     543        presets[k] = float(v) if not k.endswith('type') else v 
     544 
     545    # randomize parameters 
     546    #pars.update(set_pars)  # set value before random to control range 
     547    if opts['seed'] > -1: 
     548        pars = randomize_pars(pars, seed=opts['seed']) 
     549        print("Randomize using -random=%i"%opts['seed']) 
     550    pars.update(presets)  # set value after random to control value 
     551    constrain_pars(model_definition, pars) 
     552    constrain_new_to_old(model_definition, pars) 
     553    if opts['mono']: 
     554        pars = suppress_pd(pars) 
     555    if opts['show_pars']: 
     556        print("pars " + str(parlist(pars))) 
     557 
     558    # Create the computational engines 
     559    data, _index = make_data(opts) 
     560    if N1: 
     561        base = make_engine(model_definition, data, engines[0], opts['cutoff']) 
     562    else: 
     563        base = None 
     564    if N2: 
     565        comp = make_engine(model_definition, data, engines[1], opts['cutoff']) 
     566    else: 
     567        comp = None 
     568 
     569    # Remember it all 
     570    opts.update({ 
     571        'name'      : name, 
     572        'def'       : model_definition, 
     573        'N1'        : N1, 
     574        'N2'        : N2, 
     575        'presets'   : presets, 
     576        'pars'      : pars, 
     577        'data'      : data, 
     578        'engines'   : [base, comp], 
     579    }) 
     580 
     581    return opts 
     582 
     583def main(): 
     584    opts = parse_opts() 
     585    if opts['explore']: 
     586        explore(opts) 
     587    else: 
     588        compare(opts) 
     589 
     590def explore(opts): 
     591    import wx 
     592    from bumps.names import FitProblem 
     593    from bumps.gui.app_frame import AppFrame 
     594 
     595    problem = FitProblem(Explore(opts)) 
     596    isMac = "cocoa" in wx.version() 
     597    app = wx.App() 
     598    frame = AppFrame(parent=None, title="explore") 
     599    if not isMac: frame.Show() 
     600    frame.panel.set_model(model=problem) 
     601    frame.panel.Layout() 
     602    frame.panel.aui.Split(0, wx.TOP) 
     603    if isMac: frame.Show() 
     604    app.MainLoop() 
     605 
     606class Explore(object): 
     607    """ 
     608    Return a bumps wrapper for a SAS model comparison. 
     609    """ 
     610    def __init__(self, opts): 
     611        from bumps.cli import config_matplotlib 
     612        import bumps_model 
     613        config_matplotlib() 
     614        self.opts = opts 
     615        info = generate.make_info(opts['def']) 
     616        pars, pd_types = bumps_model.create_parameters(info, **opts['pars']) 
     617        if not opts['is2d']: 
     618            active = [base + ext 
     619                      for base in info['partype']['pd-1d'] 
     620                      for ext in ['','_pd','_pd_n','_pd_nsigma']] 
     621            active.extend(info['partype']['fixed-1d']) 
     622            for k in active: 
     623                v = pars[k] 
     624                v.range(*parameter_range(k, v.value)) 
     625        else: 
     626            for k, v in self.pars.items(): 
     627                v.range(*parameter_range(k, v.value)) 
     628 
     629        self.pars = pars 
     630        self.pd_types = pd_types 
     631 
     632    def numpoints(self): 
     633        """ 
     634        Return the number of points 
     635        """ 
     636        return len(self.pars) + 1  # so dof is 1 
     637 
     638    def parameters(self): 
     639        """ 
     640        Return a dictionary of parameters 
     641        """ 
     642        return self.pars 
     643 
     644    def nllf(self): 
     645        return 0.  # No nllf 
     646 
     647    def plot(self, view='log'): 
     648        """ 
     649        Plot the data and residuals. 
     650        """ 
     651        pars = dict((k, v.value) for k,v in self.pars.items()) 
     652        pars.update(self.pd_types) 
     653        self.opts['pars'] = pars 
     654        compare(self.opts) 
     655 
    508656 
    509657if __name__ == "__main__": 
  • sasmodels/compare_many.py

    rf4f3919 rec7e360  
    77from . import core 
    88from .kernelcl import environment 
    9 from .compare import (MODELS, randomize_model, suppress_pd, eval_sasview, 
     9from .compare import (MODELS, randomize_pars, suppress_pd, eval_sasview, 
    1010                      eval_opencl, eval_ctypes, make_data, get_demo_pars, 
    11                       columnize, constrain_pars, constrain_new_to_old) 
     11                      columnize, constrain_pars, constrain_new_to_old, 
     12                      make_engine) 
    1213 
    1314def calc_stats(target, value, index): 
     
    3435    print(','.join('"%s"'%c for c in columns)) 
    3536 
     37PRECISION = { 
     38    'fast': 1e-3, 
     39    'half': 1e-3, 
     40    'single': 5e-5, 
     41    'double': 5e-14, 
     42    'single!': 5e-5, 
     43    'double!': 5e-14, 
     44    'quad!': 5e-18, 
     45    'sasview': 5e-14, 
     46} 
    3647def compare_instance(name, data, index, N=1, mono=True, cutoff=1e-5, 
    37                      precision='double'): 
     48                     base='sasview', comp='double'): 
    3849    model_definition = core.load_model_definition(name) 
    3950    pars = get_demo_pars(model_definition) 
     
    4758    # to allow them to update values in the current scope since nonlocal 
    4859    # declarations are not available in python 2.7. 
    49     def try_model(fn, *args, **kw): 
     60    def try_model(fn, pars): 
    5061        try: 
    51             result, _ = fn(model_definition, pars_i, data, *args, **kw) 
     62            result = fn(**pars) 
    5263        except KeyboardInterrupt: 
    5364            raise 
     
    6071                result = np.NaN*data.x 
    6172        return result 
    62     def check_model(label, target, value, acceptable): 
    63         stats = calc_stats(target, value, index) 
    64         columns.extend(stats) 
    65         labels.append('GPU single') 
     73    def check_model(pars): 
     74        base_value = try_model(calc_base, pars) 
     75        comp_value = try_model(calc_comp, pars) 
     76        stats = calc_stats(base_value, comp_value, index) 
    6677        max_diff[0] = max(max_diff[0], stats[0]) 
    67         good[0] = good[0] and (stats[0] < acceptable) 
     78        good[0] = good[0] and (stats[0] < expected) 
     79        return list(stats) 
     80 
     81 
     82    calc_base = make_engine(model_definition, data, base, cutoff) 
     83    calc_comp = make_engine(model_definition, data, comp, cutoff) 
     84    expected = max(PRECISION[base], PRECISION[comp]) 
    6885 
    6986    num_good = 0 
     
    7289    for k in range(N): 
    7390        print("%s %d"%(name, k)) 
    74         pars_i, seed = randomize_model(pars) 
     91        seed = np.random.randint(1e6) 
     92        pars_i = randomize_pars(pars, seed) 
    7593        constrain_pars(model_definition, pars_i) 
    7694        constrain_new_to_old(model_definition, pars_i) 
     
    7997 
    8098        good = [True] 
    81         labels = [] 
    82         columns = [] 
    83         target = try_model(eval_sasview) 
    84         #target = try_model(eval_ctypes, dtype='double', cutoff=0.) 
    85         #target = try_model(eval_ctypes, dtype='longdouble', cutoff=0.) 
    86         if precision == 'single': 
    87             value = try_model(eval_opencl, dtype='single', cutoff=cutoff) 
    88             check_model('GPU single', target, value, 5e-5) 
    89             single_value = value  # remember for single/double comparison 
    90         elif precision == 'double': 
    91             if environment().has_type('double'): 
    92                 label = 'GPU double' 
    93                 value = try_model(eval_opencl, dtype='double', cutoff=cutoff) 
    94             else: 
    95                 label = 'CPU double' 
    96                 value = try_model(eval_ctypes, dtype='double', cutoff=cutoff) 
    97             check_model(label, target, value, 5e-14) 
    98             double_value = value  # remember for single/double comparison 
    99         elif precision == 'quad': 
    100             value = try_model(eval_opencl, dtype='longdouble', cutoff=cutoff) 
    101             check_model('CPU quad', target, value, 5e-14) 
    102         if 0: 
    103             check_model('single/double', double_value, single_value, 5e-5) 
    104  
     99        columns = check_model(pars_i) 
    105100        columns += [v for _,v in sorted(pars_i.items())] 
    106101        if first: 
     102            labels = [" vs. ".join((calc_base.engine, calc_comp.engine))] 
    107103            print_column_headers(pars_i, labels) 
    108104            first = False 
     
    110106            num_good += 1 
    111107        else: 
    112             print(("%d,"%seed)+','.join("%g"%v for v in columns)) 
     108            print(("%d,"%seed)+','.join("%s"%v for v in columns)) 
    113109    print('"good","%d/%d","max diff",%g'%(num_good, N, max_diff[0])) 
    114110 
     
    144140is set in compare.py defaults for each model. 
    145141 
    146 PRECISION is the floating point precision to use for comparisons. 
     142PRECISION is the floating point precision to use for comparisons.  If two 
     143precisions are given, then compare one to the other, ignoring sasview. 
    147144 
    148145Available models: 
     
    151148 
    152149def main(): 
    153     if len(sys.argv) != 6: 
     150    if len(sys.argv) not in (6,7): 
    154151        print_help() 
    155152        sys.exit(1) 
     
    167164        mono = sys.argv[4] == 'mono' 
    168165        cutoff = float(sys.argv[4]) if not mono else 0 
    169         precision = sys.argv[5] 
     166        base = sys.argv[5] 
     167        comp = sys.argv[6] if len(sys.argv) > 6 else "sasview" 
    170168    except: 
    171169        traceback.print_exc() 
     
    173171        sys.exit(1) 
    174172 
    175     data, index = make_data(qmax=1.0, is2D=is2D, Nq=Nq) 
     173    data, index = make_data({'qmax':1.0, 'is2d':is2D, 'nq':Nq, 'res':0., 
     174                              'accuracy': 'Low', 'view':'log'}) 
    176175    model_list = [model] if model != "all" else MODELS 
    177176    for model in model_list: 
    178177        compare_instance(model, data, index, N=count, mono=mono, 
    179                          cutoff=cutoff, precision=precision) 
     178                         cutoff=cutoff, base=base, comp=comp) 
    180179 
    181180if __name__ == "__main__": 
Note: See TracChangeset for help on using the changeset viewer.