C++ Phi函数代码示例

TString MillePedeTrees::PhiSwaps(double swapAround, const TString &tree1, const TString &tree2) const
{
  // 'true'/1 if one phi is above, one below 'swapAround'
  return 
    Parenth((Phi(tree1) += Form(">%f", swapAround)) += 
            AndL() += Phi(tree2) += Form("<%f", swapAround))
    += OrL()
    += Parenth((Phi(tree1) += Form("<%f", swapAround)) += 
               AndL() += Phi(tree2) += Form(">%f", swapAround));
}

double chisq(double z,int n)
{double s=0.,t=1.,q,h;
 int i;
 if(z<=0.) return (0.);
 if(n>3000) return Phi((exp(log(z/n)/3.)-1.+2./(9*n))/sqrt(2./(9*n)));
 h=.5*z;
 if(n&1){q=sqrt(z); t=2*exp(-.5*z-.918938533204673)/q;
    for(i=1;i<=(n-2);i+=2){ t=t*z/i; s+=t;}
    return(2*Phi(q)-1-s);
        }
 for(i=1;i<n/2;i++) { t=t*h/i; s+=t;}
 return (1.-(1+s)*exp(-h));
}

void gspmodel(void)
{
    real beta_a, mcut, vcut, vfac;
    static real *sig2 = NULL;
    real r, vmax2, sigr, sig, x, y, vr, v1, v2;
    bodyptr bp;
    vector rhat, vec1, vec2, vtmp;

    beta_a = getdparam("beta_a");
    assert(beta_a <= 1.0);
    nbody = getiparam("nbody");
    assert(nbody > 0);
    mcut = getdparam("mcut");
    assert(0.0 < mcut && mcut <= 1.0);
    vcut = getdparam("vcut");
    assert(vcut > 0.0);
    if (sig2 == NULL)
        sig2 = calc_sig2_gsp(gsp, ggsp, beta_a);
    if (btab == NULL)
	btab = (bodyptr) allocate(nbody * SizeofBody);
    vfac = rsqrt(1 - beta_a);
    for (bp = btab; bp < NthBody(btab, nbody); bp = NextBody(bp)) {
	Mass(bp) = gsp->mtot / nbody;
	r = r_mass_gsp(gsp, xrandom(0.0, mcut * gsp->mtot));
	vmax2 = -2 * rsqr(vcut) * phi_gsp(ggsp, r);
	if (vfac > 1.0)
	    vmax2 = vmax2 / rsqr(vfac);
	sigr = rsqrt(sig2_gsp(gsp, ggsp, beta_a, sig2, r));
	sig = fixsigma(sigr, rsqrt(vmax2));
	do {
	    vr = grandom(0.0, sig);
	    v1 = grandom(0.0, sig);
	    v2 = grandom(0.0, sig);
	} while (vr*vr + v1*v1 + v2*v2 > vmax2);
	picktriad(rhat, vec1, vec2);
	MULVS(Pos(bp), rhat, r);
	MULVS(Vel(bp), rhat, vr);
	MULVS(vtmp, vec1, v1 * vfac);
	ADDV(Vel(bp), Vel(bp), vtmp);
	MULVS(vtmp, vec2, v2 * vfac);
	ADDV(Vel(bp), Vel(bp), vtmp);
	Phi(bp) = phi_gsp(ggsp, r);
	Aux(bp) = Phi(bp) + 0.5 * dotvp(Vel(bp), Vel(bp));
    }
    if (getbparam("besort"))
	qsort(btab, nbody, SizeofBody, berank);
    if (getbparam("zerocm"))
	snapcenter(btab, nbody, MassField.offset);
    if (! strnull(getparam("auxvar")))
	setauxvar(btab, nbody);
}

