Changeset 181

Show
Ignore:
Timestamp:
05/15/08 10:05:12 (7 months ago)
Author:
edsuom
Message:

Unit testing model support C code

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

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

    r180 r181  
    100100                     // volatility shock v 
    101101  PyArrayObject* ri; // Inverted cross-correlation matrix, innovation shocks 
    102   PyArrayObject* c; // Scaling constants for multivariate normal pdf, given g 
     102  PyArrayObject* sc; // Scaling constants for multivariate normal pdf, given g 
    103103}; 
    104104 
    …… 
    153153 
    154154 
    155 // Matrix multiplication x*y -> +z 
    156 // 
    157 // The square matrix 'x' is a PyArrayObject and 'y' and 'z' are 1-D arrays. The 
    158 // result is added to what's already in 'z'. 
     155// Matrix multiplication x*y -> z, where the square matrix 'x' is a 
     156// PyArrayObject and 'y' and 'z' are 1-D arrays. 
    159157void matrixMultiply(PyArrayObject* x, double *y, double *z) 
    160158{ 
    161159  int kr, km; 
    162   double val, sum; 
     160  register double sum; 
    163161  const int rows = (int) x->dimensions[0]; 
    164162  // Compute the dot product for each row of x 
    …… 
    169167    for(km=0; km<rows; km++) 
    170168      sum += PP2(x, kr, km) * PA1(y, km); 
    171     PA1(z, kr) += sum; 
     169    PA1(z, kr) = sum; 
    172170  } 
    173171} 
    …… 
    180178  int k; 
    181179  const int N = (int) P->d->dimensions[0]; 
    182   for(k=0; k<N; k++) 
    183     SPA1(S, h, k) = 0; 
    184180 
    185181  // Run the VAR(1) process to get the current multivariate log-volatility 
    …… 
    188184  // Start with the matrix product for the VAR(1) term... 
    189185  matrixMultiply(P->e, S->h, S->x0); 
    190   // then compute and add the multivariate volatility shock given the proposed 
     186  // ...then compute the multivariate volatility shock given the proposed 
    191187  // normal variate... 
    192   matrixMultiply(P->rz, S->z, S->x0); 
    193   // ...plus the multivariate volatility offset 
     188  matrixMultiply(P->rz, S->z, S->x1); 
     189  // ...then add them, plus the multivariate volatility offset, to be the new h 
    194190  for(k=0; k<N; k++) 
    195     SPA1(S, x0, k) += SPP1(P, d, k); 
    196   // Now overwrite the previous log-volatility value 
    197   for(k=0; k<N; k++) 
    198     // TODO: Do this more efficiently by swapping pointers? 
    199     SPA1(S, h, k) = SPA1(S, x0, k); 
    200 } 
    201  
    202  
    203 // Computes and returns the log-likelihood of a multivariate innovation x given 
    204 // a 1-D array containing a multivariate log-volatility. 
     191    SPA1(S, h, k) = SPA1(S, x0, k) + SPA1(S, x1, k) + SPP1(P, d, k); 
     192} 
     193 
     194 
     195// Computes the log-likelihood of a multivariate innovation x given a 1-D array 
     196// containing a multivariate log-volatility. 
    205197void logLikelihood(struct sample *S, struct parameters *P, double *x) 
    206198{ 
    207   int k; 
     199  register int k; 
    208200  double *y; 
    209201  double val; 
    …… 
    235227    // exponential being done in log space by subtraction of the log of the 
    236228    // scaling coefficient 
    237     val = SPA1(S, x0, k) / SPP2(P, c, k, 0); 
    238     L -= 0.5 * val*val + SPP2(P, c, k, 1); // Fast squaring instead of pow() 
     229    val = SPA1(S, x0, k) / SPP2(P, sc, k, 0); 
     230    L -= 0.5 * val*val + SPP2(P, sc, k, 1); // Fast squaring instead of pow() 
    239231    // Since the innovation has been inverse-scaled by the square root of the 
    240232    // proposed volatility, the probability density needs to be scaled as 
    …… 
    244236 
    245237  // All done 
    246   S.L = L; 
     238  S->Lx = L; 
    247239} 
    248240 
    …… 
    250242// Returns with zero only if the proposals are close enough to finish 
    251243// likelihood computation 
    252 int notCloseEnough(struct sample *sp, int n) 
     244int notCloseEnough(struct sample *sp[], int n) 
    253245{ 
    254246  int k; 
    255247  double L; 
    256   double Lmin = +DBL_MAX;   // Anything will be less positive 
    257   double Lmax = -DBL_MAX;   // Anything will be more positive 
     248  double Lmin = +1E20;   // Anything will be less positive 
     249  double Lmax = -1E20;   // Anything will be more positive 
    258250  for(k=0; k<n; k++) { 
    259     L = PA1(sp, k).Lx; 
     251    L = sp[k]->Lx; 
    260252    if(L < Lmin) 
    261253      Lmin = L; 
    …… 
    288280// rz   Concurrent xcorr, innovation-volatility normals, [pp] 
    289281// ri   Inverse concurrent xcorr, innovation shocks, [pp] 
    290 // c    Constants for multivariate normal PDF 
     282// sc   Constants for multivariate normal PDF 
    291283 
    292284//--- Return values (tuple) --- 
    …… 
    307299// Generate parameter struct 
    308300P.d = d_array; P.e = e_array; P.g = g_array; 
    309 P.rz = rz_array; P.ri = ri_array; P.c = c_array; 
     301P.rz = rz_array; P.ri = ri_array; P.sc = sc_array; 
    310302 
    311303// Initialize various arrays 

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

    r179 r181  
    354354 
    355355 
    356 class Test_Model_logVolatilities_VAR(BaseCase): 
     356class Test_Model_logVol(BaseCase): 
    357357    N = 1000 
    358      
    359     def _draw_h(self, modelObj): 
    360         modelObj.cs = s.array([1.0, 1.0]) 
    361         modelObj.cr = [0.0] 
    362         modelObj.g = s.array([0.0]) 
    363         z = s.randn(modelObj.p, modelObj.n) 
    364         v = s.array(self._rv(modelObj.fs, modelObj.fr) * z) 
    365         return modelObj.logVolatilities( 
    366             z, v, s.zeros_like(z), modelObj.d.copy())[2] 
     358 
     359    code = """ 
     360    int k; 
     361    struct sample S; 
     362    struct parameters P; 
     363 
     364    P.d = d_array; P.e = e_array; P.rz = rz_array; 
     365    S.z = (double *) malloc(sizeof(double) * Nz[0]); 
     366    S.h = (double *) calloc(sizeof(double), Nz[0]); 
     367    S.x0 = (double *) malloc(sizeof(double) * Nz[0]); 
     368    S.x1 = (double *) malloc(sizeof(double) * Nz[0]); 
     369    for(k=0; k<Nz[0]; k++) 
     370      PA1(S.z, k) = Z1(k); 
     371 
     372    logVol(&S, &P); 
     373 
     374    for(k=0; k<Nz[0]; k++) 
     375      H1(k) = PA1(S.h, k); 
     376    """ 
     377     
     378    def _draw_h(self, z, d, e, fs, fr): 
     379        h = s.zeros_like(z) 
     380        rz = s.array(self._rv(fs, fr)) 
     381        for k, zk in enumerate(z.transpose()): 
     382            self.inline( 
     383                self.code, z=zk, 
     384                d=s.array(d), e=s.array(e), h=h[:,k], rz=rz) 
     385        return h 
    367386 
    368387    def test_univariate_LPF(self): 
    369388        e = 0.95 
    370389        # Univariate model 
    371         modelObj = self.modelFactory( 
    372             y = s.zeros((1, self.N)), 
    373             cs = [1.0], 
    374             cr = [], 
    375             d = s.array([0.0]), 
    376             e = s.array([[e]]), 
    377             g = s.array([0.0])) 
    378390        z = s.randn(1, self.N) 
    379         h = modelObj.logVolatilities( 
    380             z, z.copy(), s.zeros_like(z), s.array([0.0]))[2] 
     391        h = self._draw_h(z, d=[1], e=[[e]], fs=[1], fr=[]) 
    381392        # IIR filtering 
    382393        self._check_psd(h, e)