C++ DoubleMatrix::nc方法代码示例

void AkronItimesCTest::testAkronItimesCVal()
{
    const int m = 2;
    const int n = 3;
    const int k = 1;
    int i;
    double seq[] = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30};
    DoubleMatrix A(m,n);
    DoubleMatrix I = identity(n);
    DoubleMatrix C(n*n,k);

    std::copy(seq, seq+m*n, A.data());
    std::copy(seq, seq+n*n*k, C.data());

    DoubleMatrix AIC = AkronItimesC(A,I,C);
    CPPUNIT_ASSERT_EQUAL(m*n, AIC.nr());
    CPPUNIT_ASSERT_EQUAL(k, AIC.nc());

    DoubleMatrix ABC = AkronBtimesC(A,I,C);
    CPPUNIT_ASSERT_EQUAL(m*n, ABC.nr());
    CPPUNIT_ASSERT_EQUAL(k, ABC.nc());
    //std::cout << "ABC=" << ABC << std::endl;
    //std::cout << "AIC=" << AIC << std::endl;
    for( i = 0; i < m*n*k; i++ )
    {
        CPPUNIT_ASSERT_EQUAL(ABC.data()[i], AIC.data()[i]);
    }

}

void symmetrize(const DoubleMatrix& L, DoubleMatrix& Sym)
{
    int m = L.nr();
    int n = L.nc();
    assert(m==n);
    assert(Sym.nr() == m);
    assert(Sym.nc() == n);

    const double* pL = L.data();
    double*       pS = Sym.data();

    if( &L != &Sym )
    {
        std::copy(pL, pL+m*n, pS);
    }
    int j,i;
    for( j=0; j<n; j++ )
    {
        // Copy the lower triangle data into the resulting matrix's upper triangle
        for( i=0; i<j; i++ )
        {
            pS[i+j*m] = pL[j+i*m];
        }        
    }
}

void multiply(const DoubleMatrix &X, const DoubleMatrix &Y, DoubleMatrix& Z)
{
  assert(&X!=&Z);
  assert(&Y!=&Z);
  if( X.isEmpty() || Y.isEmpty() )
    return;

  //
  // Row Major Matrices
  //
  //    C = alpha * A * B + beta * C  --- (1)
  //
  // A : m by k         lda (stride) = k
  // B : k by n         ldb (stride) = n
  // C : m by n         ldc (stride) = n
  //
  // Column Major Matrices
  //
  //    Z = alpha * X * Y + beta * C  --- (2)
  //    Z = C^t 
  //      = alpha * B^t * A^t + beta * C^t  --- (3)
  //
  // X = B^t : n by k   ldx (stride) = n
  // Y = A^t : k by m   ldy (stride) = k
  // Z = C^t : n by m   ldz (stride) = n
  //
  int m = Y.nc();
  int k = X.nc();
  int n = X.nr();

  Z.resize( n, m );
  
  assert( X.nr() == n );
  assert( X.nc() == k );
  assert( Y.nr() == k );
  assert( Y.nc() == m );
  assert( Z.nr() == n );
  assert( Z.nc() == m );
  
  const double * pX = X.data();
  const double * pY = Y.data();
  double       * pZ = Z.data(); 
  
  int lda = n;
  int ldb = k;
  int ldc = n;

  cblas_dgemm( CblasColMajor, CblasNoTrans, CblasNoTrans, n, m, k, ALPHA, pX, lda, pY, ldb , BETA, pZ, ldc );

}

static void replaceIth( DoubleMatrix &target, int ith, const DoubleMatrix &data )
{
    assert(target.nc() == data.nc());
    assert( ith >= 0 && ith < target.nr());

    double *pTarget = target.data();
    const double *pData   = data.data();
    int j;

    for( j=0; j<target.nc(); j++ )
    {
        pTarget[j*target.nr()+ith] = pData[j];
    }
    return;
}

void lambdaTest::test(
    SpkModel<double> &model, 
    const DoubleMatrix &dvecY, 
    const DoubleMatrix &dvecAlp, 
    const DoubleMatrix &dvecB, 
    double lambdaOut, 
    DoubleMatrix &lambda_alpOut, 
    DoubleMatrix &lambda_bOut,
    const bool withD)
{
    DoubleMatrix exactTemp;

    double       dblexactTemp;
    int i;

    lambdaOut     = lambda(model, dvecY, dvecAlp, dvecB, withD);
    lambda_alpOut = lambda_alp(model, dvecY, dvecAlp, dvecB, withD);
    lambda_bOut   = lambda_b(model, dvecY, dvecAlp, dvecB, withD);

    dblexactTemp = funExactLambda(dvecAlp, dvecB, dvecY, withD);
    CPPUNIT_ASSERT_DOUBLES_EQUAL( dblexactTemp, lambdaOut, 0.0001);

    exactTemp = funExactLambda_alp(dvecAlp, dvecB, dvecY, withD);
    for( i=0; i<exactTemp.nr()*exactTemp.nc(); i++ )
    {
        CPPUNIT_ASSERT_DOUBLES_EQUAL(
            exactTemp.data()[i], 
            lambda_alpOut.data()[i], 0.0001
            );
    }
    
    exactTemp = funExactLambda_b(dvecAlp, dvecB, dvecY, withD);
    for( i=0; i<exactTemp.nr()*exactTemp.nc(); i++ )
    {
        CPPUNIT_ASSERT_DOUBLES_EQUAL( 
            exactTemp.data()[i],
            lambda_bOut.data()[i], 0.0001
            );
    }
}

