1#cython: boundscheck=False 2#cython: wraparound=False 3#cython: cdivision=False 4""" 5State Space Model - Cython tools 6 7Author: Chad Fulton 8License: Simplified-BSD 9""" 10 11# Typical imports 12cimport numpy as np 13cimport cython 14import numpy as np 15 16np.import_array() 17 18from statsmodels.src.math cimport * 19cimport scipy.linalg.cython_blas as blas 20cimport scipy.linalg.cython_lapack as lapack 21 22cdef FORTRAN = 1 23 24ctypedef fused sm_type_t: 25 np.float64_t 26 np.float32_t 27 np.complex64_t 28 np.complex128_t 29 30 31# ------------------------------------------------------------------------ 32# Array shape validation 33 34cdef validate_matrix_shape(str name, Py_ssize_t *shape, int nrows, int ncols, nobs=None): 35 if not shape[0] == nrows: 36 raise ValueError('Invalid shape for %s matrix: requires %d rows,' 37 ' got %d' % (name, nrows, shape[0])) 38 if not shape[1] == ncols: 39 raise ValueError('Invalid shape for %s matrix: requires %d columns,' 40 'got %d' % (name, shape[1], shape[1])) 41 if nobs is not None and shape[2] not in [1, nobs]: 42 raise ValueError('Invalid time-varying dimension for %s matrix:' 43 ' requires 1 or %d, got %d' % (name, nobs, shape[2])) 44 45 46cdef validate_vector_shape(str name, Py_ssize_t *shape, int nrows, nobs = None): 47 if not shape[0] == nrows: 48 raise ValueError('Invalid shape for %s vector: requires %d rows,' 49 ' got %d' % (name, nrows, shape[0])) 50 if nobs is not None and not shape[1] in [1, nobs]: 51 raise ValueError('Invalid time-varying dimension for %s vector:' 52 ' requires 1 or %d got %d' % (name, nobs, shape[1])) 53 54 55# ------------------------------------------------------------------------ 56# Blas Wrapping 57 58cdef inline copy(int size, sm_type_t* src_arr, sm_type_t* target, 59 int incx=1, int incy=1): 60 """ 61 Parameters 62 ---------- 63 size : int 64 number of elements in vectors `src_arr` and `target` 65 src_arr : sm_type_t* 66 Pointer to Array that we are copying _from_ 67 target : sm_type_t* 68 Pointer to Array that we are copying _to_ 69 incx : int, default 1 70 Specifies the increment for the elements of `src_arr` 71 incy : int, default 1 72 Specifies the increment for the elements of `target` 73 74 References 75 ---------- 76 https://software.intel.com/en-us/mkl-developer-reference-c-cblas-copy 77 """ 78 if sm_type_t is np.float32_t: 79 blas.scopy(&size, src_arr, &incx, target, &incy) 80 elif sm_type_t is np.float64_t: 81 blas.dcopy(&size, src_arr, &incx, target, &incy) 82 elif sm_type_t is np.complex64_t: 83 blas.ccopy(&size, src_arr, &incx, target, &incy) 84 elif sm_type_t is np.complex128_t: 85 blas.zcopy(&size, src_arr, &incx, target, &incy) 86 87 88cdef inline swap(int size, sm_type_t* src_arr, sm_type_t* target, 89 int incx=1, int incy=1): 90 """ 91 Parameters 92 ---------- 93 size : int 94 number of elements in vectors `src_arr` and `target` 95 src_arr : sm_type_t* 96 Pointer to first of two arrays being swapped 97 target : sm_type_t* 98 Pointer to second of two arrays being swapped 99 incx : int, default 1 100 Specifies the increment for the elements of `src_arr` 101 incy : int, default 1 102 Specifies the increment for the elements of `target` 103 104 References 105 ---------- 106 https://software.intel.com/en-us/mkl-developer-reference-c-cblas-swap 107 """ 108 109 if sm_type_t is np.float32_t: 110 blas.sswap(&size, src_arr, &incx, target, &incy) 111 elif sm_type_t is np.float64_t: 112 blas.dswap(&size, src_arr, &incx, target, &incy) 113 elif sm_type_t is np.complex64_t: 114 blas.cswap(&size, src_arr, &incx, target, &incy) 115 elif sm_type_t is np.complex128_t: 116 blas.zswap(&size, src_arr, &incx, target, &incy) 117 118 119# ------------------------------------------------------------------------ 120 121{{py: 122 123TYPES = { 124 "s": ("np.float32_t", "np.float32", "np.NPY_FLOAT32"), 125 "d": ("np.float64_t", "float", "np.NPY_FLOAT64"), 126 "c": ("np.complex64_t", "np.complex64", "np.NPY_COMPLEX64"), 127 "z": ("np.complex128_t", "complex", "np.NPY_COMPLEX128"), 128} 129 130}} 131 132{{for prefix, types in TYPES.items()}} 133{{py:cython_type, dtype, typenum = types}} 134{{py: 135combined_prefix = prefix 136combined_cython_type = cython_type 137if prefix == 'c': 138 combined_prefix = 'z' 139 combined_cython_type = 'np.complex128_t' 140if prefix == 's': 141 combined_prefix = 'd' 142 combined_cython_type = 'np.float64_t' 143}} 144 145 146cdef bint _{{prefix}}select1({{cython_type}} * a): 147 return 0 148 149cdef bint _{{prefix}}select2({{cython_type}} * a, {{cython_type}} * b): 150 return 0 151 152cdef int _{{prefix}}solve_discrete_lyapunov({{cython_type}} * a, {{cython_type}} * q, int n, int complex_step=False) except *: 153 # Note: some of this code (esp. the Sylvester solving part) cribbed from 154 # https://raw.githubusercontent.com/scipy/scipy/master/scipy/linalg/_solvers.py 155 156 # Solve an equation of the form $A'XA-X=-Q$ 157 # a: input / output 158 # q: input / output 159 cdef: 160 int i, j 161 int info 162 int inc = 1 163 int n2 = n**2 164 {{if prefix == 's' or prefix == 'c'}} 165 np.float32_t scale = 0.0 166 {{else}} 167 np.float64_t scale = 0.0 168 {{endif}} 169 {{cython_type}} tmp = 0.0 170 {{cython_type}} alpha = 1.0 171 {{cython_type}} beta = 0.0 172 {{cython_type}} delta = -2.0 173 char trans 174 cdef np.npy_intp dim[2] 175 cdef {{cython_type}} [::1,:] apI, capI, u, v 176 cdef int [::1,:] ipiv 177 # Dummy selection function, will not actually be referenced since we do not 178 # need to order the eigenvalues in the ?gees call. 179 cdef: 180 int sdim 181 int lwork = 3*n 182 bint bwork 183 cdef np.npy_intp dim1[1] 184 cdef {{cython_type}} [::1,:] work 185 cdef {{cython_type}} [:] wr 186 {{if prefix == 's' or prefix == 'c'}} 187 cdef np.float32_t [:] wi 188 {{else}} 189 cdef np.float64_t [:] wi 190 {{endif}} 191 192 # Initialize arrays 193 dim[0] = n; dim[1] = n; 194 apI = np.PyArray_ZEROS(2, dim, {{typenum}}, FORTRAN) 195 capI = np.PyArray_ZEROS(2, dim, {{typenum}}, FORTRAN) 196 u = np.PyArray_ZEROS(2, dim, {{typenum}}, FORTRAN) 197 v = np.PyArray_ZEROS(2, dim, {{typenum}}, FORTRAN) 198 ipiv = np.PyArray_ZEROS(2, dim, np.NPY_INT32, FORTRAN) 199 200 dim1[0] = n; 201 wr = np.PyArray_ZEROS(1, dim1, {{typenum}}, FORTRAN) 202 {{if prefix == 's'}} 203 wi = np.PyArray_ZEROS(1, dim1, {{typenum}}, FORTRAN) 204 {{else}} 205 wi = np.PyArray_ZEROS(1, dim1, np.NPY_FLOAT64, FORTRAN) 206 {{endif}} 207 #vs = np.PyArray_ZEROS(2, dim, {{typenum}}, FORTRAN) 208 dim[0] = lwork; dim[1] = lwork; 209 work = np.PyArray_ZEROS(2, dim, {{typenum}}, FORTRAN) 210 211 # - Solve for b.conj().transpose() -------- 212 213 # Get apI = a + I (stored in apI) 214 # = (a + eye) 215 # For: c = 2*np.dot(np.dot(inv(a + eye), q), aHI_inv) 216 copy(n2, a, &apI[0, 0]) 217 # (for loop below adds the identity) 218 219 # Get conj(a) + I (stored in capI) 220 # a^H + I -> capI 221 # For: aHI_inv = inv(aH + eye) 222 copy(n2, a, &capI[0, 0]) 223 # (for loop below adds the identity) 224 225 # Get conj(a) - I (stored in a) 226 # a^H - I -> a 227 # For: b = np.dot(aH - eye, aHI_inv) 228 # (for loop below subtracts the identity) 229 230 # Add / subtract identity matrix 231 for i in range(n): 232 apI[i,i] = apI[i,i] + 1 # apI -> a + eye 233 capI[i,i] = capI[i,i] + 1 # aH + eye 234 a[i + i*n] = a[i + i*n] - 1 # a - eye 235 236 # Solve [conj(a) + I] b' = [conj(a) - I] (result stored in a) 237 # For: b = np.dot(aH - eye, aHI_inv) 238 # Where: aHI_inv = inv(aH + eye) 239 # where b = (a^H - eye) (a^H + eye)^{-1} 240 # or b^H = (a + eye)^{-1} (a - eye) 241 # or (a + eye) b^H = (a - eye) 242 lapack.{{prefix}}getrf(&n, &n, &capI[0,0], &n, &ipiv[0,0], &info) 243 244 if not info == 0: 245 raise np.linalg.LinAlgError('LU decomposition error.') 246 247 lapack.{{prefix}}getrs("N", &n, &n, &capI[0,0], &n, &ipiv[0,0], 248 a, &n, &info) 249 250 if not info == 0: 251 raise np.linalg.LinAlgError('LU solver error.') 252 253 # Now we have b^H; we could take the conjugate transpose to get b, except 254 # that the input to the continuous Lyapunov equation is exactly 255 # b^H, so we already have the quantity we need. 256 257 # - Solve for (-c) -------- 258 259 # where c = 2*np.dot(np.dot(inv(a + eye), q), aHI_inv) 260 # = 2*(a + eye)^{-1} q (a^H + eye)^{-1} 261 # and with q Hermitian 262 # consider x = (a + eye)^{-1} q (a^H + eye)^{-1} 263 # this can be done in two linear solving steps: 264 # 1. consider y = q (a^H + eye)^{-1} 265 # or y (a^H + eye) = q 266 # or (a^H + eye)^H y^H = q^H 267 # or (a + eye) y^H = q 268 # 2. Then consider x = (a + eye)^{-1} y 269 # or (a + eye) x = y 270 271 # Solve [conj(a) + I] tmp' = q (result stored in q) 272 # For: y = q (a^H + eye)^{-1} => (a + eye) y^H = q 273 lapack.{{prefix}}getrs("N", &n, &n, &capI[0,0], &n, &ipiv[0,0], 274 q, &n, &info) 275 276 if not info == 0: 277 raise np.linalg.LinAlgError('LU solver error.') 278 279 # Replace the result (stored in q) with its (conjugate) transpose 280 for j in range(1, n): 281 for i in range(j): 282 tmp = q[i + j*n] 283 q[i + j*n] = q[j + i*n] 284 q[j + i*n] = tmp 285 286 {{if combined_prefix == 'z'}} 287 if not complex_step: 288 for i in range(n2): 289 q[i] = q[i] - q[i].imag * 2.0j 290 {{endif}} 291 292 lapack.{{prefix}}getrs("N", &n, &n, &capI[0,0], &n, &ipiv[0,0], 293 q, &n, &info) 294 295 if not info == 0: 296 raise np.linalg.LinAlgError('LU solver error.') 297 298 # q -> -2.0 * q 299 blas.{{prefix}}scal(&n2, &delta, q, &inc) 300 301 # - Solve continuous time Lyapunov -------- 302 303 # Now solve the continuous time Lyapunov equation (AX + XA^H = Q), on the 304 # transformed inputs ... 305 306 # ... which requires solving the continuous time Sylvester equation 307 # (AX + XB = Q) where B = A^H 308 309 # Compute the real Schur decomposition of a (unordered) 310 # TODO compute the optimal lwork rather than always using 3*n 311 {{if combined_prefix == 'd'}} 312 # a is now the Schur form of A; (r) 313 # u is now the unitary Schur transformation matrix for A; (u) 314 # In the usual case, we will also have: 315 # r = s, so s is also stored in a 316 # u = v, so v is also stored in u 317 # In the complex-step case, we will instead have: 318 # r = s.conj() 319 # u = v.conj() 320 lapack.{{prefix}}gees("V", "N", <lapack.{{prefix}}select2 *> &_{{prefix}}select2, &n, 321 a, &n, 322 &sdim, 323 &wr[0], &wi[0], 324 &u[0,0], &n, 325 &work[0,0], &lwork, 326 &bwork, &info) 327 {{else}} 328 lapack.{{prefix}}gees("V", "N", <lapack.{{prefix}}select1 *> &_{{prefix}}select1, &n, 329 a, &n, 330 &sdim, 331 &wr[0], 332 &u[0,0], &n, 333 &work[0,0], &lwork, 334 &wi[0], 335 &bwork, &info) 336 {{endif}} 337 338 if not info == 0: 339 raise np.linalg.LinAlgError('Schur decomposition solver error.') 340 341 # Get v (so that in the complex step case we can take the conjugate) 342 copy(n2, &u[0, 0], &v[0, 0]) 343 # If complex step, take the conjugate 344 {{if combined_prefix == 'z'}} 345 if complex_step: 346 for i in range(n): 347 for j in range(n): 348 v[i,j] = v[i,j] - v[i,j].imag * 2.0j 349 {{endif}} 350 351 # Construct f = u^H*q*u (result overwrites q) 352 # In the usual case, v = u 353 # In the complex step case, v = u.conj() 354 blas.{{prefix}}gemm("N", "N", &n, &n, &n, 355 &alpha, q, &n, 356 &v[0,0], &n, 357 &beta, &capI[0,0], &n) 358 blas.{{prefix}}gemm("C", "N", &n, &n, &n, 359 &alpha, &u[0,0], &n, 360 &capI[0,0], &n, 361 &beta, q, &n) 362 363 # DTRYSL Solve op(A)*X + X*op(B) = scale*C which is here: 364 # r*X + X*r = scale*q 365 # results overwrite q 366 copy(n2, a, &apI[0, 0]) 367 {{if combined_prefix == 'z'}} 368 if complex_step: 369 for i in range(n): 370 for j in range(n): 371 apI[j,i] = apI[j,i] - apI[j,i].imag * 2.0j 372 {{endif}} 373 lapack.{{prefix}}trsyl("N", "C", &inc, &n, &n, 374 a, &n, 375 &apI[0,0], &n, 376 q, &n, 377 &scale, &info) 378 379 # Scale q by scale 380 if not scale == 1.0: 381 blas.{{prefix}}scal(&n2, <{{cython_type}}*> &scale, q, &inc) 382 383 # Calculate the solution: u * q * v^H (results overwrite q) 384 # In the usual case, v = u 385 # In the complex step case, v = u.conj() 386 blas.{{prefix}}gemm("N", "C", &n, &n, &n, 387 &alpha, q, &n, 388 &v[0,0], &n, 389 &beta, &capI[0,0], &n) 390 blas.{{prefix}}gemm("N", "N", &n, &n, &n, 391 &alpha, &u[0,0], &n, 392 &capI[0,0], &n, 393 &beta, q, &n) 394 395 396cpdef _{{prefix}}compute_coefficients_from_multivariate_pacf({{cython_type}} [::1,:] partial_autocorrelations, 397 {{cython_type}} [::1,:] error_variance, 398 int transform_variance, 399 int order, int k_endog): 400 """ 401 Notes 402 ----- 403 404 This uses the ?trmm BLAS functions which are not available in 405 Scipy v0.11.0 406 """ 407 # Constants 408 cdef: 409 {{cython_type}} alpha = 1.0 410 {{cython_type}} beta = 0.0 411 {{cython_type}} gamma = -1.0 412 int k_endog2 = k_endog**2 413 int k_endog_order = k_endog * order 414 int k_endog_order1 = k_endog * (order+1) 415 int info, s, k 416 # Local variables 417 cdef: 418 np.npy_intp dim2[2] 419 {{cython_type}} [::1, :] initial_variance 420 {{cython_type}} [::1, :] forward_variance 421 {{cython_type}} [::1, :] backward_variance 422 {{cython_type}} [::1, :] autocovariances 423 {{cython_type}} [::1, :] forwards1 424 {{cython_type}} [::1, :] forwards2 425 {{cython_type}} [::1, :] backwards1 426 {{cython_type}} [::1, :] backwards2 427 {{cython_type}} [::1, :] forward_factors 428 {{cython_type}} [::1, :] backward_factors 429 {{cython_type}} [::1, :] tmp 430 {{cython_type}} [::1, :] tmp2 431 # Pointers 432 cdef: 433 {{cython_type}} * forwards 434 {{cython_type}} * prev_forwards 435 {{cython_type}} * backwards 436 {{cython_type}} * prev_backwards 437 # ?trmm 438 # cdef {{prefix}}trmm_t *{{prefix}}trmm = <{{prefix}}trmm_t*>Capsule_AsVoidPtr(blas.{{prefix}}trmm._cpointer) 439 440 # dim2[0] = self.k_endog; dim2[1] = storage; 441 # self.forecast = np.PyArray_ZEROS(2, dim2, {{typenum}}, FORTRAN) 442 443 # If we want to keep the provided variance but with the constrained 444 # coefficient matrices, we need to make a copy here, and then after the 445 # main loop we will transform the coefficients to match the passed variance 446 if not transform_variance: 447 initial_variance = np.asfortranarray(error_variance.copy()) 448 # Need to make the input variance large enough that the recursions 449 # do not lead to zero-matrices due to roundoff error, which would case 450 # exceptions from the Cholesky decompositions. 451 # Note that this will still not always ensure positive definiteness, 452 # and for k_endog, order large enough an exception may still be raised 453 error_variance = np.asfortranarray(np.eye(k_endog, dtype={{dtype}}) * (order + k_endog)**10) 454 455 # Initialize matrices 456 dim2[0] = k_endog; dim2[1] = k_endog; 457 forward_variance = np.PyArray_ZEROS(2, dim2, {{typenum}}, FORTRAN) 458 backward_variance = np.PyArray_ZEROS(2, dim2, {{typenum}}, FORTRAN) 459 forward_factors = np.PyArray_ZEROS(2, dim2, {{typenum}}, FORTRAN) 460 backward_factors = np.PyArray_ZEROS(2, dim2, {{typenum}}, FORTRAN) 461 tmp = np.PyArray_ZEROS(2, dim2, {{typenum}}, FORTRAN) 462 tmp2 = np.PyArray_ZEROS(2, dim2, {{typenum}}, FORTRAN) 463 464 dim2[0] = k_endog; dim2[1] = k_endog_order; 465 # \phi_{s,k}, s = 1, ..., p 466 # k = 1, ..., s+1 467 forwards1 = np.PyArray_ZEROS(2, dim2, {{typenum}}, FORTRAN) 468 forwards2 = np.PyArray_ZEROS(2, dim2, {{typenum}}, FORTRAN) 469 # \phi_{s,k}^* 470 backwards1 = np.PyArray_ZEROS(2, dim2, {{typenum}}, FORTRAN) 471 backwards2 = np.PyArray_ZEROS(2, dim2, {{typenum}}, FORTRAN) 472 473 dim2[0] = k_endog; dim2[1] = k_endog_order1; 474 autocovariances = np.PyArray_ZEROS(2, dim2, {{typenum}}, FORTRAN) 475 476 copy(k_endog2, &error_variance[0, 0], &forward_variance[0, 0]) # \Sigma_s 477 copy(k_endog2, &error_variance[0, 0], &backward_variance[0, 0]) # \Sigma_s^*, s = 0, ..., p 478 copy(k_endog2, &error_variance[0, 0], &autocovariances[0, 0]) # \Gamma_s 479 480 # error_variance_factor = linalg.cholesky(error_variance, lower=True) 481 copy(k_endog2, &error_variance[0, 0], &forward_factors[0,0]) 482 lapack.{{prefix}}potrf("L", &k_endog, &forward_factors[0,0], &k_endog, &info) 483 copy(k_endog2, &forward_factors[0, 0], &backward_factors[0, 0]) 484 485 # We fill in the entries as follows: 486 # [1,1] 487 # [2,2], [2,1] 488 # [3,3], [3,1], [3,2] 489 # ... 490 # [p,p], [p,1], ..., [p,p-1] 491 # the last row, correctly ordered, is then used as the coefficients 492 for s in range(order): # s = 0, ..., p-1 493 if s % 2 == 0: 494 forwards = &forwards1[0, 0] 495 prev_forwards = &forwards2[0, 0] 496 backwards = &backwards1[0, 0] 497 prev_backwards = &backwards2[0, 0] 498 else: 499 forwards = &forwards2[0, 0] 500 prev_forwards = &forwards1[0, 0] 501 backwards = &backwards2[0, 0] 502 prev_backwards = &backwards1[0, 0] 503 504 # Create the "last" (k = s+1) matrix 505 # Note: this is for k = s+1. However, below we then have to fill 506 # in for k = 1, ..., s in order. 507 # P L*^{-1} = x 508 # x L* = P 509 # L*' x' = P' 510 # forwards[:, s*k_endog:(s+1)*k_endog] = np.dot( 511 # forward_factors, 512 # linalg.solve_triangular( 513 # backward_factors, partial_autocorrelations[:, s*k_endog:(s+1)*k_endog].T, 514 # lower=True, trans='T').T 515 # ) 516 for k in range(k_endog): 517 copy(k_endog, &partial_autocorrelations[k, s*k_endog], &tmp[0, k], k_endog) 518 lapack.{{prefix}}trtrs("L", "T", "N", &k_endog, &k_endog, &backward_factors[0,0], &k_endog, 519 &tmp[0, 0], &k_endog, &info) 520 # {{prefix}}gemm("N", "T", &k_endog, &k_endog, &k_endog, 521 # &alpha, &forward_factors[0,0], &k_endog, 522 # &tmp[0, 0], &k_endog, 523 # &beta, &forwards[s*k_endog2], &k_endog) 524 blas.{{prefix}}trmm("R", "L", "T", "N", &k_endog, &k_endog, 525 &alpha, &forward_factors[0,0], &k_endog, 526 &tmp[0, 0], &k_endog) 527 for k in range(k_endog): 528 copy(k_endog, &tmp[k, 0], &forwards[s*k_endog2 + k*k_endog], k_endog) 529 530 # P' L^{-1} = x 531 # x L = P' 532 # L' x' = P 533 # backwards[:, s*k_endog:(s+1)*k_endog] = np.dot( 534 # backward_factors, 535 # linalg.solve_triangular( 536 # forward_factors, partial_autocorrelations[:, s*k_endog:(s+1)*k_endog], 537 # lower=True, trans='T').T 538 # ) 539 copy(k_endog2, &partial_autocorrelations[0, s*k_endog], &tmp[0, 0]) 540 lapack.{{prefix}}trtrs("L", "T", "N", &k_endog, &k_endog, &forward_factors[0, 0], &k_endog, 541 &tmp[0, 0], &k_endog, &info) 542 # {{prefix}}gemm("N", "T", &k_endog, &k_endog, &k_endog, 543 # &alpha, &backward_factors[0, 0], &k_endog, 544 # &tmp[0, 0], &k_endog, 545 # &beta, &backwards[s * k_endog2], &k_endog) 546 blas.{{prefix}}trmm("R", "L", "T", "N", &k_endog, &k_endog, 547 &alpha, &backward_factors[0,0], &k_endog, 548 &tmp[0, 0], &k_endog) 549 for k in range(k_endog): 550 copy(k_endog, &tmp[k, 0], &backwards[s*k_endog2 + k*k_endog], k_endog) 551 552 # Update the variance 553 # Note: if s >= 1, this will be further updated in the for loop 554 # below 555 # Also, this calculation will be re-used in the forward variance 556 # tmp = np.dot(forwards[:, s*k_endog:(s+1)*k_endog], backward_variance) 557 # tmpT = np.dot(backward_variance.T, forwards[:, s*k_endog:(s+1)*k_endog].T) 558 blas.{{prefix}}gemm("T", "T", &k_endog, &k_endog, &k_endog, 559 &alpha, &backward_variance[0, 0], &k_endog, 560 &forwards[s * k_endog2], &k_endog, 561 &beta, &tmp[0, 0], &k_endog) 562 # autocovariances[:, (s+1)*k_endog:(s+2)*k_endog] = tmp.copy().T 563 copy(k_endog2, &tmp[0, 0], &autocovariances[0, (s+1)*k_endog]) 564 565 # Create the remaining k = 1, ..., s matrices, 566 # only has an effect if s >= 1 567 for k in range(s): 568 # forwards[:, k*k_endog:(k+1)*k_endog] = ( 569 # prev_forwards[:, k*k_endog:(k+1)*k_endog] - 570 # np.dot( 571 # forwards[:, s*k_endog:(s+1)*k_endog], 572 # prev_backwards[:, (s-k-1)*k_endog:(s-k)*k_endog] 573 # ) 574 # ) 575 copy(k_endog2, &prev_forwards[k * k_endog2], &forwards[k * k_endog2]) 576 blas.{{prefix}}gemm("N", "N", &k_endog, &k_endog, &k_endog, 577 &gamma, &forwards[s * k_endog2], &k_endog, 578 &prev_backwards[(s - k - 1) * k_endog2], &k_endog, 579 &alpha, &forwards[k * k_endog2], &k_endog) 580 581 # backwards[:, k*k_endog:(k+1)*k_endog] = ( 582 # prev_backwards[:, k*k_endog:(k+1)*k_endog] - 583 # np.dot( 584 # backwards[:, s*k_endog:(s+1)*k_endog], 585 # prev_forwards[:, (s-k-1)*k_endog:(s-k)*k_endog] 586 # ) 587 # ) 588 copy(k_endog2, &prev_backwards[k * k_endog2], &backwards[k * k_endog2]) 589 blas.{{prefix}}gemm("N", "N", &k_endog, &k_endog, &k_endog, 590 &gamma, &backwards[s * k_endog2], &k_endog, 591 &prev_forwards[(s - k - 1) * k_endog2], &k_endog, 592 &alpha, &backwards[k * k_endog2], &k_endog) 593 594 # autocovariances[:, (s+1)*k_endog:(s+2)*k_endog] += np.dot( 595 # autocovariances[:, (k+1)*k_endog:(k+2)*k_endog], 596 # prev_forwards[:, (s-k-1)*k_endog:(s-k)*k_endog].T 597 # ) 598 blas.{{prefix}}gemm("N", "T", &k_endog, &k_endog, &k_endog, 599 &alpha, &autocovariances[0, (k+1)*k_endog], &k_endog, 600 &prev_forwards[(s - k - 1) * k_endog2], &k_endog, 601 &alpha, &autocovariances[0, (s+1)*k_endog], &k_endog) 602 603 # Create forward and backwards variances 604 # backward_variance = ( 605 # backward_variance - 606 # np.dot( 607 # np.dot(backwards[:, s*k_endog:(s+1)*k_endog], forward_variance), 608 # backwards[:, s*k_endog:(s+1)*k_endog].T 609 # ) 610 # ) 611 blas.{{prefix}}gemm("N", "N", &k_endog, &k_endog, &k_endog, 612 &alpha, &backwards[s * k_endog2], &k_endog, 613 &forward_variance[0, 0], &k_endog, 614 &beta, &tmp2[0, 0], &k_endog) 615 blas.{{prefix}}gemm("N", "T", &k_endog, &k_endog, &k_endog, 616 &gamma, &tmp2[0, 0], &k_endog, 617 &backwards[s * k_endog2], &k_endog, 618 &alpha, &backward_variance[0, 0], &k_endog) 619 # forward_variance = ( 620 # forward_variance - 621 # np.dot(tmp, forwards[:, s*k_endog:(s+1)*k_endog].T) 622 # ) 623 # forward_variance = ( 624 # forward_variance - 625 # np.dot(tmpT.T, forwards[:, s*k_endog:(s+1)*k_endog].T) 626 # ) 627 blas.{{prefix}}gemm("T", "T", &k_endog, &k_endog, &k_endog, 628 &gamma, &tmp[0, 0], &k_endog, 629 &forwards[s * k_endog2], &k_endog, 630 &alpha, &forward_variance[0, 0], &k_endog) 631 632 # Cholesky factors 633 # forward_factors = linalg.cholesky(forward_variance, lower=True) 634 # backward_factors = linalg.cholesky(backward_variance, lower=True) 635 copy(k_endog2, &forward_variance[0, 0], &forward_factors[0, 0]) 636 lapack.{{prefix}}potrf("L", &k_endog, &forward_factors[0,0], &k_endog, &info) 637 copy(k_endog2, &backward_variance[0, 0], &backward_factors[0, 0]) 638 lapack.{{prefix}}potrf("L", &k_endog, &backward_factors[0,0], &k_endog, &info) 639 640 641 # If we do not want to use the transformed variance, we need to 642 # adjust the constrained matrices, as presented in Lemma 2.3, see above 643 if not transform_variance: 644 if order % 2 == 0: 645 forwards = &forwards2[0,0] 646 else: 647 forwards = &forwards1[0,0] 648 649 # Here, we need to construct T such that: 650 # variance = T * initial_variance * T' 651 # To do that, consider the Cholesky of variance (L) and 652 # input_variance (M) to get: 653 # L L' = T M M' T' = (TM) (TM)' 654 # => L = T M 655 # => L M^{-1} = T 656 # initial_variance_factor = np.linalg.cholesky(initial_variance) 657 # L' 658 lapack.{{prefix}}potrf("U", &k_endog, &initial_variance[0,0], &k_endog, &info) 659 # transformed_variance_factor = np.linalg.cholesky(variance) 660 # M' 661 copy(k_endog2, &forward_variance[0, 0], &tmp[0, 0]) 662 lapack.{{prefix}}potrf("U", &k_endog, &tmp[0,0], &k_endog, &info) 663 # {{prefix}}potri("L", &k_endog, &tmp[0,0], &k_endog, &info) 664 665 # We need to zero out the lower triangle of L', because ?trtrs only 666 # knows that M' is upper triangular 667 for s in range(k_endog - 1): # column 668 for k in range(s+1, k_endog): # row 669 initial_variance[k, s] = 0 670 671 # Note that T is lower triangular 672 # L M^{-1} = T 673 # M' T' = L' 674 # transform = np.dot(initial_variance_factor, 675 # np.linalg.inv(transformed_variance_factor)) 676 lapack.{{prefix}}trtrs("U", "N", "N", &k_endog, &k_endog, &tmp[0,0], &k_endog, 677 &initial_variance[0, 0], &k_endog, &info) 678 # Now: 679 # initial_variance = T' 680 681 for s in range(order): 682 # forwards[:, s*k_endog:(s+1)*k_endog] = ( 683 # np.dot( 684 # np.dot(transform, forwards[:, s*k_endog:(s+1)*k_endog]), 685 # inv_transform 686 # ) 687 # ) 688 # TF T^{-1} = x 689 # TF = x T 690 # (TF)' = T' x' 691 692 # Get TF 693 copy(k_endog2, &forwards[s * k_endog2], &tmp2[0, 0]) 694 blas.{{prefix}}trmm("L", "U", "T", "N", &k_endog, &k_endog, 695 &alpha, &initial_variance[0, 0], &k_endog, 696 &tmp2[0, 0], &k_endog) 697 for k in range(k_endog): 698 copy(k_endog, &tmp2[k, 0], &tmp[0, k], k_endog) 699 # Get x' 700 lapack.{{prefix}}trtrs("U", "N", "N", &k_endog, &k_endog, &initial_variance[0,0], &k_endog, 701 &tmp[0, 0], &k_endog, &info) 702 # Get x 703 for k in range(k_endog): 704 copy(k_endog, &tmp[k, 0], &forwards[s * k_endog2 + k*k_endog], k_endog) 705 706 707 if order % 2 == 0: 708 return forwards2, forward_variance 709 else: 710 return forwards1, forward_variance 711 712cpdef _{{prefix}}constrain_sv_less_than_one({{cython_type}} [::1,:] unconstrained, int order, int k_endog): 713 """ 714 Transform arbitrary matrices to matrices with singular values less than 715 one. 716 717 Corresponds to Lemma 2.2 in Ansley and Kohn (1986). See 718 `constrain_stationary_multivariate` for more details. 719 """ 720 # Constants 721 cdef: 722 {{cython_type}} alpha = 1.0 723 int k_endog2 = k_endog**2 724 int info, i 725 # Local variables 726 cdef: 727 np.npy_intp dim2[2] 728 {{cython_type}} [::1, :] constrained 729 {{cython_type}} [::1, :] tmp 730 {{cython_type}} [::1, :] eye 731 732 dim2[0] = k_endog; dim2[1] = k_endog * order; 733 constrained = np.PyArray_ZEROS(2, dim2, {{typenum}}, FORTRAN) 734 dim2[0] = k_endog; dim2[1] = k_endog; 735 tmp = np.PyArray_ZEROS(2, dim2, {{typenum}}, FORTRAN) 736 eye = np.PyArray_ZEROS(2, dim2, {{typenum}}, FORTRAN) 737 738 eye = np.asfortranarray(np.eye(k_endog, dtype={{dtype}})) 739 for i in range(order): 740 copy(k_endog2, &eye[0, 0], &tmp[0, 0]) 741 blas.{{prefix}}gemm("N", "T", &k_endog, &k_endog, &k_endog, 742 &alpha, &unconstrained[0, i*k_endog], &k_endog, 743 &unconstrained[0, i*k_endog], &k_endog, 744 &alpha, &tmp[0, 0], &k_endog) 745 lapack.{{prefix}}potrf("L", &k_endog, &tmp[0, 0], &k_endog, &info) 746 747 copy(k_endog2, &unconstrained[0, i*k_endog], &constrained[0, i*k_endog]) 748 # constrained.append(linalg.solve_triangular(B, A, lower=lower)) 749 lapack.{{prefix}}trtrs("L", "N", "N", &k_endog, &k_endog, &tmp[0, 0], &k_endog, 750 &constrained[0, i*k_endog], &k_endog, &info) 751 return constrained 752 753cdef int _{{prefix}}ldl({{cython_type}} * A, int n) except *: 754 # See Golub and Van Loan, Algorithm 4.1.2 755 cdef: 756 int info = 0 757 int j, i, k 758 np.npy_intp dim[1] 759 np.float64_t tol = 1e-15 760 {{cython_type}} [:] v 761 762 dim[0] = n 763 v = np.PyArray_ZEROS(1, dim, {{typenum}}, FORTRAN) 764 765 for j in range(n): 766 # Compute v(1:j) 767 v[j] = A[j + j*n] 768 769 # Positive definite element: use Golub and Van Loan algorithm 770 if v[j].real < -tol: 771 info = -j 772 break 773 elif v[j].real > tol: 774 for i in range(j): 775 v[i] = A[j + i*n] * A[i + i*n] 776 v[j] = v[j] - A[j + i*n] * v[i] 777 778 # Store d(j) and compute L(j+1:n,j) 779 A[j + j*n] = v[j] 780 for i in range(j+1, n): 781 for k in range(j): 782 A[i + j*n] = A[i + j*n] - A[i + k*n] * v[k] 783 A[i + j*n] = A[i + j*n] / v[j] 784 # Positive semi-definite element: zero the appropriate column 785 else: 786 info = 1 787 for i in range(j, n): 788 A[i + j*n] 789 790 return info 791 792cpdef int {{prefix}}ldl({{cython_type}} [::1, :] A) except *: 793 _{{prefix}}ldl(&A[0,0], A.shape[0]) 794 795cdef int _{{prefix}}reorder_missing_diagonal({{cython_type}} * a, int * missing, int n): 796 """ 797 a is a pointer to an n x n diagonal array A 798 missing is a pointer to an n x 1 array 799 n is the dimension of A 800 """ 801 cdef int i, j, k, nobs 802 803 nobs = n 804 # Construct the non-missing index 805 for i in range(n): 806 nobs = nobs - missing[i] 807 808 # Perform replacement 809 k = nobs-1 810 for i in range(n-1,-1,-1): 811 if not missing[i]: 812 a[i + i*n] = a[k + k*n] 813 k = k - 1 814 else: 815 a[i + i*n] = 0 816 817cdef int _{{prefix}}reorder_missing_submatrix({{cython_type}} * a, int * missing, int n): 818 """ 819 a is a pointer to an n x n array A 820 missing is a pointer to an n x 1 array 821 n is the dimension of A 822 """ 823 cdef int i, j, k, nobs 824 825 _{{prefix}}reorder_missing_rows(a, missing, n, n) 826 _{{prefix}}reorder_missing_cols(a, missing, n, n) 827 828cdef int _{{prefix}}reorder_missing_rows({{cython_type}} * a, int * missing, int n, int m): 829 """ 830 a is a pointer to an n x m array A 831 missing is a pointer to an n x 1 array 832 n is the number of rows of A 833 m is the number of columns of A 834 """ 835 cdef int i, j, k, nobs 836 837 nobs = n 838 # Construct the non-missing index 839 for i in range(n): 840 nobs = nobs - missing[i] 841 842 # Perform replacement 843 k = nobs-1 844 for i in range(n-1,-1,-1): 845 if not missing[i]: 846 swap(m, &a[i], &a[k], n, n) 847 k = k - 1 848 849 850cdef int _{{prefix}}reorder_missing_cols({{cython_type}} * a, int * missing, int n, int m): 851 """ 852 a is a pointer to an n x m array A 853 missing is a pointer to an m x 1 array 854 n is the number of rows of A 855 m is the number of columns of A 856 """ 857 cdef int i, k, nobs 858 859 nobs = m 860 # Construct the non-missing index 861 for i in range(m): 862 nobs = nobs - missing[i] 863 864 # Perform replacement 865 k = nobs-1 866 for i in range(m-1,-1,-1): 867 if not missing[i]: 868 swap(n, &a[i*n], &a[k*n]) 869 k = k - 1 870 871 872cpdef int {{prefix}}reorder_missing_matrix({{cython_type}} [::1, :, :] A, int [::1, :] missing, int reorder_rows, int reorder_cols, int diagonal) except *: 873 cdef int n, m, T, t 874 875 n, m, T = A.shape[0:3] 876 877 if reorder_rows and reorder_cols: 878 if not n == m: 879 raise RuntimeError('Reordering a submatrix requires n = m') 880 if diagonal: 881 for t in range(T): 882 _{{prefix}}reorder_missing_diagonal(&A[0, 0, t], &missing[0, t], n) 883 else: 884 for t in range(T): 885 _{{prefix}}reorder_missing_submatrix(&A[0, 0, t], &missing[0, t], n) 886 elif diagonal: 887 raise RuntimeError('`diagonal` argument only valid with reordering a submatrix') 888 elif reorder_rows: 889 for t in range(T): 890 _{{prefix}}reorder_missing_rows(&A[0, 0, t], &missing[0, t], n, m) 891 elif reorder_cols: 892 for t in range(T): 893 _{{prefix}}reorder_missing_cols(&A[0, 0, t], &missing[0, t], n, m) 894 895 896cpdef int {{prefix}}reorder_missing_vector({{cython_type}} [::1, :] A, int [::1, :] missing) except *: 897 cdef int i, k, t, n, T, nobs 898 899 n, T = A.shape[0:2] 900 901 for t in range(T): 902 _{{prefix}}reorder_missing_rows(&A[0, t], &missing[0, t], n, 1) 903 904 905cdef int _{{prefix}}copy_missing_diagonal({{cython_type}} * a, {{cython_type}} * b, int * missing, int n): 906 """ 907 Copy the non-missing block of diagonal entries 908 909 a is a pointer to an n x n diagonal array A (copy from) 910 b is a pointer to an n x n diagonal array B (copy to) 911 missing is a pointer to an n x 1 array 912 n is the dimension of A, B 913 """ 914 cdef int i, j, k, nobs 915 916 nobs = n 917 # Construct the non-missing index 918 for i in range(n): 919 nobs = nobs - missing[i] 920 921 # Perform replacement 922 k = nobs-1 923 for i in range(nobs): 924 b[i + i*n] = a[i + i*n] 925 926 927cdef int _{{prefix}}copy_missing_submatrix({{cython_type}} * a, {{cython_type}} * b, int * missing, int n): 928 """ 929 Copy the non-missing submatrix 930 931 a is a pointer to an n x n diagonal array A (copy from) 932 b is a pointer to an n x n diagonal array B (copy to) 933 missing is a pointer to an n x 1 array 934 n is the dimension of A, B 935 """ 936 cdef int i, j, nobs 937 938 nobs = n 939 # Construct the non-missing index 940 for i in range(n): 941 nobs = nobs - missing[i] 942 943 # Perform replacement 944 for i in range(nobs): 945 copy(nobs, &a[i*n], &b[i*n]) 946 947 948cdef int _{{prefix}}copy_missing_rows({{cython_type}} * a, {{cython_type}} * b, int * missing, int n, int m): 949 """ 950 a is a pointer to an n x m array A 951 b is a pointer to an n x n diagonal array B (copy to) 952 missing is a pointer to an n x 1 array 953 n is the number of rows of A 954 m is the number of columns of A 955 """ 956 cdef int i, j, k, nobs 957 958 nobs = n 959 # Construct the non-missing index 960 for i in range(n): 961 nobs = nobs - missing[i] 962 963 # Perform replacement 964 for i in range(nobs): 965 copy(m, &a[i], &b[i], n, n) 966 967 968cdef int _{{prefix}}copy_missing_cols({{cython_type}} * a, {{cython_type}} * b, int * missing, int n, int m): 969 """ 970 a is a pointer to an n x m array A 971 b is a pointer to an n x n diagonal array B (copy to) 972 missing is a pointer to an m x 1 array 973 n is the number of rows of A 974 m is the number of columns of A 975 """ 976 cdef int i, j, nobs 977 978 nobs = m 979 # Construct the non-missing index 980 for i in range(m): 981 nobs = nobs - missing[i] 982 983 # Perform replacement 984 for i in range(nobs): 985 copy(n, &a[i*n], &b[i*n]) 986 987 988cpdef int {{prefix}}copy_missing_matrix({{cython_type}} [::1, :, :] A, {{cython_type}} [::1, :, :] B, int [::1, :] missing, int missing_rows, int missing_cols, int diagonal) except *: 989 cdef int n, m, T, t, A_T, A_t = 0, time_varying 990 991 n, m, T = B.shape[0:3] 992 A_T = A.shape[2] 993 time_varying = (A_T == T) 994 995 if missing_rows and missing_cols: 996 if not n == m: 997 raise RuntimeError('Copying a submatrix requires n = m') 998 if diagonal: 999 for t in range(T): 1000 if time_varying: 1001 A_t = t 1002 _{{prefix}}copy_missing_diagonal(&A[0, 0, A_t], &B[0, 0, t], &missing[0, t], n) 1003 else: 1004 for t in range(T): 1005 if time_varying: 1006 A_t = t 1007 _{{prefix}}copy_missing_submatrix(&A[0, 0, A_t], &B[0, 0, t], &missing[0, t], n) 1008 elif diagonal: 1009 raise RuntimeError('`diagonal` argument only valid with copying a submatrix') 1010 elif missing_rows: 1011 for t in range(T): 1012 if time_varying: 1013 A_t = t 1014 _{{prefix}}copy_missing_rows(&A[0, 0, A_t], &B[0, 0, t], &missing[0, t], n, m) 1015 elif missing_cols: 1016 for t in range(T): 1017 if time_varying: 1018 A_t = t 1019 _{{prefix}}copy_missing_cols(&A[0, 0, A_t], &B[0, 0, t], &missing[0, t], n, m) 1020 pass 1021 1022 1023cpdef int {{prefix}}copy_missing_vector({{cython_type}} [::1, :] A, {{cython_type}} [::1, :] B, int [::1, :] missing) except *: 1024 cdef int n, t, T, A_t = 0, A_T 1025 1026 n, T = B.shape[0:2] 1027 A_T = A.shape[1] 1028 time_varying = (A_T == T) 1029 1030 for t in range(T): 1031 if time_varying: 1032 A_t = t 1033 _{{prefix}}copy_missing_rows(&A[0, A_t], &B[0, t], &missing[0, t], n, 1) 1034 1035cdef int _{{prefix}}copy_index_diagonal({{cython_type}} * a, {{cython_type}} * b, int * index, int n): 1036 """ 1037 Copy the non-index block of diagonal entries 1038 1039 a is a pointer to an n x n diagonal array A (copy from) 1040 b is a pointer to an n x n diagonal array B (copy to) 1041 index is a pointer to an n x 1 array 1042 n is the dimension of A, B 1043 """ 1044 cdef int i, j, k, nobs 1045 1046 # Perform replacement 1047 for i in range(n): 1048 if index[i]: 1049 b[i + i*n] = a[i + i*n] 1050 1051 1052cdef int _{{prefix}}copy_index_submatrix({{cython_type}} * a, {{cython_type}} * b, int * index, int n): 1053 """ 1054 Copy the non-index submatrix 1055 1056 a is a pointer to an n x n diagonal array A (copy from) 1057 b is a pointer to an n x n diagonal array B (copy to) 1058 index is a pointer to an n x 1 array 1059 n is the dimension of A, B 1060 """ 1061 _{{prefix}}copy_index_rows(a, b, index, n, n) 1062 _{{prefix}}copy_index_cols(a, b, index, n, n) 1063 1064 1065cdef int _{{prefix}}copy_index_rows({{cython_type}} * a, {{cython_type}} * b, int * index, int n, int m): 1066 """ 1067 a is a pointer to an n x m array A 1068 b is a pointer to an n x n diagonal array B (copy to) 1069 index is a pointer to an n x 1 array 1070 n is the number of rows of A 1071 m is the number of columns of A 1072 """ 1073 cdef int i 1074 1075 # Perform replacement 1076 for i in range(n): 1077 if index[i]: 1078 copy(m, &a[i], &b[i], n, n) 1079 1080 1081cdef int _{{prefix}}copy_index_cols({{cython_type}} * a, {{cython_type}} * b, int * index, int n, int m): 1082 """ 1083 a is a pointer to an n x m array A 1084 b is a pointer to an n x n diagonal array B (copy to) 1085 index is a pointer to an m x 1 array 1086 n is the number of rows of A 1087 m is the number of columns of A 1088 """ 1089 cdef int i, j, k, nobs 1090 1091 # Perform replacement 1092 for i in range(m): 1093 if index[i]: 1094 copy(n, &a[i*n], &b[i*n]) 1095 1096 1097cpdef int {{prefix}}copy_index_matrix({{cython_type}} [::1, :, :] A, {{cython_type}} [::1, :, :] B, int [::1, :] index, int index_rows, int index_cols, int diagonal) except *: 1098 cdef int n, m, T, t, A_T, A_t = 0, time_varying 1099 1100 n, m, T = B.shape[0:3] 1101 A_T = A.shape[2] 1102 time_varying = (A_T == T) 1103 1104 if index_rows and index_cols: 1105 if not n == m: 1106 raise RuntimeError('Copying a submatrix requires n = m') 1107 if diagonal: 1108 for t in range(T): 1109 if time_varying: 1110 A_t = t 1111 _{{prefix}}copy_index_diagonal(&A[0, 0, A_t], &B[0, 0, t], &index[0, t], n) 1112 else: 1113 for t in range(T): 1114 if time_varying: 1115 A_t = t 1116 _{{prefix}}copy_index_submatrix(&A[0, 0, A_t], &B[0, 0, t], &index[0, t], n) 1117 elif diagonal: 1118 raise RuntimeError('`diagonal` argument only valid with copying a submatrix') 1119 elif index_rows: 1120 for t in range(T): 1121 if time_varying: 1122 A_t = t 1123 _{{prefix}}copy_index_rows(&A[0, 0, A_t], &B[0, 0, t], &index[0, t], n, m) 1124 elif index_cols: 1125 for t in range(T): 1126 if time_varying: 1127 A_t = t 1128 _{{prefix}}copy_index_cols(&A[0, 0, A_t], &B[0, 0, t], &index[0, t], n, m) 1129 1130 1131cpdef int {{prefix}}copy_index_vector({{cython_type}} [::1, :] A, {{cython_type}} [::1, :] B, int [::1, :] index) except *: 1132 cdef int n, t, T, A_t = 0, A_T 1133 1134 n, T = B.shape[0:2] 1135 A_T = A.shape[1] 1136 time_varying = (A_T == T) 1137 1138 for t in range(T): 1139 if time_varying: 1140 A_t = t 1141 _{{prefix}}copy_index_rows(&A[0, A_t], &B[0, t], &index[0, t], n, 1) 1142 1143cdef int _{{prefix}}select_cov(int k_states, int k_posdef, int k_states_total, 1144 {{cython_type}} * tmp, 1145 {{cython_type}} * selection, 1146 {{cython_type}} * cov, 1147 {{cython_type}} * selected_cov): 1148 cdef: 1149 int i, k_states2 = k_states**2 1150 {{cython_type}} alpha = 1.0 1151 {{cython_type}} beta = 0.0 1152 1153 # Only need to do something if there is a covariance matrix 1154 # (i.e k_posdof == 0) 1155 if k_posdef > 0: 1156 1157 # #### Calculate selected state covariance matrix 1158 # $Q_t^* = R_t Q_t R_t'$ 1159 # 1160 # Combine a selection matrix and a covariance matrix to get 1161 # a simplified (but possibly singular) "selected" covariance 1162 # matrix (see e.g. Durbin and Koopman p. 43) 1163 1164 # `tmp0` array used here, dimension $(m \times r)$ 1165 1166 # $\\#_0 = 1.0 * R_t Q_t$ 1167 # $(m \times r) = (m \times r) (r \times r)$ 1168 blas.{{prefix}}gemm("N", "N", &k_states, &k_posdef, &k_posdef, 1169 &alpha, selection, &k_states_total, 1170 cov, &k_posdef, 1171 &beta, tmp, &k_states) 1172 # $Q_t^* = 1.0 * \\#_0 R_t'$ 1173 # $(m \times m) = (m \times r) (m \times r)'$ 1174 blas.{{prefix}}gemm("N", "T", &k_states, &k_states, &k_posdef, 1175 &alpha, tmp, &k_states, 1176 selection, &k_states_total, 1177 &beta, selected_cov, &k_states) 1178 else: 1179 for i in range(k_states2): 1180 selected_cov[i] = 0 1181 1182{{endfor}}