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

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 StopRule::multiple (DoubleVector &vec1, DoubleVector &vec2, DoubleMatrix &proMat) {
	int row_, col_;
	proMat.resize(vec1.size());
	for (row_ = 0; row_ < vec1.size(); row_++)
		proMat[row_].resize(vec2.size());

	for (row_ = 0; row_ < vec1.size(); row_ ++)
		for (col_ = 0; col_ < vec2.size(); col_ ++)
			proMat[row_][col_] = vec1[row_] * 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 ];
}

void StopRule::readMat (char *fileName, DoubleMatrix &oriMat, int &size) {
	std::ifstream inFile_;
	inFile_.open(fileName);
	inFile_ >> size;
	oriMat.resize(size);
	for (int i = 0; i < size; i++) oriMat[i].resize(size);
	for (int row_ = 0; row_ < size; row_ ++)
		for (int col_ = 0; col_ < size; col_ ++)
			inFile_ >> oriMat[row_][col_];
	inFile_.close ();

}

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

}

long _stdcall getChartData(int bond_id, bool cluster, double* arr, unsigned char bonds[], int numbonds, long* maturities, int microWidth, int dateTolerance, int clusterLower, int clusterUpper, int history){
	SimpleRefData bbg(1);
	DoubleMatrix yieldData;
	std::vector <std::string> bond_vec;
	for (int b = 0; b < numbonds; b++){
		std::string prefix = "/isin/";
		std::string suffix;
		suffix.clear();
		for (int j = 0; j < ISIN_LENGTH; j++){
			char a = bonds[b + j*(numbonds + 1)];
			suffix += a;
		}
		bond_vec.push_back((prefix + suffix).c_str());
	}

	double** yieldArray;
	yieldArray = new double*[numbonds+1];
	for (int i = 0; i <= numbonds; i++){
		yieldArray[i] = new double[history];
	}

	int days_data = bbg.runHistData(yieldArray, bond_vec, history);

	yieldData.resize(numbonds);
	for (int i = 0; i <= numbonds; i++){
		yieldData[i].resize(days_data);
	}
	for (int i = 0; i <= numbonds; i++){
		for (int j = 0; j < days_data; j++){
			yieldData.at(i).at(j) = yieldArray[i][j];
		}
	}
	
	std::vector<double> dataLocal;
	if (cluster){
		chartData(dataLocal, cluster, bond_id, yieldData, maturities, 0, clusterLower, clusterUpper, days_data);
	}
	else{
		chartData(dataLocal, cluster, bond_id, yieldData, maturities, dateTolerance, microWidth, microWidth, days_data);
	}
	
	for (int i = 0; i < days_data; i++){
		arr[i] = dataLocal.at(i);
	}
	return days_data;
}

DoubleMatrix ParseMatrixFromText(const string& textMatrix)
{
	DoubleMatrix mat;

    //Parse the matrix
    vector<string> rows = splitString(textMatrix, "\n");
    for(int row = 0; row < rows.size(); row++)
    {
        vector<string> values = splitString(rows[row], " \t");
        for(int col = 0; col < values.size(); col++)
        {
        	if(!mat.size())
            {
                mat.resize(rows.size(), values.size());
            }

            mat(row, col) = toDouble(values[col]);
        }
    }
	return mat;
}

void StopRule::cmpLamdaMat (int k, DoubleMatrix &lamdaMat) {
	int i, j;
	lamdaMat.resize(k);
	for (i = 0; i < k; i ++)
		lamdaMat[i].resize(k);
	double muy_ = cmpMuy (k);
	for (i = 0; i < k; i ++)
		for (j = 0; j <= i; j ++) {
			/*
			lamdaMat[i][j] = (cmpGamma (2*muy_ + i + 1) * cmpGamma (muy_ + j + 1) ) /
			                 ( cmpGamma (muy_ + i + 1) * cmpGamma (j + 1) );*/

			// to fix divide by zero PROBLEM!
			lamdaMat[i][j] = cmpLnGamma (2*muy_ + i + 1) + cmpLnGamma (muy_ + j + 1) -
			                 cmpLnGamma (muy_ + i + 1) - cmpLnGamma (j + 1);

			//if (i == 98 && j == 97)
			 //     std::cout << i << "," << j << " -> " << lamdaMat[i][j] << endl;
			lamdaMat[i][j] = exp(lamdaMat[i][j]);
			lamdaMat[j][i] = lamdaMat[i][j];
		}
}

/*
 * b = [1]
 *     [1]
 */
