C++ x0函数代码示例

LinAlg::Matrix<float> SistemasLineares::GaussJacobi(LinAlg::Matrix<float> MatrizUni, unsigned MaxIterations)
{
    //Matriz Resposta
    LinAlg::Matrix<float> MatrizRes(MaxIterations, MatrizUni.getNumberOfColumns());
    LinAlg::Matrix<float> C (MatrizUni.getNumberOfRows(), MatrizUni.getNumberOfColumns() - 1);
    LinAlg::Matrix<float> g (MatrizUni.getNumberOfRows(), 1);
    LinAlg::Matrix<float> x0(C.getNumberOfColumns(), 1);

    //    //Deixa o vetor de chute inicial padronizado como vetor linha
    if(this->X0.getNumberOfColumns() < this->X0.getNumberOfRows())
        ~this->X0;

    //    //Insere o chute inicial na Matriz resposta
    for(unsigned i = 1; i < MatrizRes.getNumberOfColumns() - 1; ++i)
        x0(1,i) = this->X0(1,i);

        //Laço para contar as linhas da MatrizUni e Matriz C.
        for(unsigned i = 1; i <= MatrizUni.getNumberOfRows(); ++i)
        {   //Laço para contar as colunas da MAtrizUni e Matriz C.
            for(unsigned j = 1; j < MatrizUni.getNumberOfColumns(); ++j)
            {
                if(i != j)
                    C(i,j) = - MatrizUni(i,j)/MatrizUni(i,i);//Matriz com a diagonal zerada.
            }
            g(i,1) = MatrizUni(i,MatrizUni.getNumberOfColumns()) / MatrizUni(i,i);//Matriz dos termos independentes.
        }

        MatrizRes = ~x0;
        for(unsigned k = 1; k < MaxIterations; ++k)
        {
            x0 =  (C * x0) + g;
            MatrizRes = MatrizRes || ~x0;
        }
    return MatrizRes;
}

void CVMath::setupMatrixM1()
{
	Line r1 = Points::getInstance().getR1();
	Line r2 = Points::getInstance().getR2();
	Line r3 = Points::getInstance().getR3();
	Line r4 = Points::getInstance().getR4();
	
	l0 << r1.x, r1.y, r1.z;
	l1 << r2.x, r2.y, r2.z;
	l2 << r3.x, r3.y, r3.z;
	l3 << r4.x, r4.y, r4.z;

	x0 = l0.cross(l1);
	x0 << x0(0)/x0(2), x0(1)/x0(2), 1;
	x1 = l2.cross(l3);
	x1 << x1(0)/x1(2), x1(1)/x1(2), 1;

	l =  x0.cross(x1);

	Hp = MatrixXd(3, 3);
	Hp << 1, 0, 0,
	     0, 1, 0,
	     l(0), l(1), l(2);
	std::cout << Hp << std::endl;
	Vector3d l0;
	Vector3d l1;
	Hp_INV = MatrixXd(3, 3);
	Hp_INV = Hp.inverse().eval();
	std::cout << Hp_INV << std::endl;
}

Vector Initialize()
{
 	Vector x0(3, 0.0);

	x0(0) = 10000;
	return x0;
} 

