自学教程：C++ Gt函数代码示例

/*!  Computes the first unpolarized moment    /f[  {/cal O}_{14} = /overline{q} /left[ /gamma_1   /stackrel{/displaystyle /leftrightarrow}{D}_4  +/gamma_4   /stackrel{/displaystyle /leftrightarrow}{D}_1  /right] q  /f]    with polarized projector*/void AlgNuc3pt::calc_EnergyMomentum(const ThreeMom& mom){  OpenFile();   Fprintf(fp,"The next is: Energy momentum k4 + 4k/n");  CloseFile();  for ( int n = 0; n < num_qprop; n++ ) {    for (int i(X);i<4;i++){      DIR d = DIR(i) ;      Gamma Gx(d);      Gamma Gt(T);      Derivative Der_t(T);      Derivative Der_x(d);            Nuc3ptStru Xq_xt(mom,Gx, Der_t);      Xq_xt.Calc3pt(*u_s_prop, *q_prop[n]);      Xq_xt.Calc3pt(*d_s_prop, *q_prop[n]);            Nuc3ptStru Xq_tx(mom,Gt,Der_x);      Xq_tx.Calc3pt(*u_s_prop, *q_prop[n]);      Xq_tx.Calc3pt(*d_s_prop, *q_prop[n]);            Xq_xt += Xq_tx ;      OpenFile();       Xq_xt.Print(fp) ;       CloseFile();    }  }}

void TestTimeReversibility(Integrator const& integrator) {  Length const q_initial = 1 * Metre;  Speed const v_initial = 0 * Metre / Second;  Speed const v_amplitude = 1 * Metre / Second;  Instant const t_initial;  Instant const t_final = t_initial + 100 * Second;  Time const step = 1 * Second;  std::vector<ODE::SystemState> solution;  ODE harmonic_oscillator;  harmonic_oscillator.compute_acceleration =      std::bind(ComputeHarmonicOscillatorAcceleration,                _1, _2, _3, /*evaluations=*/nullptr);  IntegrationProblem<ODE> problem;  problem.equation = harmonic_oscillator;  ODE::SystemState const initial_state = {{q_initial}, {v_initial}, t_initial};  ODE::SystemState final_state;  problem.initial_state = &initial_state;  problem.t_final = t_final;  problem.append_state = [&final_state](ODE::SystemState const& state) {    final_state = state;  };  integrator.Solve(problem, step);  problem.initial_state = &final_state;  problem.t_final = t_initial;  integrator.Solve(problem, -step);  EXPECT_EQ(t_initial, final_state.time.value);  if (integrator.time_reversible) {    EXPECT_THAT(final_state.positions[0].value,                AlmostEquals(q_initial, 0, 8));    EXPECT_THAT(final_state.velocities[0].value,                VanishesBefore(v_amplitude, 0, 16));  } else {    EXPECT_THAT(AbsoluteError(q_initial,                              final_state.positions[0].value),                Gt(1e-4 * Metre));    EXPECT_THAT(AbsoluteError(v_initial,                              final_state.velocities[0].value),                Gt(1e-4 * Metre / Second));  }}

TEST(MatchingContainers, Foo){    MockFoo foo;    EXPECT_CALL(foo, Bar(ElementsAre(1, Gt(0), _, 5)));    EXPECT_CALL(foo, Bar(UnorderedElementsAre(2, 3)));        const vector<int> a{1, 2, 3, 5};    foo.Bar(a);        const vector<int> b{3, 2};    foo.Bar(b);}