local void hqmforces(bodyptr btab, int nbody, real M, real a, real b,
		     real tol)
{
  bodyptr bp;
  double r, mr3i, params[4], phi0, aR0, az0, abserr[3];
  static gsl_integration_workspace *wksp = NULL;
  gsl_function FPhi, F_aR, F_az;
  static double maxerr = 0.0;
  int stat[3];

  if (a == b) {					// spherical case is easy!
    for (bp = btab; bp < NthBody(btab, nbody); bp = NextBody(bp)) {
      r = absv(Pos(bp));
      Phi(bp) = - M / (a + r);
      mr3i = M * rsqr(r / (a + r)) / rqbe(r);
      MULVS(Acc(bp), Pos(bp), - mr3i);
    }
  } else {					// flattened case is harder
    if (wksp == NULL) {				// on first call, initialze
      wksp = gsl_integration_workspace_alloc(1000);
      gsl_set_error_handler_off();		// handle errors below
    }
    FPhi.function = &intPhi;
    F_aR.function = &int_aR;
    F_az.function = &int_az;
    FPhi.params = F_aR.params = F_az.params = params;
    a2(params) = rsqr(a);
    b2(params) = rsqr(b);
    for (bp = btab; bp < NthBody(btab, nbody); bp = NextBody(bp)) {
      R(params) = rsqrt(rsqr(Pos(bp)[0]) + rsqr(Pos(bp)[1]));
      z(params) = Pos(bp)[2];
      stat[0] = gsl_integration_qagiu(&FPhi, 0.0, tol, 0.0,
				      1000, wksp, &phi0, &abserr[0]);
      stat[1] = gsl_integration_qagiu(&F_aR, 0.0, tol, 0.0,
				      1000, wksp, &aR0, &abserr[1]);
      stat[2] = gsl_integration_qagiu(&F_az, 0.0, tol, 0.0,
				      1000, wksp, &az0, &abserr[2]);
      if (stat[0] || stat[1] || stat[2])	// any errors reported?
	for (int i = 0; i < 3; i++)
	  if (stat[i] != 0 && abserr[i] > maxerr) {
	    eprintf("[%s.hqmforces: warning: %s  abserr[%d] = %g]\n",
		    getprog(), gsl_strerror(stat[i]), i+1, abserr[i]);
	    maxerr = abserr[i];			// adjust reporting threshold
	  }
      Phi(bp) = - M * phi0;
      Acc(bp)[0] = - M * (Pos(bp)[0] / R(params)) * aR0;
      Acc(bp)[1] = - M * (Pos(bp)[1] / R(params)) * aR0;
      Acc(bp)[2] = - M * az0;
    }
  }
}

    void RBFInterpolation::interpolate(
        const matrix & values,
        matrix & valuesInterpolation
        )
    {
        if ( cpu && not computed )
            compute( positions, positionsInterpolation );

        assert( computed );

        if ( cpu )
        {
            matrix B, valuesLU( n_A, values.cols() ), Phi( n_B, n_A );

            if ( polynomialTerm )
            {
                Phi.resize( n_B, n_A + dimGrid + 1 );
                valuesLU.resize( n_A + dimGrid + 1, values.cols() );
            }

            valuesLU.setZero();
            valuesLU.topLeftCorner( values.rows(), values.cols() ) = values;

            B = lu.solve( valuesLU );

            evaluatePhi( positions, positionsInterpolation, Phi );

            if ( polynomialTerm )
            {
                // Include polynomial contributions in matrix Phi

                for ( int i = 0; i < Phi.rows(); i++ )
                    Phi( i, n_A ) = 1;

                Phi.topRightCorner( n_B, dimGrid ) = positionsInterpolation.block( 0, 0, n_B, dimGrid );
            }

            valuesInterpolation.noalias() = Phi * B;
        }

        if ( not cpu )
        {
            valuesInterpolation.noalias() = Hhat * values;
        }

        assert( valuesInterpolation.rows() == n_B );
        assert( values.cols() == valuesInterpolation.cols() );
    }

