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

// parallel version with arbitrary dimensions
int gridSearchParallel(DoublePyArrayVector const &controlGridArray, MaximizerCallParams &params, double &rMaxVal, DoubleVector &rArgMaxArray) {  
  int nGrids = controlGridArray.size();
  int i;
  IntVector lenArray(nGrids);
  
  for (i=0; i<nGrids; i++) {
    // check that all grids are 1-dimensional
    assert(controlGridArray[i].ndim() == 1);  
	lenArray[i] = controlGridArray[i].dims()[0];
  }

  unsigned int totalGridSize = 1;
  double totalGridSize2 = 1.0;
  for (i=0; i<nGrids; i++) {
	totalGridSize *= lenArray[i];
	totalGridSize2 *= double(lenArray[i]);
  }
  // check that total grid space isn't too large to be indexed
  assert(totalGridSize2 < double(UINT_MAX));
  
  MaxIndexFnObj2 fnObj(controlGridArray, params);  
  parallel_reduce( blocked_range<size_t>(0, totalGridSize), fnObj);  

  rMaxVal = fnObj.m_value_of_max;
  rArgMaxArray.resize(fnObj.m_argmax.size());
  copy(fnObj.m_argmax.begin(), fnObj.m_argmax.end(), rArgMaxArray.begin());
  return 1;
}

int RateGammaInvar::computePatternRates(DoubleVector &pattern_rates, IntVector &pattern_cat) {
	//cout << "Computing Gamma site rates by empirical Bayes..." << endl;

	phylo_tree->computePatternLhCat(WSL_RATECAT);

	int npattern = phylo_tree->aln->getNPattern();
	pattern_rates.resize(npattern);
	pattern_cat.resize(npattern);

	double *lh_cat = phylo_tree->_pattern_lh_cat;
	for (int i = 0; i < npattern; i++) {
		double sum_rate = 0.0, sum_lh = phylo_tree->ptn_invar[i];
		int best = 0;
        double best_lh = phylo_tree->ptn_invar[i];
		for (int c = 0; c < ncategory; c++) {
			sum_rate += rates[c] * lh_cat[c];
			sum_lh += lh_cat[c];
			if (lh_cat[c] > best_lh || (lh_cat[c] == best_lh && random_double()<0.5)) { // break tie at random
                best = c+1;
                best_lh = lh_cat[c];
            }
		}
		pattern_rates[i] = sum_rate / sum_lh;
		pattern_cat[i] = best;
		lh_cat += ncategory;
	}
    return ncategory+1;
}

int RateGamma::computePatternRates(DoubleVector &pattern_rates, IntVector &pattern_cat) {
	//cout << "Computing Gamma site rates by empirical Bayes..." << endl;

	phylo_tree->computePatternLhCat(WSL_RATECAT);

	int npattern = phylo_tree->aln->getNPattern();
	pattern_rates.resize(npattern);
	pattern_cat.resize(npattern);

	double *lh_cat = phylo_tree->_pattern_lh_cat;
	for (int i = 0; i < npattern; i++) {
		double sum_rate = 0.0, sum_lh = 0.0;
		int best = 0;
		for (int c = 0; c < ncategory; c++) {
			sum_rate += rates[c] * lh_cat[c];
			sum_lh += lh_cat[c];
			if (lh_cat[c] > lh_cat[best] || (lh_cat[c] == lh_cat[best] && random_double()<0.5))  // break tie at random
                best = c;
		}
		pattern_rates[i] = sum_rate / sum_lh;
		pattern_cat[i] = best;
		lh_cat += ncategory;
	}
    return ncategory;

//	pattern_rates.clear();
//	pattern_rates.insert(pattern_rates.begin(), ptn_rates, ptn_rates + npattern);
//	pattern_cat.resize(npattern, 0);
//	for (int i = 0; i < npattern; i++)
//		for (int j = 1; j < ncategory; j++)
//			if (fabs(rates[j] - ptn_rates[i]) < fabs(rates[pattern_cat[i]] - ptn_rates[i]))
//				pattern_cat[i] = j;
//	delete [] ptn_rates;
}

