C++ valarray::size方法代码示例

void Covariance::setCovIndPar( const valarray<double>& indParNew )
{
  // If values have been cached, then check the elements of the 
  // parameter vector that affect this covariance matrix to see
  // if the cached values are still valid. 
  if ( isCacheStatusValid() )
  {
    if ( !isDifferent( covIndPar, indParNew ) )
    { 
      return;   // Cached values are ok.
    }
    else
    {
      setCacheStatusInvalid();
    }
  }

  covIndPar.resize( indParNew.size() );
  covIndPar = indParNew;
}

  vector<double> autocovariance(const valarray<double>& x,unsigned max)
  {
    const int N = x.size();

    // specify sample size if unspecified
    if (max==0) 
      max = std::max(1+N/4,N-15);

    if (max >= N) 
      max = N-1;

    // compute mean of X
    double mean = 0;
    for(int i=0;i<N;i++)
      mean += x[i];
    mean /= N;

    // allocate covariances
    vector<double> rho(max);

    // compute each autocorrelation rho[k]
    double limit = 0.01/N;
    for(int k=0;k<max;k++) 
    {
      double total = 0;
      for(int i=0;i<N-k;i++)
	total += (x[i]-mean)*(x[i+k]-mean);

      rho[k] = total/(N-k);

      if (rho[k] < limit and k>0) {
	rho.resize(k);
	break;
      }
    }

    return rho;
  }

inline void sortArray(valarray<double>& x, vector<int> &id)
{   
	int i,j,k,t1;
	int n=x.size();
	double t2;
	for(i=0;i<n;i++)
		id.push_back(i);
	for(i=0;i<n-1;i++)
	{
		k=i;
		for(j=i+1;j<n;j++) 
			if(x[j]>x[k]) 
				k=j;
		if(k!=i)
		{
			t2=x[i]; 
			x[i]=x[k]; 
			x[k]=t2;
			t1=id[i];
			id[i]=id[k];
			id[k]=t1;
		}
	}
}

const valarray<double> elsq_x(
                    const valarray<double> &z,      // n size vector
                    const valarray<double> &h,      // n size vector
                    const valarray<double> &Q,      // n by n symmetric, positive definite
                    const valarray<double> &Qinv,   // n by n symmetric, positive definite
                    const valarray<double> &h_x,    // n by nX matrix
                    const valarray<double> &Q_x     // n^2 by nX matrix
                    )
{
  const int n = z.size();
  assert( h.size() == n );
  assert( Q.size() == n * n );
  const int nX = h_x.size() / n;
  assert( n * nX == h_x.size() );
  assert( Q_x.size() == n * n * nX );

  valarray<double> residual = z - h;
  valarray<double> residualTran = residual;   // 1 by n
  double  val;

  valarray<double> W = multiply(residualTran, n, Qinv, n);  // 1 by n
  valarray<double> rvecQinvTrans = rvec( Qinv, n );   // 1 by n^2: don't really need to transpose because it's just an array

  valarray<double> term1 = multiply(rvecQinvTrans, n*n, Q_x, nX );  // 1 by nX 

  valarray<double> term2 = multiply(W, n, h_x, nX);
  assert(term2.size()== nX);

  //
  // The loop below untangles the following statement
  // for eliminating matrix object constructions and taking advantage of
  // Q being symmetric.
  //
  //    valarray<double> term3  = AkronBtimesC( W, n, rvecQinvTrans, n, Q_x, n);
  // 
  // Let w be (z-h).
  // 
  // The orignal formula
  //    0.5 [w^T kron w^T] partial_x(Qinv) --- (eq. 1)
  //
  // is equivalent to:
  //    0.5 partial_x(k) [ w^T Qinv w ] --- (eq. 2)
  //   =0.5 SUM [ w(i) w(j) partial_x(k)(Q(i,j)) ], over 1<=i<=m and 1<=j<=m.
  //
  // For Q being symmetric, the following is true:
  //    w(i) w(j) partial_x(k)(Q(i,j)) == w(j) w(i) partial_x(k)(Q(j,i))
  //
  //
  valarray<double> term3( 0.0, nX );  // 1 by nX

  for( int k=0; k<nX; k++ )
  {
      for( int j=0; j<n; j++ )
      {
          for( int i=0; i<n; i++ )
          {
              if( i<=j )
              {
                  val = W[i] * W[j] * Q_x[ j*n+i+k*n*n ] / 2.0;
                  term3[k] += val;
                  if( i<j )
                      term3[k] += val;
              }
          }
      }
  }
  
  // Returns a 1 by nX matrix.
  return 0.5 * term1 - term2 - term3;
}

void print(const valarray<T>& v) {
   for(size_t i = 0; i < v.size(); ++i) 
     cout << '\t' << v[i];
   cout << endl;
}

void indStatistics( const valarray<double>&  indPar,
                    const valarray<double>&  dataMean_indPar,
                    const valarray<double>&  dataVariance_indPar,
                    const valarray<double>&  dataVarianceInv,
                    valarray<double>*        indParCovOut,
                    valarray<double>*        indParSEOut,                          
                    valarray<double>*        indParCorOut,
		    valarray<double>*        indParCVOut,
                    valarray<double>*        indParCIOut )
{
    using std::endl;
    using std::ends;
    //----------------------------------------------------------------
    // Preliminaries.
    //----------------------------------------------------------------
    // Return if there are no output values to compute.
    if( indParCovOut == 0 && indParSEOut == 0 && indParCorOut == 0 && 
		indParCVOut && indParCIOut == 0 ) 
       return;
  
    // Number of individual parameters
    const int nB = static_cast<int>( indPar.size() );
    assert( nB > 0 );

    // Number of data points
    const int nY = static_cast<int>( dataMean_indPar.size() ) / nB;
    assert( static_cast<int>( dataMean_indPar.size() ) == nY * nB );
    assert( nY > 0 );

    using namespace std;
    // Degree of freedom
    const int nFree = nY - nB;
    
    if( indParCIOut )
      {
	if( nFree < 1 )
	  {
	    const int max = SpkError::maxMessageLen();
	    char message[max];
	    sprintf( message, "The degree of freedom (#of measurements<%d> - #of random effects<%d>) must be positive.", nY, nB );
	    
	    throw SpkException(
			       SpkError::SPK_USER_INPUT_ERR, 
			       message,
			       __LINE__, __FILE__
			       );
	  }
      }
    //----------------------------------------------------------------
    // Calculate Covariance of individual parameter estimates 
    //----------------------------------------------------------------
    valarray<double> indParCov( nB * nB );
    
    try
      {
        indParCov = inverse( 0.5 * multiply( transpose( dataVariance_indPar, nB ), nY * nY,
					     AkronBtimesC( dataVarianceInv, nY, dataVarianceInv, 
							   nY, dataVariance_indPar, nB ), nB )
			     + multiply( transpose( dataMean_indPar, nB ), nY,
					 multiply( dataVarianceInv, nY, dataMean_indPar, nB ), nB ), 
			     nB );
      }
    catch(SpkException& e)
      {
        throw e.push( SpkError::SPK_NOT_INVERTABLE_ERR,
                      "Failed to invert information matrix",
                      __LINE__, __FILE__ );
      }
    
    if( indParCovOut != 0 )
      *indParCovOut = indParCov;
    statistics( indPar, indParCov, nFree, indParSEOut, indParCorOut, indParCVOut, indParCIOut );
}

valarray<T> rbgauss_seidel(const matrix_crs<T>& A, const valarray<T>& f, 
      const valarray<T>& v0, const T resid_tol, int& num_itr) {
   // Red-Black Gauss-Seidel method driver
   // pass in the matrix, RHS vector
   // pass in an initial guess (there is another interface which creates a
   //    random initial guess)
   // num_itr is the number of iterations to do.  defaults to -1, which will
   //   run until we reach the maximum number of iterations or reach tolerance
   //   On return, num_itr will have the number of iterations completed.

   //T resid_tol = static_cast<T>(10e-8);
   
   unsigned num_its_done;
   
   // check sizes
   if ( A.m != f.size() ) {
      cerr << "error: classical_solvers:rbgauss_seidel: dimension mismatch"
           << endl;
      exit(-1);
   }

   // it is assumed that A is odd x odd, otherwise RB GS won't work
   if ( A.m % 2 != 1 ) {
      cerr << "error: classical_solvers:rbgauss_seidel: A has even dimensions"
           << endl;
      exit(-1);
   }

   // random initial guess
   valarray<T> v = v0;
   valarray<T> resid = f-A*v;

   // number of iterations to do
   if ( num_itr == -1 ) { // just do normal solve until we reach tol or 
                          // _RBGS_MAX_ITR_
   
      // do iterations
      // norm(...,0) is infinity norm
      num_its_done = 0;
      while ( norm(resid,0) > resid_tol && num_its_done < _RBGS_MAX_ITR_) {
         
         rbgauss_seidel_it(A,f,v);
   
         resid = f-A*v;
         ++num_its_done;
      }
   }

   else if ( num_itr >= 0 ) { // do exactly num_itr iterations
      // do iterations
      // norm(...,0) is infinity norm
      unsigned num_its_todo = static_cast<unsigned>(num_itr);
      num_its_done = 0;
      while ( num_its_done < num_its_todo ) {
         
         rbgauss_seidel_it(A,f,v);
   
         ++num_its_done;
      }
   }
   
   else {
      cerr << "error: classical_solvers:rbgauss_seidel: bad input number of "
           << "iterations" << endl;
      exit(-1);
   }

   num_itr = static_cast<int>(num_its_done);
   return v;
}

 unsigned size() const {return f.size();}

 void trace(double& tr,const valarray<double>& M){
   int Dim = sqrt(double(M.size()));
   tr = 0.0;
   for(int i=0; i<Dim; ++i) tr+= M[Dim*i+i];    
 }

	void setX(valarray<double> _x){
				x.resize(_x.size());
				x=_x;
				};

	void setY(valarray<double> _y){
			y.resize(_y.size());
			y=_y;
			};

Apertures::Apertures(const Ntuple* nt,unsigned int verbose) // constructor
// 1st cleanup the aperture information directly from madx tfs and keep the result as vector<double> svec, aper_1_vec;   --- more cleanup possible later
{
  if(verbose) cout << __FILE__ << " " << __PRETTY_FUNCTION__ << " line " << __LINE__ << " constructor start verbose=" << verbose << " nt->Noent()=" << nt->Noent() << '\n';

  // copy to local const valarray's
  const valarray<double>      s_array=nt->GetVar("S");        // http://www.cplusplus.com/reference/valarray/
  const valarray<double> aper_1_array=nt->GetVar("APER_1");
  const valarray<double> aper_2_array=nt->GetVar("APER_2");
  const valarray<double> aper_3_array=nt->GetVar("APER_3");
  const valarray<double> aper_4_array=nt->GetVar("APER_4");
  bool has_betx=nt->VarExists("BETX");
  bool has_bety=nt->VarExists("BETY");
  valarray<double> betx_array(nt->Nvar()),bety_array(nt->Nvar());
  if(has_betx) betx_array=nt->GetVar("BETX");
  if(has_bety) bety_array=nt->GetVar("BETY");
  const valarray<double>      z_array=nt->GetVar("Z");
  const valarray<double>      y_array=nt->GetVar("Y");   // y directly
  const valarray<double>    x_2_array=nt->GetVar("X_2"); // x_2  is x from survey in global Euclidian coordinates

  const valarray<double> aper_x_n = aper_1_array / sqrt(betx_array); // normalized aperture in x           not working correctly in CINT
  const valarray<double> aper_y_n = aper_1_array / sqrt(bety_array); // normalized aperture in y           not working correctly in CINT

  cout << std::defaultfloat;

  bool HasX2=nt->VarExists("X_2");
  if(verbose>2)
  {
    cout
    << " aper_x_n.size()=" << aper_x_n.size()
    << " aper_1_array[0]=" << aper_1_array[0]
    << " betx_array[0]=" << betx_array[0]
    << " aper_x_n[0]=" << aper_x_n[0]
    << '\n';
    if(HasX2) cout << "X_2 exists" << '\n';
  }

  this->verbose=verbose;
  if(verbose>1)
  {
    cout << "Apertures constructor"
    << " s_array.size()=" << s_array.size()
    << " aper_1_array.size()=" << aper_1_array.size() << '\n';
  }
  double aper_x_n_min=numeric_limits<double>::max();
  double aper_y_n_min=numeric_limits<double>::max();
  isigx_min=0;
  isigy_min=0;
  vector<string> nt_StrCol_Name=nt->GetStrCol("NAME");
  vector<string> nt_StrCol_Keyw=nt->GetStrCol("KEYWORD");
  vector<string> nt_StrCol_Type=nt->GetStrCol("APERTYPE");
  if(verbose) cout
    << " nt_StrCol_Name.size()=" << nt_StrCol_Name.size()
    << " nt_StrCol_Keyw.size()=" << nt_StrCol_Keyw.size()
    << " nt_StrCol_Type.size()=" << nt_StrCol_Name.size()
    << '\n';
  // verbose=3; // CSPE
  for(unsigned int j=0;j<s_array.size();++j)
  {
    bool AperOK=aper_1_array[j]>0 && aper_1_array[j]<5;
    if( AperOK ) // reasonable between 0 and 5 m, aperture information available
    {
      AperName.push_back(nt_StrCol_Name[j]); // The name should always be defined
      if(j<nt_StrCol_Keyw.size()) AperKeyw.push_back(nt_StrCol_Keyw[j]);
      if(j<nt_StrCol_Type.size()) AperType.push_back(nt_StrCol_Type[j]);
      svec.push_back( s_array[j]);
      yvec.push_back( y_array[j]);
      aper_1_vec.push_back( aper_1_array[j]);
      if(aper_2_array.size()>0) aper_2_vec.push_back( aper_2_array[j]);
      if(aper_3_array.size()>0) aper_3_vec.push_back( aper_3_array[j]);
      if(aper_4_array.size()>0) aper_4_vec.push_back( aper_4_array[j]);
      if(HasX2)
      {
        zvec.push_back(z_array[j]);
        x_2_vec.push_back(x_2_array[j]);
      }
      else
      {
        zvec.push_back(0);
        x_2_vec.push_back(0);
      }
      if( aper_x_n[j]<aper_x_n_min )
      {
        aper_x_n_min=aper_x_n[j];
        isigx_min=j;
        if(verbose>2) cout << "new aper_x_n_min=" << aper_x_n_min << " at isigx_min=" << isigx_min << '\n';
      }
      if( aper_y_n[j]<aper_y_n_min )
      {
        aper_y_n_min=aper_y_n[j];
        isigy_min=j;
        if(verbose>2) cout << "new aper_y_n_min=" << aper_y_n_min << " at isigy_min=" << isigy_min << '\n';
      }
    }
  }
  if(verbose)
  {
    cout << "isigx_min=" << isigx_min << " isigy_min=" << isigy_min << '\n';
    cout << "aper_x_n_min=" << aper_x_n_min << "/sqrt(m) at s=" << s_array[isigx_min] << " " << nt_StrCol_Name[isigx_min] << '\n';
    cout << "aper_y_n_min=" << aper_y_n_min << "/sqrt(m) at s=" << s_array[isigy_min] << " " << nt_StrCol_Name[isigy_min] << '\n';
//.........这里部分代码省略.........

void Plot_Bend_SR_Cones(const Ntuple& nt,const Beam &Mach,const double zmin,const double zmax,const double Scale_xy,const unsigned int verbose,const bool goto_CM_units,const double sign) // see also   ~/root_git_build/tutorials/eve/jetcone.C
// --  draw tangential lines and SR cones towards z=0      w/o beam divergence, relevant for neutral tracking  --- lines ok, cones shuld be improved
// done here in TEve. Alternative could be to create each cone as a geometry which could be saved as and reloaded, to allow for individual colors, names ..
{
  if(verbose) cout << __FILE__ << " " << __PRETTY_FUNCTION__ << " line " << __LINE__ << " Scale_xy=" << Scale_xy << endl;

  unsigned int n=nt.Noent();
  vector<string> NameCol   =nt.GetStrCol("NAME");
  vector<string> KeywordCol=nt.GetStrCol("KEYWORD");

  if(gEve==NULL) cout << " *** careful *** global gEve=" << gEve << " is not defined, some general features like setting colors may cause segmentation violation" << endl;

  // use two line sets, one for start and one for end of bend
  TEveStraightLineSet* eve_line_s = new TEveStraightLineSet("BendLines"); // can display lines together as set, but only save one by one seems  http://root.cern.ch/root/html/TEveStraightLineSet.html
  TEveStraightLineSet* eve_line_e = new TEveStraightLineSet("BendLines"); // can display lines together as set, but only save one by one seems  http://root.cern.ch/root/html/TEveStraightLineSet.html
  eve_line_s->SetLineWidth(2); eve_line_s->SetLineColor(kRed);   // draw a   red line from start of bend
  eve_line_e->SetLineWidth(2); eve_line_e->SetLineColor(kGreen); // draw a green line from end   of bend

  double length=1;
  Plot_axis_arrow("x",length,Scale_xy,verbose);
  Plot_axis_arrow("y",length,Scale_xy,verbose);
  if(zmax>20) length=100;
  Plot_axis_arrow("z",length,Scale_xy,verbose);

  // PlotCone
  TEveStraightLineSet* eve_line_m = new TEveStraightLineSet("BendLines"); // middle vector
  eve_line_m->SetLineWidth(1); eve_line_m->SetLineColor(kGray); eve_line_m->SetLineStyle(7); // https://root.cern.ch/root/html/TAttLine.html   use 7 to get dashed

  const vector<const char*> koord =GetKoordNames(nt);;
  const valarray<double> xvec=nt.GetVar(koord[0]);
  const valarray<double> yvec=nt.GetVar(koord[1]);
  const valarray<double> zvec=nt.GetVar(koord[2]);
  const valarray<double> angle=nt.GetVar("ANGLE");
  const valarray<double> theta=nt.GetVar("THETA");
  const valarray<double>   phi=nt.GetVar("PHI");
  const valarray<double>   psi=nt.GetVar("PSI");
  const valarray<double> EcritBend=nt.GetVar("EcritBend");
  const valarray<double> ngamBend=nt.GetVar("ngamBend");

  Vec3 z_unit_vec(0,0,1);
  if(sign<0) z_unit_vec=-1.*z_unit_vec;  // (0,0,-1)

  if(verbose) cout << __FILE__ << " " << __FUNCTION__  << " line " << __LINE__ << setprecision(3) << " n=" << n << " yvec.size()=" << yvec.size() << '\n';
  for(unsigned int i=1; i<n; ++i) // loop over elements
  {
    if(KeywordCol[i].find("BEND") !=string::npos) // sbend or rbend
    {
      Vec3 p0_s(xvec[i-1],yvec[i-1],zvec[i-1]); // Start of beand from previous element, ordered by increasing s
      Vec3 p0_e(  xvec[i],  yvec[i],  zvec[i]); // End   of bend

      if(verbose>1) cout << __FILE__ << " " << __FUNCTION__  << " line " << __LINE__ << " i=" << i << " zvec[i]=" << zvec[i] << '\n';
      
      if(zvec[i] > zmin && zvec[i] < zmax) // in z within range
      {
        // power
        double U0=ngamBend[i]*MeanSyn*EcritBend[i];
        double PowBend=Mach.ibeam*1.e9*U0;

        //---- cone declaration --- here in loop to get separate cones -- change color..
        TEveBoxSet* cones = new TEveBoxSet(NameCol[i].c_str());  // see  http://root.cern.ch/root/html/TEveBoxSet.html    $EDITOR $ROOTSYS/tutorials/eve/boxset_cones.C
        // cones->SetPickable(kTRUE);

        cones->UseSingleColor(); // seems needed for Transparency, do not make shape too flat, then always grey
        cones->SetMainColor(kGreen);    // see  $EDITOR ~/root_git/include/root/TEveDigitSet.h ~/root_git/src/graf3d/eve/src/TEveDigitSet.cxx         void DigitColor(Color_t ci, Char_t transparency);    https://root.cern.ch/root/html/TColor.html
        Char_t MaxTransp=99;  // seen in boxset.C  boxset_single_color   50 is half transparent, 90 is very transparent   https://root.cern.ch/root/html/TEveElement.html
        Char_t DeltaTransp=1;
        cones->SetMainTransparency(MaxTransp);
        if(PowBend>1.e1) { cones->SetMainTransparency(MaxTransp-1*DeltaTransp); }
        if(PowBend>1.e2) { cones->SetMainTransparency(MaxTransp-2*DeltaTransp); }
        if(PowBend>1.e3) { cones->SetMainTransparency(MaxTransp-3*DeltaTransp); cones->SetMainColor(kRed-2); }
        if(PowBend>1.e4) { cones->SetMainTransparency(MaxTransp-4*DeltaTransp); cones->SetMainColor(kRed-3); }
        if(PowBend>1.e5) { cones->SetMainTransparency(MaxTransp-5*DeltaTransp); cones->SetMainColor(kRed-4); }
        cones->Reset(TEveBoxSet::kBT_EllipticCone, kFALSE, 64); // EllipticCone
        //-----

        Mat3x3 WCS_s(WCS_mat3(theta[i-1],phi[i-1],psi[i-1])); // matrix, bend start
        Mat3x3 WCS_e(WCS_mat3(theta[i],    phi[i],  psi[i])); // matrix, bend end
        Vec3 dirvec_s=WCS_s*z_unit_vec; // tanget bend start
        Vec3 dirvec_e=WCS_e*z_unit_vec; // tanget bend end
        if(verbose>1) cout << left << setw(12) << NameCol[i] << setw(6) << " zvec[i]=" << setw(6) << zvec[i] << " dirvec_s " << dirvec_s.Print() << '\n' << WCS_s.Print();
        if(verbose>1) cout << left << setw(12) << NameCol[i] << setw(6) << " zvec[i]=" << setw(6) << zvec[i] << " dirvec_e " << dirvec_e.Print() << '\n' << WCS_e.Print();
        Vec3 p1_s=p0_s+dirvec_s; // start bend vector + direction unit vector
        Vec3 p1_e=p0_e+dirvec_e; // end   bend vector + direction unit vector
        double len_s=1;
        double len_e=1;
        if(fabs(p1_e.r[2])<fabs(p0_e.r[2])) // get length to reach to 0 in z
        {
          len_s=fabs(p0_s.r[2]/(p1_s.r[2]-p0_s.r[2]));
          len_e=fabs(p0_e.r[2]/(p1_e.r[2]-p0_e.r[2]));
          if(verbose>1) cout << " pointing in z to 0  len_s=" << setw(10) << len_s << " p1_s.r[2]-p0_s.r[2]= " << setw(10) << p1_s.r[2]-p0_s.r[2] << '\n';
          if(verbose>1) cout << " pointing in z to 0  len_e=" << setw(10) << len_e << " p1_e.r[2]-p0_e.r[2]= " << setw(10) << p1_e.r[2]-p0_e.r[2] << '\n';
        }
        p1_s=p0_s+len_s*dirvec_s;
        p1_e=p0_e+len_e*dirvec_e;
        if(verbose) cout << setw(5) << i << setw(15) << NameCol[i-1] << " line from " << p0_s.Print() << " to " << p1_s.Print() << " theta=" << setw(12) << theta[i-1] << '\n';
        if(verbose) cout << setw(5) << i << setw(15) << NameCol[i]   << " line from " << p0_e.Print() << " to " << p1_e.Print() << " theta=" << setw(12) << theta[i] << '\n';
        if(fabs(p1_e.r[0])> zmax)
        {
          if(verbose) cout << " x1=p1_e.r[0]=" << p1_e.r[0] << " larger than zmax=" << zmax << " skip" << '\n';
          continue;
//.........这里部分代码省略.........

	unsigned dataSize(){return data.size();};

 double fraction(const valarray<bool>& v) {
   return fraction(count(v),v.size(),0);
 }