void polymodel(void)
{
  gsl_interp_accel *pmsplacc = gsl_interp_accel_alloc();
  bodyptr p;
  real rad, phi, vel, psi, vr, vp, a, E, J;
  vector rhat, vtmp, vper;

  for (p = btab; p < NthBody(btab, nbody); p = NextBody(p)) {
    rad = rad_m(xrandom(0.0, mtot));
    phi = gsl_spline_eval(pmspline, (double) rad, pmsplacc);
    vel = pick_v(phi);
    psi = pick_psi();
    vr = vel * rcos(psi);
    vp = vel * rsin(psi);
    Mass(p) = mtot / nbody;
    pickshell(rhat, NDIM, 1.0);
    MULVS(Pos(p), rhat, rad);
    pickshell(vtmp, NDIM, 1.0);
    a = dotvp(vtmp, rhat);
    MULVS(vper, rhat, - a);
    ADDV(vper, vper, vtmp);
    a = absv(vper);
    MULVS(vper, vper, vp / a);
    MULVS(Vel(p), rhat, vr);
    ADDV(Vel(p), Vel(p), vper);
    Phi(p) = phi;
    E = phi + 0.5 * rsqr(vel);
    J = rad * ABS(vp);
    Aux(p) = Kprime * rpow(phi1 - E, npol - 1.5) * rpow(J, 2 * mpol);
  }
  gsl_interp_accel_free(pmsplacc);
}

local void sum1force(bodyptr btab, int nbody, int i0, real eps2)
{
    bodyptr bp0, bp;
    double phi0, acc0[NDIM];
    vector pos0, dr;
    real dr2, drab, mri, mr3i;
    int j;

    bp0 = NthBody(btab, i0);
    phi0 = 0.0;
    CLRV(acc0);
    SETV(pos0, Pos(bp0));
    for (j = 0; j < nbody; j++)
	if (j != i0) {
	    bp = NthBody(btab, j);
	    DOTPSUBV(dr2, dr, Pos(bp), pos0);
	    dr2 += eps2;
	    drab = rsqrt(dr2);
	    mri = Mass(bp) / drab;
	    phi0 -= mri;
	    mr3i = mri / dr2;
	    ADDMULVS(acc0, dr, mr3i);
	}
    Phi(bp0) = phi0;
    SETV(Acc(bp0), acc0);
}

/********************************************************************
* The executable "e" runs as:  ./e <File> <Size>, where:
* File: is the text file with the network (graph) info
* Size: is the number of Monte Carlo trials.
* Seed: is the seed for the random numbers generator
********************************************************************/
int main(int argc, char *argv[]){
  int i,j,k,S,size;
  double X,V,Q,t;
  pt_net n;
  clock_t  ti;


  if(argc<4){
    printf("\n Some input data is missing! Run as:\n");
    printf("\n ./executable <File> <Size> <Seed>\n\n");
    exit(1);
  }
  if(argc>4){
    printf("\n You've entered more data than necessary! Run as:\n\n");
    printf("\n ./executable <File> <Size> <Seed>\n\n");
    exit(1);
  }
  if((size=atoi(argv[2]))<1){
    printf("\n The number of trials can not be less than 1! Run as:\n");
    printf("\n ./executable <File> <Size> <Seed>\n\n");
    exit(1);
  }

  if((seed=atoi(argv[3]))<0){
    printf("\n The seed can not be negative! Run as:\n");
    printf("\n ./executable <File> <Size> <Seed>\n\n");
    exit(1);
  }

  n=Initialize(argv[1]);

  //---------------- THE CRUDE MONTE CARLO ALGORITHM ----------------
  ti = clock();

  S=0;
  X=0.0;
  V=0.0;
  for(k=0;k<size;k++){
    Fail(n);
    X=(1-Phi(n));
    S+=X;
    V+=X*X;
  }
  Q=(double)S/size;
  V=(V/size-Q*Q)/(size-1);

  t=(double)(clock()-ti)/CLOCKS_PER_SEC;
 //------------------ PRINT OF OUTPUT AND RESULTS ------------------
 printf("\n  CRUDE MONTE CARLO:");
 printf("\n  Network: %s   Replications: %d   ExecTime=%f",argv[1],size,t);
 printf("\n\n  Q =%1.16f = %1.2e  (Unreliability)",Q,Q);
 printf("\n  V =%1.16f = %1.2e  (Variance)",V,V);
 printf("\n  SD=%1.16f = %1.2e  (Standard Deviation)",sqrt(V),sqrt(V));
 printf("\n  RE=%1.16f = %1.2f%%     (Relative Error)",sqrt(V)/Q,100*sqrt(V)/Q);
 printf("\n  ------------------------------------------------------------\n\n");
 //-----------------------------------------------------------------

  return 0;

}