void RRSchedule::init2DEmpty(DoubleVector &_invect, int row, int col)
{
    _invect.clear();
    _invect.resize(row);
    for (int i = 0; i < row; i++)
        allocate1D(_invect[i], col);
}

void RateMeyerHaeseler::getRates(DoubleVector &rates) {
	rates.clear();
	if (empty()) {
		rates.resize(phylo_tree->aln->size(), 1.0);
	} else {
		rates.insert(rates.begin(), begin(), end());
	} 
}

void StopRule::multiple (DoubleVector &vec1, DoubleMatrix &mat2, DoubleVector &proVec) {
	int row_, col_;
	proVec.resize(mat2[0].size());
	for (col_ = 0; col_ < mat2[0].size(); col_ ++) {
		proVec[col_] = 0.0;
		for (row_ = 0; row_ < mat2.size(); row_ ++)
			proVec[col_] += vec1[row_] * mat2[row_][col_];
	}
}

void HeatControl2Delta::initialize()
{
    U.clear();
    U.resize(N2+1, N1+1);

    DoubleVector x = O;

#ifdef POWER_OPTIMIZE
    x.resize( 2*L + (M+1)*L );
#else
    x.resize( 2*L );
#endif

#ifdef POWER_OPTIMIZE
    for (unsigned int k=0; k<=M; k++)
    {
        x[2*L + 0*(M+1) + k] = g1(k*ht);
        x[2*L + 1*(M+1) + k] = g2(k*ht);
        x[2*L + 2*(M+1) + k] = g3(k*ht);
    }
#endif

    px = &x;
    IParabolicEquation2D::calculateMVD(U, h1, h2, ht, N1, N2, M, a1, a2);

    puts("+------------------------------------------------------------------------------------------------------------------------------------------------------------------+");
    IPrinter::printMatrix(U, 10, 10);
    printf("eo: [%12.8f, %12.8f] [%12.8f, %12.8f] [%12.8f, %12.8f]\n", x[0], x[1], x[2], x[3], x[4], x[5]);
    puts("+------------------------------------------------------------------------------------------------------------------------------------------------------------------+");

    write("optimal.txt", U);

    FILE* f = fopen("optimal1.txt", "w");
    IPrinter::printMatrix(U, N1, N2, NULL, f);
    fclose(f);
}

DoubleVector& get_cleaved_masses ( const string& protein, const IntVector& cleavageIndex )
{
	static DoubleVector cleavedMassArray ( 36000 );
	StringSizeType numAA = protein.length ();
	if ( numAA > cleavedMassArray.size () ) cleavedMassArray.resize ( numAA );

	for ( StringSizeType i = 0, j = 0 ; i < numAA ; ) {
		double mass = 0.0;
		StringSizeType index = cleavageIndex [j];
		while ( i <= index ) {
			mass += amino_acid_wt [protein[i++]];
		}
		cleavedMassArray [j++] = mass;
	}
	return cleavedMassArray;
}

void Example1::runge_kutta(OrdDifEquation *ode, double x0, double y0, unsigned int N, DoubleVector &y, double h)
{
    double k1 = 0.0;
    double k2 = 0.0;
    double k3 = 0.0;
    double k4 = 0.0;

    y.clear();
    y.resize(N+1,0.0);

    y[0] = y0;
    for (unsigned int i=0; i<N; i++)
    {
        double _x = x0 + i*h;
        double _y = y[i];

        k1 = ode->fx(_x, _y);
        k2 = ode->fx(_x+h/2.0, _y+(h/2.0)*k1);
        k3 = ode->fx(_x+h/2.0, _y+(h/2.0)*k2);
        k4 = ode->fx(_x+h, _y+h*k3);
        y[i+1] = y[i] + (h/6.0) * (k1 + 2.0*k2 + 2.0*k3 + k4);;
    }
}