/*!  Computes the conserved vector current using the  /f[  {/cal O}_/mu = /overline{q} /gamma_/mu q  /f]    all possible mometa are inserted*/void AlgNuc3pt::calc_Cons_Vector(int Nmom, ThreeMom* mom){       Gamma Gt(T);  Nuc3ptCons VectCurr(Gt);  const int MaxNmom=50;  if(Nmom>MaxNmom)    ERR.General(cname,"calc_Cons_Vector","Nmom(%d)>MaxNmom(%d)",Nmom,MaxNmom);  Nuc3ptCons *VectCurrp[MaxNmom][4];  for(int ip(0);ip<Nmom;ip++) for (int i(X);i<4;i++){      DIR d = DIR(i) ;      Gamma G(d);      VectCurrp[ip][i] = new Nuc3ptCons(mom[ip],G) ;    }    for ( int n = 0; n < num_qprop; n++ ) {    QPropW* quark = new QPropWGaussSrc(*q_prop[n]);        VectCurr.Calc3pt(*u_s_prop,*quark);    for(int ip(0);ip<Nmom;ip++)       for (int i(X);i<4;i++)	VectCurrp[ip][i]->Calc3pt(*u_s_prop,*quark);        u_s_prop->DeleteQPropLs();        if(Nuc3pt_arg->DoConserved == 1) {      char dummy[30];      d_s_prop->RestoreQPropLs_ftom(dummy);    }        VectCurr.Calc3pt(*d_s_prop,*quark);    for(int ip(0);ip<Nmom;ip++) for (int i(X);i<4;i++)				  VectCurrp[ip][i]->Calc3pt(*d_s_prop,*quark);        d_s_prop->DeleteQPropLs();        OpenFile();     VectCurr.Print(fp) ;    for(int ip(0);ip<Nmom;ip++) for (int i(X);i<4;i++)				  VectCurrp[ip][i]->Print(fp) ;    CloseFile();        delete quark;  }  for(int ip(0);ip<Nmom;ip++)     for (int i(X);i<4;i++)      delete VectCurrp[ip][i];}

void TestTermination(    Integrator const& integrator) {  Length const q_initial = 1 * Metre;  Speed const v_initial = 0 * Metre / Second;  Instant const t_initial;  Instant const t_final = t_initial + 163 * Second;  Time const step = 42 * Second;  int const steps = static_cast<int>(std::floor((t_final - t_initial) / step));  int evaluations = 0;  std::vector<ODE::SystemState> solution;  ODE harmonic_oscillator;  harmonic_oscillator.compute_acceleration =      std::bind(ComputeHarmonicOscillatorAcceleration,                _1, _2, _3, &evaluations);  IntegrationProblem<ODE> problem;  problem.equation = harmonic_oscillator;  ODE::SystemState const initial_state = {{q_initial}, {v_initial}, t_initial};  problem.initial_state = &initial_state;  problem.t_final = t_final;  problem.append_state = [&solution](ODE::SystemState const& state) {    solution.push_back(state);  };  integrator.Solve(problem, step);  EXPECT_EQ(steps, solution.size());  EXPECT_THAT(solution.back().time.value,              AllOf(Gt(t_final - step), Le(t_final)));  switch (integrator.composition) {    case BA:    case ABA:      EXPECT_EQ(steps * integrator.evaluations, evaluations);      break;    case BAB:      EXPECT_EQ(steps * integrator.evaluations + 1, evaluations);      break;    default:      LOG(FATAL) << "Invalid composition";  }  Length q_error;  Speed v_error;  for (int i = 0; i < steps; ++i) {    Time const t = solution[i].time.value - t_initial;    EXPECT_THAT(t, AlmostEquals((i + 1) * step, 0));  }}

void TestTermination(Integrator const& integrator) {  Length const q_initial = 1 * Metre;  Speed const v_initial = 0 * Metre / Second;  Instant const t_initial;  Instant const t_final = t_initial + 1630 * Second;  Time const step = 42 * Second;  int const steps = static_cast<int>(std::floor((t_final - t_initial) / step));  int evaluations = 0;  std::vector<ODE::SystemState> solution;  ODE harmonic_oscillator;  harmonic_oscillator.compute_acceleration =      std::bind(ComputeHarmonicOscillatorAcceleration,                _1, _2, _3, &evaluations);  IntegrationProblem<ODE> problem;  problem.equation = harmonic_oscillator;  ODE::SystemState const initial_state = {{q_initial}, {v_initial}, t_initial};  problem.initial_state = &initial_state;  auto append_state = [&solution](ODE::SystemState const& state) {    solution.push_back(state);  };  auto const instance =      integrator.NewInstance(problem, std::move(append_state), step);  integrator.Solve(t_final, *instance);  EXPECT_EQ(steps, solution.size());  EXPECT_THAT(solution.back().time.value,              AllOf(Gt(t_final - step), Le(t_final)));  Length q_error;  Speed v_error;  for (int i = 0; i < steps; ++i) {    Time const t = solution[i].time.value - t_initial;    EXPECT_THAT(t, AlmostEquals((i + 1) * step, 0));  }}

