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

// do partitions on each variable selected in the redundent group
void RddtClust::doPartitions(){
    DoubleVector dimData;

    this->m_partitionPerVariable.clear();
    //m_OkcData->GetDimensionData(&dimData,);
    for(unsigned i = 0; i<this->m_nRddtDims.size(); i++){

        dimData.clear();
        this->getDimData(dimData,m_nRddtDims[i]);
        //this->m_rddtOp->getOkcData()->
                //GetDimensionData(dimData,m_nRddtDims[i]);

        DataPartitions* partitions = new DataPartitions();
        partitions->setRddtClust(this);
        partitions->setVarIdx(m_nRddtDims[i]);
        partitions->cutData(dimData);
        m_partitionPerVariable.push_back(partitions);
    }
    for(unsigned i = 0; i<this->m_rddtDims.size(); i++){
        dimData.clear();
        //this->m_rddtOp->getOkcData()->
                //GetDimensionData(dimData,m_rddtDims[i]);
        //dimData.clear();
        this->getDimData(dimData,m_rddtDims[i]);
        DataPartitions* partitions = new DataPartitions();
        partitions->setRddtClust(this);
        partitions->setVarIdx(m_rddtDims[i]);
        partitions->cutData(dimData);
        m_partitionPerVariable.push_back(partitions);
    }
}

TEST(String,NumFractionalDigits) {
  DoubleVector values;
  values.push_back(128196.198);
  values.push_back(19671.281);
  values.push_back(218528.28);
  values.push_back(192.186);

  std::pair<unsigned,unsigned> result = numFractionalDigits(values,3u);
  EXPECT_EQ(0u,result.first); EXPECT_EQ(0u,result.second);

  result = numFractionalDigits(values,4u);
  EXPECT_EQ(0u,result.first); EXPECT_EQ(1u,result.second);

  result = numFractionalDigits(values,7u);
  EXPECT_EQ(1u,result.first); EXPECT_EQ(4u,result.second);

  values.clear();
  values.push_back(0.189678);
  values.push_back(0.001869168);
  values.push_back(0.7198);

  result = numFractionalDigits(values,2u);
  EXPECT_EQ(2u,result.first); EXPECT_EQ(4u,result.second);

  result = numFractionalDigits(values,8u);
  EXPECT_EQ(8u,result.first); EXPECT_EQ(10u,result.second);

  values.clear();
  values.push_back(0.07592);
  values.push_back(198.82);
  values.push_back(210.28);
  values.push_back(0.628);

  result = numFractionalDigits(values,2u);
  EXPECT_EQ(0u,result.first); EXPECT_EQ(3u,result.second);

  result = numFractionalDigits(values,3u);
  EXPECT_EQ(0u,result.first); EXPECT_EQ(4u,result.second);

  result = numFractionalDigits(values,5u);
  EXPECT_EQ(2u,result.first); EXPECT_EQ(6u,result.second);

  values.clear();
  values.push_back(0.0);

  result = numFractionalDigits(values,2u);
  EXPECT_EQ(1u,result.first); EXPECT_EQ(1u,result.second);

  result = numFractionalDigits(values,5u);
  EXPECT_EQ(4u,result.first); EXPECT_EQ(4u,result.second);
}

int main(int argc, char **argv)
{
  std::string cornerString="   0.1     0.1      0.1   20 20 10";
  
  typedef boost::tokenizer<boost::char_separator<char> >            Tokenizer;
  typedef std::vector<double>                                       DoubleVector;
  typedef boost::tokenizer<boost::char_separator<char> >::iterator  TokIterator;

  boost::char_separator<char> sep(", ");
  Tokenizer tok(cornerString, sep);
  DoubleVector corners;
  corners.clear();

  for (TokIterator it=tok.begin(); it!=tok.end(); ++it ) {
    double q;
    printf("");
    sscanf(it->c_str(), "%lf", &q);
    //DPrintf(DebugReaction," %lf  ['%s'];  ", q, it->c_str());
    corners.push_back( q );
  }
  printf("the corners are:\n   ");
  printf("%e", corners[0]);
  for(int i=1; i<corners.size(); ++i ) {
    printf(", %e", corners[i]);
  }
  printf("\n");
}

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 Fresnel::spiral (DoubleVector &spiralX, DoubleVector &spiralY) const
{
    spiralX.clear ();
    spiralY.clear ();
    spiralX.push_back (0.0);
    spiralY.push_back (0.0);

    Complex sp;

    int    fn     = currentFresnelNumber;
    double R      = holeRadius;
    double b      = observerDistance;
    double b2     = b * b;
    double innerR = 0.0;
    double outerR = 0.0;
    double innerD = b;
    double outerD = 0.0;
    double dd     = 0.0;

    for (int n = 0; n < fn + 1; ++n) {
        outerR = zoneOuterRadius (n);
        outerD = sqrt (outerR*outerR + b2);
        dd     = (outerD - innerD) / accuracySpiral;
        bool stop = false;

        for (int i = 0; i < accuracySpiral && !stop; ++i) {
            innerR = sqrt (innerD*innerD - b2);

            innerD += dd;
            outerR  = sqrt (innerD*innerD - b2);

            if (outerR >= R) {
                outerR = R;
                stop = true;
            }

            sp += amplitude (innerR, outerR);
            spiralX.push_back (sp.im);
            spiralY.push_back (sp.re);
        }
    }
}

