C++ DoubleMatrix类代码示例

void MatrixHeatImaging(const DoubleMatrix &m, double minimum, double maximum, QPixmap &pixmap, int width, int height)
{
    pixmap = QPixmap(QSize(width, height));
    pixmap.fill(Qt::white);
    QPainter painter(&pixmap);

    unsigned int rows = m.rows();
    unsigned int cols = m.cols();

    for (unsigned int j=0; j<rows; j++)
    {
        for (unsigned int i=0; i<cols; i++)
        {
            double u = m.at(j,i);
            double ratio = 0.0;
            if (minimum!=maximum)
                ratio = 2.0 * (u-minimum) / (maximum - minimum);
            int b = int(MAX(0, 255*(1 - ratio)));
            int r = int(MAX(0, 255*(ratio - 1)));
            int g = 255 - b - r;
            QColor c(r, g, b);
            painter.setPen(c);
            painter.drawPoint(i,height-j-1);
        }
    }
}

static DoubleMatrix funExactLambda_b(const DoubleMatrix &dvecAlp, 
                                     const DoubleMatrix &dvecB, 
                                     const DoubleMatrix &dvecY,
                                     const bool withD){
	
	// Updated 1-8-01 Alyssa
	// fixed for const correctness
    DoubleMatrix term1(1, dvecB.nr()), term2(1, dvecB.nr());
    const double *pdAlp   = dvecAlp.data();
    const double *pdB     = dvecB.data();
    const double *pdY     = dvecY.data();
    double *pdTerm1 = term1.data();
    double *pdTerm2 = term2.data();

    term1.fill(0.0);
    term2.fill(0.0);

    pdTerm1[0] = 1.0/pdB[0]  
                -0.5 * (pow(pdB[0],-2.0)*pow((pdY[0]-pdAlp[1]-pdB[1]),2.0)
                        +pow(pdB[0],-2.0)*pow((pdY[1]-pdAlp[1]-pdB[1]),2.0));
    pdTerm1[1] = -(pdY[0]+pdY[1]-2.0*pdAlp[1]-2.0*pdB[1])/pdB[0] ;

    if( withD ){
        //pdAns[0] = 1.0/pdB[0] - pow((1.0 - pdAlp[1] - pdB[1]*pdB[1]), 2.0) + pdB[0]/pdAlp[0];
        //pdAns[1] = -2.0 * (1.0 - pdAlp[1] - pdB[1]) / pdB[0] + pdB[1]/pdAlp[0];
        pdTerm2[0] = pdB[0]/pdAlp[0];
        pdTerm2[1] = pdB[1]/pdAlp[0];
        return add( term1, term2 );
    }
    else{
        return term1;
    }
}

/*!
 * \brief assertSizeEqual Asserts that the two matrices
 * given as parameters are of equal size.
 * Throws an invalid argument exception if that is not the case.
 * \param A First matrix.
 * \param B Second matrix.
 */
void assertSizeEqual(const DoubleMatrix &A, const DoubleMatrix &B)
{
    if(A.size() != B.size())
    {
        throw std::invalid_argument("Matrices are not of equal size!");
    }
}

MatlabEngine::DoubleMatrix
MatlabEngine::GetMatrix( const string& inExp )
{
    DoubleMatrix result;
    mxArray* ans = GetMxArray( inExp );
    if( ans )
    {
        int numDims = mxGetNumberOfDimensions( ans );
        const int* dims = mxGetDimensions( ans );
        if( numDims != 2 )
            bcierr << "Can only handle two dimensions" << endl;
        result.resize( dims[ 0 ], vector<double>( dims[ 1 ] ) );
        double* value = mxGetPr( ans );
        if( value )
        {
            int indices[] = { 0, 0 };
            for( size_t i = 0; i < result.size(); ++i )
            {
                indices[ 0 ] = i;
                for( size_t j = 0; j < result[ i ].size(); ++j )
                {
                    indices[ 1 ] = j;
                    result[ i ][ j ] = value[ mxCalcSingleSubscript( ans, 2, indices ) ];
                }
            }
        }
        mxDestroyArray( ans );
    }
    return result;
}

void TimeWarp::createCombinedMatrix(DoubleMatrix myChromaMatrix, DoubleVector energyVector, DoubleMatrix* chromaEnergyMatrix){
	chromaEnergyMatrix->clear();
	int sizeRatio = energyVector.size() / myChromaMatrix.size();//
	printf("COMBINE: size of my chroma is %i\n", (int) myChromaMatrix.size());// energyVector.size() / myChromaMatrix.size();
	printf("COMBINED: size ratio of energy to chroma is %i \n", sizeRatio);
	int chromaSize = myChromaMatrix.size();	
//	printf("index is %i\n", index);
	
	for (int i = 0;i < energyVector.size();i++){
		DoubleVector d;
		int index = min(chromaSize-1, (int) floor(i/sizeRatio));

		for (int y = 0;y < 12;y++){
			d.push_back(myChromaMatrix[index][y]);//
		}

		
		d.push_back(energyVector[i]);
					(*chromaEnergyMatrix).push_back(d);
	}
	printf("COMBINED: size of chroma energy is %i\n", (int)(*chromaEnergyMatrix).size());

//	int x = (int)(*chromaEnergyMatrix).size()/3;
//	printf("energy[%i] %f \n", x, energyVector[x]);
/*
 for (int y = 0;y < 13;y++){
		printf("chroma[%i][%i] %f \n", x, y, myChromaMatrix[x/sizeRatio][y]);
		printf("c[%i][%i] %f \n", x, y, (*chromaEnergyMatrix)[x][y]);
	}
		printf("\n");
*/	

}

