Changeset ac98886 in sasmodels


Ignore:
Timestamp:
Aug 14, 2016 11:15:15 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:
725ee36
Parents:
bde38b5
Message:

restore support for mixture models

File:
1 edited

Legend:

Unmodified
Added
Removed
  • sasmodels/mixture.py

    r7ae2b7f rac98886  
    1111*ProductModel(P, S)*. 
    1212""" 
     13from __future__ import print_function 
     14 
    1315from copy import copy 
    1416import numpy as np  # type: ignore 
     
    1618from .modelinfo import Parameter, ParameterTable, ModelInfo 
    1719from .kernel import KernelModel, Kernel 
     20from .details import make_details 
    1821 
    1922try: 
    2023    from typing import List 
    21     from .details import CallDetails 
    2224except ImportError: 
    2325    pass 
     
    4345        # to support vector parameters 
    4446        prefix = chr(ord('A')+k) + '_' 
    45         combined_pars.append(Parameter(prefix+'scale')) 
     47        scale =  Parameter(prefix+'scale', default=1.0, 
     48                           description="model intensity for " + part.name) 
     49        combined_pars.append(scale) 
    4650        for p in part.parameters.kernel_parameters: 
    4751            p = copy(p) 
     
    5155                p.length_control = prefix + p.length_control 
    5256            combined_pars.append(p) 
     57    #print("pars",combined_pars) 
    5358    parameters = ParameterTable(combined_pars) 
     59    parameters.max_pd = sum(part.parameters.max_pd for part in parts) 
    5460 
    5561    model_info = ModelInfo() 
     
    7076    # Remember the component info blocks so we can build the model 
    7177    model_info.composition = ('mixture', parts) 
     78    model_info.demo = {} 
     79    return model_info 
    7280 
    7381 
     
    7886        self.parts = parts 
    7987 
    80     def __call__(self, q_vectors): 
     88    def make_kernel(self, q_vectors): 
    8189        # type: (List[np.ndarray]) -> MixtureKernel 
    8290        # Note: may be sending the q_vectors to the n times even though they 
     
    104112        self.info =  model_info 
    105113        self.kernels = kernels 
    106  
    107     def __call__(self, call_details, value, weight, cutoff): 
     114        self.dtype = self.kernels[0].dtype 
     115 
     116    def __call__(self, call_details, values, cutoff, magnetic): 
    108117        # type: (CallDetails, np.ndarray, np.ndarry, float) -> np.ndarray 
    109         scale, background = value[0:2] 
     118        scale, background = values[0:2] 
    110119        total = 0.0 
    111120        # remember the parts for plotting later 
    112121        self.results = [] 
    113         for kernel, kernel_details in zip(self.kernels, call_details.parts): 
    114             part_result = kernel(kernel_details, value, weight, cutoff) 
    115             total += part_result 
    116             self.results.append(part_result) 
     122        offset = 2 # skip scale & background 
     123        parts = MixtureParts(self.info, self.kernels, call_details, values) 
     124        for kernel, kernel_details, kernel_values in parts: 
     125            #print("calling kernel", kernel.info.name) 
     126            result = kernel(kernel_details, kernel_values, cutoff, magnetic) 
     127            #print(kernel.info.name, result) 
     128            total += result 
     129            self.results.append(result) 
    117130 
    118131        return scale*total + background 
     
    123136            k.release() 
    124137 
     138 
     139class MixtureParts(object): 
     140    def __init__(self, model_info, kernels, call_details, values): 
     141        self.model_info = model_info 
     142        self.parts = model_info.composition[1] 
     143        self.kernels = kernels 
     144        self.call_details = call_details 
     145        self.values = values 
     146        self.spin_index = model_info.parameters.npars + 2 
     147        #call_details.show(values) 
     148 
     149    def __iter__(self): 
     150        # type: () -> PartIterable 
     151        self.part_num = 0 
     152        self.par_index = 2 
     153        self.mag_index = self.spin_index + 3 
     154        return self 
     155 
     156    def next(self): 
     157        # type: () -> Tuple[List[Callable], CallDetails, np.ndarray] 
     158        if self.part_num >= len(self.parts): 
     159            raise StopIteration() 
     160        info = self.parts[self.part_num] 
     161        kernel = self.kernels[self.part_num] 
     162        call_details = self._part_details(info, self.par_index) 
     163        values = self._part_values(info, self.par_index, self.mag_index) 
     164        values = values.astype(kernel.dtype) 
     165        #call_details.show(values) 
     166 
     167        self.part_num += 1 
     168        self.par_index += info.parameters.npars + 1 
     169        self.mag_index += 3 * len(info.parameters.magnetism_index) 
     170 
     171        return kernel, call_details, values 
     172 
     173    def _part_details(self, info, par_index): 
     174        # type: (ModelInfo, int) -> CallDetails 
     175        full = self.call_details 
     176        # par_index is index into values array of the current parameter, 
     177        # which includes the initial scale and background parameters. 
     178        # We want the index into the weight length/offset for each parameter. 
     179        # Exclude the initial scale and background, so subtract two, but each 
     180        # component has its own scale factor which we need to skip when 
     181        # constructing the details for the kernel, so add one, giving a 
     182        # net subtract one. 
     183        index = slice(par_index - 1, par_index - 1 + info.parameters.npars) 
     184        length = full.length[index] 
     185        offset = full.offset[index] 
     186        # The complete weight vector is being sent to each part so that 
     187        # offsets don't need to be adjusted. 
     188        part = make_details(info, length, offset, full.num_weights) 
     189        return part 
     190 
     191    def _part_values(self, info, par_index, mag_index): 
     192        # type: (ModelInfo, int, int) -> np.ndarray 
     193        #print(info.name, par_index, self.values[par_index:par_index + info.parameters.npars + 1]) 
     194        scale = self.values[par_index] 
     195        pars = self.values[par_index + 1:par_index + info.parameters.npars + 1] 
     196        nmagnetic = len(info.parameters.magnetism_index) 
     197        if nmagnetic: 
     198            spin_state = self.values[self.spin_index:self.spin_index + 3] 
     199            mag_index = self.values[mag_index:mag_index + 3 * nmagnetic] 
     200        else: 
     201            spin_state = [] 
     202            mag_index = [] 
     203        nvalues = self.model_info.parameters.nvalues 
     204        nweights = self.call_details.num_weights 
     205        weights = self.values[nvalues:nvalues+2*nweights] 
     206        zero = self.values.dtype.type(0.) 
     207        values = [[scale, zero], pars, spin_state, mag_index, weights] 
     208        # Pad value array to a 32 value boundary 
     209        spacer = (32 - sum(len(v) for v in values)%32)%32 
     210        values.append([zero]*spacer) 
     211        values = np.hstack(values).astype(self.kernels[0].dtype) 
     212        return values 
Note: See TracChangeset for help on using the changeset viewer.