inline ext::shared_ptr<OneFactorModel::ShortRateDynamics>
 ExtendedCoxIngersollRoss::dynamics() const {
     return ext::shared_ptr<ShortRateDynamics>(
                         new Dynamics(phi_, theta(), k() , sigma(), x0()));
 }

 // A function that prevents us putting too much stock in small sample
 // sets.  Returns a number between 2.0 and 1.0, depending on the number
 // of samples.  5 or more samples yields one; fewer scales linearly from
 // 2.0 at 1 sample to 1.0 at 5.
 double confidence_factor(int samples) {
   if (samples > 4) return 1.0;
   else return  1.0 + sigma() * ((double)(5 - samples))/2.0;
 }

 double get_new_prediction(TruncatedSeq* seq) {
   return MAX2(seq->davg() + sigma() * seq->dsd(),
               seq->davg() * confidence_factor(seq->num()));
 }

//.........这里部分代码省略.........
    // 4-vector DOT metric 1 -1 -1 -1
    RooFormulaVar p1Dph2("p1Dph2","@0*@[email protected]",RooArgList(E1,Eph2,*p1v3Dph2));

    // 3-vector DOT
    RooFormulaVar* p2v3Dph2 = new RooFormulaVar("p2v3Dph2",
         "@0*@1*( (TMath::Cos(@2)*TMath::Cos(@3))/(TMath::Sin(@2)*TMath::Sin(@3))+TMath::Cos(@[email protected]))",
         RooArgList(*pT2,*pTph2,*theta2,*thetaph2,*phi2,*phiph2));
    // 4-vector DOT metric 1 -1 -1 -1
    RooFormulaVar p2Dph2("p2Dph2","@0*@[email protected]",RooArgList(E2,Eph2,*p2v3Dph2));

    // fsrPhoton1 DOT fsrPhoton2

    // 3-vector DOT
    RooFormulaVar* ph1v3Dph2 = new RooFormulaVar("ph1v3Dph2",
         "@0*@1*( (TMath::Cos(@2)*TMath::Cos(@3))/(TMath::Sin(@2)*TMath::Sin(@3))+TMath::Cos(@[email protected]))",
         RooArgList(*pTph1,*pTph2,*thetaph1,*thetaph2,*phiph1,*phiph2));    
    // 4-vector DOT metric 1 -1 -1 -1
    RooFormulaVar ph1Dph2("ph1Dph2","@0*@[email protected]",RooArgList(Eph1,Eph2,*ph1v3Dph2));

    // mZ1

    RooFormulaVar* mZ1;
    mZ1 = new RooFormulaVar("mZ1","TMath::Sqrt(2*@[email protected]*@[email protected]*@2)",RooArgList(p1D2,*m1,*m2));
    if(p4sZ1ph_.size()==1)
      mZ1 = new RooFormulaVar("mZ1","TMath::Sqrt(2*@0+2*@1+2*@[email protected]*@[email protected]*@4)",
                                    RooArgList(p1D2, p1Dph1, p2Dph1, *m1,*m2));
    if(p4sZ1ph_.size()==2)
      mZ1 = new RooFormulaVar("mZ1","TMath::Sqrt(2*@0+2*@1+2*@2+2*@3+2*@4+2*@[email protected]*@[email protected]*@7)",
                              RooArgList(p1D2,p1Dph1,p2Dph1,p1Dph2,p2Dph2,ph1Dph2, *m1,*m2));

    if(debug_) cout<<"mZ1 is "<<mZ1->getVal()<<endl;

    // pTerrs, 1,2,ph1,ph2
    RooRealVar sigmaZ1_1("sigmaZ1_1", "sigmaZ1_1", pTerrZ1_1);
    RooRealVar sigmaZ1_2("sigmaZ1_2", "sigmaZ1_2", pTerrZ1_2);

    RooRealVar sigmaZ1_ph1("sigmaZ1_ph1", "sigmaZ1_ph1", pTerrZ1_ph1);
    RooRealVar sigmaZ1_ph2("sigmaZ1_ph2", "sigmaZ1_ph2", pTerrZ1_ph2);

    // resolution for decay products
    RooGaussian gauss1("gauss1","gaussian PDF", *pT1RECO, *pT1, sigmaZ1_1);
    RooGaussian gauss2("gauss2","gaussian PDF", *pT2RECO, *pT2, sigmaZ1_2);

    RooGaussian gaussph1("gaussph1","gaussian PDF", *pTph1RECO, *pTph1, sigmaZ1_ph1);
    RooGaussian gaussph2("gaussph2","gaussian PDF", *pTph2RECO, *pTph2, sigmaZ1_ph2);

    RooRealVar bwMean("bwMean", "m_{Z^{0}}", 91.187);
    RooRealVar bwGamma("bwGamma", "#Gamma", 2.5);

    RooRealVar sg("sg", "sg", sgVal_);
    RooRealVar a("a", "a", aVal_);
    RooRealVar n("n", "n", nVal_);

    RooCBShape CB("CB","CB",*mZ1,bwMean,sg,a,n);
    RooRealVar f("f","f", fVal_);

    RooRealVar mean("mean","mean",meanVal_);
    RooRealVar sigma("sigma","sigma",sigmaVal_);
    RooRealVar f1("f1","f1",f1Val_);

    RooGenericPdf RelBW("RelBW","1/( pow(mZ1*mZ1-bwMean*bwMean,2)+pow(mZ1,4)*pow(bwGamma/bwMean,2) )", RooArgSet(*mZ1,bwMean,bwGamma) );

    RooAddPdf RelBWxCB("RelBWxCB","RelBWxCB", RelBW, CB, f);
    RooGaussian gauss("gauss","gauss",*mZ1,mean,sigma);
    RooAddPdf RelBWxCBxgauss("RelBWxCBxgauss","RelBWxCBxgauss", RelBWxCB, gauss, f1);