void TelluriumData::setData(const DoubleMatrix& theData)
{
    mTheData = theData;
	mColumnNames.clear();
	for (int k = 0; k < theData.getColNames().size(); ++k) {
		mColumnNames.add(theData.getColNames().at(k));
	}
    RRPLOG(Logger::LOG_DEBUG) << "Simulation Data =========== \n" << mTheData;
    check();
}

void EXPECT_MATRIX_DOUBLE_EQ(DoubleMatrix A, DoubleMatrix B, string message = ""){
	EXPECT_EQ(A.size(), B.size()) << "Number of lines is different. " << message;
	EXPECT_EQ(A[0].size(), B[0].size()) << "Number of columns is different. " << message;
	unsigned int N = A.size();
	unsigned int M = A[0].size();
	for (unsigned int i = 0; i < N; i++){
		for (unsigned int j = 0; j < M; j++)
			EXPECT_FLOAT_EQ(A[i][j], B[i][j]) << message;
	}
}

void StopRule::multiple (DoubleMatrix &mat1, DoubleVector &vec2, DoubleVector &proVec) {
	int row_, col_;
	proVec.resize(mat1.size());

	for (row_ = 0; row_ < mat1.size (); row_ ++) {
		proVec[row_] = 0.0;
		for (col_ = 0; col_ < mat1[0].size(); col_ ++)
			proVec[row_] += mat1[row_][col_] * vec2[col_];
	}
}

double
MatlabEngine::GetScalar( const string& inExp )
{
    DoubleMatrix result = GetMatrix( inExp );
    if( result.empty() || result[ 0 ].empty() )
    {
        bcierr << "Could not get scalar value \"" << inExp << "\"" << endl;
        result.resize( 1, vector<double>( 1 ) );
    }
    return result[ 0 ][ 0 ];
}

// T^ijk = T^ijl M_lk
Tensor Tensor::operator*(const DoubleMatrix & M) const {
  if (M.displayRows() != 3 || M.displayCols() !=3) {
    ostringstream ii;
    ii << "Tensor " << *this << " * matrix" << M << "of wrong size.\n";
    throw ii.str();
  }
  Tensor T;
  int i; 
  for (i=1; i<=3; i++) T(i) = display(i) * M;
  return T;
}

void lmWorker::calculateCovariance()
{
    //Check if Hessain is invertible
    DoubleMatrix mat = mTheHost.mHessian.getValue();

    ls::ComplexMatrix temp(mat); //Get a complex matrix from a double one. Imag part is zero
    ls::ComplexMatrix Inv = GetInverse(temp);

    DoubleMatrix temp2(mat.RSize(), mat.CSize());
    temp2 = Inv;

    mTheHost.mCovarianceMatrix.setValue(temp2);
}

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]);
    }

}

// Does T^ijk = T^ljk M_li
Tensor Tensor::product(const DoubleMatrix & M) const {
  if (M.displayRows() != 3 || M.displayCols() !=3) {
    ostringstream ii;
    ii << "Tensor " << *this << " * matrix" << M << "of wrong size.\n";
    throw ii.str();
  }
  Tensor T;
  int i,j,k,l; 
  for (i=1; i<=3; i++)
    for (j=1; j<=3; j++)
      for (k=1; k<=3; k++) 
	for (l=1; l<=3; l++) 
	  T(i, j, k) = T(i, j, k) + display(l, j, k) * M.display(l, i);
  return T;
}

// T^kij = M_il T^klj
Tensor operator*(const DoubleMatrix &M, const Tensor &T) {
  if (M.displayRows() !=3 || M.displayCols() != 3) {
    ostringstream ii;
    ii << "DoubleMatrix " << M << " * Tensor" << T << "of wrong size\n";
    throw ii.str();
  }
  Tensor temp;
  int i,j,k,l; 
  for(k=1; k<=3; k++)
    for(i=1; i<=3; i++) 
      for(j=1; j<=3; j++)
	for(l=1; l<=3; l++)
	  temp(k, i, j) += M.display(i, l) * T.display(k, l, j);
  return temp;  
}

void chartData(std::vector<double>& results, bool cluster, int subjectBond, DoubleMatrix data, long* maturities, int date_tolerance, int lower_date, int upper_date, int length){
	DoubleMatrix lowerBonds;
	DoubleMatrix upperBonds;
	std::vector<double>timeSeries;
	upperAndLowerBonds(cluster, subjectBond, data, maturities, date_tolerance, lower_date, upper_date, upperBonds, lowerBonds);
	flyTimeSeries(results, lowerBonds, data.at(subjectBond), upperBonds, length);
}