Changeset 179

Show
Ignore:
Timestamp:
05/14/08 17:30:57 (7 months ago)
Author:
edsuom
Message:

Getting support C code for new model working and tested

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

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

    r178 r179  
    8686// ---------------------------------------------------------------------------- 
    8787 
    88 // Normal (zero mean, unit variance) probability density 
    89 double npdf(double x) 
    90 { 
    91   // 1/sqrt(2*pi) = 0.39894228... 
    92   return 0.398942280401 * exp(-0.5 * pow(x, 2)); 
    93 } 
    94  
    95  
    96 // Normal (zero mean, unit variance) cumulative density 
    97 double ncdf(double x) 
    98 { 
    99   double ax, t, tp, sum; 
    100   // The cdf is computed from the absolute value of x and corrected at the end 
    101   // if necessary 
    102   ax = fabs(x); 
    103   // Compute t, the variable used in the series summation 
    104   t = 1.0 / (1.0 + 0.2316419*ax); 
    105   // Efficiently compute the series summation 
    106   tp = t; 
    107   sum  = tp * 0.31938;  // t 
    108   tp *= t; 
    109   sum -= tp * 0.35656;  // t^2 
    110   tp *= t; 
    111   sum += tp * 1.78148;  // t^3 
    112   tp *= t; 
    113   sum -= tp * 1.82125;  // t^4 
    114   tp *= t; 
    115   sum += tp * 1.33027;  // t^5 
    116   // Compute the cdf, assuming x to be positive 
    117   sum *= npdf(ax); 
    118   // Correct if x < 0 
    119   if(x > 0) 
    120     sum = 1 - sum; 
     88 
     89// Model parameters 
     90struct parameters 
     91{ 
     92  PyArrayObject* d;     // Multivariate log-volatility offset 
     93  PyArrayObject* r_zv;  // Cross-correlation matrix, multivariate normal z to 
     94                        // volatility shock v 
     95  PyArrayObject* e;     // VAR(1) coefficient, volatility shock to log-volatility 
     96  PyArrayObject* ri_y;  // Inverted cross-correlation matrix, innovation shocks 
     97  PyArrayObject* g;     // Correlation, multivariate normal z vs. multivariate 
     98                        // normal for innovation shocks 
     99  PyArrayObject* c;     // Scaling constants for multivariate normal pdf 
     100}; 
     101 
     102 
     103// A single normal variate and the log-volatility sample derived from it 
     104struct sample 
     105{ 
     106  double *z;    // Multivariate normal proposal 
     107  double *h;    // Multivariate log-volatility value 
     108}; 
     109 
     110 
     111// Uniformly distributed random variate in (a,b) 
     112double uniform(double a, double b) 
     113{ 
     114  // Seed the PRNG on the first run only 
     115  static int seeded = 0; 
     116  if (!seeded) { 
     117    int i; 
     118    int j = 0; 
     119    unsigned short seed[3]; 
     120    for(i=0;i<3;i++) { 
     121      do { 
     122        j++; 
     123      } while (j < (i+1)*10000); 
     124      seed[i] = clock()*clock() % 65537; 
     125    } 
     126    seed48(seed); 
     127    seeded = 1; 
     128  } 
     129  // Compute and return a variate from U(a,b) 
     130  return (b - a) * drand48() + a; 
     131} 
     132 
     133 
     134// Matrix multiplication x*y -> +z 
     135// 
     136// The square matrix 'x' is a PyArrayObject and 'y' and 'z' are 1-D arrays. The 
     137// result is added to what's already in 'z'. 
     138void matrixMultiply(PyArrayObject* x, double *y, double *z) 
     139{ 
     140  int kr, km; 
     141  double val, sum; 
     142  const int rows = (int) x->dimensions[0]; 
     143  // Compute the dot product for each row of x 
     144  for(kr=0; kr<rows; kr++) { 
     145    // For this row of y and z... 
     146    sum = 0; 
     147    // Matrix multiplication dot product 
     148    for(km=0; km<rows; km++) 
     149      sum += PP2(x, kr, km) * PA1(y, km); 
     150    PA1(z, kr) += sum; 
     151  } 
     152} 
     153 
     154 
     155// Computes the log-volatility for a sample, given a 1-D array containing a 
     156// multivariate normal z for the sample, a like-dimensioned array h0 containing 
     157// a log-volatility value previous to the sample, and the model parameters. 
     158// 
     159// The z array is in a supplied sample struct, which the function modifies with 
     160// the log-volatility value corresponding to the proposal. 
     161void logVol(struct sample *S, struct parameters *P, double *h0) 
     162{ 
     163  int k; 
     164  const int N = (int) P->d->dimensions[0]; 
     165  for(k=0; k<N; k++) 
     166    SPA1(S, h, k) = 0; 
     167 
     168  // Run the VAR(1) process to get the current multivariate log-volatility 
     169  // value 
     170 
     171  // Start with the matrix product for the VAR(1) term... 
     172  matrixMultiply(P->e, h0, S->h); 
     173  // then compute and add the multivariate volatility shock given the proposed 
     174  // normal variate... 
     175  matrixMultiply(P->r_zv, S->z, S->h); 
     176  // ...plus the multivariate volatility offset 
     177  for(k=0; k<N; k++) 
     178    SPA1(S, h, k) += SPP1(P, d, k); 
     179} 
     180 
     181 
     182// Computes and returns the log-likelihood of a multivariate innovation x given 
     183// a 1-D array containing a multivariate log-volatility. 
     184double logLikelihood(struct sample *S, struct parameters *P, double *x) 
     185{ 
     186  int k; 
     187  double *y; 
     188  double val; 
     189  double L = 0; 
     190  const int N = (int) P->d->dimensions[0]; 
     191  y = (double *) malloc(2*N * sizeof(double)); 
     192 
     193  // Inverse-scale the innovation by the square root of the volatilities in 
     194  // preparation for decorrelation 
     195  for(k=0; k<N; k++) 
     196    PA1(y, k) = exp(-0.5 * SPA1(S, h, k)) * PA1(x, k); 
     197     
     198  // Now decorrelate the equi-variance innovations in preparation for the 
     199  // log-likelihood computation for this multivariate log-volatility proposal 
     200  matrixMultiply(P->ri_y, y, y+N); 
     201 
     202  // Now compute log Pr(x[t]|h[t]), the log-likelihood of the decorrelated 
     203  // multivariate innovation for this sample, conditional on the sample's 
     204  // log-volatility. 
     205   
     206  // Now the log-likelihoood. For each time series... 
     207  for(k=0; k<N; k++) { 
     208    // Offset the decorrelated innovation by the mean, which is the correlation 
     209    // between the normal variates for innovation and proposed volatility 
     210    // shocks, times the value of the normal variate value underlying the 
     211    // volatility shock for this sample 
     212    PA1(y, k) = PA1(y, k+N) - SPP1(P, g, k) * SPA1(S, z, k); 
     213    // The log-likelihood for this sample is simply the argument of the 
     214    // exponential of the normal pdf, with the reciprocal scaling of the 
     215    // exponential being done in log space by subtraction of the log of the 
     216    // scaling coefficient 
     217    val = PA1(y, k) / SPP2(P, c, k, 0); 
     218    L -= 0.5 * val*val + SPP2(P, c, k, 1); // Fast squaring instead of pow() 
     219    // Since the innovation has been inverse-scaled by the square root of the 
     220    // proposed volatility, the probability density needs to be scaled as 
     221    // well. In log space, this is a subtraction of half the log-volatility. 
     222    L -= 0.5 * SPA1(S, h, k); 
     223  } 
     224 
    121225  // All done 
    122   return sum; 
    123 } 
    124  
    125  
    126 // Truncated normal (zero mean, unit variance) probability density 
    127 double tnpdf(double x, double a, double b) 
    128 { 
    129   double y; 
    130   if(x < a || x > b) { 
    131     y = 0; 
    132   } else { 
    133     y = npdf(x) / (ncdf(b) - ncdf(a)); 
    134   } 
    135   return y; 
    136 } 
    137   
    138  
    139  
    140 // Model.zProposal 
    141 // 
    142 // Supplied variables 
    143 // ---------------------------------------------------------------------------- 
    144 // z    Independent normal variates, [pn] array (updated) 
    145 // wiggle = half-width of independent uniform proposal 
    146  
    147 #define  w  (double)wiggle 
    148  
    149 int k0, k1; 
    150 double x, zk, c, cp, L; 
    151  
    152 // Seed the PRNG on the first run only 
    153 static int seeded = 0; 
    154 if (!seeded) { 
    155   int i; 
    156   int j = 0; 
    157   unsigned short seed[3]; 
    158   for(i=0;i<3;i++) { 
    159     do { 
    160       j++; 
    161     } while (j < (i+1)*10000); 
    162     seed[i] = clock()*clock() % 65537; 
    163   } 
    164   seed48(seed); 
    165   seeded = 1; 
    166 } 
    167  
    168 // The joint log-density of all the jumps 
    169 L = 0; 
    170 // for each multivariate sample... 
    171 for(k0=0; k0<Nz[1]; k0++) { 
    172  
    173   // Do a random walk (a single uniformly distributed step) from this column of 
    174   // independent random-walk normal variates to produce a multivariate iid 
    175   // normal proposal. It is this on which we build a multivariate 
    176   // log-volatility proposal for this multivariate sample. 
    177    
    178   // For each time series... 
    179   for(k1=0; k1<Nz[0]; k1++) { 
    180     // Compute the scale value c as a bit more* than the maximum pdf (regular 
    181     // normal distribution) at the endpoints and middle of the truncated range 
    182     c = 0; 
    183     zk = Z2(k1, k0); 
    184     for(x = -w; x < w; x += w) { 
    185       cp = 1.12 * npdf(zk + x); // *Scaling by 1.12 works for w <= 1.0 
    186       if(cp > c) 
    187         c = cp; 
    188     } 
    189     // Accept-reject sampling for the truncated normal distribution 
    190     do { 
    191       // Generate a proposal within the truncated range 
    192       x = zk + 2*w*(drand48() - 0.5); 
    193       // Compute the pdf of the proposal, using the regular normal distribution 
    194       cp = npdf(x); 
    195       // Accept (and exit the loop) when the scaled uniform variate is less 
    196       // than the proposal's pdf 
    197     } while(c*drand48() > cp); 
    198     // Update this element of the IID normal array in place with the accepted 
    199     // proposal 
    200     Z2(k1, k0) = x; 
    201     // Add the log-density of the accepted proposal 
    202     L += log(tnpdf(x, zk-w, zk+w)); 
    203   } 
    204 } 
    205  
    206 return_val = L; 
    207  
    208  
    209  
    210 // Model.logVolatilities 
     226  free(y); 
     227  return L; 
     228} 
     229 
     230 
     231// Model.hybridGibbs 
    211232//  
    212233// Supplied variables 
    …… 
    214235 
    215236//--- Supplied arrays --- 
    216 // z    Independent normal variates for log-volatilities, [pn] array 
    217 // v    Log-volatility shocks, [pn] 
     237// z    Independent normal variates, [pn] array (updated) 
     238// h    Simulated last-sample multivariate log-volatility, [p] (overwritten) 
    218239// x    Innovation values, [pn] array 
    219 // h    Log-volatility values, [p(n+1)] array (last n samples overwritten) 
    220240// ri   Inverse, concurrent cross-correlations between innovation shocks, [pp] 
    221 // y    Miscellaneous multivariate scratchpad, [p8] array (overwritten) 
     241// c    Constants for multivariate normal PDF 
    222242 
    223243//--- Model parameters --- 
    …… 
    226246// g    Volatility/innovation shock correlations, [p] vector 
    227247 
    228  
    229 // The multivariate scratchpad y is used as follows, by column: 
    230 // ---------------------------------------------------------------------------- 
    231 // 0    Bivariate norm PDF: sqrt(1-sigma**2) 
    232 // 1    Bivariate norm PDF, log-scaling: log(sqrt(2*pi)*sqrt(1-sigma**2)) 
    233 // 2    Innovation values after inverse-scaling by log-volatility 
    234 // 3    Inverse-scaled innovation values after decorrelation 
    235  
    236  
    237 int k0, k1, km; 
    238 double sum; 
    239 double L = 0; 
     248//--- Return values (tuple) --- 
     249// 0    Lx: Log-likelihood of innovations given simulated log-volatilities 
     250// 1    Lh: Log-likelihood of simulated log-volatilities 
     251 
     252 
     253int k0, k1; 
     254double sum, L; 
    240255 
    241256py::tuple result(2); 
    242  
    243257 
    244258// For each sample... 
    245259for(k0=0; k0<Nz[1]; k0++) { 
    246260 
    247   // Run the VAR(1) process to get the current multivariate log-volatility 
    248   // value 
    249  
    250   // For each time series... 
    251   for(k1=0; k1<Nz[0]; k1++) { 
    252     // The volatility offset plus the shock... 
    253     sum = D1(k1) + V2(k1, k0); 
    254     // ...plus the dot product for the VAR(1) term... 
    255     for(km=0; km<Nz[0]; km++) 
    256       sum += E2(k1, km) * H2(km, k0); 
    257     // ...is the current log-volatility value 
    258     H2(k1, k0+1) = sum; 
    259   } 
    260    
    261   // Now inverse-scale the innovations for this sample by the square root of 
    262   // the volatilities in preparation for decorrelating the innovations 
    263   for(k1=0; k1<Nz[0]; k1++) 
    264     Y2(k1, 2) = exp(-0.5 * H2(k1, k0+1)) * X2(k1, k0); 
    265      
    266   // Now decorrelate the equi-variance innovations in preparation for the 
    267   // log-likelihood computation for this multivariate log-volatility proposal 
    268      
    269   // For each time series... 
    270   for(k1=0; k1<Nz[0]; k1++) { 
    271     // The product of the inverted cross-correlation matrix and the correlated 
    272     // variates... 
    273     sum = 0; 
    274     for(km=0; km<Nz[0]; km++) 
    275       sum += RI2(k1, km) * Y2(km, 2); 
    276     // ...is the decorrelated multivariate innovation for this sample, given 
    277     // its log-volatility 
    278     Y2(k1, 3) = sum; 
    279   } 
    280      
    281   // Now compute log Pr(x[t]|h[t]), the log-likelihood of the decorrelated 
    282   // multivariate innovation for this sample, conditional on the just-computed 
    283   // log-volatility h[t]. 
    284    
    285   // For each time series... 
    286   for(k1=0; k1<Nz[0]; k1++) { 
    287     // Offset the decorrelated innovation by the mean, which is the correlation 
    288     // between the normal variates for innovation and proposed volatility 
    289     // shocks, times the value of the normal variate value underlying the 
    290     // volatility shock for this sample 
    291     Y2(k1, 3) = Y2(k1, 3)  -  G1(k1) * Z2(k1, k0); 
    292     // The log-likelihood for this sample is simply the argument of the 
    293     // exponential of the normal pdf, with the reciprocal scaling of the 
    294     // exponential being done in log space by subtraction of the log of the 
    295     // scaling coefficient 
    296     L -= 0.5 * pow(Y2(k1, 3) / Y2(k1, 0), 2) + Y2(k1, 1); 
    297     // Since the innovation has been inverse-scaled by the square root of the 
    298     // proposed volatility, the probability density needs to be scaled as 
    299     // well. In log space, this is a subtraction of half the log-volatility. 
    300     L -= 0.5 * H2(k1, k0+1); 
    301   } 
    302  
    303   // The probability density, on its own, of the first log-volatility sample is 
    304   // proportional to the probability density of the innovation for that sample 
    305   // given the log-volatility. So we'll be returning that value, too. 
    306   if(k0 == 0) 
    307     result[0] = L; 
    308 } 
    309  
    310 // Return the first log-likelihood along with the total of the log-likelihoods 
    311 result[1] = L; 
     261  // TODO 
     262  ; 
     263} 
     264 
     265// All done 
     266free(P->d); 
     267free(P->r_zv); 
     268free(P->e); 
     269free(P->ri_y); 
     270free(P->g); 
     271free(P->c0); 
     272free(P->c1); 
    312273return_val = result; 

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

    r178 r179  
    180180 
    181181    """ 
    182     _logU_index = -1 
    183  
    184182    keyAttrs = {'y':None, 'Mv':1} 
    185183 
    …… 
    224222        return cv 
    225223 
    226     def decaySamples(self, tol): 
    227         """ 
    228         Returns as an int the number of samples needed for the effect of an 
    229         impulse to the log-volatility VAR(1) process to decay below I{tol}. 
    230         """ 
    231         if 'ds' not in self.cache: 
    232             p = self.e.shape[0] 
    233             # The impulse 
    234             x = s.column_stack([s.ones(p), s.zeros((p, self.n))]) 
    235             # Simulate the impulse response with a VAR object 
    236             self.cache['ds'] = VAR(x).forward(x, s.zeros(p), self.e).max(0) 
    237         I = s.less(self.cache['ds'], tol).nonzero()[0] 
    238         if len(I): 
    239             return I.min() 
    240         return self.n 
    241      
    242     def logUniform(self): 
    243         """ 
    244         Returns a log-uniform variate taken from a large array that is 
    245         generated and refreshed periodically. Generating the uniform variates 
    246         and taking their log in a single large array operation is more 
    247         efficient than doing so one value at a time. 
    248         """ 
    249         logU = getattr(self, '_logU_array', []) 
    250         if self._logU_index >= len(logU)-1: 
    251             self._logU_index = -1 
    252             # The efficient array operation 
    253             logU = self._logU_array = s.log(s.rand(10000)) 
    254         self._logU_index += 1 
    255         return logU[self._logU_index] 
    256  
    257     def zProposal(self, z, wiggle): 
    258         """ 
    259         Does a random walk step from the supplied 2-D array I{z} of zero-mean, 
    260         unit-variance normal random variates with a jump from a uniform 
    261         distribution having width +/- I{wiggle}. 
    262          
    263         Returns a copy of the array with all elements changed from the original 
    264         values, independently, along with the total log-density of the jumps. 
    265         """ 
    266         z = z.copy() 
    267         L = self.inline('z', 'wiggle') 
    268         return z, L 
    269      
    270224    def logVolatilities(self, z, v, x, h0): 
    271225        """ 
    …… 
    313267        Updates the IID normal variates in place. Returns the likelihood of the 
    314268        innovations in I{x}, given the simulated log-volatilities, along with 
    315         the last sample's multivariate log-volatility value, the 
    316         log-probability of the simulated result, and the fractional acceptance 
    317         rate of the MCMC steps. 
     269        the last sample's multivariate log-volatility value, and the 
     270        log-probability of the simulated result. 
    318271         
    319272        The method will account for the effect of log-volatility proposals on 
    …… 
    323276        """ 
    324277        L_total = 0 
    325         acceptances = 0 
    326278        N = self.decaySamples(tol) 
    327279        # Initialize volatility shocks from cross-correlation of the IID 
    …… 
    364316                z[:,k] = zp[:,k] 
    365317                v[:,k] = vp[:,k] 
    366                 acceptances += 1 
    367318            # Update stuff for the next sample 
    368319            h0 = LV[2][:,0] 
    …… 
    372323            if L_total < -1E6: 
    373324                return -s.inf, s.zeros_like(h0), 0.0, 0.0 
    374         return L_total, h0, L_result, float(acceptances)/self.n 
     325        return L_total, h0, L_result 
    375326     
    376327 

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

    r178 r179  
    2020""" 
    2121 
    22 import os.path 
     22import os.path, textwrap 
    2323import scipy as s 
    2424from scipy import stats, random, linalg, signal 
    …… 
    161161        return model.Model(**kw) 
    162162 
     163    def inline(self, code, **kw): 
     164        """ 
     165        Do an inline call with the model's Weaver API and support code. 
     166        """ 
     167        return self.modelFactory().inlineDirect(textwrap.dedent(code), **kw) 
     168     
    163169    def _covarMatrix(self, cs, cr): 
    164170        i = 0 
    …… 
    316322            [[9.0, 1.2, 3.0], [1.2, 16.0, 6.0], [3.0, 6.0, 25.0]]) 
    317323 
    318     def test_logUniform(self): 
    319         N = 35000 
    320         import timeit 
    321         tObserved = timeit.Timer( 
    322             "M.logUniform()", 
    323             "import model\nM = model.Model()").timeit(N) 
    324         tBaseline = timeit.Timer( 
    325             "s.log(s.rand())", 
    326             "import scipy as s").timeit(N) 
    327         percent = 100 * tObserved/tBaseline 
    328         if VERBOSE: 
    329             print "The efficient method took %2d%% as long as the stupid way" \ 
    330                   % percent 
    331         self.failUnless(percent < 40) 
    332  
    333     def test_zProposal_dist(self): 
    334         N = 10000 
    335         z = s.zeros((2,3)) 
    336         modelObj = self.modelFactory() 
    337         z1 = s.empty(N) 
    338         z2 = s.empty(N) 
    339         for k in xrange(N): 
    340             z, L_result = modelObj.zProposal(z, 0.8) 
    341             z1[k] = z[0,0] 
    342             z2[k] = z[1,1] 
    343         self.failUnlessNormal(z1) 
    344         self.failUnlessNormal(z2) 
    345         self._check_uncorr(z1, z2) 
    346  
    347     def test_zProposal_L(self): 
    348         N = 1000 
    349         z0 = s.ones((1,1)) 
    350         modelObj = self.modelFactory() 
    351         z = s.empty(N) 
    352         L = s.empty(N) 
    353         for k in xrange(N): 
    354             zk, L_result = modelObj.zProposal(z0, 0.8) 
    355             z[k] = zk[0,0] 
    356             L[k] = L_result 
    357         fig = self.fig 
    358         sp = fig.add_subplot(211) 
    359         sp.hist(z) 
    360         sp = fig.add_subplot(212) 
    361         sp.plot(z, s.exp(L), '.') 
    362      
    363     def test_decaySamples(self): 
    364         x = s.zeros((5, 1000)) 
    365         x[2,0] = 1.0 # The impulse 
    366         var = model.VAR(x) 
    367         a = s.zeros(5) 
    368         modelObj = self.modelFactory() 
    369         for j in xrange(20): 
    370             # Keep all coefficients small enough to ensure stability 
    371             b = modelObj.e = 0.3*s.rand(5,5) 
    372             y = var.forward(x, a, b) 
    373             for tol in s.logspace(-5, -1, 10): 
    374                 N = modelObj.decaySamples(tol) 
    375                 print tol, N 
    376                 failed = abs(y[:,N].max()) > tol 
    377                 if (N > 100 or failed) and VERBOSE: 
    378                     N_plot = min([3*N, 30]) 
    379                     sp = self.fig.add_subplot(111) 
    380                     sp.semilogy(y[:,:N_plot].transpose()) 
    381                     sp.axvline(N) 
    382                     sp.axhline(tol) 
    383                     self.failIf( 
    384                         failed, 
    385                         "Estimate of %d samples " % N +\ 
    386                         "to decay within %f tolerance was insufficient" % tol) 
    387          
     324    def test_support_uniform(self): 
     325        code = """ 
     326        return_val = uniform((double) a, (double) b); 
     327        """ 
     328        N = 100 
     329        for x, y in ((1, 2), (-1, +3)): 
     330            z = s.array([self.inline(code, a=x, b=y) for k in xrange(N)]) 
     331            self.failUnless(len(z), N) 
     332            self.failUnless(z.min() > x) 
     333            self.failUnless(z.max() < y) 
     334            self.failUnlessEqual(len(s.unique(z)), N) 
     335 
     336    def test_support_matrixMultiply(self): 
     337        code = """ 
     338        int k; 
     339        double b[8]; 
     340        double c[8]; 
     341        for(k=0; k<8; k++) { 
     342          b[k] = 1.5*k; 
     343          c[k] = 0; 
     344        } 
     345        matrixMultiply(a_array, b, c); 
     346        for(k=0; k<8; k++) 
     347          R1(k) = c[k]; 
     348        """ 
     349        a = s.arange(64, dtype='d').reshape(8,8) 
     350        b = 1.5 * s.arange(8).reshape((8,1)) 
     351        r = s.empty(8) 
     352        self.inline(code, a=a, r=r) 
     353        self.failUnlessElementsEqual(r, s.matrix(a)*b) 
     354 
    388355 
    389356class Test_Model_logVolatilities_VAR(BaseCase): 

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

    r175 r179  
    3434 
    3535//---------------------------------------- 
     36// TestableWeaver.support 
     37 
     38struct vectorPackage 
     39{ 
     40  PyArrayObject* x; 
     41  PyArrayObject* y; 
     42  double result; 
     43}; 
     44 
     45void dotProduct(struct vectorPackage *vp) 
     46{ 
     47  int k; 
     48  int N = (int) vp->x->dimensions[0]; 
     49  vp->result = 0; 
     50  for(k=0; k<N; k++) 
     51    vp->result += SPP1(vp, x, k) * SPP1(vp, y, k); 
     52} 
     53 
     54//---------------------------------------- 
    3655// TestableWeaver.foo 
    3756 
    …… 
    4362X1(0) += 100; 
    4463//---------------------------------------- 
     64 
    4565 
    4666// Bogus.something 
    …… 
    6484        return x[0] 
    6585 
     86    def dotProduct(self, x, y): 
     87        code = "\n".join([ 
     88            line.strip() for line in [ 
     89            "struct vectorPackage vectors;", 
     90            "vectors.x = x_array;", 
     91            "vectors.y = y_array;", 
     92            "dotProduct(&vectors);", 
     93            "return_val = vectors.result;"]]) 
     94        return self.inlineDirect(code, x=x, y=y) 
     95 
    6696 
    6797class Test_Weaver(util.TestCase): 
    …… 
    6999        self.weaver = TestableWeaver() 
    70100 
    71     def test_inline_alone(self): 
    72         self.weaver._code = {'support':None, 'foo':CODE} 
    73         self.failUnlessEqual(self.weaver.foo(s.array([1])), 1164) 
    74  
     101    def test_support(self): 
     102        x = s.randn(100) 
     103        y = s.randn(100) 
     104        self.weaver._code = self.weaver._parseCode(CODE) 
     105        self.failUnlessEqual(self.weaver.dotProduct(x, y), s.dot(x, y)) 
     106     
    75107    def test_parseCode(self): 
    76108        code = self.weaver._parseCode(CODE) 

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

    r175 r179  
    2626from pkg_resources import resource_string 
    2727from scipy.weave import * 
     28from scipy.weave.base_info import custom_info 
     29 
     30 
     31class my_info(custom_info): 
     32    _extra_compile_args = ['-w'] 
     33    __macros = [ 
     34        # Access an element i of a 1-D array whose pointer is pa 
     35        'PA1(pa,i)', 
     36        "*(pa + i)", 
     37 
     38        # Access an element i of a 1-D array whose pointer is pa in a structure 
     39        # whose pointer is s 
     40        'SPA1(s,pa,i)', 
     41        "*(s->pa + i)", 
     42 
     43        # Access an element i, j of a 2-D Numpy array object whose pointer is 
     44        # ppa 
     45        'PP2(ppa,i,j)', 
     46        "(*((double*)(ppa->data + (i)*ppa->strides[0] + (j)*ppa->strides[1])))", 
     47 
     48        # Access an element i of a 1-D Numpy array object whose pointer is ppa 
     49        # in a structure whose pointer is s 
     50        'SPP1(s,ppa,i)', 
     51        "(*((double*)(s->ppa->data + (i)*s->ppa->strides[0])))", 
     52         
     53        # Access an element i, j of a 2-D Numpy array object whose pointer is 
     54        # ppa in a structure whose pointer is s 
     55        'SPP2(s,ppa,i,j)', 
     56        "(*((double*)(s->ppa->data + "+\ 
     57        "(i)*s->ppa->strides[0] + (j)*s->ppa->strides[1])))", 
     58        ] 
     59 
     60    def __init__(self): 
     61        k = 0 
     62        macros = ["'%s'" % (x,) for x in self.__macros] 
     63        while(k < len(macros)): 
     64            name, value = macros[k:k+2] 
     65            self.add_define_macro((name, value)) 
     66            self.add_undefine_macro(name) 
     67            k += 2 
    2868 
    2969 
    …… 
    4585    line. 
    4686    """ 
     87    infoObj = my_info() 
     88     
    4789    def _parseCode(self, text): 
    4890        re_delimiter = re.compile( 
    …… 
    111153        return inline( 
    112154            self.code[methodName], args, local_dict, 
    113             extra_compile_args=['-w'], support_code=self.code['support']) 
     155            customize=self.infoObj, support_code=self.code['support']) 
    114156 
    115              
    116          
     157    def inlineDirect(self, code, **kw): 
     158        """ 
     159        Call this method with a string of C I{code} and keywords defining the 
     160        names and values of variables to be passed to the code function. 
     161 
     162        Using this method entails less CPU overhead (but also less convenience 
     163        and flexibility) than using L{inline}. 
     164        """ 
     165        return inline( 
     166            code, kw.keys(), kw, 
     167            customize=self.infoObj, support_code=self.code['support'])