//.........这里部分代码省略.........
		  Number   p_solid = 0.;

		  grad_u_mat(0) = grad_u_mat(1) = grad_u_mat(2) = 0;
		
		  for (unsigned int d = 0; d < dim; ++d) {
			std::vector<Number> u_undefo;
			std::vector<Number> u_undefo_ref;

			//Fills the vector di with the global degree of freedom indices for the element. :dof_indicies
			
			

			
			Last_non_linear_soln.get_dof_map().dof_indices(elem, undefo_index,d);
			Last_non_linear_soln.current_local_solution->get(undefo_index, u_undefo);
			reference.current_local_solution->get(undefo_index, u_undefo_ref);

			for (unsigned int l = 0; l != n_u_dofs; l++){
			   grad_u_mat(d).add_scaled(dphi[l][qp], u_undefo[l]+u_undefo_ref[l]); 
			}
		  }
          
		  for (unsigned int l=0; l<n_p_dofs; l++)
		  {
			p_solid += psi[l][qp]*Last_non_linear_soln.current_local_solution->el(dof_indices_p[l]);
		}
		
		Point rX;
		material.init_for_qp(rX,grad_u_mat, p_solid, qp,0, p_solid,es);
		Real J=material.J;
		 Real I_1=material.I_1;
		 Real I_2=material.I_2;
		 Real I_3=material.I_3;
		 RealTensor sigma=material.sigma;
		 
		 av_pressure=av_pressure + p_solid*JxW[qp];
		 
		 /*
		 		 		std::cout<<"grad_u_mat(0)" << grad_u_mat(0) <<std::endl;

		 		std::cout<<" J " << J <<std::endl;

		std::cout<<" sigma " << sigma <<std::endl;
		*/
		 Real sigma_sum_sq=pow(sigma(0,0)*sigma(0,0)+sigma(0,1)*sigma(0,1)+sigma(0,2)*sigma(0,2)+sigma(1,0)*sigma(1,0)+sigma(1,1)*sigma(1,1)+sigma(1,2)*sigma(1,2)+sigma(2,0)*sigma(2,0)+sigma(2,1)*sigma(2,1)+sigma(2,2)*sigma(2,2),0.5);
		 
		// std::cout<<" J " << J <<std::endl;


		 sum_jac_postvars=sum_jac_postvars+JxW[qp];
 

		
		for (unsigned int i=0; i<n_u_dofs; i++){
          Fu(i) += I_1*JxW[qp]*phi[i][qp];
          Fv(i) += I_2*JxW[qp]*phi[i][qp];
          Fw(i) += I_3*JxW[qp]*phi[i][qp];

	        Fx(i) += sigma_sum_sq*JxW[qp]*phi[i][qp];
          Fy(i) += J*JxW[qp]*phi[i][qp];
          Fz(i) += 0*JxW[qp]*phi[i][qp];
    }
    
    
 
		for (unsigned int i=0; i<n_p_dofs; i++){

int main(int argc, char **argv) {
	clock_t t1, t2;
	t1 = clock();

	IloEnv env;
	IloModel model(env);
	IloCplex cplex(model);

	/************************** Defining the parameters ************************************/
	IloInt N;							//No. of nodes
	IloInt M;							//No. of calls
	Num2DMatrix links(env);				//Defines the topology of the network
	IloNumArray call_demand(env);		//link bandwidth requirement of each call
	IloNumArray call_revenue(env);		//revenue generated from the call
	IloNumArray call_origin(env);		//origin node index of each call
	IloNumArray call_destination(env);	//destination node index of each call
	Num2DMatrix Q(env);					//Bandwidth capacity of each link
	Num2DMatrix sigma(env);				//Standard deviation of service times on link (i,j)
	Num2DMatrix cv(env);				//coefficient of variation of service times on the link (i,j)
	IloNum C;							//Unit queueing delay cost per unit time
	IloNumArray R_approx_init(env);

	ifstream fin;
	const char* filename = "BPP_data_sample - Copy.txt";
	if (argc > 1)
		filename = argv[1];
		
	fin.open(filename);
	//fin.open("BPP_10node_navneet.txt");
	fin >> links >> call_origin >> call_destination >> call_demand >>
		call_revenue >> Q >> cv >> R_approx_init >> C ;
	cout << "Reading Data from the file - "<<filename<<endl;

	N = links.getSize();
	M = call_origin.getSize();

	IloInt H = R_approx_init.getSize();
	Num3DMatrix R_approx(env, N);			//The tangential linear function approximation to R.
	for (IloInt i=0; i<N; i++) {
		R_approx[i] = Num2DMatrix(env, N);
		for (IloInt j=0; j<N; j++) {
			R_approx[i][j] = IloNumArray(env, H);
			for (IloInt h=0; h<H; h++)
				R_approx[i][j][h] = R_approx_init[h];
		}
	}

	/************************** Defining the parameters ENDS ************************************/

	/************* Defining the variables defined in the model formulation **********************/
	IloNumVarArray Y(env, M, 0, 1, ILOINT); //Variable to define whether a call m is routed or not
	IloNumArray Y_sol(env, M); //Solution values

	NumVar3DMatrix X(env, N); //Variable to define whether a call m is routed along path i-j
	Num3DMatrix X_sol(env, N);
	for (IloInt i=0; i<N; i++) {
		X[i] = NumVar2DMatrix(env, N);
		X_sol[i] = Num2DMatrix(env, N);
		for (IloInt j=0; j<N; j++) {
			X[i][j] = IloNumVarArray(env, M, 0, 1, ILOINT);
			X_sol[i][j] = IloNumArray(env, M);
		}
	}

	NumVar3DMatrix W(env, N); //Variable to define whether a call m is routed along path i-j
	for (IloInt i=0; i<N; i++) {
		W[i] = NumVar2DMatrix(env, N);
		for (IloInt j=0; j<N; j++)
			W[i][j] = IloNumVarArray(env, M, 0, 1, ILOINT);
	}

	NumVar2DMatrix R(env, (IloInt)N); //The linearization Variable
	for (IloInt i=0; i<N; i++)
		R[i] = IloNumVarArray(env, (IloInt)N, 0, IloInfinity, ILOFLOAT);
	
	/************* Defining the variables defined in the model formulation ENDS *****************/

	/**************************** Defining the Constraints *******************************/
	// Constraint #1 : Flow Conservation Constraint
	for (IloInt m=0; m<M; m++) {
		for (IloInt i=0; i<N; i++) {
			IloExpr constraint1(env);
			for (IloInt j=0; j<N; j++) {
				if (links[i][j] == 1)
					constraint1 += W[i][j][m];
			}
			for (IloInt j=0; j<N; j++) {
				if (links[j][i] == 1)
					constraint1 += -W[j][i][m];
			}
			
			if (i == call_origin[m])
				model.add(constraint1 == Y[m]);
			else if (i == call_destination[m])
				model.add(constraint1 == -Y[m]);
			else 
				model.add(constraint1 == 0);

			constraint1.end();
		}
//.........这里部分代码省略.........

const Vector&  TwentyEightNodeBrickUP::getResistingForce( )

{

	int i, j, jk, k, k1;

	double xsj;

	static Matrix B(6, 3);

	double volume = 0.;



//	printf("calling getResistingForce()\n");

	resid.Zero();



	//compute basis vectors and local nodal coordinates

	computeBasis( ) ;

	//gauss loop to compute and save shape functions

	for( i = 0; i < nintu; i++ ) {

		// compute Jacobian and global shape functions

		Jacobian3d(i, xsj, 0);

		//volume element to also be saved

		dvolu[i] = wu[i] * xsj ;

		volume += dvolu[i];

	} // end for i

	//printf("volume = %f\n", volume);

	volume = 0.;

	for( i = 0; i < nintp; i++ ) {

		// compute Jacobian and global shape functions

		Jacobian3d(i, xsj, 1);

		//volume element to also be saved

		dvolp[i] = wp[i] * xsj ;

		volume += dvolp[i];

	} // end for i

	//printf("volume = %f\n", volume);



	// Loop over the integration points

	for (i = 0; i < nintu; i++) {



		// Get material stress response

		const Vector &sigma = materialPointers[i]->getStress();



		// Perform numerical integration on internal force

		//P = P + (B^ sigma) * intWt(i)*intWt(j) * detJ;

		//P.addMatrixTransposeVector(1.0, B, sigma, intWt(i)*intWt(j)*detJ);

		for (j = 0; j < nenu; j++) {

			if (j<nenp)

				jk = j*4;

			else

				jk = nenp*4 + (j-nenp)*3;



			B(0,0) = shgu[0][j][i];

			B(0,1) = 0.;

			B(0,2) = 0.;

			B(1,0) = 0.;

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

void bug_3()
{
    INT res;
    cout << "should be :  division by 0 \n";
    copy(rectangle(Pt2di(-10,-13),Pt2di(103,105)),FY/FX,sigma(res));
}

void bug_4()
{
    REAL res;
    cout << "should be :  division by 0 \n";
    copy(rectangle(Pt2di(-10,-13),Pt2di(103,105)),FY/(FX+1.0),sigma(res));
}

     cBenchCorrel (Pt2di aSz) :
         mIm1(aSz.x,aSz.y),
         mDup1(aSz.x,aSz.y),
         mPds1(aSz.x,aSz.y),
         mIm2(aSz.x,aSz.y),
         mDup2(aSz.x,aSz.y),
         mPds2(aSz.x,aSz.y)
     {


           ELISE_COPY(mIm1.all_pts(),frandr(),mIm1.out()|mDup1.out());
           ELISE_COPY(mIm2.all_pts(),frandr(),mIm2.out()|mDup2.out());


           ELISE_COPY(mIm1.all_pts(),frandr(),mPds1.out());
           ELISE_COPY(mIm1.all_pts(),frandr(),mPds2.out());



           ElFFTCorrelCirc(mDup1,mDup2);
           Im2D_REAL8 aCPad = ElFFTCorrelPadded(mIm1, mIm2);


           REAL anEps = (1+10*NRrandom3()) * 1e-2;

           REAL aRatioSurf = aSz.x * aSz.y * (1+NRrandom3()) / 6.0;

           Im2D_REAL8 aCNC = ElFFTCorrelNCPadded(mIm1, mIm2,anEps,aRatioSurf);
           Im2D_REAL8 aPdsCNC = ElFFTPonderedCorrelNCPadded
                          (
                                mIm1.in(),
                                mIm2.in(),
                                aSz,
                                mPds1.in(),
                                mPds2.in(),
                                anEps,
                                aRatioSurf
                          );

           for (INT x =-1 ; x<= aSz.x ; x++)
           {
               for (INT y =-1 ; y<= aSz.y ; y++)
               {
                    REAL aSElise;
                    ELISE_COPY
                    (
	                mIm1.all_pts(),
		          mIm1.in()[Virgule(mod(FX+x,aSz.x),mod(FY+y,aSz.y))] 
		        * mIm2.in(),
		        sigma(aSElise)
                    );

	            REAL aSFFT = mDup1.data()[mod(y,aSz.y)][mod(x,aSz.x)];

	            BENCH_ASSERT(std::abs(aSElise-aSFFT)<epsilon);

                    Pt2di aDec(x,y);
                    VerifCorrelCNC(aDec,false,aCPad, anEps,aCNC,aRatioSurf);
                    VerifCorrelCNC(aDec,true,aCPad, anEps,aPdsCNC,aRatioSurf);

               }
           }
     }

void bug_2()
{
    REAL res;
    cout << "should be :  incompatible type in Out Reduction \n";
    copy(rectangle(Pt2di(-10,-13),Pt2di(103,105)),FY+FX*3,sigma(res));
}

void main(void) {
	cout << "\nTesting blank constructor";
	tabloid t;

	cout << "\nTesting constructor of a tabloid of shape 311";
	Partition lambda;
	lambda.Add(3);
	lambda.Add(2);
	lambda.Add(1);
	tabloid s(lambda);
	
	cout << "\nGetting the standard tabloid of shape 311";
	s = standardoid(lambda);

	cout << "\nTesting output";
	cout << s;

	cout << "\nTesting Sn action";
	Perm sigma(6);
	sigma[1]=6;
	sigma[2]=3;
	sigma[3]=1;
	sigma[4]=2;
	sigma[5]=5;
	sigma[6]=4;
		cout << "\nApplying ";
	    cout << sigma;
		cout << sigma * s;

	
	cout << "\nNote that ";
	cout << s;
	cout << " has " << s.RowNumber() << " of rows";

	cout << " and the second row is \n";
	rows temp = s.Row(2);
	for (rows::const_iterator iter2 = temp.begin(); iter2 != temp.end(); ++iter2) {
		cout << " " << *iter2;
	}


// reorder the rows of a tabloid Note it no longer have the shape of a paritition

	Perm tau(3);
	tau[1]=2;
	tau[2]=3;
	tau[3]=1;
	
	cout << "\nReordering rows of ";
	cout << s;
	cout << " with ";
	cout << tau;
	tabloid q = reorderrows(tau, s);
	cout << q;

//	void insert(const vector<int> row);
	cout << "\nNow we insert a row of 7, 8, 9";
	vector<int> rowtemp;
	rowtemp.push_back(7);
	rowtemp.push_back(8);
	rowtemp.push_back(9);
	q.insert(rowtemp);
	cout << q;

//	void Add(const int row, const int value);
	cout << "\nNow we add a 10 to row 2";
	q.Add(2,10);
	cout << q;
	int crap;
	cin >> crap;
}

/*--------------------------------------------------------*/
int main(int argc, char *argv[])
{	
        char    *parfname, *inpfname, *outfname;
        int     i, nobj, nterm, ik, ik0, delta, niter;
        float   *x, *nx, *nzx, *zx,
                *y, *ny, *nzy, *zy,
                dx, dy, rr, thresh, sx, sy;
        double  *coeffx, *coeffy;
        FILE    *outf;
        PARAMS  par;

/* IO stuff */

  if (argc != 4)
  {
    printf("\n\tUSAGE: xygrid  parameter_file  input_list  coeff_file\n\n");
    exit(1);
  }

  parfname = argv[1];
  inpfname = argv[2];
  outfname = argv[3];

  readpar(parfname, &par);

  nterm = (par.ndeg+1)*(par.ndeg+2)/2;

  nobj=readcoor(inpfname, &x, &y, &zx, &zy);

  if (!(coeffx=(double *)calloc(nterm, sizeof(double))))
    errmess("calloc(coeffx)");
  if (!(coeffy=(double *)calloc(nterm, sizeof(double))))
    errmess("calloc(coeffy)");

/* make initial fit */
  fit(x, y, zx, par.ndeg, nobj, coeffx);
  fit(x, y, zy, par.ndeg, nobj, coeffy);

  sx = sigma(x, y, zx, par.ndeg, nobj, coeffx);
  sy = sigma(x, y, zy, par.ndeg, nobj, coeffy);

  thresh = par.sigmaf*par.sigmaf*(sx*sx + sy*sy);

  ik = nobj;
  delta = 1;
  niter = 0;

  if (!(nx=(float *)calloc(nobj, sizeof(float))))
    errmess("readcoor: calloc(nx)");
  if (!(ny=(float *)calloc(nobj, sizeof(float))))
    errmess("readcoor: calloc(ny)");
  if (!(nzx=(float *)calloc(nobj, sizeof(float))))
    errmess("readcoor: calloc(nzx)");
  if (!(nzy=(float *)calloc(nobj, sizeof(float))))
    errmess("readcoor: calloc(nzy)");

/* sigma clipping of the fit until MAX_NITER reached or nothing changes */

  while ((delta > 0) && (niter < par.maxniter))
  {
    ik0 = ik;
    ik  = 0;

    for (i=0; i<nobj; i++)
    {
      dx  = poly(x[i], y[i], par.ndeg, coeffx) - zx[i];
      dy  = poly(x[i], y[i], par.ndeg, coeffy) - zy[i];

      rr  = dx*dx + dy*dy;

      nx[ik]  = x[i];
      ny[ik]  = y[i];
      nzx[ik] = zx[i];
      nzy[ik] = zy[i];

      if (rr < thresh) ++ik;
    }

    delta = ik0 - ik;

    fit(nx, ny, nzx, par.ndeg, ik, coeffx);
    fit(nx, ny, nzy, par.ndeg, ik, coeffy);

    sx = sigma(nx, ny, nzx, par.ndeg, ik, coeffx);
    sy = sigma(nx, ny, nzy, par.ndeg, ik, coeffy);

    niter++;
  }

  free(x);
  free(y);
  free(zx);
  free(zy);
  free(nx);
  free(ny);
  free(nzx);
  free(nzy);

/* print results and store coefficients in binary file */
//.........这里部分代码省略.........

 inline void ExtendedCoxIngersollRoss::generateArguments() {
     phi_ = FittingParameter(termStructure(), theta(), k(), sigma(), x0());
 }

bool Elasticity::evalSol (Vector& s, const Vectors& eV, const FiniteElement& fe,
                          const Vec3& X, bool toLocal, Vec3* pdir) const
{
  if (eV.empty())
  {
    std::cerr <<" *** Elasticity::evalSol: No solutions vector."<< std::endl;
    return false;
  }
  else if (!eV.front().empty() && eV.front().size() != fe.dNdX.rows()*nsd)
  {
    std::cerr <<" *** Elasticity::evalSol: Invalid displacement vector."
	      <<"\n     size(eV) = "<< eV.front().size() <<"   size(dNdX) = "
	      << fe.dNdX.rows() <<","<< fe.dNdX.cols() << std::endl;
    return false;
  }

  // Evaluate the deformation gradient, dUdX, and/or the strain tensor, eps
  Matrix Bmat;
  Tensor dUdX(nDF);
  SymmTensor eps(nsd,axiSymmetry);
  if (!this->kinematics(eV.front(),fe.N,fe.dNdX,X.x,Bmat,dUdX,eps))
    return false;

  // Add strains due to temperature expansion, if any
  double epsT = this->getThermalStrain(eV.back(),fe.N,X);
  if (epsT != 0.0) eps -= epsT;

  // Calculate the stress tensor through the constitutive relation
  Matrix Cmat;
  SymmTensor sigma(nsd, axiSymmetry || material->isPlaneStrain()); double U;
  if (!material->evaluate(Cmat,sigma,U,fe,X,dUdX,eps))
    return false;
  else if (epsT != 0.0 && nsd == 2 && material->isPlaneStrain())
    sigma(3,3) -= material->getStiffness(X)*epsT;

  Vec3 p;
  bool havePval = false;
  if (toLocal && wantPrincipalStress)
  {
    // Calculate principal stresses and associated direction vectors
    if (sigma.size() == 4)
    {
      SymmTensor tmp(2); tmp = sigma; // discard the sigma_zz component
      havePval = pdir ? tmp.principal(p,pdir,2) : tmp.principal(p);
    }
    else
      havePval = pdir ? sigma.principal(p,pdir,2) : sigma.principal(p);

    // Congruence transformation to local coordinate system at current point
    if (locSys) sigma.transform(locSys->getTmat(X));
  }

  s = sigma;

  if (toLocal)
    s.push_back(sigma.vonMises());

  if (havePval)
  {
    s.push_back(p.x);
    s.push_back(p.y);
    if (sigma.dim() == 3)
      s.push_back(p.z);
  }

  return true;
}

void bug_6()
{
    INT res;
    cout << "should be :  incompatible dimension in  fonction coordinate \n";
    copy(rectangle(Pt2di(-10,-13),Pt2di(103,105)),FZ,sigma(res));
}

int main(){
    int i, N;
    char c, file[100];
    
    // check that gnuplot is present on the system
    int gnupl = control();
    if(gnupl == 1){
        printf("\nYou need gnuplot to graph the results.");
        printf("\nInstall it with: sudo apt-get install gnuplot\n\n");
        exit(2);
    }
    
    printf("Enter the name or path of file: ");
    fgets(file,sizeof(file),stdin);
    file[strlen(file)-1] = '\0';
    
    // the file's lines number is the number of points to be saved
    N = linesFile(file);
    if(N <= 2){
        printf("\nError: insufficient data number.\n");
        exit(2);
    }
    
    // creating data's arrays
    double *x = calloc(N,sizeof(double));
    double *y = calloc(N,sizeof(double));
    double *errors = calloc(N,sizeof(double));
    if(x == NULL || y == NULL || errors == NULL){
            perror("\nerror");
            printf("\n");
            exit(1);
    }
    
    // reading from file
    FILE *inputFile = fopen(file,"r");
    if(inputFile == NULL){
        perror("\nError");
        exit(1);
    }
    
    for(i=0; i<N; i++){
        fscanf(inputFile,"%lf %lf %lf\n",&x[i],&y[i],&errors[i]);
    }
    
    fclose(inputFile);
    
    // determine linear coefficients
    double M = Mbest(N,x,y,errors);
    double Q = Qbest(N,x,y,errors,M);
    double sigmaM = fabs(uM(N,x,y,errors,M,Q));
    double sigmaQ = fabs(uQ(N,x,errors,sigmaM));
    
    // defining best sigma(Y) and correlation coefficient
    double sigmaY = fabs(bestSigma(N,x,y,M,Q)); // <-- residuals analysis
    double cov = covariance(N,x,y);
    double cor = correlation(N,x,y);
    double lCov = linearParamCovariance(N,mean(x,N),sigma(x,N,mean(x,N)),sigmaY,M,Q);
    double lCor = linearParamCorrelation(N,x);
    
    // Chi-square test
    int freedomDegrees = N - 2; // infer 2 parameters (): M and Q
    double chi2 = 0, rChi2;
    
    for(i=0; i<N; i++){
        chi2 += pow(y[i] - ((M * x[i]) + Q),2) / pow(errors[i],2);
    }
    
    rChi2 = chi2 / freedomDegrees;
    
    printf("\nThe best linear fit Y = mX + q is:");
    printf("\nm = %.3lf\tsigma(m) = %.3lf\nq = %.3lf\tsigma(q) = %.3lf",M,sigmaM,Q,sigmaQ);
    printf("\n\nBest sigma(Y) = %.3lf",sigmaY);
    printf("\nCov(X,Y) = %.3lf",cov);
    printf("\nCor(X,Y) = %.3lf",cor);
    printf("\nCov(m,c) = %.3lf",lCov);
    printf("\nCor(m,c) = %.3lf",lCor);
    //printf("\nChi square = %.3lf",chi2);
    printf("\nReduced Chi square = %.3lf",rChi2);
    
    // interpolation and extrapolation
    int choice;
    double pointX, pointY, sigmaPointY, alpha;
    printf("\n\nDo you want to extrapolate a point with the calculated linear regression? (1 = YES | 0 = NO): ");
    scanf("%d",&choice);
    if(choice == 1){
        extrapolation(N,x,errors,sigmaY,M,Q);
    }
    
    // creating fit
    printf("\nPlotting fit...\n");
    FILE *data = fopen("data.dat","w");
    if(data == NULL){
        perror("\nError");
        exit(1);
    }
    
    // writing experimental datas
    for(i=0; i<N; i++){
        fprintf(data,"%lf %lf %lf\n",x[i],y[i],errors[i]);
    }
//.........这里部分代码省略.........

void cMMTP::ContAndBoucheTrou()
{
   int aDist32Close = 6;
   int aNbErod = 6;

   // 1- Quelques fitre morpho de base, pour calculer les points eligibles au bouche-trou
   int aLabelOut = 0;
   //int aLabelIn = 1;
   int aLabelClose = 2;
   int aLabelFront = 3;

   ELISE_COPY(mImMasqInit.all_pts(),mImMasqInit.in(),mImLabel.out());
   ELISE_COPY(mImLabel.border(2),aLabelOut,mImLabel.out());

      // 1.1 calcul des point dans le fermeture
   ELISE_COPY
   (
          select
          (
             mImLabel.all_pts(),
             close_32(mImLabel.in(0),aDist32Close) && (mImLabel.in()==aLabelOut)
          ),
          aLabelClose,
          mImLabel.out()
   );
   ELISE_COPY(mImLabel.border(2),aLabelOut,mImLabel.out());


      // 1.2 erosion de ces points
   Neighbourhood V4 = Neighbourhood::v4();
   Neighbourhood V8 = Neighbourhood::v8();
   Neigh_Rel aRelV4(V4);

   Liste_Pts_U_INT2 aLFront(2);
   ELISE_COPY
   (
          select
          (
             mImLabel.all_pts(),
             (mImLabel.in(0)==aLabelClose) &&  (aRelV4.red_max(mImLabel.in(0)==aLabelOut))
          ),
          aLabelFront,
          mImLabel.out() | aLFront
    );
    for (int aK=0 ; aK<aNbErod ; aK++)
    {
        Liste_Pts_U_INT2 aLNew(2);
        ELISE_COPY
        (
               dilate
               (
                  aLFront.all_pts(),
                  mImLabel.neigh_test_and_set(Neighbourhood::v4(),2,3,20)
               ),
               aLabelFront,
               aLNew
         );
         aLFront = aLNew;
    }
    ELISE_COPY(select(mImLabel.all_pts(),mImLabel.in()==aLabelFront),0,mImLabel.out());

    // Au cas ou on ferait un export premature
    ELISE_COPY(mImMasqFinal.all_pts(),mImLabel.in()!=0,mImMasqFinal.out());

    int aSomMaskF;
    ELISE_COPY(mImMasqFinal.all_pts(),mImLabel.in()==1,sigma(aSomMaskF));
    if (aSomMaskF < 100) return;
    // std::cout << "aSomMaskFaSomMaskF " << aSomMaskF << "\n";
   // 2- Dequantifiication, adaptee au image a trou

       Im2D_REAL4 aProfCont(mSzTiep.x,mSzTiep.y,0.0);
       {
           Im2D_INT2 aPPV = BouchePPV(mImProf,mImLabel.in()==1);

           ElImplemDequantifier aDeq(mSzTiep);
           aDeq.DoDequantif(mSzTiep,aPPV.in());
           ELISE_COPY(aProfCont.all_pts(),aDeq.ImDeqReelle(),aProfCont.out());

           ELISE_COPY(select(aProfCont.all_pts(),mImLabel.in()!=1),0,aProfCont.out());
       }

       
    //Im2D_REAL4 aImInterp(mSzTiep.x,mSzTiep.y);
    TIm2D<REAL4,REAL8> aTInterp(mContBT);

   // 3- Bouchage "fin" des trour par moinde L2
          // 3.1 Valeur initial

                 // Filtrage gaussien
    Fonc_Num aFMasq = (mImLabel.in(0)==1);
    Fonc_Num aFProf = (aProfCont.in(0) * aFMasq);
    for (int aK=0 ; aK<3 ; aK++)
    {
          aFMasq = rect_som(aFMasq,1) /9.0;
          aFProf = rect_som(aFProf,1) /9.0;
    }

    ELISE_COPY
    (
         mContBT.all_pts(),
//.........这里部分代码省略.........

CFX_PtrArray* CBC_ReedSolomonDecoder::RunEuclideanAlgorithm(
    CBC_ReedSolomonGF256Poly* a,
    CBC_ReedSolomonGF256Poly* b,
    int32_t R,
    int32_t& e) {
  if (a->GetDegree() < b->GetDegree()) {
    CBC_ReedSolomonGF256Poly* temp = a;
    a = b;
    b = temp;
  }
  CBC_ReedSolomonGF256Poly* rsg1 = a->Clone(e);
  BC_EXCEPTION_CHECK_ReturnValue(e, NULL);
  CBC_AutoPtr<CBC_ReedSolomonGF256Poly> rLast(rsg1);
  CBC_ReedSolomonGF256Poly* rsg2 = b->Clone(e);
  BC_EXCEPTION_CHECK_ReturnValue(e, NULL);
  CBC_AutoPtr<CBC_ReedSolomonGF256Poly> r(rsg2);
  CBC_ReedSolomonGF256Poly* rsg3 = m_field->GetOne()->Clone(e);
  BC_EXCEPTION_CHECK_ReturnValue(e, NULL);
  CBC_AutoPtr<CBC_ReedSolomonGF256Poly> sLast(rsg3);
  CBC_ReedSolomonGF256Poly* rsg4 = m_field->GetZero()->Clone(e);
  BC_EXCEPTION_CHECK_ReturnValue(e, NULL);
  CBC_AutoPtr<CBC_ReedSolomonGF256Poly> s(rsg4);
  CBC_ReedSolomonGF256Poly* rsg5 = m_field->GetZero()->Clone(e);
  BC_EXCEPTION_CHECK_ReturnValue(e, NULL);
  CBC_AutoPtr<CBC_ReedSolomonGF256Poly> tLast(rsg5);
  CBC_ReedSolomonGF256Poly* rsg6 = m_field->GetOne()->Clone(e);
  BC_EXCEPTION_CHECK_ReturnValue(e, NULL);
  CBC_AutoPtr<CBC_ReedSolomonGF256Poly> t(rsg6);
  while (r->GetDegree() >= R / 2) {
    CBC_AutoPtr<CBC_ReedSolomonGF256Poly> rLastLast = rLast;
    CBC_AutoPtr<CBC_ReedSolomonGF256Poly> sLastLast = sLast;
    CBC_AutoPtr<CBC_ReedSolomonGF256Poly> tLastlast = tLast;
    rLast = r;
    sLast = s;
    tLast = t;
    if (rLast->IsZero()) {
      e = BCExceptionR_I_1IsZero;
      BC_EXCEPTION_CHECK_ReturnValue(e, NULL);
    }
    CBC_ReedSolomonGF256Poly* rsg7 = rLastLast->Clone(e);
    BC_EXCEPTION_CHECK_ReturnValue(e, NULL);
    CBC_AutoPtr<CBC_ReedSolomonGF256Poly> rTemp(rsg7);
    r = rTemp;
    CBC_ReedSolomonGF256Poly* rsg8 = m_field->GetZero()->Clone(e);
    BC_EXCEPTION_CHECK_ReturnValue(e, NULL);
    CBC_AutoPtr<CBC_ReedSolomonGF256Poly> q(rsg8);
    int32_t denominatorLeadingTerm = rLast->GetCoefficients(rLast->GetDegree());
    int32_t dltInverse = m_field->Inverse(denominatorLeadingTerm, e);
    BC_EXCEPTION_CHECK_ReturnValue(e, NULL);
    while (r->GetDegree() >= rLast->GetDegree() && !(r->IsZero())) {
      int32_t degreeDiff = r->GetDegree() - rLast->GetDegree();
      int32_t scale =
          m_field->Multiply(r->GetCoefficients(r->GetDegree()), dltInverse);
      CBC_ReedSolomonGF256Poly* rsgp1 =
          m_field->BuildMonomial(degreeDiff, scale, e);
      BC_EXCEPTION_CHECK_ReturnValue(e, NULL);
      CBC_AutoPtr<CBC_ReedSolomonGF256Poly> build(rsgp1);
      CBC_ReedSolomonGF256Poly* rsgp2 = q->AddOrSubtract(build.get(), e);
      BC_EXCEPTION_CHECK_ReturnValue(e, NULL);
      CBC_AutoPtr<CBC_ReedSolomonGF256Poly> temp(rsgp2);
      q = temp;
      CBC_ReedSolomonGF256Poly* rsgp3 =
          rLast->MultiplyByMonomial(degreeDiff, scale, e);
      BC_EXCEPTION_CHECK_ReturnValue(e, NULL);
      CBC_AutoPtr<CBC_ReedSolomonGF256Poly> multiply(rsgp3);
      CBC_ReedSolomonGF256Poly* rsgp4 = r->AddOrSubtract(multiply.get(), e);
      BC_EXCEPTION_CHECK_ReturnValue(e, NULL);
      CBC_AutoPtr<CBC_ReedSolomonGF256Poly> temp3(rsgp4);
      r = temp3;
    }
    CBC_ReedSolomonGF256Poly* rsg9 = q->Multiply(sLast.get(), e);
    BC_EXCEPTION_CHECK_ReturnValue(e, NULL);
    CBC_AutoPtr<CBC_ReedSolomonGF256Poly> temp1(rsg9);
    CBC_ReedSolomonGF256Poly* rsg10 = temp1->AddOrSubtract(sLastLast.get(), e);
    BC_EXCEPTION_CHECK_ReturnValue(e, NULL);
    CBC_AutoPtr<CBC_ReedSolomonGF256Poly> temp2(rsg10);
    s = temp2;
    CBC_ReedSolomonGF256Poly* rsg11 = q->Multiply(tLast.get(), e);
    BC_EXCEPTION_CHECK_ReturnValue(e, NULL);
    CBC_AutoPtr<CBC_ReedSolomonGF256Poly> temp5(rsg11);
    CBC_ReedSolomonGF256Poly* rsg12 = temp5->AddOrSubtract(tLastlast.get(), e);
    BC_EXCEPTION_CHECK_ReturnValue(e, NULL);
    CBC_AutoPtr<CBC_ReedSolomonGF256Poly> temp6(rsg12);
    t = temp6;
  }
  int32_t sigmaTildeAtZero = t->GetCoefficients(0);
  if (sigmaTildeAtZero == 0) {
    e = BCExceptionIsZero;
    BC_EXCEPTION_CHECK_ReturnValue(e, NULL);
  }
  int32_t inverse = m_field->Inverse(sigmaTildeAtZero, e);
  BC_EXCEPTION_CHECK_ReturnValue(e, NULL);
  CBC_ReedSolomonGF256Poly* rsg13 = t->Multiply(inverse, e);
  BC_EXCEPTION_CHECK_ReturnValue(e, NULL);
  CBC_AutoPtr<CBC_ReedSolomonGF256Poly> sigma(rsg13);
  CBC_ReedSolomonGF256Poly* rsg14 = r->Multiply(inverse, e);
  BC_EXCEPTION_CHECK_ReturnValue(e, NULL);
  CBC_AutoPtr<CBC_ReedSolomonGF256Poly> omega(rsg14);
  CFX_PtrArray* temp = new CFX_PtrArray;
  temp->Add(sigma.release());
//.........这里部分代码省略.........

int

TwentyEightNodeBrickUP::getResponse(int responseID, Information &eleInfo)

{

	static Vector stresses(nintu*6);



	if (responseID == 1)

		return eleInfo.setVector(this->getResistingForce());



	else if (responseID == 2)

		return eleInfo.setMatrix(this->getTangentStiff());



    else if (responseID == 3)

        return eleInfo.setMatrix(this->getMass());



    else if (responseID == 4)

        return eleInfo.setMatrix(this->getDamp());



	else if (responseID == 5) {



		// Loop over the integration points

		int cnt = 0;

		for (int i = 0; i < nintu; i++) {

			// Get material stress response

			const Vector &sigma = materialPointers[i]->getStress();

			stresses(cnt++) = sigma(0);

			stresses(cnt++) = sigma(1);

			stresses(cnt++) = sigma(2);

			stresses(cnt++) = sigma(3);

			stresses(cnt++) = sigma(4);

			stresses(cnt++) = sigma(5);

		}

		return eleInfo.setVector(stresses);



	}

	else



		return -1;

}