static double compareAgainst(
            const DoubleMatrix &dvecZ,     // column vector
            const DoubleMatrix &dvecH,     // column vector
            const DoubleMatrix &dmatQ,     // symmetric matrix Q(x)
            const DoubleMatrix &dmatInvQ   // inverse of Q
            )
{

    assert( dvecZ.nc() == 1 );
    assert( dvecH.nc() == 1 );
    assert( dmatQ.nr() == dmatQ.nc() );
    assert( hasPosDet( dmatQ ) );

    DoubleMatrix dmatQ2PI( dmatQ.nr(), dmatQ.nc() );
    
	const double *pQ    = dmatQ.data();
    double *pQ2PI = dmatQ2PI.data();

    for( int i=0; i<dmatQ.nr()*dmatQ.nc(); i++ )
        pQ2PI[i] = pQ[i] * 2.0 * PI;

    // Compute det(Q2PI) = b * 2^c .
    double b;
    long c;
    det( dmatQ2PI , &b, &c );

    // Compute log(det(Q2PI)).
    double dTerm1 = log( b ) + c * log( 2.0 );

    DoubleMatrix  dmatR = subtract( dvecZ, dvecH );

    DoubleMatrix  term2 = multiply(multiply(transpose(dmatR),dmatInvQ), dmatR);
    double        dTerm2 = oneByOneToScalar(term2);

    return (dTerm1 + dTerm2) / 2.0;
}

static bool check(DoubleMatrix dmatB,
                  int rows,
                  int cols,
                  int d1, int d2, int d3, 
                  int d4, int d5, int d6,
                  int d7, int d8, int d9)
{
    using namespace std;
    bool isOkay = true;
    int m = dmatB.nr();
    int n = dmatB.nc();
    double *pdB = dmatB.data();

    if( m != rows || n != n ){
        isOkay = false;
        return isOkay;
    }

    int d[9];
    d[0] = d1;
    d[1] = d2;
    d[2] = d3;
    d[3] = d4;
    d[4] = d5;
    d[5] = d6;
    d[6] = d7;
    d[7] = d8;
    d[8] = d9;
    
    for( int i=0; i<rows*cols && isOkay; i++ ){
        if( d[i] != pdB[i] ){
            cout << "d[" << i << "] received was " << d[i] << endl;
            cout << "mat[" << i << "] was " << pdB[i] << endl;
            dmatB.print();
            isOkay = false;
        }
    }
    return isOkay;
}