void MyWidget::process()
{
    ui->plainTextEdit->clear();
    variant = ui->comboBox->currentIndex();
    print(QString("Variant %1").arg(variant));
    print(QString("[%1; %2]").arg(x0(), 0, 'g', 5).arg(xFinish(), 0, 'g', 3));

    //calc

    qreal x = x0(), y = y0(), z = z0();
    for (int i = 1; x < xFinish() || qFuzzyCompare(x, xFinish()); ++i)
    {
        qreal k1 = h() * f(x, y, z);
        qreal q1 = h() * g(x, y, z);

        qreal k2 = h() * f(x + h() / 2.0, y + k1 / 2.0, z + q1 / 2.0);
        qreal q2 = h() * g(x + h() / 2.0, y + k1 / 2.0, z + q1 / 2.0);

        qreal k3 = h() * f(x + h() / 2.0, y + k2 / 2.0, z + q2 / 2.0);
        qreal q3 = h() * g(x + h() / 2.0, y + k2 / 2.0, z + q2 / 2.0);

        qreal k4 = h() * f(x + h(), y + k3, z + q3);
        qreal q4 = h() * g(x + h(), y + k3, z + q3);

        print(QString("#%1").arg(i));
        print(QString("y(%1) = %2").arg(x, 0, 'g', 5).arg(y, 0, 'g', 5));
        print(QString("y_ex(%1) = %2").arg(x, 0, 'g',5).arg(yt(x), 0, 'g', 5));
        print(QString("error(%1) = %2").arg(x, 0, 'g',5).arg(ypo(y, x), 0, 'g', 5));

        x += h();
        y += (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0;
        z += (q1 + 2.0 * q2 + 2.0 * q3 + q4) / 6.0;

    }
}

         void FoncATrou_OPB_Comp<Type>::calc_buf
         (
               INT **     output,
               INT ***    input
         )
{
    for (INT d= 0; d<dim_out() ; d++)
    {
        convert
        (
             output[d]+x0(),
             input[std::min(d,dim_in()-1)][0]+x0(),
             tx()
        );
    }

     INT a0 = _cpt[ycur()];
     INT a1 = _cpt[ycur()+1];
     for (INT a=a0 ; a<a1 ; a++)
     {
         INT x = _tx[a];
         for (INT d= 0; d<dim_out() ; d++)
             output[d][x] = _v[d][a];
     }
}

void xRungeKutta4 :: halfTest (void)
{
	gpstk::Matrix<double> x0(2,1), truncError(2,1);
	x0(0,0) = 0.001; // Initial angle in radians
	x0(1,0) = 0.0; // Initial angular velocity in radians/second
	
	PendulumIntegrator pModel(x0,0.);
	
	double g = 9.81, L = 1.0;
	pModel.setPhysics(g,L);
	
	double deltaT = .00001;  // Step size in seconds for integrator
	
	double time = 0;
	double Nper = 2.5; // number of periods
	
	double addError = 0; //Total Error for angle
	double addDotError = 0; //Total Error for rate of change in angle
	
	long count = 0;
	
	while (pModel.getTime() < Nper * (2*3.14159265/sqrt(g/L)))
	{
		pModel.integrateTo((count++)*deltaT,truncError);
		addError += fabs(truncError(0,0));
		addDotError += fabs(truncError(1,0));
	}
	
	CPPUNIT_ASSERT_DOUBLES_EQUAL(-x0(0,0),pModel.getState()(0,0),addError*2);
	CPPUNIT_ASSERT_DOUBLES_EQUAL(x0(1,0),pModel.getState()(1,0),addDotError*2);
	
}

mat nnls_solver(const mat & H, mat mu, const umat & mask, int max_iter, double rel_tol, int n_threads)
{
	/****************************************************************************************************
	 * Description: sequential Coordinate-wise algorithm for non-negative least square regression problem
	 * 		A x = b, s.t. x[!m] >= 0, x[m] == 0
	 * Arguments:
	 * 	H         : A^T * A
	 * 	mu        : -A^T * b
	 * 	mask      : a mask matrix (m) of same dim of x
	 * 	max_iter  : maximum number of iterations
	 * 	rel_tol   : stop criterion, minimum change on x between two successive iteration
	 * 	n_threads : number of threads
	 * Return:
	 * 	x : solution to argmin_{x, x>=0} ||Ax - b||_F^2
	 * Reference:
	 * 	http://cmp.felk.cvut.cz/ftp/articles/franc/Franc-TR-2005-06.pdf
	 * Author:
	 * 	Eric Xihui Lin <[email protected]>
	 * Version:
	 * 	2015-11-16
	 ****************************************************************************************************/

	mat x(H.n_cols, mu.n_cols, fill::zeros);
	if (n_threads < 0) n_threads = 0;
	bool is_masked = !mask.empty();


	#pragma omp parallel for num_threads(n_threads) schedule(dynamic)
	for (int j = 0; j < mu.n_cols; j++)
	{
		if (is_masked && arma::all(mask.col(j))) 
			continue;
		vec x0(H.n_cols);
		// x0.fill(-9999);
		double tmp;
		int i = 0;
		double err1, err2 = 9999;
		do {
			// break if all entries of col_j are masked
			x0 = x.col(j);
			err1 = err2;
			err2 = 0;
			for (int k = 0; k < H.n_cols; k++)
			{
				if (is_masked && mask(k,j) > 0) continue;
				tmp = x(k,j) - mu(k,j) / H(k,k);
				if (tmp < 0) tmp = 0;
				if (tmp != x(k,j))
				{
					mu.col(j) += (tmp - x(k,j)) * H.col(k);
				}
				x(k,j) = tmp;
				tmp = std::abs(x(k,j) - x0(k));
				if (tmp > err2) err2 = tmp;
			}
		} while(++i < max_iter && std::abs(err1 - err2) / (err1 + 1e-9) > rel_tol);
	}
	return x;
}

Vector Initialize()
{
 	Vector x0(3, 0.0);

	x0(0) = 400;
	x0(1) = 798; 
	return x0;
} 

void MinimizerThread::run()
{
  boost::numeric::ublas::vector< double > x0(2);
  x0(0) = x0_;
  x0(1) = y0_;
  
  if(!solver_->isExternalSolver())
  {
    NativeSolver *nativeSolver = dynamic_cast< NativeSolver * >(solver_);
    
    if(C_ != NULL)
      solver_->setup(*objFunc_, x0, Solver::DefaultSetup(), *C_);
    else
      solver_->setup(*objFunc_, x0, Solver::DefaultSetup(), NoConstraints());
    results_.iterates.clear();
    results_.iterates.push_back(x0);
    if(nativeSolver->getM() > 1)
    {
      results_.iteratePointSets.clear();
      results_.iteratePointSets.push_back(nativeSolver->getXArray());
    }
    
    unsigned int nIter = 0;
    NativeSolver::IterationStatus status = NativeSolver::ITERATION_CONTINUE;
    while(status == NativeSolver::ITERATION_CONTINUE && 
          (nativeSolver->hasBuiltInStoppingCriterion() || 
           !stopCrit_->test(*nativeSolver)) && 
          nIter < 10000 && !terminationFlag_)
    {
      status = nativeSolver->iterate();
      results_.iterates.push_back(nativeSolver->getX());
      if(nativeSolver->getM() > 1)
        results_.iteratePointSets.push_back(nativeSolver->getXArray());
      nIter++;
    }
    results_.xMin = nativeSolver->getX()[0];
    results_.yMin = nativeSolver->getX()[1];
    results_.fMin = nativeSolver->getFVal();
    results_.numIter = nativeSolver->getNumIter();
    results_.numFuncEval = nativeSolver->getNumFuncEval();
    results_.numGradEval = nativeSolver->getNumGradEval();
  }
  else
  {
    // TODO: stopping criterion for LMBM?
    boost::shared_ptr< Solver::Results > results = 
      solver_->solve(*objFunc_, x0, GradNormTest(1e-8));
    
    results_.xMin        = results->xMin[0];
    results_.yMin        = results->xMin[1];
    results_.fMin        = results->fMin;
    results_.numIter     = results->numIter;
    results_.numFuncEval = results->numFuncEval;
    results_.numGradEval = results->numGradEval;
    //results_.iterates    = results->iterates; // TODO: iterates?
  }
}

//..
// Our verification program simply instantiate several 'MyGenericContainer'
// templates with the two test types above, and checks that the allocator
// slot is as expected:
//..
    int usageExample()
    {
        bslma::TestAllocator ta0;
        bslma::TestAllocator ta1;
//..
// With 'MyTestTypeWithNoBslmaAllocatorTraits', the slot should never be set.
//..
        MyTestTypeWithNoBslmaAllocatorTraits x;

        allocSlot = &ta0;
        MyGenericContainer<MyTestTypeWithNoBslmaAllocatorTraits> x0(x);
        ASSERT(&ta0 == allocSlot);

        allocSlot = &ta0;
        MyGenericContainer<MyTestTypeWithNoBslmaAllocatorTraits> x1(x, &ta1);
        ASSERT(&ta0 == allocSlot);
//..
// With 'MyTestTypeWithBslmaAllocatorTraits', the slot should be set to the
// allocator argument, or to 0 if not specified:
//..
        MyTestTypeWithBslmaAllocatorTraits y;

        allocSlot = &ta0;
        MyGenericContainer<MyTestTypeWithBslmaAllocatorTraits> y0(y);
        ASSERT(0 == allocSlot);

        allocSlot = &ta0;
        MyGenericContainer<MyTestTypeWithBslmaAllocatorTraits> y1(y, &ta1);
        ASSERT(&ta1 == allocSlot);

        return 0;
    }

    void Bidomain::DoOdeRhs(
            const Array<OneD, const  Array<OneD, NekDouble> >&inarray,
            Array<OneD,        Array<OneD, NekDouble> >&outarray,
            const NekDouble time)
    {
        int nq = m_fields[0]->GetNpoints();
        m_cell->TimeIntegrate(inarray, outarray, time);
        if (m_stimDuration > 0 && time < m_stimDuration)
        {
            Array<OneD,NekDouble> x0(nq);
            Array<OneD,NekDouble> x1(nq);
            Array<OneD,NekDouble> x2(nq);
            Array<OneD,NekDouble> result(nq);

            // get the coordinates
            m_fields[0]->GetCoords(x0,x1,x2);

            LibUtilities::EquationSharedPtr ifunc
                    = m_session->GetFunction("Stimulus", "u");
            ifunc->Evaluate(x0,x1,x2,time, result);

            Vmath::Vadd(nq, outarray[0], 1, result, 1, outarray[0], 1);
        }        
        Vmath::Smul(nq, 1.0/m_capMembrane, outarray[0], 1, outarray[0], 1);
    }

void MainWindow::plot3(QCustomPlot *customPlot)
{
  QVector<double> x0(pts), y0(pts);

  update_w(w,1.0);
  double y_inc = 1.0; 
  for (int i=0; i<pts; ++i)   {
    x0[i] = (double)0.5*i/pts;
    y0[i] = y_inc*(w[i]);
  }

  customPlot->graph()->setData(x0, y0);
  customPlot->graph()->setLineStyle(QCPGraph::lsLine); 
  customPlot->graph()->setScatterStyle(QCPScatterStyle::ssNone); // (QCP::ScatterStyle)(rand()%9+1));
  QPen graphPen;
  graphPen.setColor(QColor(Qt::GlobalColor(5+graph_counter)));
  graphPen.setWidthF(2);
  customPlot->graph()->setPen(graphPen);
  customPlot->xAxis->setRange(0,0.5);
  customPlot->yAxis->setRange(-90,10);
  customPlot->legend->setVisible(true);
  //  customPlot->setInteraction(QCustomPlot::iSelectPlottables);
  customPlot->graph()->setName(QString(shape.c_str()));
  customPlot->replot();
}

QgsLineString *QgsMapToolAddRectangle::rectangleToLinestring( const bool isOriented ) const
{
  std::unique_ptr<QgsLineString> ext( new QgsLineString() );
  if ( mRectangle.isEmpty() )
  {
    return ext.release();
  }

  QgsPoint x0( mRectangle.xMinimum(), mRectangle.yMinimum() );

  QgsPoint x1, x2, x3;
  if ( isOriented )
  {
    const double perpendicular = 90.0 * mSide;
    x1 = x0.project( mDistance1, mAzimuth );
    x3 = x0.project( mDistance2, mAzimuth + perpendicular );
    x2 = x1.project( mDistance2, mAzimuth + perpendicular );
  }
  else
  {
    x1 = QgsPoint( mRectangle.xMinimum(), mRectangle.yMaximum() );
    x2 = QgsPoint( mRectangle.xMaximum(), mRectangle.yMaximum() );
    x3 = QgsPoint( mRectangle.xMaximum(), mRectangle.yMinimum() );
  }
  ext->addVertex( x0 );
  ext->addVertex( x1 );
  ext->addVertex( x2 );
  ext->addVertex( x3 );
  ext->addVertex( x0 );

  return ext.release();
}

RANDOM_HILL_CLIMBING::Parameter RANDOM_HILL_CLIMBING::move(RANDOM_HILL_CLIMBING::Parameter par)
{
	boost::variate_generator<RNG&, boost::uniform_real<Float_t> > r(_rng, boost::uniform_real<Float_t>(-1,1));
	Parameter result;
	
	Float_t factor=4;
	
	Float_t maxpos=0.1*factor;
	Float_t maxang=osg::PI/45*factor;
	do
	{
		result._pos=Vec3_t(r(), r(), r())*maxpos;
	} while (result._pos.length()>maxpos);
	result._pos=par._pos+result._pos;
	
	boost::variate_generator<RNG&, boost::uniform_real<Float_t> > rsimple(_rng, boost::uniform_real<Float_t>(0,1));
	Float_t x0(rsimple());
	Float_t t1(2*osg::PI*rsimple()), t2(2*osg::PI*rsimple());
	Float_t s1(sin(t1)), s2(sin(t2)), c1(cos(t1)), c2(cos(t2));
	Float_t r1(sqrt(1-x0)), r2(sqrt(x0));
	Quat_t randomori=Quat_t(s1*r1, c1*r1, s2*r2, c2*r2);
	result._ori.slerp(maxang, Quat_t(), randomori);
	Quat_t temp=par._ori*result._ori;
	result._ori=temp;
	
	// bond rotatio
	result._bond=vector<Float_t>(_fit->RotCount());
	for (unsigned int i=0; i!=_fit->RotCount(); ++i)
	{
		result._bond[i]=par._bond[i]+r()*maxang;
	}
	return result;
}

int gauss(dvec& x, dvec& w) {

    int n = x.size();
    double dist = 1;
    double tol = 1e-15;
    int iter = 0;

    // Use Chebyshev-Gauss nodes and initial guess
    initguess(x);
//    chebnodes(x);

    dvec x0(x);
    dvec L(n,0.0);
    dvec Lp(n,0.0);
    dvec a(n,0.0);
    dvec b(n,0.0);

    rec_legendre(a,b);

    // Iteratively correct nodes with Newton's method
    while(dist>tol) {     
        newton_step(a,b,x,L,Lp);
        dist = dist_max(x,x0); 
        ++iter;
        x0 = x;
    } 

    // Compute Weights
    for(int i=0;i<n;++i){
        w[i] = 2.0/((1-x[i]*x[i])*(Lp[i]*Lp[i]));
    }

    return iter;
}