SCM

SCM Repository

[matrix] Annotation of /pkg/src/pdIdent.c
ViewVC logotype

Annotation of /pkg/src/pdIdent.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 10 - (view) (download) (as text)

1 : bates 10 #include "Mutils.h"
2 :    
3 :     SEXP pdIdent_gradient(SEXP x, SEXP Ain,
4 :     SEXP nlev)
5 :     {
6 :     int nlevVal = asInteger((SEXP)nlev);
7 :     SEXP param = GET_SLOT((SEXP)x, install("param"));
8 :     int* dims = INTEGER(getAttrib((SEXP)Ain, R_DimSymbol));
9 :     int m = dims[0];
10 :     int n = dims[1];
11 :     int q = asInteger(GET_SLOT((SEXP)x, install("Ncol")));
12 :     double* Amat = REAL((TYPEOF((SEXP)Ain) == REALSXP) ? (SEXP)Ain :
13 :     coerceVector((SEXP)Ain, REALSXP));
14 :     int i, j;
15 :     double sum;
16 :    
17 :     if (q <= 0) {
18 :     error("Uninitialized pdIdent object");
19 :     }
20 :     if (m != n || m != q) {
21 :     error("A must be a %d by %d matrix", q, q);
22 :     }
23 :     if (nlevVal <= 0) {
24 :     error("nlev must by > 0");
25 :     }
26 :     sum = 0.;
27 :     for (j = 0; j < q; j++) {
28 :     for (i = 0; i < q; i++) {
29 :     sum += Amat[i + j * q] * Amat[i + j * q];
30 :     }
31 :     }
32 :     return ScalarReal((double)(q * nlevVal) -
33 :     sum * exp(REAL(param)[0] * 2.));
34 :     }
35 :    
36 :     SEXP pdIdent_EMupdate(SEXP x, SEXP nlev, SEXP Ain)
37 :     {
38 :     int nlevVal = asInteger(nlev);
39 :     SEXP param = GET_SLOT(x, install("param"));
40 :     int* dims = INTEGER(getAttrib(Ain, R_DimSymbol));
41 :     int m = dims[0];
42 :     int n = dims[1];
43 :     int q = asInteger(GET_SLOT(x, install("Ncol")));
44 :     double* Amat = REAL((TYPEOF(Ain) == REALSXP)?duplicate(Ain):
45 :     coerceVector(Ain, REALSXP));
46 :     int i, j;
47 :     double sum;
48 :    
49 :     if (q <= 0) {
50 :     error("Uninitialized pdIdent object");
51 :     }
52 :     if (m != n || m != q) {
53 :     error("A must be a %d by %d matrix", q, q);
54 :     }
55 :     if (nlevVal <= 0) {
56 :     error("nlev must by > 0");
57 :     }
58 :     sum = 0.;
59 :     for (j = 0; j < q; j++) {
60 :     for (i = 0; i < q; i++) {
61 :     sum += Amat[i + j * q] * Amat[i + j * q];
62 :     }
63 :     }
64 :     REAL(param)[0] = -log(sum/((double) nlevVal*q))/2;
65 :     return x;
66 :     }

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge