C++ valarray类代码示例

inline void f_harmonic(valarray<double> &r, valarray<double> &k, const int dim)

{
    //constants
    const double m = 1.0;
    const double omega = 1.0;

    //get number of particles
    const int N_p = r.size() / (2 * dim);


    //reset k vector
    for (unsigned int i = 0; i < k.size(); i++)
    {
        k[i] = 0.0;
    }



    for (int d = 0; d < dim; d++)
    {

        for (int i = 0; i < N_p; i++)

        {
            k[(2 * i + 1)*dim + d] = (-1.) * m * omega * omega * r[i * 2 * dim + d];
            k[i * 2 * dim + d] = r[(2 * i + 1) * dim + d] / m;
        }

    }


}

T norm(const valarray<T>& v, const unsigned p) {
   T res = 0, scale = 0., absvi, ssq = 1., tmp;

   if ( v.size() < 1 ) {
      return static_cast<T>(0.);
   }

   else if ( v.size() == 1 ) {
      return abs(v[0]);
   }

   else {
      switch (p) {
         case 0:
            //for (size_t i = 0; i < v.size(); ++i) {
            //   if ( res < abs(v[i]) )
            //      res = abs(v[i]);
            //}
            //return res;
            for (auto it = begin(v); it != end(v); ++it) {
               if ( res < abs(*it) )
                  res = abs(*it);
            }
            return res;
            //return MAX(abs(v.min()), abs(v.max()));

         case 1:
            for (size_t i = 0; i < v.size(); ++i) {
               res += abs(v[i]);
            }
            return res;
            //for ( auto it = begin(v); it != end(v); ++it)
            //   res += abs(*it);
            //return res;

         case 2:
            for (size_t i = 0; i < v.size(); ++i) {
               if ( v[i] != 0. ) {
                  absvi = abs(v[i]);
                  if ( scale < absvi ) {
                     tmp = scale/absvi;
                     ssq = 1. + ssq*tmp*tmp;
                     scale = absvi;
                  }
                  else {
                     tmp = absvi/scale;
                     ssq += tmp*tmp;
                  }
               }
            }
            res = scale*sqrt(ssq);
            return res;

         default:
            cerr << "error: utils:norm(vector): unsupported p-norm: p = " << p 
                 << endl;
            exit(-1);
      }  
   }
}

	void SetUpCL(valarray<cl_float> & data, valarray<cl_float> & data_err, valarray<cl_float> & model, valarray<cl_float> & output)
	{
		unsigned int test_size = data.size();

		assert(data_err.size() == test_size);
		assert(model.size() == test_size);
		assert(output.size() == test_size);

		// Init OpenCL and the routine
		cl = new COpenCL(CL_DEVICE_TYPE_GPU);
		zero = new CRoutine_Zero(cl->GetDevice(), cl->GetContext(), cl->GetQueue());
		zero->SetSourcePath(LIBOI_KERNEL_PATH);
		zero->Init();
		square = new CRoutine_Square(cl->GetDevice(), cl->GetContext(), cl->GetQueue());
		square->SetSourcePath(LIBOI_KERNEL_PATH);
		square->Init();
		r = new CRoutine_Chi(cl->GetDevice(), cl->GetContext(), cl->GetQueue(), zero, square);
		r->SetSourcePath(LIBOI_KERNEL_PATH);
		r->Init(test_size);

		// Make OpenCL buffers for the data, data_err, model, and output.
		data_cl = clCreateBuffer(cl->GetContext(), CL_MEM_READ_WRITE, sizeof(cl_float) * test_size, NULL, NULL);
		data_err_cl = clCreateBuffer(cl->GetContext(), CL_MEM_READ_WRITE, sizeof(cl_float) * test_size, NULL, NULL);
		model_cl = clCreateBuffer(cl->GetContext(), CL_MEM_READ_WRITE, sizeof(cl_float) * test_size, NULL, NULL);
		output_cl = clCreateBuffer(cl->GetContext(), CL_MEM_READ_WRITE, sizeof(cl_float) * test_size, NULL, NULL);

		// Fill the input buffer
		int err = CL_SUCCESS;
		err = clEnqueueWriteBuffer(cl->GetQueue(), data_cl, CL_TRUE, 0, sizeof(cl_float) * test_size, &data[0], 0, NULL, NULL);
		err = clEnqueueWriteBuffer(cl->GetQueue(), data_err_cl, CL_TRUE, 0, sizeof(cl_float) * test_size, &data_err[0], 0, NULL, NULL);
		err = clEnqueueWriteBuffer(cl->GetQueue(), model_cl, CL_TRUE, 0, sizeof(cl_float) * test_size, &model[0], 0, NULL, NULL);
	}