void IFunctional::toVector(const Parameter &prm, DoubleVector &pv) const
{
    unsigned int Lc = prm.Lc;
    unsigned int Lo = prm.Lo;

    //unsigned int length = 0;
    //if (optimizeK) length += Lc*Lo;
    //if (optimizeZ) length += Lc*Lo;
    //if (optimizeC) length += 2*Lc;
    //if (optimizeO) length += 2*Lo;

    pv.clear();
    pv.resize(2*Lc*Lo+2*Lo+2*Lc);

    //if (optimizeK)
    {
        for (unsigned int i=0; i<Lc; i++)
        {
            for (unsigned int j=0; j<Lo; j++)
            {
                pv[i*Lo + j] = prm.k[i][j];
            }
        }
    }

    //if (optimizeZ)
    {
        for (unsigned int i=0; i<Lc; i++)
        {
            for (unsigned int j=0; j<Lo; j++)
            {
                pv[i*Lo + j + Lc*Lo] = prm.z[i][j];
            }
        }
    }

    //if (optimizeO)
    {
        for (unsigned int j=0; j<Lo; j++)
        {
            pv[2*j + 0 + 2*Lc*Lo] = prm.xi[j].x;
            pv[2*j + 1 + 2*Lc*Lo] = prm.xi[j].y;
        }
    }

    //if (optimizeC)
    {
        for (unsigned int i=0; i<Lc; i++)
        {
            pv[2*i + 0 + 2*Lo + 2*Lc*Lo] = prm.eta[i].x;
            pv[2*i + 1 + 2*Lo + 2*Lc*Lo] = prm.eta[i].y;
        }
    }

    //    for (unsigned int i=0; i<Lc; i++)
    //    {
    //        x[2*i + 0 + 2*Lc*Lo] = prm.eta[i].x;
    //        x[2*i + 1 + 2*Lc*Lo] = prm.eta[i].y;

    //        for (unsigned int j=0; j<Lo; j++)
    //        {
    //            x[i*Lo + j] = prm.k[i][j];
    //            x[i*Lo + j + Lc*Lo] = prm.z[i][j];

    //            x[2*j + 0 + 2*Lc + 2*Lc*Lo] = prm.xi[j].x;
    //            x[2*j + 1 + 2*Lc + 2*Lc*Lo] = prm.xi[j].y;
    //        }
    //    }
}

