Changeset 161

Show
Ignore:
Timestamp:
04/28/08 21:37:01 (7 months ago)
Author:
edsuom
Message:

Working on prior and parameter containers, starting to tie everything together

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

  • projects/AsynCluster/trunk/svpmc/model.py

    r160 r161  
    9292        """ 
    9393        Returns a list of L{tseries.TimeSeries} objects loaded and transformed 
    94         as defined in my specification. 
     94        as defined in my specification along with a dict of descriptions for 
     95        the time series, keyed by name. 
    9596        """ 
    9697        tsList = [] 
     98        descriptionDict = {} 
    9799        for name, filePath, transform, description in self.tables['series']: 
    98100            filePath = os.path.join(self.specDir, filePath) + ".dat" 
    …… 
    101103                getattr(tsObject, methodName)() 
    102104            tsList.append(tsObject) 
    103         return tsList 
     105            descriptionDict[name] = description 
     106        return tsList, descriptionDict 
    104107 
    105108    def priorContainer(self): 
    …… 
    107110        Returns an instance of L{params.PriorContainer} with 
    108111        L{params.FlexArray} objects loaded with a L{params.Prior} object for 
    109         each array element of each of my parameters. 
    110         """ 
     112        each array element of each of my parameters, along with a dict of 
     113        descriptions for the parameters, keyed by name. 
     114        """ 
     115        descriptionDict = {} 
    111116        p = len(self.tables['series']) 
    112117        container = params.PriorContainer() 
    …… 
    138143                                    continue 
    139144                            pa[j,k] = params.Prior(**kw) 
    140         return container 
    141              
     145            descriptionDict[name] = description 
     146        return container, descriptionDict 
     147 
    142148 
    143149class ModelManager(object): 
    …… 
    224230                # An error from the remote likelihood method call is treated as 
    225231                # infinitely low likelihood 
    226                 paramContainer.pLogLikelihood = -s.inf 
     232                paramContainer.Lx = -s.inf 
    227233            up = pack.Unpacker(result) 
    228234            L = paramContainer.L = up() 
    …… 
    427433        return L0, L, h[:,1:] 
    428434 
    429     def hybridGibbs(self, z, x, sigma, rv, tol=1E-4): 
     435    def hybridGibbs(self, z, x, rv, sigma, tol=1E-4): 
    430436        """ 
    431437        Does one iteration of a hybrid Gibbs sampler for the log-volatilities, 

  • projects/AsynCluster/trunk/svpmc/params.py

    r160 r161  
    6868    array. 
    6969 
     70    You can add two like-dimensioned instances of me. Each element of the sum 
     71    will be the sum of the respective objects of the added L{FlexArray} 
     72    instances. 
     73 
     74    You can call the L{call} method on an instance of me with a method name and 
     75    zero or more arguments. The result will be a numeric array of the same 
     76    shape as the instance, each element of which will have the numeric result 
     77    of a call of the named method on each of the corresponding objects in the 
     78    instance. Any argument that is an array (or L{FlexArray} instance) of the 
     79    same shape as my instance will have each the value of each element 
     80    corresponding to the called object's element extracted and used as the 
     81    method argument instead of the original array argument. 
     82 
     83    You can get and set an attribute of an instance's objects in array fashion 
     84    by simply setting or getting the named attribute of the instance. 
     85 
    7086    B{Caveat}: Unlike the behavior of SciPy arrays, Setting elements of a 
    71       subset or raveled version of an instance of me does I{not} set the 
    72       corresponding element of the original instance. 
    73      
     87    subset or raveled version of an instance of me does I{not} set the 
     88    corresponding element of the original instance. 
    7489    """ 
    7590    def __init__(self, *shape): 
    76         self.shape = shape 
    77         self.N = 1 
    78         for dim in shape: 
    79             self.N *= dim 
     91        self._shape = tuple([int(x) for x in shape]) 
     92        self._N = 1 
     93        for dim in self._shape: 
     94            self._N *= dim 
    8095        # Object list 
    81         self.O = [None] * self.N 
     96        self._O = [None] * self._N 
    8297        # Index array 
    83         self.I = s.arange(self.N).reshape(*shape) 
     98        self._I = s.arange(self._N).reshape(*self._shape) 
    8499 
    85100    def __call__(self, I): 
    86101        if isinstance(I, int): 
    87             return self.O[I] 
     102            return self._O[I] 
    88103        newVersion = self.__class__(*I.shape) 
    89104        for j, k in enumerate(I.ravel()): 
    90             newVersion.O[j] = self.O[k] 
     105            newVersion._O[j] = self._O[k] 
    91106        return newVersion 
    92107 
    93108    def __len__(self): 
    94         return self.N 
    95  
     109        return self.shape[0] 
     110     
    96111    def __iter__(self): 
    97         self.k = -1 
     112        self._k = -1 
    98113        return self 
    99114 
    100115    def next(self): 
    101         self.k += 1 
    102         if self.k < self.N: 
    103             return self(self.I[self.k]) 
     116        self._k += 1 
     117        if self._k < self._N: 
     118            return self(self._I[self._k]) 
    104119        raise StopIteration 
    105120     
    106121    def __getitem__(self, key): 
    107         return self(self.I[key]) 
     122        return self(self._I[key]) 
    108123 
    109124    def __setitem__(self, key, value): 
    110         I = self.I[key] 
     125        I = self._I[key] 
    111126        if isinstance(I, int): 
    112             self.O[I] = value 
     127            self._O[I] = value 
    113128        elif I.shape == s.shape(value): 
    114129            if not isinstance(value, (list, tuple)): 
    115130                value = value.ravel() 
    116131            for j, k in enumerate(I.ravel()): 
    117                 self.O[k] = value[j] 
     132                self._O[k] = value[j] 
    118133        else: 
    119134            raise SHAPE_MISMATCH 
    120135 
    121136    def ravel(self): 
    122         return self(self.I.ravel()) 
     137        return self(self._I.ravel()) 
    123138 
    124139    def reshape(self, *shape): 
    125         return self(self.I.reshape(*shape)) 
     140        return self(self._I.reshape(*shape)) 
    126141 
    127142    def __add__(self, other): 
    128143        if other.shape != self.shape: 
    129144            raise SHAPE_MISMATCH 
    130         newVersion = self(self.I) 
    131         for k, thisObj in enumerate(newVersion.O): 
    132             thisObj += other.O[k] 
     145        newVersion = self(self._I) 
     146        for k, thisObj in enumerate(newVersion._O): 
     147            thisObj += other._O[k] 
    133148        return newVersion 
    134149 
    135     def call(self, methodName, *args, **kw): 
     150    def __getattr__(self, name): 
     151        """ 
     152        For each of my objects, get the specified numeric attribute. 
     153 
     154        Returns the result (must be a numeric scalar) of each attribute value 
     155        in an array of the same shape as my object array. 
     156        """ 
     157        if name == 'shape': 
     158            return self._shape 
     159        x = s.empty(self.shape) 
     160        for k in self._I.ravel(): 
     161            thisValue = getattr(self._O[k], name) 
     162            if not isinstance(thisValue, (int, long, float)): 
     163                raise ValueError( 
     164                    "Attribute must have a numeric scalar value") 
     165            x.ravel()[k] = thisValue 
     166        return x 
     167 
     168    def __setattr__(self, name, value): 
     169        """ 
     170        For each of my objects, sets the specified numeric attribute to the 
     171        supplied value. 
     172 
     173        If the value is an array (or L{FlexArray}) with the same shape as me, 
     174        each object's attribute value will be set to the value of the 
     175        corresponding array element. 
     176        """ 
     177        if name.startswith('_'): 
     178            object.__setattr__(self, name, value) 
     179        else: 
     180            if s.shape(value) == self.shape: 
     181                values = value.ravel() 
     182            else: 
     183                values = [value] * self._N 
     184            for k, thisValue in enumerate(values): 
     185                setattr(self._O[k], name, thisValue) 
     186     
     187    def call(self, methodName, *args): 
    136188        """ 
    137189        For each of my objects, call the method specified by I{methodName} with 
    138         any further args and keywords supplied. 
    139  
    140         Returns the result (must be a numeric scalar) of each call in an array 
    141         of the same shape as my object array. 
    142         """ 
     190        any further args and keywords supplied. For any argument that is an 
     191        array (or L{FlexArray}) with the same shape as me, only the 
     192        corresponding array element will be supplied as the argument to each 
     193        object's method call. 
     194 
     195        Returns the result of each call in an array of the same shape as my 
     196        object array. The results must all be numeric scalars. 
     197 
     198        TODO: If the result of each call is not a numeric scalar, the result 
     199        will be another instance of me with the results (of whatever type) as 
     200        its objects. 
     201        """ 
     202        indicesOfArrayArgs = [ 
     203            k for k, arg in enumerate(args) if s.shape(arg) == self.shape] 
     204        theseArgs = [arg for arg in args] 
    143205        x = s.empty(self.shape) 
    144         for k in self.I.ravel(): 
    145             thisValue = getattr(self.O[k], methodName)(*args, **kw) 
     206        for j in self._I.ravel(): 
     207            for k in indicesOfArrayArgs: 
     208                theseArgs[k] = args[k].ravel()[j] 
     209            thisValue = getattr(self._O[j], methodName)(*theseArgs) 
    146210            if not isinstance(thisValue, (int, long, float)): 
    147211                raise ValueError( 
    148212                    "Called method must return a numeric scalar value") 
    149             x.ravel()[k] = thisValue 
     213            x.ravel()[j] = thisValue 
    150214        return x 
    151215 
    152216 
    153 class PriorContainer(object): 
    154     """ 
    155     I am a container for array-like instances of L{FlexArray} that contain 
    156     L{Prior} objects for the array elements of my model's parameters. 
    157     """ 
     217class ParameterContainer(Parameterized): 
     218    """ 
     219    I am a container for arrays of model parameters. 
     220 
     221    In addition to the parameter arrays, I house as key attributes: 
     222 
     223        - A latent parameter array I{z}. 
     224 
     225        - A scalar I{Lx} with the log-likelihood of the model's data given my 
     226          parameters and the latent parameter I{z}. 
     227 
     228        - A scalar I{Lp} with the prior log-probability density of the 
     229          parameters. 
     230 
     231        - A scalar I{Lj} with the log-probability density of the proposal jump 
     232          resulting in my parameters. 
     233 
     234    """ 
     235    keyAttrs = {'z':[], 'Lx':None, 'Lp':None, 'Lj':None} 
    158236    paramNames = PARAM_NAMES 
    159  
    160     def priorArray(self, paramName, *shape): 
    161         array = getattr(self, paramName, None) 
    162         if array is None: 
    163             array = FlexArray(*shape) 
    164             setattr(self, paramName, array) 
    165         return array 
    166237 
    167238 
    …… 
    218289    def pdf(self, x): 
    219290        """ 
    220         Returns probability values of the supplied parameter value I{x}. 
    221         """ 
    222         return self.distObj.pdf(x) 
    223  
    224     def rvs(self, N): 
    225         """ 
    226         Returns I{N} parameter values drawn from my distribution. 
    227         """ 
    228         return self.distObj.rvs(N) 
    229  
    230     def rds(self, N, wiggle): 
    231         """ 
    232         Returns I{N} values that deviate from zero by random amounts drawn from 
    233         a normal distribution with a pre-scaled standard deviation that is set 
    234         to approximate the 70% probability width of my distribution. (That is 
    235         the +/- 1 standard deviation width of a normal distribution.) 
    236  
    237         The deviates are then scaled by a I{wiggle} value. 
     291        Returns probability density of the supplied parameter value I{x}. 
     292        """ 
     293        return float(self.distObj.pdf(x)) 
     294 
     295    def rvs(self): 
     296        """ 
     297        Returns a parameter value drawn from my distribution. 
     298        """ 
     299        return self.distObj.rvs(1)[0] 
     300 
     301    def rJump(self, wiggle): 
     302        """ 
     303        Returns a value that jumps from zero by a random amount that, before 
     304        scaling by the supplied I{wiggle}, is drawn from a normal distribution 
     305        with a pre-scaled standard deviation that is set to approximate the 70% 
     306        probability width of my distribution. (That is the +/- 1 standard 
     307        deviation width of a normal distribution.) 
     308        """ 
     309        return self.rdsObj.rvs(1)[0] 
     310 
     311    def pJump(self, x, wiggle): 
     312        """ 
     313        Returns the probability density of the supplied jump value I{x}, given 
     314        that it was generated by a call to my L{rJump} method with the 
     315        specified I{wiggle}. 
     316 
     317        For my Gaussian-PDF deviate generator, the probability density of the 
     318        jump is inversely proportional to the scaling to its standard deviation 
     319        imparted by I{wiggle}. Thus, the unit-normal PDF is divided by the 
     320        I{wiggle} in the result. 
     321        """ 
     322        return self.rdsObj.pdf(x/wiggle) / wiggle 
     323 
     324 
     325class PriorContainer(object): 
     326    """ 
     327    I am a container for array-like instances of L{FlexArray} that contain 
     328    L{Prior} objects for the array elements of my model's parameters. 
     329    """ 
     330    paramNames = PARAM_NAMES 
     331 
     332    def priorArray(self, paramName, *shape): 
     333        array = getattr(self, paramName, None) 
     334        if array is None: 
     335            array = FlexArray(*shape) 
     336            setattr(self, paramName, array) 
     337        return array 
     338 
     339    def new(self): 
     340        """ 
     341        Returns a new parameter container with values intialized from a random 
     342        draw of my priors. 
     343        """ 
     344        Lp = 0 
     345        paramContainer = ParameterContainer() 
     346        for name in self.paramNames: 
     347            priorFlexArray = getattr(self, name) 
     348            paramArray = priorFlexArray.call('rvs') 
     349            Lp += s.log(priorFlexArray.call('pdf', paramArray)).sum() 
     350            setattr(paramContainer, name, paramArray) 
     351        paramContainer.Lp = Lp 
     352        return paramContainer 
     353 
     354    def proposal(self, paramContainer, wiggle): 
     355        """ 
     356        Returns a new instance of L{ParameterContainer} with a random-walk 
     357        proposal based on the supplied I{paramContainer}, my priors, and the 
     358        supplied I{wiggle}. 
     359 
     360        The returned parameter container object will have a copy of the 
     361        supplied instance's latent parameter array I{z} and new values of I{Lp} and 
     362        I{Lj} corresponding to the proposal. 
    238363        """ 
    239364        if wiggle < 0.0: 
    …… 
    241366                "Specify a positive wiggle value, '%s' obj invalid" \ 
    242367                % str(wiggle)) 
    243         return wiggle * self.rdsObj.rvs(N) 
    244  
    245     def pds(self, X, wiggle): 
    246         """ 
    247         Returns the probability densities of the supplied zero-mean deviates, 
    248         under the assumption that they were generated from my Gaussian-PDF 
    249         L{rds} method with the specified I{wiggle}. 
    250  
    251         For my Gaussian-PDF deviate generator, the probability density of a 
    252         given deviate is inversely proportional to the scaling to its standard 
    253         deviation imparted by I{wiggle}. Thus, the unit-normal PDF is divided 
    254         by the I{wiggle} in the result. 
    255         """ 
    256         return self.rdsObj.pdf(X/wiggle) / wiggle 
    257  
    258  
    259 class ParameterContainer(Parameterized): 
    260     """ 
    261     I am a container for arrays of model parameters. 
    262     """ 
    263     keyAttrs = {'z':[], 'L':None} 
    264     paramNames = PARAM_NAMES 
    265  
    266     def __add__(self, other): 
    267         newVersion = self.__class__(z=self.z) 
    268         newVersion.setParamSequence([ 
    269             getattr(self, name) + getattr(other, name) 
    270             for name in self.paramNames]) 
     368        if not isinstance(paramContainer, ParameterContainer): 
     369            raise TypeError( 
     370                "You must supply a ParameterContainer as the first "+\ 
     371                "argument, not a '%s'" % type(paramContainer)) 
     372        newVersion = ParameterContainer(z=s.array(paramContainer.z).copy()) 
     373        Lp, Lj = 0, 0 
     374        for name in self.paramNames: 
     375            priorFlexArray = getattr(self, name) 
     376            jumps = priorArray.call('rJump', wiggle) 
     377            paramArray = getattr(paramContainer, name) + jumps 
     378            Lp += s.log(priorFlexArray.call('pdf', paramArray)).sum() 
     379            Lj += s.log(priorFlexArray.call('pJump', jumps)).sum() 
     380            setattr(newVersion, name, paramArray) 
     381        newVersion.Lp = Lp 
     382        newVersion.Lj = Lj 
    271383        return newVersion 
    272384 
     385 
     386 

  • projects/AsynCluster/trunk/svpmc/pmc.py

    r151 r161  
    2626import scipy as s 
    2727from scipy import random 
    28 from zope.interface import implements 
    29 from twisted.internet import defer, interfaces 
     28from Scientific.IO import NetCDF 
     29from twisted.internet import defer 
    3030 
    3131import asynqueue 
    …… 
    3737 
    3838 
    39 class IterationProducer(object): 
     39class Recorder(object): 
    4040    """ 
    41     For each iteration of an instance of L{MCMC}, I produce a 2-D array of rows 
    42     that each contain the elements of one chain's vector of parameters. 
    43  
    44     @ivar consumer: A list of the consumers for whom I'm producing 
    45       information. 
    46      
     41    For each iteration of an instance of L{PMC}, I update an open NetCDF file 
     42    with the posterior draws represented by a flex-array of parameter 
     43    containers. 
    4744    """ 
    48     implements(interfaces.IPushProducer) 
    49      
    50     def __init__(self, mcmcObj): 
    51         self.mcmcObj = mcmcObj 
    52         self.producing = False 
    53         self.consumers = [] 
     45    def __init__(self, filePath): 
     46        self.cdf = NetCDF.NetCDFFile(filePath, 'w') 
     47         
    5448 
    5549    def restart(self): 
    …… 
    9387class PMC(object): 
    9488    """ 
    95     Population MCMC with per Cappe et al. 
     89    Population Monte Carlo with per Cappe et al. 
    9690    """ 
    97     wiggle = 1.0 
    9891    chunkSize = 5000 
    99     V = [1.0, 0.5, 0.1, 0.05, 0.001] 
    100  
     92    V = [1.0, 0.5, 0.1, 0.05, 0.01] 
     93     
    10194    def __init__(self, modelManager, **kw): 
    10295        self.mgr = modelManager 
    103         self.iterationProducer = IterationProducer(self) 
    104         self.consumers = [] 
     96        self.recorder = Recorder(self.mgr) 
    10597        for name, value in kw.iteritems(): 
    10698            setattr(self, name, value) 
    10799 
    108     def setupRun(self, N_chains, **kw): 
    109         """ 
    110         All subclasses setup their runs with certain common parameters, and 
    111         with a restarted iteration producer. 
    112         """ 
    113         self.N_chains = N_chains 
    114         N_iter = kw.get('N_iter', getattr(self, 'N_iter', 1000)) 
    115         N_bi = kw.get('N_bi', getattr(self, 'N_bi', N_iter/2)) 
    116         if 'wiggle' in kw: 
    117             self.wiggle = kw['wiggle'] 
     100    def setupRun(self, N_iter, N_members): 
     101        """ 
     102        Set up a new run with various parameters and a restarted iteration 
     103        producer. 
     104        """ 
     105        self.N_iter = N_iter 
     106        self.N_members = N_members 
    118107        self.iterationProducer.restart() 
    119         return N_iter, N_bi 
    120108 
    121109    def subscribe(self, consumer): 
    …… 
    148136        """ 
    149137        Returns a 3-tuple of info for jump proposals on the parameter 
    150         containers in the supplied parameter array I{X}, given a jump variance 
    151         I{v}. Proposals with zero probability of occuring, given the 
     138        containers in the supplied 1-D parameter array I{X}, given a jump 
     139        deviation I{v}. Proposals with zero probability of occuring, given the 
    152140        parameters' prior distributions, are censored from the result. (There's 
    153141        no point considering the likelihood of proposals that cannot occur.) 
    154142         
    155         The 3-tuple contains: 
    156  
    157             1. a parameter array of valid proposed parameter 
    158                sequences (if any, given the priors), 
    159  
    160             2. the log-probability density of the proposals given the 
    161                parameters' prior distributions, and 
    162  
    163             3. the log-probability density of the proposal jumps. 
     143        The 3-tuple contains like-dimensioned 1-D arrays, with the following 
     144        respective elements: 
     145 
     146            1. A valid proposed parameter container. 
     147             
     148            2. The log-probability density of the proposal given the 
     149               parameters' prior distributions. 
     150 
     151            3. The log-probability density of the proposal jumps. 
    164152         
    165153        """ 
    …… 
    222210     
    223211    @defer.deferredGenerator 
    224     def run(self, N_chains, **kw): 
    225         """ 
    226         Does a PMC run with I{N_chains} population members and jump variances 
     212    def run(self, N_members, **kw): 
     213        """ 
     214        Does a PMC run with I{N_members} population members and jump variances 
    227215        in the supplied 1-D array I{V}. The number of population members for 
    228216        each jump variance will go up or down, depending on the performance for 
    …… 
    245233        @param N_bi: The number of burn-in iterations to discard before 
    246234          producing any. 
    247  
    248         @param wiggle: A value between 0 and 1 that globally scales all jump 
    249           variances to smaller values from those supplied in I{V}. Default 
    250           value is 1.0 for no global scaling. 
    251235 
    252236        @param V: A sequence of variance values to use instead of the default. 
    …… 
    272256            # minimum size of two apiece. 
    273257            R = s.clip( 
    274                 s.round_(N_chains*R), rMin, N_chains-rMin*(P-1)).astype(int) 
    275             # Twiddle the biggest one as needed to keep sum = N_chains 
    276             R[s.argmax(R)] += N_chains - sum(R) 
     258                s.round_(N_members*R), rMin, N_members-rMin*(P-1)).astype(int) 
     259            # Twiddle the biggest one as needed to keep sum = N_members 
     260            R[s.argmax(R)] += N_members - sum(R) 
    277261            # Replace the old list and subset index 
    278262            self.R = R 
    279             self.Is = s.arange(N_chains) 
    280  
    281         N_iter, N_bi = self.setupRun(N_chains, **kw) 
     263            self.Is = s.arange(N_members) 
     264 
     265        N_iter, N_bi = self.setupRun(N_members, **kw) 
    282266        # The variance settings (default of 6) 
    283267        V = kw.pop('V', self.V) 
    284268        # Some constants 
    285269        P = len(V) 
    286         rMin = max([2, int(round(0.01*N_chains))]) 
     270        rMin = max([2, int(round(0.01*N_members))]) 
    287271        # Initialize some arrays 
    288         X = self.mgr.rPriors(N_chains) 
     272        X = self.mgr.rPriors(N_members) 
    289273        normalize(s.ones(P)) 
    290         V = self.wiggle * s.array(V) 
     274        V = s.array(V) 
    291275        # Iteration 
    292276        for i in xrange(N_iter + N_bi): 
    …… 
    311295            Wnz = [w for w in W if len(w)] 
    312296            if Wnz: 
    313                 I = resample(s.concatenate(Wnz), N_chains, logWeights=True) 
     297                I = resample(s.concatenate(Wnz), N_members, logWeights=True) 
    314298                if len(I): 
    315299                    X = s.row_stack([XP[k] for k in xrange(P) if len(W[k])])[I] 
    …… 
    318302            else: 
    319303                R = s.ones(P) 
    320             # Normalize R to maintain a total of N_chains population members 
     304            # Normalize R to maintain a total of N_members population members 
    321305            normalize(R) 
    322306            # Provide some info 
    …… 
    326310                ", ".join(["%9.3g" % x for x in X.mean(0)])) 
    327311            print msg 
    328         self.iterationProducer.stopProducing() 
     312        self.recorder.done() 

  • projects/AsynCluster/trunk/svpmc/test/project-spec.txt

    r160 r161  
    4141 
    4242 
     43# Specify priors of each parameter array by element(s) with one or more rows in 
     44# the parameter tables below. Subsequent lines in a given table override any 
     45# conflicting specifications of earlier rows of that table. 
     46#  
     47# 1-D array elements are specified by integer or ':' for all elements. 
     48# 
     49# 2-D array elements are specified by x,y integer pairs (no spaces!), ':' for 
     50# all elements, or 'I' for elements along the leading diagonal. 
     51 
     52 
    4353a           Distribution    Loc         Scale       a       b 
    4454------------------------------------------------------------------------------- 

  • projects/AsynCluster/trunk/svpmc/test/test_model.py

    r160 r161  
    5252 
    5353    def test_timeSeries(self): 
    54         tsList = self.sp.timeSeries() 
     54        tsList = self.sp.timeSeries()[0] 
    5555        self.failUnlessEqual(len(tsList), 3) 
    5656 
    5757    def test_priorContainer(self): 
    58         container = self.sp.priorContainer() 
    59         self.fail("TODO") 
     58        container = self.sp.priorContainer()[0] 
     59        self.failUnlessEqual(container.a.shape, (3,)) 
     60        self.failUnlessEqual(container.b.shape, (3,3)) 
     61        self.failUnlessEqual(container.cr.shape, (6,)) 
    6062 
    6163 
    …… 
    614616        z = s.randn(1, self.N) 
    615617        rv = self.rv(self.model.fs, self.model.fr) 
     618        print "RV\n", rv, "\n" 
    616619        for k in xrange(N): 
    617620            L_this = self.model.hybridGibbs(z, self.x, rv, sigma) 
    …… 
    619622            print "%03d: L=%6.1f" % (k, L_this) 
    620623        h = self.model.logVolatilities( 
    621             z, z.copy(), self.x, self.model.d.copy())[-1] 
     624            z, s.array(rv*z), self.x, self.model.d.copy())[-1] 
    622625        if VERBOSE: 
    623626            fig = self.fig 
    624             vectors = (self.x[0,:], self.h, h[0,:]) 
     627            vectors = (self.x[0,:], (self.h, h[0,:])) 
    625628            for k, vector in enumerate(vectors): 
    626629                spNumber = 100*(len(vectors)+1) + k + 11 
    627630                sp = fig.add_subplot(spNumber) 
    628                 sp.plot(vector) 
     631                if not isinstance(vector, (tuple, list)): 
     632                    vector = [vector] 
     633                for v in vector: 
     634                    sp.plot(v) 
    629635            fig.add_subplot(spNumber+1).plot(L) 
    630636            self._check_corr(self.h, h[0,:]) 
    …… 
    642648        self._setupModel(e=[[0.0]]) 
    643649        self._runHG(20, 0.5) 
     650 
     651    def test_hybridGibbs_highSigma(self): 
     652        self._setupModel(e=[[0.5]], fs=[[2.0]]) 
     653        self._runHG(20, 0.5) 

  • projects/AsynCluster/trunk/svpmc/test/test_params.py

    r151 r161  
    109109        def scale(self, x): 
    110110            return self.c * x 
     111        def addThese(self, x, y): 
     112            return self.c + x + y 
     113 
     114    class Scaler(object): 
     115        def __init__(self, scale): 
     116            self.scale = scale 
     117        def scaleThingy(self, thingy): 
     118            return self.scale * thingy.c 
    111119 
    112120    def setUp(self): 
    …… 
    125133                self.failUnlessEqual(z[j, k].c, 100*j+k+y[j, k]) 
    126134 
    127     def test_call(self): 
     135    def test_getattr(self): 
     136        y = self.x.c 
     137        for j in xrange(4): 
     138            for k in xrange(4): 
     139                self.failUnlessEqual(y[j, k], 100*j+k) 
     140 
     141    def test_setattr_scalar(self): 
     142        self.x.c = 0 
     143        for j in xrange(4): 
     144            for k in xrange(4): 
     145                self.failUnlessEqual(self.x[j, k].c, 0) 
     146 
     147    def test_setattr_array(self): 
     148        self.x.c = s.arange(16).reshape(4,4) 
     149        for j in xrange(4): 
     150            for k in xrange(4): 
     151                self.failUnlessEqual(self.x[j, k].c, 4*j+k) 
     152     
     153    def test_call_scalarArgs(self): 
    128154        y = self.x.call('scale', 10) 
    129155        for j in xrange(4): 
    130156            for k in xrange(4): 
    131157                self.failUnlessEqual(y[j, k], 1000*j+10*k) 
     158 
     159    def test_call_arrayArg(self): 
     160        y = self.x.call('addThese', s.arange(16).reshape(4,4), 10000) 
     161        for j in xrange(4): 
     162            for k in xrange(4): 
     163                self.failUnlessEqual( 
     164                    y[j, k], 100*j+k + 4*j+k + 10000, 
     165                    "Error doing addThese method for j=%d, k=%d" % (j,k)) 
     166 
     167    def test_call_flexArrayArg(self): 
     168        arrayScale = params.FlexArray(4, 4) 
     169        for j in xrange(4): 
     170            for k in xrange(4): 
     171                arrayScale[j, k] = self.Scaler(j*k) 
     172        y = arrayScale.call('scaleThingy', self.x) 
     173        for j in xrange(4): 
     174            for k in xrange(4): 
     175                self.failUnlessEqual( 
     176                    y[j, k], j*k * (100*j + k), 
     177                    "Error doing addThese method for j=%d, k=%d" % (j,k)) 
     178 
     179 
     180class Test_ParameterContainer(util.TestCase): 
     181    pass 
     182 
    132183 
    133184 
    …… 
    140191        self.failUnless(s.less(X, 3).all()) 
    141192        self.failUnless(s.equal(p.pdf(X), 0.5).all()) 
     193 
     194 
     195class Test_PriorContainer(util.TestCase): 
     196    p = 3 
     197    shapes = [(p,), (p,p), (0.5*(p**2+p),)  ] 
     198    scales = [ 0.1,  0.2,   0.3             ] 
     199     
     200    def setUp(self): 
     201        self.priorContainer = params.PriorContainer() 
     202        self.priorContainer.paramNames = ['a', 'b', 'c'] 
     203        kw = {'dname':'norm', 'loc':0} 
     204        for i , name in enumerate(self.priorContainer.paramNames): 
     205            shape = self.shapes[i] 
     206            kw['scale'] = self.scales[i] 
     207            priorArray = self.priorContainer.priorArray(name, *shape) 
     208            for j in xrange(int(shape[0])): 
     209                if len(shape) == 1: 
     210                    priorArray[j] = params.Prior(**kw) 
     211                else: 
     212                    for k in xrange(int(shape[1])): 
     213                        priorArray[j,k] = params.Prior(**kw) 
     214 
     215    def test_init(self): 
     216        for k, name in enumerate(self.priorContainer.paramNames): 
     217            priorArray = getattr(self.priorContainer, name) 
     218            self.failUnlessEqual(priorArray.shape, self.shapes[k]) 
     219 
     220    def test_new(self): 
     221        paramContainer = self.priorContainer.new() 
     222        self.failUnless(isinstance(paramContainer, params.ParameterContainer)) 
     223        self.failUnless(isinstance(paramContainer.Lp, float)) 
     224        for k, name in enumerate(self.priorContainer.paramNames): 
     225            paramArray = getattr(paramContainer, name) 
     226            self.failUnlessEqual(paramArray.shape, self.shapes[k]) 
     227 
     228    def test_proposal(self): 
     229        pass 
     230