EXTERN_ENV
#define global extern

#include "stdinc.h"

/*
 * HACKGRAV: evaluate grav field at a given particle.
 */

void hackgrav(bodyptr p, long ProcessId)
{
   Local[ProcessId].pskip = p;
   SETV(Local[ProcessId].pos0, Pos(p));
   Local[ProcessId].phi0 = 0.0;
   CLRV(Local[ProcessId].acc0);
   Local[ProcessId].myn2bterm = 0;
   Local[ProcessId].mynbcterm = 0;
   Local[ProcessId].skipself = FALSE;
   hackwalk(ProcessId);
   Phi(p) = Local[ProcessId].phi0;
   SETV(Acc(p), Local[ProcessId].acc0);
#ifdef QUADPOLE
   Cost(p) = Local[ProcessId].myn2bterm + NDIM * Local[ProcessId].mynbcterm;
#else
   Cost(p) = Local[ProcessId].myn2bterm + Local[ProcessId].mynbcterm;
#endif
}

//function reads in the real and imaginary data, and rotates them to output signal into R and noise into I.
void rotate_data(double * R, double * I, const int size)
{

double phi;
double sumR=0, sumI=0;
std::complex<double> m(0,-1);

for(unsigned int h=0; h<30;h++)
{
sumR=sumR+R[h];
sumI=sumI+I[h];
}
phi = atan2(sumI,sumR);

std::complex<double> Phi(phi,0);

for(unsigned int ii=0;ii<size;ii++)
{
std::complex<double> M(R[ii],I[ii]);
M=M*exp(m*Phi);
R[ii]=real(M);
I[ii]=imag(M);
}


}

inline void trajectory(double* act) {
    /*trajectory tau, t, d, p*/
    int k, j = 1;
    p_str[1].t = p_str[1].p = 0.;
    while ( p_str[j].t <= Tmax ) {
        p_str[j].tau = log( 1/(RAND() + 0.0000001 ) ) / kappa;
        p_str[j].d = randd();
        j++;
        p_str[j].t = p_str[j-1].t + p_str[j-1].tau;
        p_str[j].p = p_str[j-1].p + p_str[j-1].d;
    }


    if (j>1) p_str[j-1].tau = Tmax - p_str[j-1].t;
    else p_str[j].t = Tmax;
    /*end trajectory*/
    j--;

    if (j == 1) act[1] = act[2] = act[3] = act[4] = 0.;
    else {
        act[1] = j-1;
        act[2] = p_str[j].p;
        act[3] = act[4] = 0.;
        for (k=1; k<=j-1; k++)  {
            act[3] += p_str[k].tau * (p_str[k].p + E0*(p_str[k+1].t-p_str[k].tau/2.));
            act[4] += Phi(p_str[k].d) +
                      ( M_PI + p_str[k].tau*( p_str[k].p + E0*( p_str[k].t-p_str[k].tau/2. )*( p_str[k].t-p_str[k].tau/2. )
                                              + E0*E0 * p_str[k].tau * p_str[k].tau * p_str[k].tau / 12. )
                      ) / 2.;
        }
    }
}

void XZHydrostatic_GeoPotential<EvalT, Traits>::
evaluateFields(typename Traits::EvalData workset)
{
#ifndef ALBANY_KOKKOS_UNDER_DEVELOPMENT
  for (int cell=0; cell < workset.numCells; ++cell) {
    for (int node=0; node < numNodes; ++node) {
      for (int level=0; level < numLevels; ++level) {
        ScalarT sum =
        PhiSurf(cell,node) +
        0.5 * Pi(cell,node,level) * E.delta(level) / density(cell,node,level);
        for (int j=level+1; j < numLevels; ++j) sum += Pi(cell,node,j)     * E.delta(j)     / density(cell,node,j);

        Phi(cell,node,level) = sum;
        
        //std::cout <<"Inside GeoP, cell, node, PhiSurf(cell,node)="<<cell<<
        //", "<<node<<", "<<PhiSurf(cell,node) <<std::endl;
      }
    }
  }

  /* OG Debugging statements
  std::cout << "Printing PHI at level 0 ----------------------------------------- \n";
  //for(int level=0; level < numLevels; ++level){
  for (int node=0; node < numNodes; ++node) {
	 //std::cout << "lev= " << level <<  ", phi = " << Phi(23,0,level) <<"\n";
	std::cout << "node = " << node <<  ", phi = " << Phi(23,node,0) <<"\n";
  }
  //}*/

#else
  Kokkos::parallel_for(XZHydrostatic_GeoPotential_Policy(0,workset.numCells),*this);
  cudaCheckError();

#endif
}