bool cellGreater(Cell cell, double val) {    return cellRelOp(Gt(), cell, val);}

//.........这里部分代码省略.........	    Nuc3pt_arg->mt[nt]+=Nuc3pt_arg->t_sink; // locations of the proton sinks            Nuc3pt_arg->mt[nt]=Nuc3pt_arg->mt[nt]%(GJP.Tnodes()*GJP.TnodeSites());	    q_prop[n]->GaussSmearSinkProp(Nuc3pt_arg->mt[nt],q_prop[n]->GaussArg());          }	  	} //end smeared sink      }      //Now do the 3pt functions      // If doing coherent sinks, don't calculate the 3pt functions until all the       // forward propagators have been calculated. --MFL      int do_seq = 0;//      int t_sink;      if(Nuc3pt_arg->calc_seqQ != MULT_SEQ 	&& Nuc3pt_arg->calc_seqQ != WRITE_MULT_SEQ 	&& Nuc3pt_arg->calc_seqQ != READ_MULT_SEQ) {	do_seq = 1;        t_sink = ts + Nuc3pt_arg->t_sink;      }      // once all the forward propagators have been calculated,      // do the sequential propagaotrs. --MFL      else if (i_source == num_qprop-1) {	do_seq = 1;	t_sink = Nuc3pt_arg->t_source + Nuc3pt_arg->t_sink;      }      if(Nuc3pt_arg->DoUnPolarized && do_seq)	{	  // for conserved vector currents	  Gamma Gt(T);	  Nuc3ptCons VectCurr(Gt);	  Nuc3ptCons *VectCurrp[4];	  if(Nuc3pt_arg->DoConserved) {	  for (int i(X);i<4;i++){	    DIR d = DIR(i) ;	    Gamma G(d);	    VectCurrp[i] = new Nuc3ptCons(Nmom,G) ;	  }	  }	  OpenFile();	  Fprintf(fp,"UnPolarized Zero mom./n");	  CloseFile();  	  //first do the zero momentum un-polarized stuff	  //up-quark	  GetTheSeqPropagator(t_sink,qmass,			      PROT_U_SEQ,ZeroMom,PPAR);	  //down-quark	  GetTheSeqPropagator(t_sink,qmass,			      PROT_D_SEQ,ZeroMom,PPAR);	  calc_Scalar();         //needed for the sigma term	  calc_Vector();	 //Vector current	  calc_X_q_b();          //<x>_q (b) does not need non-zero momentum	  for(int i(0) ; i<Nmom ; i++) calc_Vector(sink_mom[i]);	  //conserved current	  if(Nuc3pt_arg->DoConserved) {	    if(GJP.Snodes()==2) u_s_prop->SwapQPropLs();

    virtual void init(Voxel& voxel)    {        half_odf_size = voxel.ti.half_vertices_count;        unsigned int b_count = voxel.bvalues.size();        icosa_data.resize(half_odf_size*3);        for (unsigned int index = 0; index < half_odf_size; ++index)            std::copy(voxel.ti.vertices[index].begin(),voxel.ti.vertices[index].end(),icosa_data.begin()+index*3);        float interop_angle = voxel.param[0]/180.0*M_PI;        float smoothing_angle = voxel.param[1]/180.0*M_PI;        Ht.resize(half_odf_size*b_count); // n * m        // H=phi(acos(QtV))        for (unsigned int n = 0,index = 0; n < half_odf_size; ++n)            for (unsigned int m = 0; m < b_count; ++m,++index)            {                float value = std::abs(                                  voxel.bvectors[m]*image::vector<3,float>(voxel.ti.vertices[n]));                Ht[index] = spherical_guassian(value,interop_angle);            }        iHtH.resize(half_odf_size*half_odf_size);        iHtH_pivot.resize(half_odf_size);        image::matrix::square(Ht.begin(),iHtH.begin(),image::dyndim(half_odf_size,b_count));        image::matrix::lu_decomposition(iHtH.begin(),iHtH_pivot.begin(),image::dyndim(half_odf_size,half_odf_size));        // vector of angles        std::vector<float> C(3*k); // 3 by k matrix        for (unsigned int index = 0; index < k; ++index)        {            C[index] = std::cos(2.0*M_PI*((float)index+1)/((float)k));            C[k+index] = std::sin(2.0*M_PI*((float)index+1)/((float)k));            C[k+k+index] = 0.0;        }        // RC        std::vector<float> G(half_odf_size*half_odf_size);        std::vector<float> icosa_data_r(half_odf_size*3);        for (unsigned int gi = 0; gi < half_odf_size; ++gi)        {            image::vector<3,float> u(voxel.ti.vertices[gi]);            float r[9];// a 3-by-3 matrix            rotation_matrix(r,u.begin());            std::vector<float> Gt(half_odf_size*k); // a half_odf_size-by-k matrix            // 	Gt = icosa_data*r*C;            image::matrix::product(icosa_data.begin(),r,icosa_data_r.begin(),image::dyndim(half_odf_size,3),image::dim<3,3>());            image::matrix::product(icosa_data_r.begin(),C.begin(),Gt.begin(),image::dyndim(half_odf_size,3),image::dyndim(3,k));            for (unsigned int i = 0; i < Gt.size(); ++i)                Gt[i] = spherical_guassian(std::abs(Gt[i]),interop_angle);            unsigned int posgi = gi*half_odf_size;            for (unsigned int i = 0,posi = 0; i < half_odf_size; ++i,posi+=k)                G[posgi+i] = std::accumulate(Gt.begin()+posi,Gt.begin()+posi+k,0.0);        }        // add smoothing to G        std::vector<float> S(half_odf_size*half_odf_size);        for (unsigned int i = 0; i < half_odf_size; ++i)        {            float sum = 0.0;            for (unsigned int j = 0,index = i*half_odf_size; j < half_odf_size; ++j,++index)                sum +=                    S[index] = spherical_guassian(std::abs(voxel.ti.vertices_cos(i,j)),smoothing_angle);            for (unsigned int j = 0,index = i*half_odf_size; j < half_odf_size; ++j,++index)                S[index] /= sum;        }        sG.resize(half_odf_size*half_odf_size);        //sG = S*G;        image::matrix::product(S.begin(),G.begin(),sG.begin(),image::dyndim(half_odf_size,half_odf_size),image::dyndim(half_odf_size,half_odf_size));    }

void TestConvergence(Integrator const& integrator,                     Time const& beginning_of_convergence) {  Length const q_initial = 1 * Metre;  Speed const v_initial = 0 * Metre / Second;  Speed const v_amplitude = 1 * Metre / Second;  AngularFrequency const ω = 1 * Radian / Second;  Instant const t_initial;  Instant const t_final = t_initial + 100 * Second;  Time step = beginning_of_convergence;  int const step_sizes = 50;  double const step_reduction = 1.1;  std::vector<double> log_step_sizes;  log_step_sizes.reserve(step_sizes);  std::vector<double> log_q_errors;  log_step_sizes.reserve(step_sizes);  std::vector<double> log_p_errors;  log_step_sizes.reserve(step_sizes);  std::vector<ODE::SystemState> solution;  ODE harmonic_oscillator;  harmonic_oscillator.compute_acceleration =      std::bind(ComputeHarmonicOscillatorAcceleration,                _1, _2, _3, /*evaluations=*/nullptr);  IntegrationProblem<ODE> problem;  problem.equation = harmonic_oscillator;  ODE::SystemState const initial_state = {{q_initial}, {v_initial}, t_initial};  problem.initial_state = &initial_state;  ODE::SystemState final_state;  auto const append_state = [&final_state](ODE::SystemState const& state) {    final_state = state;  };  for (int i = 0; i < step_sizes; ++i, step /= step_reduction) {    auto const instance = integrator.NewInstance(problem, append_state, step);    integrator.Solve(t_final, *instance);    Time const t = final_state.time.value - t_initial;    Length const& q = final_state.positions[0].value;    Speed const& v = final_state.velocities[0].value;    double const log_q_error = std::log10(        AbsoluteError(q / q_initial, Cos(ω * t)));    double const log_p_error = std::log10(        AbsoluteError(v / v_amplitude, -Sin(ω * t)));    if (log_q_error <= -13 || log_p_error <= -13) {      // If we keep going the effects of finite precision will drown out      // convergence.      break;    }    log_step_sizes.push_back(std::log10(step / Second));    log_q_errors.push_back(log_q_error);    log_p_errors.push_back(log_p_error);  }  double const q_convergence_order = Slope(log_step_sizes, log_q_errors);  double const q_correlation =      PearsonProductMomentCorrelationCoefficient(log_step_sizes, log_q_errors);  LOG(INFO) << "Convergence order in q : " << q_convergence_order;  LOG(INFO) << "Correlation            : " << q_correlation;#if !defined(_DEBUG)  EXPECT_THAT(RelativeError(integrator.order, q_convergence_order),              Lt(0.05));  EXPECT_THAT(q_correlation, AllOf(Gt(0.99), Lt(1.01)));#endif  double const v_convergence_order = Slope(log_step_sizes, log_p_errors);  double const v_correlation =      PearsonProductMomentCorrelationCoefficient(log_step_sizes, log_p_errors);  LOG(INFO) << "Convergence order in p : " << v_convergence_order;  LOG(INFO) << "Correlation            : " << v_correlation;#if !defined(_DEBUG)  // SPRKs with odd convergence order have a higher convergence order in p.  EXPECT_THAT(      RelativeError(integrator.order + (integrator.order % 2),                    v_convergence_order),      Lt(0.03));  EXPECT_THAT(v_correlation, AllOf(Gt(0.99), Lt(1.01)));#endif}

Gt operator>(const FieldType& fld, const FieldType& o2) {    return Gt(fld, o2);}

 static bool Eq(Value lhs, Value rhs) { return !Lt(lhs, rhs) && !Gt(lhs, rhs); }

void Parser::executeAction(int production) try {    if (d_token__ != _UNDETERMINED_)        pushToken__(d_token__); // save an already available token    // $insert defaultactionreturn    // save default non-nested block $$    if (int size = s_productionInfo[production].d_size)        d_val__ = d_vsp__[1 - size];    switch (production) {        // $insert actioncases        case 1:#line 43 "parser.yy"        {            d_val__.get<Tag__::basic>() = d_vsp__[0].data<Tag__::basic>();            res = d_val__.get<Tag__::basic>();        } break;        case 2:#line 51 "parser.yy"        {            d_val__.get<Tag__::basic>() = add(d_vsp__[-2].data<Tag__::basic>(),                                              d_vsp__[0].data<Tag__::basic>());        } break;        case 3:#line 54 "parser.yy"        {            d_val__.get<Tag__::basic>() = sub(d_vsp__[-2].data<Tag__::basic>(),                                              d_vsp__[0].data<Tag__::basic>());        } break;        case 4:#line 57 "parser.yy"        {            d_val__.get<Tag__::basic>() = mul(d_vsp__[-2].data<Tag__::basic>(),                                              d_vsp__[0].data<Tag__::basic>());        } break;        case 5:#line 60 "parser.yy"        {            d_val__.get<Tag__::basic>() = div(d_vsp__[-2].data<Tag__::basic>(),                                              d_vsp__[0].data<Tag__::basic>());        } break;        case 6:#line 63 "parser.yy"        {            auto tup = parse_implicit_mul(d_vsp__[-2].data<Tag__::string>());            if (neq(*std::get<1>(tup), *one)) {                d_val__.get<Tag__::basic>() = mul(                    std::get<0>(tup),                    pow(std::get<1>(tup), d_vsp__[0].data<Tag__::basic>()));            } else {                d_val__.get<Tag__::basic>()                    = pow(std::get<0>(tup), d_vsp__[0].data<Tag__::basic>());            }        } break;        case 7:#line 73 "parser.yy"        {            d_val__.get<Tag__::basic>() = pow(d_vsp__[-2].data<Tag__::basic>(),                                              d_vsp__[0].data<Tag__::basic>());        } break;        case 8:#line 76 "parser.yy"        {            d_val__.get<Tag__::basic>() = rcp_static_cast<const Basic>(                Lt(d_vsp__[-2].data<Tag__::basic>(),                   d_vsp__[0].data<Tag__::basic>()));        } break;        case 9:#line 79 "parser.yy"        {            d_val__.get<Tag__::basic>() = rcp_static_cast<const Basic>(                Gt(d_vsp__[-2].data<Tag__::basic>(),                   d_vsp__[0].data<Tag__::basic>()));        } break;        case 10:#line 82 "parser.yy"        {            d_val__.get<Tag__::basic>() = rcp_static_cast<const Basic>(                Le(d_vsp__[-2].data<Tag__::basic>(),                   d_vsp__[0].data<Tag__::basic>()));        } break;        case 11:#line 85 "parser.yy"        {            d_val__.get<Tag__::basic>() = rcp_static_cast<const Basic>(                Ge(d_vsp__[-2].data<Tag__::basic>(),                   d_vsp__[0].data<Tag__::basic>()));        } break;//.........这里部分代码省略.........

TEST_P(SimpleHarmonicMotionTest, Convergence) {  parameters_.initial.positions.emplace_back(SIUnit<Length>());  parameters_.initial.momenta.emplace_back(Speed());  parameters_.initial.time = Time();#if defined(_DEBUG)  parameters_.tmax = 1 * SIUnit<Time>();#else  parameters_.tmax = 100 * SIUnit<Time>();#endif  parameters_.sampling_period = 0;  parameters_.Δt = GetParam().beginning_of_convergence;  int const step_sizes = 50;  double const step_reduction = 1.1;  std::vector<double> log_step_sizes;  log_step_sizes.reserve(step_sizes);  std::vector<double> log_q_errors;  log_step_sizes.reserve(step_sizes);  std::vector<double> log_p_errors;  log_step_sizes.reserve(step_sizes);  for (int i = 0; i < step_sizes; ++i, parameters_.Δt /= step_reduction) {    integrator_->SolveTrivialKineticEnergyIncrement<Length>(        &ComputeHarmonicOscillatorAcceleration,        parameters_,        &solution_);    double const log_q_error = std::log10(        std::abs(solution_[0].positions[0].value / SIUnit<Length>() -                 Cos(solution_[0].time.value *                     SIUnit<AngularFrequency>())));    double const log_p_error = std::log10(        std::abs(solution_[0].momenta[0].value / SIUnit<Speed>() +                 Sin(solution_[0].time.value *                     SIUnit<AngularFrequency>())));    if (log_q_error <= -13 || log_p_error <= -13) {      // If we keep going the effects of finite precision will drown out      // convergence.      break;    }    log_step_sizes.push_back(std::log10(parameters_.Δt / SIUnit<Time>()));    log_q_errors.push_back(log_q_error);    log_p_errors.push_back(log_p_error);  }  double const q_convergence_order = Slope(log_step_sizes, log_q_errors);  double const q_correlation =      PearsonProductMomentCorrelationCoefficient(log_step_sizes, log_q_errors);  LOG(INFO) << GetParam();  LOG(INFO) << "Convergence order in q : " << q_convergence_order;  LOG(INFO) << "Correlation            : " << q_correlation;#if 0  LOG(INFO) << "Convergence data for q :/n" <<      BidimensionalDatasetMathematicaInput(log_step_sizes, log_q_errors);#endif#if !defined(_DEBUG)  EXPECT_THAT(RelativeError(GetParam().convergence_order, q_convergence_order),              Lt(0.02));  EXPECT_THAT(q_correlation, AllOf(Gt(0.99), Lt(1.01)));#endif  double const v_convergence_order = Slope(log_step_sizes, log_p_errors);  double const v_correlation =      PearsonProductMomentCorrelationCoefficient(log_step_sizes, log_p_errors);  LOG(INFO) << "Convergence order in p : " << v_convergence_order;  LOG(INFO) << "Correlation            : " << v_correlation;#if 0  LOG(INFO) << "Convergence data for p :/n" <<      BidimensionalDatasetMathematicaInput(log_step_sizes, log_q_errors);#endif#if !defined(_DEBUG)  // SPRKs with odd convergence order have a higher convergence order in p.  EXPECT_THAT(     RelativeError(((GetParam().convergence_order + 1) / 2) * 2,                   v_convergence_order),     Lt(0.02));  EXPECT_THAT(v_correlation, AllOf(Gt(0.99), Lt(1.01)));#endif}

TEST_F(Man
C++ Guard函数代码示例
C++ Grid函数代码示例