const DoubleMatrix backDiv(const DoubleMatrix &dmatA, const DoubleMatrix &dmatB)
{
  // A is assumed to be square.
  int m = dmatA.nr();
  int n = dmatA.nc();
  assert( m == n );

  // B is m by l matrix, where l is the number of right hand sides.
  int l = dmatB.nc();
  assert( dmatB.nr() == m );

  if( dmatA.isEmpty() || dmatB.isEmpty() )
    return DoubleMatrix( 0, 0 );

  //==============================================================
  // First decompose A into LU such that A = P * L * U, 
  // where P is the permutation matrix,
  // L is the lower triangle and the U the upper triangle.
  //
  // We use CLAPACK's DGETRF() which does LU decomposition
  // with partial (ie. row interchanges only) pivoting.
  //==============================================================
  
  // enum CBLAS_ORDER order =: (CblasColMajor | CblasRowMajor)
  //
  // If order = CblasColMajor, the array, a, is assumed to 
  // hold each matrix A's column in the contiguous manner
  // in memory (ie. A is said to be in the column major order).
  // If order = CblasRowMajor, the array, a, is assumed to
  // hold each matrix A's row in the contiguous manner
  // in memory (ie. A is said to be in the row major order).
  enum CBLAS_ORDER order = CblasColMajor;

  // double *a
  //
  // (on entry) a points to the elements of matrix A(m,n) 
  // in the column major order if "order" = CblasColMajor, 
  // or in the row major order if "order" = CblasRowMajor.
  //
  // (on exit) The lower triangle (j<=i) is replaced by L
  // and the upper triangle (j>i) is replaced by U.
  double a[m*n];
  copy( dmatA.data(), dmatA.data()+m*n, a );

  // int lda
  // 
  // The leading dimension of A.  
  // If A is in the column major order, lda = m.
  // If A is in the row major order, lda = n.
  int lda = m;

  // int ipiv(m)
  //
  // (on exit) The i-th row in A was interchanged with the row
  // indicated by the value in ipiv[i].
  int ipiv[m];

  int info = clapack_dgetrf( order, m, n, a, lda, ipiv );
  if( info < 0 )
    {
      char mess[ SpkError::maxMessageLen() ];
      snprintf( mess, SpkError::maxMessageLen(), "Solution of a system of linear equations using the LU decomposition failed: \n the %s argument to the function that performs the LU decomposition  had an illegal value.", 
               intToOrdinalString( -info, ONE_IS_FIRST_INT ).c_str() );
      throw SpkException( SpkError::SPK_UNKNOWN_ERR, mess, __LINE__, __FILE__ );
    }
  else if( info > 0 )
    {
      char mess[ SpkError::maxMessageLen() ];
      snprintf( mess, SpkError::maxMessageLen(), "Solution of a system of linear equations using the LU decomposition failed: \nthe %s diagonal element of U is exactly zero.", 
               intToOrdinalString( info, ONE_IS_FIRST_INT ).c_str() );
      throw SpkException( SpkError::SPK_NOT_POS_DEF_ERR, mess, __LINE__, __FILE__ );
    }

  //==============================================================
  // Solve A x = B for x using the LU computed in the previous
  // step.
  // Note that A is now assumed to be square: m = n.
  //==============================================================
  
  // int rhs
  //
  // The number of right hand sides (ie. the number of columns of B).
  int nrhs = l;

  // int ldb
  // The leading dimension of B.
  // If B is in the column major order, ldb = m.
  // If B is in the row major order, ldb = l.
  int ldb = m;

  // double *x
  //
  // (on entry) x points to the elements of B in the column major
  // order if "order" = CblasColMajor or in the row major otherwise.
  // (on exit) x points to the solution matrix, x (m=n by l).
  DoubleMatrix X( dmatB );
  double * x = X.data();// This points to b on entry and contains the solution x upon exit

  info = clapack_dgetrs( order, CblasNoTrans, n, nrhs, a, lda, ipiv, x, ldb );
  if( info < 0 )
//.........这里部分代码省略.........

extern void spk_non_par(
	size_t                             level       ,
	SpkModel< CppAD::AD<double> >     &admodel     ,
	SpkModel<double>                  &model       ,
	const DoubleMatrix                &N           ,
	const DoubleMatrix                &y           ,
	const DoubleMatrix                &max_itr     ,
	const DoubleMatrix                &epsilon     ,
	const DoubleMatrix                &blow        ,
	const DoubleMatrix                &bup         ,
	const DoubleMatrix                &Bin         ,
	DoubleMatrix                      &Bout        ,
	DoubleMatrix                      &lamout      ,
	DoubleMatrix                      &Pout  )
{
	// temporary indices
	size_t i, j, k;

	// temporary double pointer
	double       *ptr;
	const double *ptr_c;

	// number of discrete measure points
	size_t J = Bin.nc();

	// number of random effects
	size_t n = blow.nr();
	
	// ------------ Arguments to non_par::opt_measure --------------------
	assert( max_itr.nr()  == 2 );
	mat2cpp::matrix<size_t> maxitr(2, 1);
	maxitr(0, 0) = size_t( *(max_itr.data() + 0) );
	maxitr(1, 0) = size_t( *(max_itr.data() + 1) );

	assert( epsilon.nr()  == 5 );
	mat2cpp::matrix<double> eps(5, 1);
	eps(0, 0)    = *(epsilon.data() + 0);
	eps(1, 0)    = *(epsilon.data() + 1);
	eps(2, 0)    = *(epsilon.data() + 2);
	eps(3, 0)    = *(epsilon.data() + 3);
	eps(4, 0)    = *(epsilon.data() + 4);

	// input likelihood function
	Like like(admodel, model, N, y, n);

	// input number of individuals in the population
	size_t M = N.nr();

	// input lower limit for the random effects
	mat2cpp::matrix<double> xLow(1, n);
	ptr_c = blow.data();
	for(k = 0; k < n; k++)
		xLow(0, k) = ptr_c[k];

	// input upper limit for the random effects
	mat2cpp::matrix<double> xUp(1, n);
	ptr_c = bup.data();
	for(k = 0; k < n; k++)
		xUp(0, k) = ptr_c[k];

	// input and return discrete measure points
	mat2cpp::matrix<double> X(J, n);
	ptr_c = Bin.data();
	for(j = 0; j < J; j++)
	{	for(k = 0; k < n; k++)
			X(j, k) = ptr_c[k + j * n];
	}

	// return weight corresponding to each measure oint
	mat2cpp::matrix<double> lambda(J, 1);

	// return convergence information
	mat2cpp::matrix<double> info; 

	// return status message
	const char *msg;

	// -----------------------------------------------------------------
	msg = non_par::opt_measure( level, maxitr, eps, 
		&like, M, xLow, xUp, X, lambda, info
	);
	// -----------------------------------------------------------------

	if( strcmp(msg, "ok") != 0 )
	{	throw SpkException(
      			SpkError::SPK_NON_PAR_ERR,
			msg,
      			__LINE__,
      			__FILE__ 
		);
	}
	// determine number of discrete measure points
	assert( n == X.size2() );
	J = X.size1();

	// dimension the return matrices
	Bout.resize(n, J);
	lamout.resize(J, 1);
	Pout.resize(M, J);

//.........这里部分代码省略.........

static double oneByOneToScalar( const DoubleMatrix &dmatA ){
    assert( dmatA.nr() == 1 );
    assert( dmatA.nc() == 1 );

    return (dmatA.data())[0];
}

const DoubleMatrix elsq_x(const DoubleMatrix &dvecR,      // m size vector, z - h
                    const DoubleMatrix &dmatQ,      // m by m symmetric, positive definite
                    const DoubleMatrix &dmatQinv,   // m by m symmetric, positive definite
                    const DoubleMatrix &dmatH_x,    // m by n matrix
                    const DoubleMatrix &dmatQ_x     // m*m by n matrix
                    )
{
    using namespace std;
    // x:   m by 1
    // z:   m by 1
    // h:   m by 1
    // Q:   m by m
    // invQ:m by m
    // h_x: m by n
    // Q_x: m*m by n
    const int m = dvecR.nr();
    const int n = dmatQ_x.nc();
    double  val;

    assert( dmatQ.nr() == m );
    assert( dmatQ.nc() == m );
    assert( dmatQinv.nr() == m );
    assert( dmatQinv.nc() == m );
    assert( dmatH_x.nr() == m );
    assert( dmatH_x.nc() == n );
    assert( dmatQ_x.nr() == m*m );
    assert( dmatQ_x.nc() == n );

    transpose(dvecR, dvecRTrans);
    multiply(dvecRTrans, dmatQinv, drowW);
    transpose( rvec( dmatQinv ), rvecQinvTrans );

    multiply(rvecQinvTrans, dmatQ_x, drowTerm1);
    assert(drowTerm1.nr()== 1);
    assert(drowTerm1.nc()== n);

    multiply(drowW, dmatH_x, drowTerm2);
    assert(drowTerm2.nr()== 1);
    assert(drowTerm2.nc()== n);

    //
    // The loop below untangles the following statement
    // for eliminating matrix object constructions and taking advantage of
    // Q being symmetric.
    //
    //    DoubleMatrix dmatTerm3  = AkronBtimesC( drowW, rvecQinvTrans, dmatQ_x );
    // 
    // Let w be (z-h).
    // 
    // The orignal formula
    //    0.5 [w^T kron w^T] partial_x(Qinv) --- (eq. 1)
    //
    // is equivalent to:
    //    0.5 partial_x(k) [ w^T Qinv w ] --- (eq. 2)
    //   =0.5 SUM [ w(i) w(j) partial_x(k)(Q(i,j)) ], over 1<=i<=m and 1<=j<=m.
    //
    // For Q being symmetric, the following is true:
    //    w(i) w(j) partial_x(k)(Q(i,j)) == w(j) w(i) partial_x(k)(Q(j,i))
    //
    //
    const double* pW   = drowW.data();
    const double* pQ_x = dmatQ_x.data();
    drowTerm3.resize(1,n);
    drowTerm3.fill(0.0);
    double* pTerm3 = drowTerm3.data();

    for( int k=0; k<n; k++ )
    {
        for( int j=0; j<m; j++ )
        {
            for( int i=0; i<m; i++ )
            {
                if( i<=j )
                {
                    val = pW[i] * pW[j] * pQ_x[ j*m+i+k*m*m ] / 2.0;
                    drowTerm3.data()[k] += val;
                    if( i<j )
                        pTerm3[k] += val;
                }
            }
        }
    }
    
    subtract( subtract(mulByScalar(drowTerm1, 0.5), drowTerm2 )
                                        ,drowTerm3, drowAns );


    // elsq_x returns a 1 by n matrix
    assert(drowAns.nr()==1);
    assert(drowAns.nc()==n);
    return drowAns;

}