void apply_op(valarray<T>& lhs, valarray<T> const& x, valarray<T> const& y, Op op = Op{}) {
	uint64_t size = std::min(x.size(), y.size());
	size = std::min(size, lhs.size()); // probably not needed
	for (uint64_t k = 0; k < size; k += 1) {
		lhs[k] = op(x[k], y[k]);
	}
}

	void MakeChiZeroBuffers(valarray<cl_float> & data, valarray<cl_float> & data_err, valarray<cl_float> & model, valarray<cl_float> & output, unsigned int test_size)
	{
		unsigned int n = 2*test_size;
		// Create buffers
		data.resize(n);
		data_err.resize(n);
		model.resize(n);
		output.resize(n);

		valarray<cl_float2> temp = CModel::GenerateUVSpiral_CL(test_size);

		// Set data = model to produce a zero chi result.
		for(int i = 0; i < test_size; i++)
		{
			data[i] = temp[i].s0;
			data[test_size + i] = temp[i].s1;

			model[i] = temp[i].s0;
			model[test_size + i] = temp[i].s1;

			// 1% error on amplitudes, 10% error on phases
			data_err[i] = amp_err * data[i];
			data_err[test_size + i] = phi_err * data[i];
		}
	}

static double compareAgainst(
            const valarray<double> &z,     // column vector
            const valarray<double> &h,     // column vector
            const valarray<double> &Q,     // symmetric matrix Q(x)
            const valarray<double> &invQ   // inverse of Q
            )
{
    const int n = z.size();

    assert( Q.size() == n * n );
    assert( h.size() == n );

    //assert( hasPosDet( dmatQ ) );

    valarray<double> Q2PI( n * n );
    
    Q2PI = Q * 2.0 * PI;

    // Compute det(Q2PI) = b * 2^c .
    double b;
    long c;
    det( Q2PI, n, &b, &c );

    // Compute log(det(Q2PI)).
    double dTerm1 = log( b ) + c * log( 2.0 );

    valarray<double>  residual = z - h;

    valarray<double>  term2 = multiply( multiply( transpose( residual, n ), n, invQ, n ), n, residual, 1 );
    double        dTerm2 = term2[0];

    return (dTerm1 + dTerm2) / 2.0;
}

// return dy/dx for each particular equation
valarray<double> dydx( const double x, const valarray<double>& y ) {

    if ( y.size() != 4 ) throw;

    valarray<double> tmpdydx(0.0,y.size());
    try {
        // slices would be easier, but oh well...
        double z = y[0];
        double p = y[1];
        double zbar = y[2];
        double pbar = y[3];
        
        double zdot = ( 1 + z*zbar ) * ( pbar + (pbar*zbar)*z );
        double pdot = 2 * zbar * ( p*pbar + (p*z)*(pbar*zbar) ) +
            ( 1 + z*zbar ) * ( (pbar*zbar)*z );
        double zbardot = ( 1 + z*zbar ) * ( p + (p*z)*zbar );
        double pbardot = 2 * z * ( p*pbar + (p*z)*(pbar*zbar) ) +
            ( 1 + z*zbar ) * ( (p*z)*zbar );

        // recombine into tmpdydx... again, slices would be easier
        tmpdydx[0] = zdot;
        tmpdydx[1] = pdot;
        tmpdydx[2] = zbardot;
        tmpdydx[3] = pbardot;

    }
    catch(out_of_range) {
        cerr << "oops" << endl;
        exit(1);
    }

    return tmpdydx;

}

inline void rk4(valarray<double> &r, double h, double &t, const int dim, void (*f)(valarray<double> &, valarray<double> &, const int))

{
    //define k, r_step, and temp arrays to hold eom evaluations
    valarray<double> k(r.size());
    valarray<double> r_step(r.size());
    valarray<double> temp(r.size());

    const double half_h = h / 2.;

    //1st rk4 step
    f(r, k, dim);
    r_step = h * (1. / 6.) * k;
    temp = r + half_h * k;

    //2nd
    f(temp, k, dim);
    r_step += h * (1. / 3.) * k;
    temp = r + half_h * k;

    //3rd
    f(temp, k, dim);
    r_step += h * (1. / 3.) * k;
    temp = r + h * k;

    //4th
    f(temp, k, dim);

    //advance r in time
    r += r_step + h * (1. / 6.) * k;


    //advance time
    t += h;
}

BroadEstepper1::BroadEstepper1(CanaryOptions& opts,
                               valarray<double> vals_in,
                               CanaryPrior & prior_in,
                               vector<int> & clusters_in)
  : _opts(opts)
  {
  N = vals_in.size();
  G = clusters_in.size();
  _vals.resize(vals_in.size()); _vals = vals_in;
  _prior = prior_in;
  _cvec.resize(G); _cvec = clusters_in;
  _prop.resize(G,1.0/G);
  _mean.resize(G);
  _var.resize(G);

  for (unsigned int k=0; k<G; k++)
    {
    int cpos = _cvec[k];
//  The algorithm converges to a cycle dependent on initial values.  If it
//  converged to a point - as it should - it would make sense to uncomment.
//    _prop[k] = _prior.prop()[cpos];
    _mean[k] = _prior.mean()[cpos];
    _var[k] = _prior.var()[cpos];
    }

  _prob.ReSize(N,G);
  _prob = 0.0;

  update_prob();
  }

