SCM

SCM Repository

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

View of /pkg/src/pdIdent.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 10 - (download) (as text) (annotate)
Mon Mar 22 20:20:05 2004 UTC (15 years, 2 months ago) by bates
File size: 1756 byte(s)
Initial import
#include "Mutils.h"

SEXP pdIdent_gradient(SEXP x, SEXP Ain,
		      SEXP nlev)
{
    int nlevVal = asInteger((SEXP)nlev);
    SEXP param = GET_SLOT((SEXP)x, install("param"));
    int* dims = INTEGER(getAttrib((SEXP)Ain, R_DimSymbol));
    int m = dims[0];
    int n = dims[1];
    int q = asInteger(GET_SLOT((SEXP)x, install("Ncol")));
    double* Amat = REAL((TYPEOF((SEXP)Ain) == REALSXP) ? (SEXP)Ain :
                        coerceVector((SEXP)Ain, REALSXP));
    int i, j;
    double sum;
    
    if (q <= 0) {
	error("Uninitialized pdIdent object");
    }
    if (m != n || m != q) {
	error("A must be a %d by %d matrix", q, q);
    }
    if (nlevVal <= 0) {
	error("nlev must by > 0");
    }
    sum = 0.;
    for (j = 0; j < q; j++) {
	for (i = 0; i < q; i++) {
	    sum += Amat[i + j * q] * Amat[i + j * q];
	}
    }
    return ScalarReal((double)(q * nlevVal) -
		      sum * exp(REAL(param)[0] * 2.));
}

SEXP pdIdent_EMupdate(SEXP x, SEXP nlev, SEXP Ain)
{
    int nlevVal = asInteger(nlev);
    SEXP param = GET_SLOT(x, install("param"));
    int* dims = INTEGER(getAttrib(Ain, R_DimSymbol));
    int m = dims[0];
    int n = dims[1];
    int q = asInteger(GET_SLOT(x, install("Ncol")));
    double* Amat = REAL((TYPEOF(Ain) == REALSXP)?duplicate(Ain):
                        coerceVector(Ain, REALSXP));
    int i, j;
    double sum;
    
    if (q <= 0) {
	error("Uninitialized pdIdent object");
    }
    if (m != n || m != q) {
	error("A must be a %d by %d matrix", q, q);
    }
    if (nlevVal <= 0) {
	error("nlev must by > 0");
    }
    sum = 0.;
    for (j = 0; j < q; j++) {
	for (i = 0; i < q; i++) {
	    sum += Amat[i + j * q] * Amat[i + j * q];
	}
    }
    REAL(param)[0] = -log(sum/((double) nlevVal*q))/2;
    return x;
}

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