void RddtClust::getDimData(DoubleVector& dimData, int idx){
    using namespace Rcpp;
    XmdvToolMainWnd* m_mwd = this->m_rddtOp->m_pm->getMainWnd();
    RInside r = m_mwd->getRinstance();
    //qDebug()<<"idx qt "<<idx;
    r["currentDim"] = idx;

    std::string cmd = "currentDim <- as.integer(currentDim)+1;\n"
            //"print(currentDim);\n"
            "dimData <- M[,currentDim];\n"
            "dimData";
    SEXP evals = r.parseEval(cmd);
    NumericVector rDimData(evals);

    dimData.clear();
    for(int i = 0; i< rDimData.length(); i++){
        dimData.push_back(rDimData[i]);
    }
}

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 NucleusGridFunction::
readParameterFile( ParameterReader &params )
{
  //read:
  // nucleus type=none, file, box, sphere
  // nucleus file= <file.cwn>
  // nucleus corners= x0 y0 z0 x1 y1 z1
  // nucleus center= x0 y0 z0
  // nucleus radius= r0
  // nucleus boundary thickness= t0

  const std::string keyBoundaryThickness="nucleus boundary thickness";
  double thickness;
  params.get( "nucleus boundary thickness", thickness, -1. );
  if ( thickness>=0. ) nucleusBoundaryThickness=thickness;

  typedef boost::tokenizer<>::iterator TokIterator;
  typedef std::vector<double>          DoubleVector;

  bool noNucleus=false;
  std::string nucleusType;
  params.get( "nucleus type", nucleusType, "");
  if ( (nucleusType == "box" )   || (nucleusType == "cubic")) {
    Nucleus::NucleusShape nucleusShape=Nucleus::BoxNucleus;
    std::string doubleList="";
    params.get("nucleus corners",  doubleList ); 
    DPrintf(DebugSolver,"..NucleusGridFunction: nucleus type=box.\n");
    DPrintf(DebugSolver,"....  nucleus corners = '%s'\n", doubleList.c_str());

    boost::tokenizer<> tok(doubleList);
    DoubleVector       corners;

    DPrintf(DebugSolver,"    corners are = ");
    corners.clear();
    for (TokIterator it=tok.begin(); it!=tok.end(); ++it) {
      double q;
      sscanf(it->c_str(), "%lf", &q);
      DPrintf(DebugSolver," %lf ", q);
      corners.push_back( q );
    }
    DPrintf(DebugSolver,"\n");
    //..create the nucleic info
    if( corners.size() >5 ) { 
      Nucleus thisNuc;
      thisNuc.setID(1);
      thisNuc.setBoundaryThickness( nucleusBoundaryThickness );
      thisNuc.setCorners(corners[0], corners[1], corners[2],
			 corners[3], corners[4], corners[5]);
      thisNuc.setShape( nucleusShape );

      nucleus.push_back( thisNuc );
    }
  }
  else if ( (nucleusType == "sphere") || ( nucleusType == "spherical")) {
    Nucleus::NucleusShape nucleusShape=Nucleus::SphericalNucleus;
    double radius;
    params.get("nucleus radius",radius, 20.);
    std::string doubleList="";
    params.get("nucleus center",  doubleList ); 
    DPrintf(DebugSolver,"..NucleusGridFunction: nucleus type=sphere not supported yet:\n");
    DPrintf(DebugSolver,"..  nucleus center = '%s'\n", doubleList.c_str());

    boost::tokenizer<> tok(doubleList);
    DoubleVector       center;

    DPrintf(DebugSolver,"    center is = ");
    center.clear();
    for (TokIterator it=tok.begin(); it!=tok.end(); ++it) {
      double q;
      sscanf(it->c_str(), "%lf", &q);
      DPrintf(DebugSolver," %lf ", q);
      center.push_back( q );
    }
    DPrintf(DebugSolver,"\n");
    //..create the nucleic info
    double x0=0., y0=0., z0=0.;
    if( center.size() >1) {
      x0=center[0];
      y0=center[1];
    }
    if( center.size()>2) {
      z0=center[2];
    }

    Nucleus thisNuc;
    thisNuc.setID(1);
    
    thisNuc.setBoundaryThickness( nucleusBoundaryThickness );
    thisNuc.setCenter( x0, y0, z0);
    thisNuc.setRadius( radius);
    thisNuc.setShape( nucleusShape );
    nucleus.push_back( thisNuc );

  }
  else if ( nucleusType == "file" ) {
    std::string nucleusFileName="";
    params.get("nucleus file",nucleusFileName);
    DPrintf(DebugSolver,"nuclei read in from file='%s'\n", nucleusFileName.c_str());
    readCellNucleusFile( nucleusFileName );
  }
//.........这里部分代码省略.........

bool
ExemplarClient::getRandomNumbers(const int      howMany,
                                 DoubleVector & result)
{
    ODL_OBJENTER(); //####
    ODL_I1("howMany = ", howMany); //####
    ODL_P1("result = ", &result); //####
    bool okSoFar = false;

    try
    {
        if (0 < howMany)
        {
            yarp::os::Bottle parameters;
            ServiceResponse  response;

            parameters.addInt(howMany);
            reconnectIfDisconnected();
            if (send(MpM_SIMPLE_REQUEST_, parameters, response))
            {
                if (howMany == response.count())
                {
                    result.clear();
                    okSoFar = true;
                    for (int ii = 0; ii < howMany; ++ii)
                    {
                        yarp::os::Value retrieved(response.element(ii));

                        if (retrieved.isDouble())
                        {
                            result.push_back(retrieved.asDouble());
                        }
                        else
                        {
                            ODL_LOG("! (retrieved.isDouble())"); //####
                            okSoFar = false;
                            break;
                        }

                    }
                }
                else
                {
                    ODL_LOG("! (howMany == response.count())"); //####
                    ODL_S1s("response = ", response.asString()); //####
                }
            }
            else
            {
                ODL_LOG("! (send(MpM_SIMPLE_REQUEST_, parameters, response))"); //####
            }
        }
        else
        {
            ODL_LOG("! (0 < howMany)"); //####
        }
    }
    catch (...)
    {
        ODL_LOG("Exception caught"); //####
        throw;
    }
    ODL_OBJEXIT_B(okSoFar); //####
    return okSoFar;
} // ExemplarClient::getRandomNumbers

void BackwardParabolicIBVP::gridMethod(DoubleVector &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(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[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[n] - ht * f(isn, tn);
        }

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

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

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

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

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

        layerInfo(p, m);
    }

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

void TimeWarp::calculatePartJointSimilarityMatrix(DoubleVector* firstEnergyVector, DoubleVector* secondEnergyVector, DoubleMatrix* chromaSimMatrix, DoubleMatrix* simMatrix, const int& startX, const int& startY, int endX, int endY){
//	printf("PART SIM CALC Calculate similarity : pointers : startX %i startY %i, endX %i endY %i  ", startX, startY, endX, endY);
//	printf("PART SIM ENERGY VEC LENGTHS %i and %i\n", (int) (*firstEnergyVector).size(), (int) (*secondEnergyVector).size());
	
	conversionFactor = (int) round((*firstEnergyVector).size() / (*chromaSimMatrix).size() );
//	printf("conversion fac %i", conversionFactor);
	
	simMatrix->clear();
	
	double energyProportion = 0.2;
	double chromaProportion = 1 - energyProportion;
	double distance, firstSum, secondSum;
	
	endX = min (endX, (int)(*firstEnergyVector).size()-1);//in case out of size
	endY = min( endY+1, (int)(*secondEnergyVector).size());
//	endY = min( endY, (int)(*secondEnergyVector).size()-1);//added to switch 

	int lastChromaYvalue = 0;
	int chromaIndexY = 0;
	double chromaComponent = 0;
	double 	chromaContribution = 0;

	//int  chromaXsize = (int)(*chromaSimMatrix).size()-1;
	//int chromaYsize;
	
	DoubleVector d;	
	
	for (int x = startX;x <= endX;x++){
		d.clear();
		
		//now need the chroma part
		int chromaIndexX = min(x / conversionFactor, (int)(*chromaSimMatrix).size()-1);
		//chromaYsize = (int)(*chromaSimMatrix)[chromaIndexX].size()-1;
		
		for (int y = startY;y < endY;y++){
			chromaIndexY = min(y / conversionFactor, (int)(*chromaSimMatrix)[chromaIndexX].size()-1);
			
			//was thinking to restrict the y part too, but not working yet
			if (chromaIndexY != lastChromaYvalue){			
			
			chromaComponent = (*chromaSimMatrix)[chromaIndexX][chromaIndexY];
				lastChromaYvalue = chromaIndexY;
				chromaContribution  = chromaProportion * chromaComponent;
			}
			
			distance = (*firstEnergyVector)[x] * (*secondEnergyVector)[y];//energy similarity
			distance *= energyProportion;
			distance += chromaContribution;
			
			distance = chromaComponent;
			
			d.push_back( distance);	
		}	//end for y
		
		(*simMatrix).push_back(d);
		
	}//end for x
	
	
//	printf("..part JOINT sim size: %i, height: %i \n", (int) (*simMatrix).size(), (int) (*simMatrix)[0].size());
	
}//end self sim

void TimeWarp::calculateRestrictedCausalChromaSimilarityMatrix(DoubleMatrix& firstChromaMatrix, DoubleMatrix& secondChromaMatrix, DoubleMatrix& simMatrix, const int& Xstart, const int& Ystart, const int& Xlimit, const int& yLimit){
//calculates the chroma only similarity matrix - used to reduce computation later when doing the joint energy and chroma matrix
	DoubleVector d;
	
	int size = 0;
	if (simMatrix.size() > 0){
		size = simMatrix[0].size();
	}
	
	if (secondChromaMatrix.size() > size ){
	//	printf("sim matrix size %i and second matrix size %i Xstart %i\n", simMatrix.size(), secondChromaMatrix.size(), Xstart);
		for (int x = 0;x < Xstart;x++){
			simMatrix[x].push_back(0);
		}
		
		int endX = min(Xlimit, (int)firstChromaMatrix.size());
		
		
		for (int x = Xstart;x < endX;x++){
	
			if (x < simMatrix.size()){
				//i.e. entries already exist
				for (int y = (int)simMatrix[x].size();y < secondChromaMatrix.size();y++){
					double distance;
					if (useDotProduct)
						distance = getChromaSimilarity(x, y, &firstChromaMatrix, &secondChromaMatrix);
					else
						distance = getEuclideanDistance(x, y, &firstChromaMatrix, &secondChromaMatrix);
					
					//	printf("X %i Y %i dist %f\n", x, y, distance);
					simMatrix[x].push_back(distance);
				}
			}
			else {
				DoubleVector d;
				for (int y = 0;y < secondChromaMatrix.size();y++){
					double distance;
					if (useDotProduct)
						distance = getChromaSimilarity(x, y, &firstChromaMatrix, &secondChromaMatrix);
					else
						distance = getEuclideanDistance(x, y, &firstChromaMatrix, &secondChromaMatrix);
					
					//	printf("new row X %i Y %i dist %f\n", x, y, distance);
					d.clear();
					d.push_back( distance);	
				}
				simMatrix.push_back(d);
			}
			
		}//end for good section
		
		for (int x = endX;x < firstChromaMatrix.size(); x++){
			if (x < simMatrix.size()){
				simMatrix[x].push_back(0);
			}else{
				d.clear();
				d.push_back(0);
				simMatrix.push_back(d);
			}
		}//end for x to end
		
	}
	//	if (size > 0)
	//		printf("Causial CHROMA ONLY SIM SIZE %i x %i; ", (int)simMatrix.size(), (int) simMatrix[0].size());
	//printf("First matrix SIZE %i , SECOND size %i\n", (int)firstChromaMatrix.size(), (int) secondChromaMatrix.size());	
	

}