double FullCovariance::calcLogdetFromCholesky( const valarray<double>& chol ) const
{
  //------------------------------------------------------------
  // Preliminaries.
  //------------------------------------------------------------

  int nChol = static_cast<int>(sqrt(static_cast<double>(chol.size())));


  //------------------------------------------------------------
  // Validate the current state (debug version only).
  //------------------------------------------------------------

  assert( chol.size() == nChol * nChol );


  //------------------------------------------------------------
  // Compute the logdet.
  //------------------------------------------------------------

  // Sum the logs of the diagonals.
  double logsum = 0.0;
  for ( int i = 0; i < nChol; i++ )
  {
    logsum += log( choleskyCached[i + i * nChol] );
  }

  // Return
  //                      T                       2         ---
  //     logdet( Chol Chol  )  =  log[ det( Chol )  ]  =  2 >    log( chol   )  .
  //                                                        ---           ii
  return 2.0 * logsum;
}

vector<double> convert(const valarray<double>& v1)
{
  vector<double> v2(v1.size());
  for(int i=0;i<v1.size();i++)
    v2[i] = v1[i];
  return v2;
}

double GradientProjection::computeSteepestDescentVector(
        valarray<double> const &b,
        valarray<double> const &x,
        valarray<double> &g) const {
    // find steepest descent direction
    //  g = 2 ( b - A x )
    //    where: A = denseQ + sparseQ
    //  g = 2 ( b - denseQ x) - 2 sparseQ x
    //
    //  except the 2s don't matter because we compute 
    //  the optimal stepsize anyway
    COLA_ASSERT(x.size()==b.size() && b.size()==g.size());
    g = b;
    for (unsigned i=0; i<denseSize; i++) {
        for (unsigned j=0; j<denseSize; j++) {
            g[i] -= (*denseQ)[i*denseSize+j]*x[j];
        }
    }
    // sparse part:
    if(sparseQ) {
        valarray<double> r(x.size());
        sparseQ->rightMultiply(x,r);
        g-=r;
    }
    return computeStepSize(g,g);
}

int drawCairo(const string& fname,
        const valarray<double>& Xin, const valarray<double>& Yin, 
        const Hull& hull) {
#ifdef CAIRO_HAS_SVG_SURFACE
    unsigned n=Xin.size();
    assert(Yin.size()==n);

    // normalise coords to range 0-1
    valarray<double> X=Xin, Y=Yin;
    X-=X.min();
    Y-=Y.min();
    X/=X.max();
    Y/=Y.max();

    Cairo::RefPtr<Cairo::SvgSurface> surface =
        Cairo::SvgSurface::create(fname, width+2*border, height+2*border);

    Cairo::RefPtr<Cairo::Context> cr = Cairo::Context::create(surface);

    cr->save(); // save the state of the context
    cr->set_source_rgba(0.0, 0.0, 0.0, 0.7);
    // draw a circle at each coordinate
    for(unsigned i=0;i<n;i++) {
        dot(cr,xcoord(X[i]),ycoord(Y[i]));
    }

    cr->set_source_rgba(0.0, 0.0, 0.0, 0.3);
    cr->move_to(xcoord(X[hull[0]]),ycoord(Y[hull[0]]));
    for(unsigned i=1;i<hull.size();i++) {
        cr->line_to(xcoord(X[hull[i]]),ycoord(Y[hull[i]]));
    }
    cr->line_to(xcoord(X[hull[0]]),ycoord(Y[hull[0]]));
    cr->stroke();
    cr->set_source_rgba(0.0, 0.0, 0.0, 1.);
    for(vector<unsigned>::const_iterator i=hull.begin();i!=hull.end();++i) {
        unsigned j=*i;
        stringstream ss;
        ss<<j;
        printf("p[%d]=(%f,%f)\n",j,X[j],Y[j]);
        cr->move_to(xcoord(X[j]),ycoord(Y[j]));
        cr->show_text(ss.str());
        cr->stroke();
    }
    cr->restore();

    cr->show_page();

    cout << "Wrote SVG file \"" << fname << "\"" << endl;
    return 0;

#else

    cout << "You must compile cairo with SVG support for this example to work."
        << endl;
    return 1;

#endif

}

  /// FIXME - this could be faster by a factor of 2 - we are sorting TWICE!
  pair<double,double> confidence_interval(const valarray<double>& values,double P) 
  {
    vector<double> values2(values.size());
    for(int i=0;i<values.size();i++)
      values2[i] = values[i];

    return confidence_interval(values2,P);
  }

valarray<int> iterate_zs(valarray<cmplx>& zs, const complex<float>& c,
                         size_t max_iters) {
    valarray<int> ns(zs.size());
#pragma omp parallel for schedule(runtime)
    for (size_t i = 0; i < zs.size(); i++)
        ns[i] = iterate_z(zs[i], c, max_iters);
    return ns;
}