void setauxvar(bodyptr btab, int nbody)
{
    bodyptr bp;
    vector jvec;
    real jtot, etot, r0, r1, r;

    if (streq(getparam("auxvar"), "mass"))
	for (bp = btab; bp < NthBody(btab, nbody); bp = NextBody(bp))
	    Aux(bp) = mass_gsp(ggsp, absv(Pos(bp)));
    else if (streq(getparam("auxvar"), "rperi"))
	for (bp = btab; bp < NthBody(btab, nbody); bp = NextBody(bp)) {
	    CROSSVP(jvec, Pos(bp), Vel(bp));
	    jtot = absv(jvec);
	    etot = 0.5 * dotvp(Vel(bp), Vel(bp)) + Phi(bp);
	    r0 = 0.0;
	    r1 = absv(Pos(bp));
	    r = 0.5 * (r0 + r1);
	    while ((r1 - r0) > TOL * r) {
		if (rsqrt(2 * (etot - phi_gsp(ggsp, r))) > jtot/r)
		    r1 = r;
		else
		    r0 = r;
		r = 0.5 * (r0 + r1);
	    }
	    Aux(bp) = r;
	}
    else
	error("%s: unknown auxvar option %s\n",
	      getargv0(), getparam("auxvar"));
}

local void sumforces(bodyptr btab, int nbody, bodyptr gtab, int ngrav,
		     real eps2)
{
  bodyptr bp, gp;
  double phi0, acc0[NDIM];
  vector pos0, dr;
  real dr2, dr2i, dr1i, mdr1i, mdr3i;

  for (bp = btab; bp < NthBody(btab, nbody); bp = NextBody(bp)) {
    phi0 = 0.0;
    CLRV(acc0);
    SETV(pos0, Pos(bp));
    for (gp = gtab; gp < NthBody(gtab, ngrav); gp = NextBody(gp)) {
      DOTPSUBV(dr2, dr, Pos(gp), pos0);
      dr2i = ((real) 1.0) / (dr2 + eps2);
      dr1i = rsqrt(dr2i);
      mdr1i = Mass(gp) * dr1i;
      mdr3i = mdr1i * dr2i;
      phi0 -= mdr1i;
      ADDMULVS(acc0, dr, mdr3i);
    }
    Phi(bp) = phi0;
    SETV(Acc(bp), acc0);
  }
}

    void RBFInterpolation::buildPhi(
        const matrix & positions,
        const matrix & positionsInterpolation
        )
    {
        n_A = positions.rows();
        n_B = positionsInterpolation.rows();
        int phiColsOld = Phi.cols();

        if ( polynomialTerm )
            Phi.conservativeResize( n_B, n_A + dimGrid + 1 );
        else
            Phi.conservativeResize( n_B, n_A );

        int nNewPoints = Phi.cols() - phiColsOld;

        if ( nNewPoints == Phi.cols() )
            nNewPoints = n_A;

        scalar r = 0;

        for ( int i = 0; i < nNewPoints; i++ )
        {
            int index = Phi.cols() - (i + 1);

            if ( polynomialTerm )
                index = Phi.cols() - 1 - dimGrid - (i + 1);

            for ( int j = 0; j < n_B; j++ )
            {
                r = ( positions.row( index ) - positionsInterpolation.row( j ) ).norm();
                Phi( j, index ) = rbfFunction->evaluate( r );
            }
        }
    }