void lambdaTest::setB(DoubleMatrix &b)
{
    b.resize(2,1);
    b.fill(rand()/(pow(2.0,32.0)) * 100.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);

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

/*Inputs: Destinaion Excel grid, array of ISINS, number of bonds, array of maturities, date tolerance for flies,  */
long _stdcall getFlyData(double* arr, unsigned char bonds[], int numbonds, long* maturities, int microWidth, int dateTolerance, int clusterLower, int clusterUpper, int history){
	SimpleRefData bbg(1);
	DoubleMatrix yieldData;
	DoubleMatrix micro;
	DoubleMatrix cluster;
	std::vector <std::string> bond_vec;
	std::wstring stemp2 = std::to_wstring(numbonds);
	LPCWSTR sw2 = stemp2.c_str();
	OutputDebugString(L"numbonds ");

	OutputDebugString(sw2);

	for (int b = 0; b <= numbonds; b++){
		std::string prefix = "/isin/";
		std::string suffix;
		suffix.clear();
		for (int j = 0; j < ISIN_LENGTH; j++){
			char a = bonds[b + j*(numbonds + 1)];
			suffix += a;
		}
		std::wstring stemp = std::wstring(suffix.begin(), suffix.end());
		LPCWSTR sw = stemp.c_str();
		OutputDebugString(sw);
		OutputDebugString(L"\n");

		bond_vec.push_back((prefix + suffix).c_str());
	}

	double** yieldArray;
	yieldArray = new double*[numbonds+1];
	for (int i = 0; i <= numbonds; i++){
		yieldArray[i] = new double[history];
	}

	int days_data = bbg.runHistData(yieldArray, bond_vec, history);
	std::wstring stemp = std::to_wstring(days_data);
	LPCWSTR sw = stemp.c_str();
	OutputDebugString(sw);

	yieldData.resize(numbonds+1);
	for (int i = 0; i <= numbonds; i++){
		yieldData[i].resize(days_data);
	}
	for (int i = 0; i <= numbonds; i++){
		for (int j = 0; j < days_data; j++){
			yieldData.at(i).at(j) = yieldArray[i][j];
		}
	}

	DoubleMatrix microMatrix;
	DoubleMatrix clusterMatrix;
	for (int i = 0; i < numbonds; i++){
		double microLocal[8];
		double clusterLocal[8];
		getFlyMetrics(microLocal, false, i, yieldData, maturities, dateTolerance, microWidth, microWidth, days_data);
		getFlyMetrics(clusterLocal, true, i, yieldData, maturities, 0, clusterLower, clusterUpper, days_data);
		double flyLocal[16];
		for (int j = 0; j < 8; j++){
			flyLocal[j] = microLocal[j];
			flyLocal[8 + j] = clusterLocal[j];
		}

		for (int k = 0; k < 16; k++){
			arr[i + k*(numbonds + 1)] = flyLocal[k];
		}
	}
	delete[] yieldArray;
	return 0;
}

void BackwardParabolicIBVP::gridMethod(DoubleMatrix &p, SweepMethodDirection direction) const
{
    typedef void (*t_algorithm)(const double*, const double*, const double*, const double*, double*, unsigned int);
    t_algorithm algorithm = &tomasAlgorithm;
    if (direction == ForwardSweep) algorithm = &tomasAlgorithmL2R;
    if (direction == BackwardSweep) algorithm = &tomasAlgorithmR2L;

    Dimension time = mtimeDimension;
    Dimension dim1 = mspaceDimension.at(0);

    double ht = time.step();
    unsigned int minM = time.min();
    unsigned int maxM = time.max();
    unsigned int M = maxM-minM;

    double hx = dim1.step();
    unsigned int minN = dim1.min();
    unsigned int maxN = dim1.max();
    unsigned int N = maxN-minN;

    double h = ht/(hx*hx);

    p.clear();
    p.resize(M+1, N+1);

    double *ka = (double*) malloc(sizeof(double)*(N-1));
    double *kb = (double*) malloc(sizeof(double)*(N-1));
    double *kc = (double*) malloc(sizeof(double)*(N-1));
    double *kd = (double*) malloc(sizeof(double)*(N-1));
    double *rx = (double*) malloc(sizeof(double)*(N-1));

    /* initial condition */
    SpaceNodePDE isn;
    for (unsigned int n=0; n<=N; n++)
    {
        isn.i = n+minN;
        isn.x = isn.i*hx;
        p[M][n] = initial(isn);
    }
    layerInfo(p, M);

    SpaceNodePDE lsn;
    lsn.i = minN;
    lsn.x = minN*hx;

    SpaceNodePDE rsn;
    rsn.i = maxN;
    rsn.x = maxN*hx;

    TimeNodePDE tn;
    for (unsigned int m=M-1; m!=UINT32_MAX; m--)
    {
        tn.i = m+minM;
        tn.t = tn.i*ht;

        for (unsigned int n=1; n<=N-1; n++)
        {
            isn.i = n+minN;
            isn.x = isn.i*hx;

            double alpha = -a(isn,tn)*h;
            double betta = 1.0 - 2.0*alpha;

            ka[n-1] = alpha;
            kb[n-1] = betta;
            kc[n-1] = alpha;
            kd[n-1] = p[m+1][n] - ht * f(isn, tn);
        }

        ka[0]   = 0.0;
        kc[N-2] = 0.0;

        /* border conditions */
        p[m][0] = boundary(lsn, tn);
        p[m][N] = boundary(rsn, tn);

        kd[0]   += a(lsn,tn) * h * p[m][0];
        kd[N-2] += a(rsn,tn) * h * p[m][N];

        (*algorithm)(ka, kb, kc, kd, rx, N-1);

        for (unsigned int n=1; n<=N-1; n++) p[m][n] = rx[n-1];

        layerInfo(p, m);
    }

    free(ka);
    free(kb);
    free(kc);
    free(kd);
    free(rx);
}

void ppedOpt( SpdModel&             model,
              Optimizer&            xOptInfo,
              const DoubleMatrix&   dvecXLow,
              const DoubleMatrix&   dvecXUp,
              const DoubleMatrix&   dvecXIn,
              DoubleMatrix*         pdvecXOut,
              const DoubleMatrix&   dvecXStep,
              Optimizer&            alpOptInfo,
              const DoubleMatrix&   dvecAlpLow,
              const DoubleMatrix&   dvecAlpUp,
              const DoubleMatrix&   dvecAlpIn,
              DoubleMatrix*         pdvecAlpOut,
              const DoubleMatrix&   dvecAlpStep,
              double*               pdPhiTildeOut,
              DoubleMatrix*         pdrowPhiTilde_xOut,
              DoubleMatrix*         pdmatPhiTilde_x_xOut )
{
  //------------------------------------------------------------
  // Preliminaries.
  //------------------------------------------------------------

  using namespace std;

  // If no evaluation is requested, return immediately.
  if ( !pdvecXOut && 
       !pdPhiTildeOut &&
       !pdrowPhiTilde_xOut &&
       !pdmatPhiTilde_x_xOut )
  {
    return;
  }

  // Get the number of design and fixed population parameters.
  int nX   = dvecXIn.nr();
  int nAlp = dvecAlpIn.nr();


  //------------------------------------------------------------
  // Validate the inputs (Debug mode).
  //------------------------------------------------------------

  // Check the length of the design parameter vector.
  assert( nX == model.getNDesPar() );


  //------------------------------------------------------------
  // Validate the inputs (All modes).
  //------------------------------------------------------------

  // [Revisit - Warm Start Disabled - Mitch]
  // In order to simplify the testing of this function,
  // warm starts have been disabled for now.
  if ( xOptInfo.getIsWarmStart() || alpOptInfo.getIsWarmStart() )
  {
    throw SpkException(
      SpkError::SPK_USER_INPUT_ERR,
      "One of the input Optimizer objects requested a warm start, which have been disabled for now.",
      __LINE__,
      __FILE__ );
  }


  //------------------------------------------------------------
  // Prepare the values required to calculate the output values.
  //------------------------------------------------------------

  // If this value is required to calculate the output values, 
  // initialize it so that it wil be calculated.  Otherwise, 
  // set its pointer to zero so that it won't be calculated.
  DoubleMatrix dvecXOutTemp;
  DoubleMatrix* pdvecXOutTemp = &dvecXOutTemp;
  if ( pdvecXOut )
  {
    dvecXOutTemp.resize( nX, 1 );
  }
  else
  {
    pdvecXOutTemp = 0;
  }

  // If this value is required to calculate the output values, 
  // initialize it so that it wil be calculated.  Otherwise, 
  // set its pointer to zero so that it won't be calculated.
  DoubleMatrix dvecAlpOutTemp;
  DoubleMatrix* pdvecAlpOutTemp = &dvecAlpOutTemp;
  if ( pdvecAlpOut )
  {
    dvecAlpOutTemp.resize( nAlp, 1 );
  }
  else
  {
    pdvecAlpOutTemp = 0;
  }

  // If this value is required to calculate the output values, 
  // initialize it so that it wil be calculated.  Otherwise, 
  // set its pointer to zero so that it won't be calculated.
  double dLTildeOut;
  double* pdLTildeOut = &dLTildeOut;
  if ( pdPhiTildeOut || pdrowPhiTilde_xOut || pdmatPhiTilde_x_xOut )
//.........这里部分代码省略.........

/*
 * y = [1]
 *     [1]
 */
void lambdaTest::setY(DoubleMatrix &y)
{
    y.resize(2,1);
    y.fill(rand()/(pow(2.0,32.0)) * 100.0 );
}

/*
 * alpha = [1]
 *         [1]
 */
void lambdaTest::setAlp(DoubleMatrix &alp)
{
    alp.resize(2,1);
    alp.fill(rand()/(pow(2.0,32.0)) * 100.0);
}