void IFunctional::gradient(const DoubleVector &pv, DoubleVector &g) const
{
    IFunctional *ifunc = const_cast<IFunctional*>(this);
    unsigned int L = mTimeDimension.sizeN();
    double ht = mTimeDimension.step();

    ifunc->fromVector(pv, ifunc->mParameter);
    forward->setParameter(mParameter);
    backward->setParameter(mParameter);

    DoubleMatrix u;
    DoubleMatrix p;

    vector<ExtendedSpaceNode2D> u_info;
    forward->calculateMVD(u, u_info, true);

    backward->u = &u;
    backward->U = &ifunc->U;
    backward->info = &u_info;
    vector<ExtendedSpaceNode2D> p_info;
    backward->calculateMVD(p, p_info, true);

    g.clear();
    g.resize(pv.length(), 0.0);
    unsigned int gi = 0;

    // k
    if (optimizeK)
    {
        for (unsigned int i=0; i<mParameter.Lc; i++)
        {
            ExtendedSpaceNode2D &pi = p_info[i];

            for (unsigned int j=0; j<mParameter.Lo; j++)
            {
                ExtendedSpaceNode2D &uj = u_info[j];

                double grad_Kij = 0.0;

                grad_Kij += 0.5 * (pi.value(0)+2.0*r*gpi(i,0,u_info)*sgn(g0i(i,0,u_info))) * (uj.value(0) - mParameter.z.at(i,j));
                for (unsigned int m=1; m<=L-1; m++)
                {
                    grad_Kij += (pi.value(m)+2.0*r*gpi(i,m,u_info)*sgn(g0i(i,m,u_info))) * (uj.value(m) - mParameter.z.at(i,j));
                }
                grad_Kij += 0.5 * (pi.value(L)+2.0*r*gpi(i,L,u_info)*sgn(g0i(i,L,u_info))) * (uj.value(L) - mParameter.z.at(i,j));
                grad_Kij *= -ht;

                g[gi++] = grad_Kij + 2.0*regEpsilon*(mParameter.k.at(i,j) - mParameter0.k.at(i,j));
            }
        }
    }
    else
    {
        for (unsigned int i=0; i<mParameter.Lc; i++)
        {
            for (unsigned int j=0; j<mParameter.Lo; j++)
            {
                g[gi++] = 0.0;
            }
        }
    }

    // z
    if (optimizeZ)
    {
        for (unsigned int i=0; i<mParameter.Lc; i++)
        {
            ExtendedSpaceNode2D &pi = p_info[i];

            for (unsigned int j=0; j<mParameter.Lo; j++)
            {
                double grad_Zij = 0.0;

                grad_Zij += 0.5 * (pi.value(0)+2.0*r*gpi(i,0,u_info)*sgn(g0i(i,0,u_info))) * mParameter.k.at(i,j);
                for (unsigned int m=1; m<=L-1; m++)
                {
                    grad_Zij += (pi.value(m)+2.0*r*gpi(i,m,u_info)*sgn(g0i(i,m,u_info)))  * mParameter.k.at(i,j);
                }
                grad_Zij += 0.5 * (pi.value(L)+2.0*r*gpi(i,L,u_info)*sgn(g0i(i,L,u_info))) * mParameter.k.at(i,j);
                grad_Zij *= ht;

                g[gi++] = grad_Zij + 2.0*regEpsilon*(mParameter.z[i][j] - mParameter0.z[i][j]);
            }
        }
    }
    else
    {
        for (unsigned int i=0; i<mParameter.Lc; i++)
        {
            for (unsigned int j=0; j<mParameter.Lo; j++)
            {
                g[gi++] = 0.0;
            }
        }
    }

    // xi
    if (optimizeO)
    {
        for (unsigned int j=0; j<mParameter.Lo; j++)
//.........这里部分代码省略.........

void CplexSolver::rc(DoubleVector & result)const{
	result.resize(ncols());
	CPXgetdj(_env, _prob, result.data(), 0, ncols() - 1);
}

void CplexSolver::solution(int i, DoubleVector & x) const {
	x.resize(ncols());
	CPXgetsolnpoolx(_env, _prob, (int) i, x.data(), 0, (int) (x.size() - 1));
}

// single-threaded grid search with an arbitrary number of dimensions
int gridSearch(DoublePyArrayVector const &controlGridArray, MaximizerCallParams &params, double &rMaxVal, DoubleVector &rArgMaxArray) {  
  int nGrids = controlGridArray.size();
  int i;
  IntVector lenArray(nGrids), strideArray(nGrids), dataIndexArray(nGrids);
  
  for (i=0; i<nGrids; i++) {
    // check that all grids are 1-dimensional
    assert(controlGridArray[i].ndim() == 1);  
	lenArray[i] = controlGridArray[i].dims()[0];
	strideArray[i] = controlGridArray[i].strides()[0];
	dataIndexArray[i] = 0;
  }

  DoubleVector argArray(nGrids);  
  int nIter = 0;
  int nMaxMultiplicity = 0;  
  bool bDone = false;
  
  // check for zero size, if one grid has zero size then there's nothing to do
  for (i=0; i<nGrids; i++) {
    if (lenArray[i] == 0) {
	  bDone = true;
	  break;
	}
  }
  double max = -DBL_MAX;
  while (!bDone) {
    double result;
	// get double array args from char* data
	for (i=0; i<nGrids; i++) {
	  char *pData = (char*) (controlGridArray[i].array().data());
	  char *pArg = pData + (strideArray[i] * dataIndexArray[i]);
	  argArray[i] = *(double*) pArg;
	}
	result = params.objectiveFunction(argArray);
	if (nIter == 0) {			// first iteration
      max = result;
	  rArgMaxArray.resize(argArray.size());
	  copy(argArray.begin(), argArray.end(), rArgMaxArray.begin());		
	  nMaxMultiplicity = 1;
    } else {
      if (result > max) {
        max = result;
		rArgMaxArray.resize(argArray.size());
		copy(argArray.begin(), argArray.end(), rArgMaxArray.begin());		
	    nMaxMultiplicity = 1;
	  } else if (result == max) {
	    nMaxMultiplicity++;
	  }
	}
	nIter++;

    // increment indices for next iteration
	bDone = true;
	for (i=nGrids-1; i>=0; i--) {
	  // increment index i
	  dataIndexArray[i] += 1;
	  if (dataIndexArray[i] == lenArray[i]) {
	    // cycle this index
		dataIndexArray[i] = 0;
	  } else {
	    bDone = false;
		break;
	  }
	  // if we make it here, all indices have cycled, therefore we are done
	}
	
  }
  rMaxVal = max;
  return nMaxMultiplicity;
}

//////////////////////////////////////////////////////////////////////
// CInterpolationEngine Interface
//////////////////////////////////////////////////////////////////////
void CInterpolationEngine::RecalculateVector() throw(CInterpolationException)
{
	// Sort point vector
	std::sort( m_pts.begin(), m_pts.end() );

	// Spline interpolation is used for more than 2 points
	long i, size = m_pts.size();
	if( size < MIN_SPLINE_POINTS )
		return;

	// Calculate points offsets
	m_h.resize( size );
	for( i = 0; i < size - 1; i++ )
		m_h[i] = m_pts[i+1].x - m_pts[i].x;

	// Fill auxiliary arrays
	DoubleVector lam, mu, d;
	lam.resize( size );
	mu.resize( size );
	d.resize( size );
	for( i = 1; i < size - 1; i++ )
	{
		lam[i] = m_h[i]/( m_h[i-1] + m_h[i] );
		mu[i]  = 1 - lam[i];
		d[i]   = 2 * (3 + m_dr) * (lam[i] * (m_pts[i].f - m_pts[i-1].f)/m_h[i-1] +
			     mu[i] * (m_pts[i+1].f - m_pts[i].f)/m_h[i]);
	}
	d[0] = (3 + m_dr) * (m_pts[1].f - m_pts[0].f)/m_h[0] -
		   m_h[0]/m_dR * (m_pts[0].f/m_h[1]/(m_h[1]+m_h[0]) - m_pts[1].f/m_h[1]/m_h[0] + m_pts[2].f/m_h[0]/(m_h[1] + m_h[0]));
	d[size-1] = (3 + m_dr) * (m_pts[size-1].f - m_pts[size-2].f)/m_h[size-2] + 
				m_h[size-2]/m_dR * (m_pts[size-1].f/m_h[size-3]/(m_h[size-2] + m_h[size-3]) - m_pts[size-2].f/m_h[size-2]/m_h[size-3] +
				m_pts[size-3].f/m_h[size-2]/(m_h[size-2] + m_h[size-3]));

	// Forward
	DoubleVector P, Q;
	P.resize( size );
	Q.resize( size );
	P[0] = -1/( 2 + m_dr );
	Q[0] = d[0]/( 2 + m_dr );
	for( i = 1; i < size - 1; i++ )
	{
		P[i] = -mu[i]/(lam[i] * P[i-1] + 2 + m_dr);
		Q[i] = P[i] * (2 * lam[i] * Q[i-1] - d[i]);
	}

	// Backward
	DoubleVector m;
	m.resize( size );
	m[size-1] = (d[size-1] - Q[size-2])/(2 + m_dr + P[size-2]);
	for( i = size - 2; i >= 0; i-- )
		m[i] = P[i] * m[i+1] + Q[i];

	// Factor arrays
	m_A.resize( size - 1 );
	m_B.resize( size - 1 );
	m_C.resize( size - 1 );
	m_D.resize( size - 1 );
	for( i = 0; i < size - 1; i++ )
	{
		m_C[i] = (-(3 + m_dr) * (m_pts[i+1].f - m_pts[i].f) + m_h[i] * m[i]   + (2 + m_dr) * m_h[i] * m[i+1])/(3 + 4 * m_dr + m_dr * m_dr);
		m_D[i] = ( (3 + m_dr) * (m_pts[i+1].f - m_pts[i].f) - m_h[i] * m[i+1] - (2 + m_dr) * m_h[i] * m[i])/(3 + 4 * m_dr + m_dr * m_dr);
		m_A[i] = m_pts[i+1].f - m_C[i];
		m_B[i] = m_pts[i].f - m_D[i];
	}
}