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

/*------------------------------------------------------------------*/void AzOptOnTree::updateTreeWeights(AzRgfTreeEnsemble *ens) const{  int dtree_num = ens->size();   int tx;   for (tx = 0; tx < dtree_num; ++tx) {    ens->tree_u(tx)->resetWeights();   }  const double *weight = weights()->point();   double const_val = constant();   ens->set_constant(const_val);   int num = tree_feat->featNum();   int fx;     for (fx = 0; fx < num; ++fx) {    if (weight[fx] != 0) {      const AzTrTreeFeatInfo *fp = tree_feat->featInfo(fx);       ens->tree_u(fp->tx)->setWeight(fp->nx, weight[fx]);     }  }}

// This is inherently inefficient as this function will be called n times (up to 4 times)// for every overlapping pixel. I'm doing this as to reduce memory use. Since it will// only be done once, it may not matter too much.float IdMap::calculateWeights(int x, int y, FrameStack* frame, 			      std::vector<FrameStack*> stacks, float pixelFactor){  std::vector<float> weights(stacks.size());  std::vector<int> border_distances(stacks.size());  int maxBorderDistance = 0;  for(uint i=0; i < stacks.size(); ++i){    border_distances[i] = stacks[i]->nearestBorderGlobal(x, y);    maxBorderDistance = border_distances[i] > maxBorderDistance ? border_distances[i] : maxBorderDistance;  }  // override pixelFactor with MaxBorderDistance  pixelFactor = maxBorderDistance ? (float)maxBorderDistance : pixelFactor;  int offset = -1;  float w_sum = 0;  for(uint i=0; i < stacks.size(); ++i){    if(stacks[i] == frame)      offset = i;    int border_distance = stacks[i]->nearestBorderGlobal(x, y);    weights[i] = border_distance > (int)pixelFactor ? 1.0 : (float)(1 + border_distance) / (1.00 + pixelFactor);    weights[i] = weights[i] < 0 ? 0 : weights[i];    w_sum += weights[i];  }  // float w_sum = 0;  // for(uint i=0; i < weights.size(); ++i)  //   w_sum += weights[i];  if(!w_sum){    std::cerr << "idMap::calculateWeights: w_sum is 0" << std::endl;    std::cerr << "/tx,y: " << x << "," << y << " offset : " << offset	      << "/tstacks.size() " << stacks.size() << std::endl;    exit(1);  }  for(uint i=0; i < weights.size(); ++i)    weights[i] /= w_sum;  if(offset < 0){    std::cerr << "idMap::calculateWeights: offset not defined returning 0 weighting" << std::endl;    return(0);  }  return(weights[offset]);}

	//! Static function to calculate the weights.	static std::vector<double> calculate_weights(const std::vector< std::vector<double> >& acceptance_probabilities_optimization)	{	  // Abbreviation for the number of optimization steps	  unsigned int optimization_steps = acceptance_probabilities_optimization.size();	  // Calculate the weights	  std::vector<double> weights(optimization_steps, 0.0);	  double weights_sum = 0.0;	  for (unsigned int r = 0; r < optimization_steps; ++r)	  {	    double sigma_squared = 0.0;	    for (unsigned int i = 0; i < acceptance_probabilities_optimization[r].size(); ++i)	      sigma_squared += 1.0/pow(acceptance_probabilities_optimization[r][i], 2);	  	    weights[r] = 1.0 / sqrt(sigma_squared);	    weights_sum += weights[r];	  }	  	  // Normalize the weights	  for (unsigned int r = 0; r < optimization_steps; ++r) weights[r] /= weights_sum;	  return weights;	}

std::vector<double> MultiSequence::ComputePositionBasedSequenceWeights() const{    std::vector<double> weights(sequences.size(), 0.0);    std::vector<int> counts(256);    for (int i = 1; i <= sequences[0]->GetLength(); i++)    {        int diversity = 0;        std::fill(counts.begin(), counts.end(), 0);                for (size_t j = 0; j < sequences.size(); j++)        {            if (counts[BYTE(sequences[j]->GetData()[i])] == 0) diversity++;            ++(counts[BYTE(sequences[j]->GetData()[i])]);        }                for (size_t j = 0; j < sequences.size(); j++)            weights[j] += 1.0 / (diversity * counts[BYTE(sequences[j]->GetData()[i])]);    }    weights /= Sum(weights);    return weights;}

void LandmarkSelectionBase::selectLandmarks( MultiReferenceBase* myframes ) {  // Select landmarks  myframes->clearFrames(); select( myframes );  plumed_assert( myframes->getNumberOfReferenceFrames()==nlandmarks );  // Now calculate voronoi weights  if( !novoronoi ) {    unsigned rank=action->comm.Get_rank();    unsigned size=action->comm.Get_size();    std::vector<double> weights( nlandmarks, 0.0 );    for(unsigned i=rank; i<action->data.size(); i+=size) {      unsigned closest=0;      double mindist=distance( action->getPbc(), action->getArguments(), action->data[i], myframes->getFrame(0), false );      for(unsigned j=1; j<nlandmarks; ++j) {        double dist=distance( action->getPbc(), action->getArguments(), action->data[i], myframes->getFrame(j), false );        if( dist<mindist ) { mindist=dist; closest=j; }      }      weights[closest] += getWeightOfFrame(i);    }    action->comm.Sum( &weights[0], weights.size() );    myframes->setWeights( weights );  }}

voidtestTranslationRotationMatrix (const IMATH_INTERNAL_NAMESPACE::M44d& mat){    std::cout << "Testing known translate/rotate matrix:/n " << mat;    typedef IMATH_INTERNAL_NAMESPACE::Vec3<T> Vec;    static IMATH_INTERNAL_NAMESPACE::Rand48 rand (2047);    size_t numPoints = 7;    std::vector<Vec> from;  from.reserve (numPoints);    std::vector<Vec> to;    to.reserve (numPoints);    for (size_t i = 0; i < numPoints; ++i)    {        IMATH_INTERNAL_NAMESPACE::V3d a (rand.nextf(), rand.nextf(), rand.nextf());        IMATH_INTERNAL_NAMESPACE::V3d b = a * mat;        from.push_back (Vec(a));        to.push_back (Vec(b));    }    std::vector<T> weights (numPoints, T(1));    const IMATH_INTERNAL_NAMESPACE::M44d m1 = procrustesRotationAndTranslation (&from[0], &to[0], &weights[0], numPoints);    const IMATH_INTERNAL_NAMESPACE::M44d m2 = procrustesRotationAndTranslation (&from[0], &to[0], numPoints);    const T eps = sizeof(T) == 8 ? 1e-8 : 1e-4;    for (size_t i = 0; i < numPoints; ++i)    {        const IMATH_INTERNAL_NAMESPACE::V3d a = from[i];        const IMATH_INTERNAL_NAMESPACE::V3d b = to[i];        const IMATH_INTERNAL_NAMESPACE::V3d b1 = a * m1;        const IMATH_INTERNAL_NAMESPACE::V3d b2 = a * m2;        assert ((b - b1).length() < eps);        assert ((b - b2).length() < eps);    }    std::cout << "  OK/n";}

void move_particles (std::vector <Particle> *particles,        std::vector <double> *field){    auto num_particles = particles->size();    std::vector <double> weights (2);    std::vector <int> points (2);    for (auto& particle : *particles)    {        if (CIC){            weighing (&particle, &weights[0]);        }        else if (ZERO_ORDER){            zero_order_weighing (&particle, &weights[0]);        }        adjacent_points (&particle, &points[0]);        auto left_field = field->at (points [0]) * weights [0];        auto right_field = field->at (points [1]) * weights [1];        auto accel_0 = find_accel(left_field + right_field, &particle);        particle.inc_pos (accel_0);        if (CIC){            weighing (&particle, &weights[0]);        }        else if (ZERO_ORDER){            zero_order_weighing (&particle, &weights[0]);        }        adjacent_points (&particle, &points[0]);        left_field = field->at (points [0]) * weights [0];        right_field = field->at (points [1]) * weights [1];        auto accel_1 = find_accel(left_field + right_field, &particle);        particle.inc_vel (accel_0, accel_1);    }}

void Foam::domainDecomposition::distributeCells(){    Info<< "/nCalculating distribution of cells" << endl;    cpuTime decompositionTime;    autoPtr<decompositionMethod> decomposePtr = decompositionMethod::New    (        decompositionDict_    );    scalarField cellWeights;    if (decompositionDict_.found("weightField"))    {        word weightName = decompositionDict_.lookup("weightField");        volScalarField weights        (            IOobject            (                weightName,                time().timeName(),                *this,                IOobject::MUST_READ,                IOobject::NO_WRITE            ),            *this        );        cellWeights = weights.primitiveField();    }    cellToProc_ = decomposePtr().decompose(*this, cellWeights);    Info<< "/nFinished decomposition in "        << decompositionTime.elapsedCpuTime()        << " s" << endl;}

void TestWeightMatrix_Forward() {  int n_outputs = 2;  InputLayer input_layer(3);  OutputLayer output_layer(n_outputs);  WeightMatrix weights(input_layer, output_layer);  weights.set(0, 0, 0.5);  weights.set(1, 0, -2.0);  weights.set(2, 0, 1.5);  weights.set(0, 1, 1.0);  weights.set(1, 1, 0.7);  weights.set(2, 1, -1.0);  weights.setBias(0, 0.8);  weights.setBias(1, -0.3);  std::vector<double> inputs;  inputs.push_back(-2.0);  inputs.push_back(1.0);  inputs.push_back(3.0);  input_layer.receiveInput(inputs);  std::vector<double> transition = weights.fire(input_layer);  assert(transition.size() == 2);  assert(transition[0] == 2.3);  assert(transition[1] == -4.6);  output_layer.receiveInput(transition);  assert(output_layer.getInput(0) == 2.3);  assert(output_layer.getInput(1) == -4.6);  assert(output_layer.getOutput(0) == 2.3);  assert(output_layer.getOutput(1) == -4.6);  printPass("TestWeightMatrix_Forward()");}

Double KMAlgo::ComputeMssc(IPartition const & x, KMInstance const & instance) {	RectMatrix centers(x.maxNbLabels(), instance.nbAtt());	centers.assign(0);	DoubleVector weights(x.maxNbLabels(), 0);	for (auto const & l : x.usedLabels()) {		for (auto const & i : x.observations(l)) {			weights[l] += instance.weight(i);			for (size_t d(0); d < instance.nbAtt(); ++d)				centers.plus(l, d, instance.get(i, d) * instance.weight(i));		}	}	Double result(0);	for (size_t i(0); i < instance.nbObs(); ++i) {		size_t const l(x.label(i));		for (size_t d(0); d < instance.nbAtt(); ++d)			result += instance.weight(i)					* std::pow(							instance.get(i, d) - centers.get(l, d) / weights[l],							2);	}	return result;}

    double eval(const std::vector<double> &at, bool &valid){        //Allow for the weights to be part of the optimization		//by assuming theyre tacked on to the end of the at vector		//NOTE: the last weight is the one implied by the 1-sum(others)		//Determine where the weights start		int offset = this->P->dimDesign;		//Separate the set of "Problem" variables		const Homotopy::designVars_t vars(at.begin(), at.begin() + offset);        Homotopy::designVars_t weights( at.begin() + offset, at.end() );        //Evaluate the objectives with the designVars        Homotopy::objVars_t result( this->e.eval( vars, valid ) );            //Return immediately if result is unavailable        if(!valid) return 0.0;        else {            //Calculate the remaining weight and add it            weights.push_back( 1.0 - std::accumulate(weights.begin(), weights.end(), 0.0) );            //Return the weighted sum            return std::inner_product( weights.begin(), weights.end(), result.begin(), 0.0 );        }	}

std::vector<double>  computeInterpolationWeights2d(const GenericPoint&  thepoint,                                                    const Container&     pt_v) {   if (pt_v.size() !=3)   {      std::cout << " Compute interpolation weights.. error.. wrong number of points: " << pt_v.size() << std::endl;      return std::vector<double>();   }   double  c1[2], c2[2];   c1[0] = pt_v[1][0] - pt_v[0][0];   c1[1] = pt_v[1][1] - pt_v[0][1];           c2[0] = pt_v[2][0] - pt_v[0][0];   c2[1] = pt_v[2][1] - pt_v[0][1];                   double  det(det2x2(c1, c2));   if (det == 0.0)   {      throw gsse::numerical_calculation_error(":: interpolation :: 2D :: coefficients.. determinant to small.. ");   }           double  rhs[2];   std::vector<double> weights(3);           rhs[0] = thepoint[0] - pt_v[0][0];   rhs[1] = thepoint[1] - pt_v[0][1];   weights[1] = det2x2(rhs, c2)  / det;   weights[2] = det2x2(c1,  rhs) / det;   weights[0] = 1.0 - weights[1] - weights[2];   return weights;   // here .. copy constructor // think about it.. [RH]}

Layer::Matrix Layer::getFlattenedWeights() const{	Matrix weights(1, totalWeights());		size_t position = 0;	for(auto matrix = begin(); matrix != end(); ++matrix)	{		std::memcpy(&weights.data()[position], &matrix->data()[0],			matrix->size() * sizeof(float));		position += matrix->size();	}	for(auto matrix = begin_bias(); matrix != end_bias(); ++matrix)	{		std::memcpy(&weights.data()[position], &matrix->data()[0],			matrix->size() * sizeof(float));		position += matrix->size();	}		return weights;}

void ScoreIndexManager::Debug_PrintLabeledScores(std::ostream& os, const ScoreComponentCollection& scc) const{	std::vector<float> weights(scc.m_scores.size(), 1.0f);	Debug_PrintLabeledWeightedScores(os, scc, weights);}

void _jit_avx512_core_fp32_wino_conv_4x3_t<is_fwd>::_execute_data_W_S_G_D(        float *inp_ptr, float *out_ptr, float *wei_ptr, float *bias_ptr,        const memory_tracking::grantor_t &scratchpad) const {    const auto &jcp = kernel_->jcp;    const auto &p_ops = attr_->post_ops_;    const int inph = is_fwd ? jcp.ih : jcp.oh;    const int inpw = is_fwd ? jcp.iw : jcp.ow;    const int outh = is_fwd ? jcp.oh : jcp.ih;    const int outw = is_fwd ? jcp.ow : jcp.iw;    /* Notation:       FWD: dimM:oc, dimN:ntiles, dimK:ic,       BWD: dimM:ic, dimN:ntiles, dimK:oc,       FWD/BWD: V: src/diff_dst transform, U:weight transform,                M:dst/diff_src transform  */    array_offset_calculator<float, 5> input(inp_ptr,            jcp.mb, jcp.dimK/jcp.dimK_reg_block, inph, inpw,            jcp.dimK_reg_block);    array_offset_calculator<float, 5> output(out_ptr,            jcp.mb, jcp.dimM/jcp.dimM_simd_block, outh, outw,            jcp.dimM_simd_block);    array_offset_calculator<float, 6> weights(wei_ptr,            jcp.oc/jcp.oc_simd_block, jcp.ic/jcp.ic_simd_block, jcp.kh, jcp.kw,            jcp.ic_simd_block, jcp.oc_simd_block);    array_offset_calculator<float, 2> bias(bias_ptr,            jcp.dimM/jcp.dimM_simd_block, jcp.dimM_simd_block);    array_offset_calculator<float, 8> M(is_fwd            ? scratchpad.template get<float>(key_wino_M)            : scratchpad.template get<float>(key_wino_V),            jcp.dimN_nb_block, jcp.dimM_nb_block,            alpha, alpha,            jcp.dimN_block, jcp.dimM_block * jcp.dimM_reg_block,            jcp.dimN_reg_block, jcp.dimM_simd_block);    auto wino_wei = (jcp.prop_kind == prop_kind::forward_inference)            ? wei_ptr            : scratchpad.template get<float>(key_wino_U);    array_offset_calculator<float, 8> U(wino_wei,            jcp.dimM_nb_block,            alpha, alpha,            jcp.dimK_nb_block,            jcp.dimM_block * jcp.dimM_reg_block, jcp.dimK_block,            jcp.dimK_reg_block, jcp.dimM_simd_block);    array_offset_calculator<float, 8> V(is_fwd            ? scratchpad.template get<float>(key_wino_V)            : scratchpad.template get<float>(key_wino_M),            jcp.dimN_nb_block, alpha, alpha,            jcp.dimN_block, jcp.dimK_nb_block,            jcp.dimK_block, jcp.dimN_reg_block, jcp.dimK_reg_block);    const bool wants_padded_bias = jcp.with_bias        && jcp.oc_without_padding != jcp.oc;    float last_slice_bias[simd_w] = {0};    if (wants_padded_bias) {        for (int oc = 0; oc < jcp.oc_without_padding % jcp.oc_simd_block; ++oc)            last_slice_bias[oc] = bias(jcp.dimM / jcp.dimM_simd_block - 1, oc);    }PRAGMA_OMP(parallel)    {        parallel_nd_in_omp(jcp.mb, jcp.dimK_nb_block, jcp.dimK_block,                [&](int img, int K_blk1, int K_blk2) {                input_transform_data(img, jcp,                    &(input(img, K_blk1 * jcp.dimK_block + K_blk2,                            0, 0, 0)),                        &(V(0, 0, 0, 0, K_blk1, K_blk2, 0, 0)));                });        if (jcp.prop_kind != prop_kind::forward_inference) {            parallel_nd_in_omp(jcp.nb_oc, jcp.nb_ic, (jcp.oc_block * jcp.oc_reg_block),                (jcp.ic_block * jcp.ic_reg_block),                [&](int ofm1, int ifm1, int ofm2, int ifm2) {                    float *U_base_ptr = is_fwd                        ? &(U(ofm1, 0, 0, ifm1, ofm2, ifm2, 0, 0))                        : &(U(ifm1, 0, 0, ofm1, ifm2, ofm2, 0, 0));                    weight_transform_data(jcp,                        &(weights(                                ofm1 * jcp.oc_block * jcp.oc_reg_block + ofm2,                                ifm1 * jcp.ic_block * jcp.ic_reg_block + ifm2,                                0, 0, 0, 0)),                        U_base_ptr);            });        }PRAGMA_OMP(barrier)        parallel_nd_in_omp(jcp.dimN_nb_block, alpha, alpha, jcp.dimM_nb_block,            [&](int N_blk1, int oj, int oi, int M_blk1) {            for (int K_blk1 = 0; K_blk1 < jcp.dimK_nb_block;                 K_blk1++)            for (int N_blk2 = 0; N_blk2 < jcp.dimN_block; N_blk2++)                kernel_->gemm_loop_ker(                        (float *)&(M(N_blk1, M_blk1, oj, oi,                            N_blk2, 0, 0, 0)),                        (const float *)&(U(M_blk1, oj, oi,                            K_blk1, 0, 0, 0, 0)),//.........这里部分代码省略.........

int Dmc_method::calcBranch() {   int totwalkers=mpi_info.nprocs*nconfig;  Array1 <doublevar> weights(totwalkers);  Array1 <doublevar> my_weights(nconfig);    for(int walker=0; walker < nconfig; walker++)    my_weights(walker)=pts(walker).weight;#ifdef USE_MPI  MPI_Allgather(my_weights.v,nconfig, MPI_DOUBLE, weights.v,nconfig,MPI_DOUBLE, MPI_Comm_grp);#else  weights=my_weights;#endif  Array1 <int> my_branch(nconfig);  Array1 <int> nwalkers(mpi_info.nprocs);  nwalkers=0;  int root=0;  if(mpi_info.node==root) {  //this if/else clause may be refactored out    Array1 <int> branch(totwalkers);    //----Find which walkers branch/die    //we do it on one node since otherwise we'll have different random numbers!    //we'll assign the weight for each copy that will be produced    //this is the core of the branching algorithm..    //my homegrown algo, based on Umrigar, Nightingale, and Runge    branch=-1;    long int time_a=clock();    match_walkers(weights,branch);    long int time_b=clock();    single_write(cout,"matching: ",double(time_b-time_a)/CLOCKS_PER_SEC,"/n");    for(int w=0; w< totwalkers; w++) {      if(branch(w)==-1) branch(w)=1;    }    //----end homegrown algo    //count how many walkers each node will have    //without balancing    int walk=0;    for(int n=0; n< mpi_info.nprocs; n++) {       for(int i=0; i< nconfig; i++) {        nwalkers(n)+=branch(walk);        walk++;      }      //cout << "nwalkers " << n << " " << nwalkers(n) << endl;    }    //now send nwalkers and which to branch to all processors    for(int i=0; i< nconfig; i++) {       my_branch(i)=branch(i);      my_weights(i)=weights(i);    }    time_a=clock();#ifdef USE_MPI    MPI_Bcast(nwalkers.v, mpi_info.nprocs, MPI_INT, mpi_info.node, MPI_Comm_grp);    for(int i=1; i< mpi_info.nprocs; i++) {      MPI_Send(branch.v+i*nconfig,nconfig,MPI_INT,i,0,MPI_Comm_grp);      MPI_Send(weights.v+i*nconfig, nconfig, MPI_DOUBLE, i,0,MPI_Comm_grp);    }#endif    time_b=clock();    single_write(cout,"sending branch: ",double(time_b-time_a)/CLOCKS_PER_SEC,"/n");                 }  else { #ifdef USE_MPI    MPI_Bcast(nwalkers.v, mpi_info.nprocs, MPI_INT, root, MPI_Comm_grp);    MPI_Status status;    MPI_Recv(my_branch.v,nconfig, MPI_INT,root,0,MPI_Comm_grp, &status);    MPI_Recv(my_weights.v,nconfig, MPI_DOUBLE,root,0,MPI_Comm_grp, &status);#endif	  }  //--end if/else clause  long int time_a=clock();  for(int i=0; i< nconfig; i++) {     pts(i).weight=my_weights(i);  }    //Now we all have my_branch and nwalkers..we need to figure out who   //needs to send walkers to whom--after this, nwalkers should be a flat array equal to   //nconfig(so don't try to use it for anything useful afterwards)  vector <Queue_element> send_queue;  int curr_needs_walker=0;  int nnwalkers=nwalkers(mpi_info.node); //remember how many total we should have  for(int i=0; i< mpi_info.nprocs; i++) {     while(nwalkers(i) > nconfig) {      if(nwalkers(curr_needs_walker) < nconfig) {         nwalkers(curr_needs_walker)++;        nwalkers(i)--;        send_queue.push_back(Queue_element(i,curr_needs_walker));        //cout << mpi_info.node << ":nwalkers " << nwalkers(i) << "  " << nwalkers(curr_needs_walker) << endl;        //cout << mpi_info.node << ":send " << i << "  " << curr_needs_walker << endl;      }      else {         curr_needs_walker++;      }    }  }    for(int i=0; i< mpi_info.nprocs; i++) assert(nwalkers(i)==nconfig);  int killsize=0;//.........这里部分代码省略.........

//.........这里部分代码省略.........  }  else if (responseID == 9) {    return eleInfo.setVector(q);  }  else if (responseID == 19) {    static Matrix kb(3,3);    this->getBasicStiff(kb);    return eleInfo.setMatrix(kb);  }  // Chord rotation  else if (responseID == 3) {    return eleInfo.setVector(crdTransf->getBasicTrialDisp());  }  // Plastic rotation  else if (responseID == 4) {    static Vector vp(3);    static Vector ve(3);    const Matrix &kb = this->getInitialBasicStiff();    kb.Solve(q, ve);    vp = crdTransf->getBasicTrialDisp();    vp -= ve;    return eleInfo.setVector(vp);  }  // Curvature sensitivity  else if (responseID == 5) {    /*      Vector curv(numSections);      const Vector &v = crdTransf->getBasicDispGradient(1);            double L = crdTransf->getInitialLength();      double oneOverL = 1.0/L;      //const Matrix &pts = quadRule.getIntegrPointCoords(numSections);      double pts[2];      pts[0] = 0.0;      pts[1] = 1.0;            // Loop over the integration points      for (int i = 0; i < numSections; i++) {	int order = theSections[i]->getOrder();	const ID &code = theSections[i]->getType();	//double xi6 = 6.0*pts(i,0);	double xi6 = 6.0*pts[i];	curv(i) = oneOverL*((xi6-4.0)*v(1) + (xi6-2.0)*v(2));      }            return eleInfo.setVector(curv);    */    Vector curv(numSections);    /*    // Loop over the integration points    for (int i = 0; i < numSections; i++) {      int order = theSections[i]->getOrder();      const ID &code = theSections[i]->getType();      const Vector &dedh = theSections[i]->getdedh();      for (int j = 0; j < order; j++) {	if (code(j) == SECTION_RESPONSE_MZ)	  curv(i) = dedh(j);      }    }    */    return eleInfo.setVector(curv);  }  // Basic deformation sensitivity  else if (responseID == 6) {      const Vector &dvdh = crdTransf->getBasicDisplSensitivity(1);    return eleInfo.setVector(dvdh);  }  else if (responseID == 7) {    //const Matrix &pts = quadRule.getIntegrPointCoords(numSections);    double xi[maxNumSections];    beamInt->getSectionLocations(numSections, L, xi);    Vector locs(numSections);    for (int i = 0; i < numSections; i++)      locs(i) = xi[i]*L;    return eleInfo.setVector(locs);  }  else if (responseID == 8) {    //const Vector &wts = quadRule.getIntegrPointWeights(numSections);    double wt[maxNumSections];    beamInt->getSectionWeights(numSections, L, wt);    Vector weights(numSections);    for (int i = 0; i < numSections; i++)      weights(i) = wt[i]*L;    return eleInfo.setVector(weights);  }  else    return Element::getResponse(responseID, eleInfo);}

	void withoutReplacementImpl(withoutReplacementArgs& args)	{		boost::numeric::ublas::matrix<double>& matrix = args.matrix;		int n = args.n;		int seed = args.seed;		double alpha = args.alpha;		if(matrix.size1() != matrix.size2())		{			throw std::runtime_error("Matrix must be square");		}		if(n < 1)		{			throw std::runtime_error("Input n must be at least 1");		}		int dimension = matrix.size1();		boost::mt19937 randomSource;		randomSource.seed(seed);		std::vector<double> columnSums(dimension,0);		for(int row = 0; row < dimension; row++)		{			for(int column = 0; column < dimension; column++)			{				columnSums[column] += std::fabs(matrix(row, column));			}		}				std::vector<double> currentColumnSums(dimension), newCurrentColumnSums;		std::vector<int> availableColumns(dimension), newAvailableColumns;		std::vector<int> availableRows(dimension), newAvailableRows;		std::vector<int> previousEntry(1), newPreviousEntry;		std::vector<bool> usedRows(dimension), newUsedRows;		std::vector<mpfr_class> weights(1, 1), newWeights;		//Get out the choices at the start. 		std::vector<possibility> possibilities;		for(int i = 0; i < dimension; i++)		{			for(int j = 0; j < dimension; j++)			{				possibility nextPos;				nextPos.parent = 0;				nextPos.previousEntry = i;				nextPos.newEntry = j;				possibilities.push_back(nextPos);			}		}		//These choices are all derived from a single particle at the start		std::copy(columnSums.begin(), columnSums.end(), currentColumnSums.begin());		for(int i = 0; i < dimension; i++)		{			availableColumns[i] = i;			availableRows[i] = i;		}		std::fill(usedRows.begin(), usedRows.end(), false);		sampling::sampfordFromParetoNaiveArgs samplingArgs;		samplingArgs.n = n;		samplingArgs.weights.resize(dimension*dimension);		for(int i = 0; i < (int)possibilities.size(); i++)		{			possibility& pos = possibilities[i];			if(pos.previousEntry == pos.newEntry)			{				samplingArgs.weights[i] = alpha * std::fabs(matrix(pos.previousEntry, pos.newEntry)) / (dimension*(dimension - 1));			}			else 			{				samplingArgs.weights[i] = std::fabs(matrix(pos.previousEntry, pos.newEntry)) / dimension;			}		}		int currentSamples = 1;		for(int permutationCounter = 0; permutationCounter < dimension; permutationCounter++)		{			//Draw sample			if((int)possibilities.size() <= n)			{				newCurrentColumnSums.resize(possibilities.size()*dimension);				newAvailableColumns.resize(possibilities.size()*(dimension - permutationCounter-1));				newAvailableRows.resize(possibilities.size()*(dimension - permutationCounter-1));				newUsedRows.resize(possibilities.size()*dimension);				newPreviousEntry.resize(possibilities.size());				int newAvailableColumnsIndex = 0;				int newAvailableRowsIndex = 0;				for(int i = 0; i < (int)possibilities.size(); i++)				{					possibility& pos = possibilities[i];					std::copy(currentColumnSums.begin() + pos.parent * dimension, currentColumnSums.begin() + (pos.parent + 1) * dimension, newCurrentColumnSums.begin() + dimension * i);					newCurrentColumnSums[dimension*i + pos.newEntry] -= std::fabs(matrix(pos.previousEntry, pos.newEntry));					std::copy(usedRows.begin() + pos.parent*dimension, usedRows.begin() + (pos.parent + 1) * dimension, newUsedRows.begin() + dimension*i);					newUsedRows[dimension*i + pos.previousEntry] = true;					for(int j = 0; j < dimension - permutationCounter; j++)					{						if(availableColumns[j + pos.parent*(dimension - permutationCounter)] != pos.newEntry)						{							newAvailableColumns.at(newAvailableColumnsIndex) = availableColumns.at(j + pos.parent*(dimension - permutationCounter));							newAvailableColumnsIndex++;						}						if(availableRows[j + pos.parent*(dimension - permutationCounter)] != pos.previousEntry)						{//.........这里部分代码省略.........

//.........这里部分代码省略.........			project.readAndResize("edges", beEdges);			std::set<int> nodes;			for (int n : beNodes)				nodes.insert(n);			std::set<std::pair<int, int>> edges;			for (int i = 0; i < beEdges.shape(1); i++)				edges.insert(						std::make_pair(							std::min(beEdges(i, 0), beEdges(i, 1)),							std::max(beEdges(i, 0), beEdges(i, 1))));			for (Crag::NodeIt n(crag); n != lemon::INVALID; ++n)				bestEffort->node[n] = nodes.count(crag.id(n));			for (Crag::EdgeIt e(crag); e != lemon::INVALID; ++e)				bestEffort->edge[e] = edges.count(						std::make_pair(								std::min(crag.id(crag.u(e)), crag.id(crag.v(e))),								std::max(crag.id(crag.u(e)), crag.id(crag.v(e)))));		} else {			LOG_USER(logger::out) << "reading ground-truth" << std::endl;			ExplicitVolume<int> groundTruth;			volumeStore.retrieveGroundTruth(groundTruth);			LOG_USER(logger::out) << "finding best-effort solution" << std::endl;			overlapLoss = new OverlapLoss(crag, groundTruth);			bestEffort = new BestEffort(crag, *overlapLoss);		}		Loss* loss = 0;		bool  destructLoss = false;		if (optionLoss.as<std::string>() == "hamming") {			LOG_USER(logger::out) << "using Hamming loss" << std::endl;			loss = new HammingLoss(crag, *bestEffort);			destructLoss = true;		} else if (optionLoss.as<std::string>() == "overlap") {			LOG_USER(logger::out) << "using overlap loss" << std::endl;			if (!overlapLoss) {				LOG_USER(logger::out) << "reading ground-truth" << std::endl;				ExplicitVolume<int> groundTruth;				volumeStore.retrieveGroundTruth(groundTruth);				LOG_USER(logger::out) << "finding best-effort solution" << std::endl;				overlapLoss = new OverlapLoss(crag, groundTruth);			}			loss = overlapLoss;		} else {			LOG_ERROR(logger::out)					<< "unknown loss: "					<< optionLoss.as<std::string>()					<< std::endl;			return 1;		}		if (optionNormalizeLoss) {			LOG_USER(logger::out) << "normalizing loss..." << std::endl;			loss->normalize(crag, MultiCut::Parameters());		}		Oracle oracle(				crag,				nodeFeatures,				edgeFeatures,				*loss,				*bestEffort);		std::vector<double> weights(nodeFeatures.dims() + edgeFeatures.dims(), 0);		optimizer.optimize(oracle, weights);		storeVector(weights, optionFeatureWeights);		if (destructLoss && loss != 0)			delete loss;		if (overlapLoss)			delete overlapLoss;		if (bestEffort)			delete bestEffort;	} catch (boost::exception& e) {		handleException(e, std::cerr);	}}

void QCAD::ResponseFieldIntegral<EvalT, Traits>::evaluateFields(typename Traits::EvalData workset){  // Zero out local response  for (typename PHX::MDField<ScalarT>::size_type i=0;        i<this->local_response.size(); i++)    this->local_response[i] = 0.0;  typename std::vector<PHX::MDField<ScalarT,Cell,QuadPoint> >::const_iterator it;  if(opRegion->elementBlockIsInRegion(workset.EBName)) {    ScalarT term, val; //, dbI = 0.0;    std::size_t n, max, nExtraMinuses, nOneBits, nBits = fields.size();    //DEBUG    //std::size_t nContrib1 = 0, nContrib2 = 0;    //ScalarT dbMaxRe[10], dbMaxIm[10];    //for(std::size_t i=0; i<10; i++) dbMaxRe[i] = dbMaxIm[i] = 0.0;    for (std::size_t cell=0; cell < workset.numCells; ++cell) {      if(!opRegion->cellIsInRegion(cell)) continue;      for (std::size_t qp=0; qp < numQPs; ++qp) {	val = 0.0;	//Loop over all possible combinations of Re/Im parts which form product terms and 	// add the relevant ones (depending on whether we're returning the overall real or	// imaginary part of the integral) to get the integrand value for this (cell,qp).	// We do this by mapping the Re/Im choice onto a string of N bits, where N is the 	// number of fields being multiplied together. (0 = RePart, 1 = ImPart)	//nContrib1++; //DEBUG	//for(it = fields.begin(); it != fields.end(); ++it)	max = (std::size_t)std::pow(2.,(int)nBits);//	max = pow(2.0,static_cast<int>(nBits));	for(std::size_t i=0; i<max; i++) {	  // Count the number of 1 bits, and exit early if 	  //  there's a 1 bit for a field that is not complex	  nOneBits = nExtraMinuses = 0;	  for(n=0; n<nBits; n++) {	    if( (0x1 << n) & i ) { // if n-th bit of i is set (use Im part of n-th field)	      if(!fieldIsComplex[n]) break;	      if(conjugateFieldFlag[n]) nExtraMinuses++;	      nOneBits++;	    }		  }	  if(n < nBits) continue;  // we exited early, signaling this product can't contribute	  //check if this combination of Re/Im parts contributes to the overall Re or Im part we return	  if( (bReturnImagPart && nOneBits % 2) || (!bReturnImagPart && nOneBits % 2 == 0)) {	    term = (nOneBits % 4 >= 2) ? -1.0 : 1.0; //apply minus sign if nOneBits % 4 == 2 (-1) or == 3 (-i)	    if(nExtraMinuses % 2) term *= -1.0;      //apply minus sign due to conjugations	    //nContrib2++;	    //multiply fields together	    for(std::size_t m=0; m<nBits; m++) { 	      if( (0x1 << m) & i ) {		term *= fields_Imag[m](cell,qp);		//if( abs(fields_Imag[m](cell,qp)) > dbMaxIm[m]) dbMaxIm[m] = abs(fields_Imag[m](cell,qp));	      }		      else {		term *= fields[m](cell,qp);		//if( abs(fields[m](cell,qp)) > dbMaxRe[m]) dbMaxRe[m] = abs(fields[m](cell,qp));	      }	    }	    val += term;  //add term to overall integrand	  }	}	val *= weights(cell,qp) * scaling; //multiply integrand by volume	//dbI += val; //DEBUG        this->local_response(cell) += val;	this->global_response(0) += val;      }    }    //DEBUG    /*if(fieldNames.size() > 1) {      std::cout << "DB: " << (bReturnImagPart == true ? "Im" : "Re") << " Field Integral - int(";      for(std::size_t i=0; i<fieldNames.size(); i++) 	std::cout << fieldNames[i] << "," << (conjugateFieldFlag[i] ? "-" : "") << (fieldIsComplex[i] ? fieldNames_Imag[i] : "X") << " * ";      std::cout << " dV) -- I += " << dbI << "  (ebName = " << workset.EBName << 	" contrib1=" << nContrib1 << " contrib2=" << nContrib2 << ")" << std::endl;      std::cout << "DB MAX of Fields Re: " << dbMaxRe[0] << "," << dbMaxRe[1] << "," << dbMaxRe[2] << "," << dbMaxRe[3] << std::endl;      std::cout << "DB MAX of Fields Im: " << dbMaxIm[0] << "," << dbMaxIm[1] << "," << dbMaxIm[2] << "," << dbMaxIm[3] << std::endl;      }*/  }  // Do any local-scattering necessary  PHAL::SeparableScatterScalarResponse<EvalT,Traits>::evaluateFields(workset);}

//.........这里部分代码省略.........      {        //MIP        for (idxY = 0; idxY <= maxY; idxY++)        {          //useYMin = ((idxY + outExt[2]) <= wholeExtent[2]) ? 0 : -inIncs[1];          //useYMax = ((idxY + outExt[2]) >= wholeExtent[3]) ? 0 : inIncs[1];          for (idxX = 0; idxX <= maxX; idxX++)          {            //useXMin = ((idxX + outExt[0]) <= wholeExtent[0]) ? 0 : -inIncs[0];            //useXMax = ((idxX + outExt[0]) >= wholeExtent[1]) ? 0 : inIncs[0];            double sum = 0;            for(int z = _minZ; z<= _maxZ;z++)            {              T value = inPtr[z*inIncs[2]];              sum += value;            }            // do X axis            *outPtr = static_cast<T>(invNum*sum);            outPtr++;            inPtr++;          }          outPtr += outIncY;          inPtr += inIncY;        }      }      break;  case vtkMitkThickSlicesFilter::WEIGHTED:    {      const int size = _maxZ-_minZ;      std::vector<double> weights(size);      double mean = 0.5 * double(_minZ + _maxZ);      double sigma_sq = double(size) / 6.0;      sigma_sq *= sigma_sq;      double sum = 0;      int i=0;      for(int z = _minZ+1; z<= _maxZ;z++)      {        double val = exp(-(((double)z-mean)/sigma_sq));        weights[i++] = val;        sum += val;      }      for(i=0; i<size; i++)      {        weights[i] /= sum;      }      for (idxY = 0; idxY <= maxY; idxY++)      {        //useYMin = ((idxY + outExt[2]) <= wholeExtent[2]) ? 0 : -inIncs[1];        //useYMax = ((idxY + outExt[2]) >= wholeExtent[3]) ? 0 : inIncs[1];        for (idxX = 0; idxX <= maxX; idxX++)        {          //useXMin = ((idxX + outExt[0]) <= wholeExtent[0]) ? 0 : -inIncs[0];          //useXMax = ((idxX + outExt[0]) >= wholeExtent[1]) ? 0 : inIncs[0];          T mip = inPtr[_minZ*inIncs[2]];          i=0;          double mymip = 0;          for(int z = _minZ+1; z<= _maxZ;z++)          {            double value = inPtr[z*inIncs[2]];            mymip+=value*weights[i++];

//.........这里部分代码省略.........            AutoPtr<ErrorEstimator> error_estimator;            // To solve to a tolerance in this problem we            // need a better estimator than Kelly            if (global_tolerance != 0.)            {                // We can't adapt to both a tolerance and a mesh                // size at once                libmesh_assert_equal_to (nelem_target, 0);                UniformRefinementEstimator *u =                    new UniformRefinementEstimator;                // The lid-driven cavity problem isn't in H1, so                // lets estimate L2 error                u->error_norm = L2;                error_estimator.reset(u);            }            else            {                // If we aren't adapting to a tolerance we need a                // target mesh size                libmesh_assert_greater (nelem_target, 0);                // Kelly is a lousy estimator to use for a problem                // not in H1 - if we were doing more than a few                // timesteps we'd need to turn off or limit the                // maximum level of our adaptivity eventually                error_estimator.reset(new KellyErrorEstimator);            }            // Calculate error            std::vector<Real> weights(9,1.0);  // based on u, v, p, c, their adjoints, and source parameter            // Keep the same default norm type.            std::vector<FEMNormType>            norms(1, error_estimator->error_norm.type(0));            error_estimator->error_norm = SystemNorm(norms, weights);            error_estimator->estimate_error(system, error);            // Print out status at each adaptive step.            Real global_error = error.l2_norm();            std::cout << "Adaptive step " << a_step << ": " << std::endl;            if (global_tolerance != 0.)                std::cout << "Global_error = " << global_error                          << std::endl;            if (global_tolerance != 0.)                std::cout << "Worst element error = " << error.maximum()                          << ", mean = " << error.mean() << std::endl;            if (global_tolerance != 0.)            {                // If we've reached our desired tolerance, we                // don't need any more adaptive steps                if (global_error < global_tolerance)                    break;                mesh_refinement.flag_elements_by_error_tolerance(error);            }            else            {                // If flag_elements_by_nelem_target returns true, this                // should be our last adaptive step.                if (mesh_refinement.flag_elements_by_nelem_target(error))                {

	void potential_and_gradient(const Eigen::VectorXd& parameters, const Eigen::VectorXd& hyperparameters, View& view, double& potential, Eigen::VectorXd& gradient)	{		// Loop over layers to calculate weights part of potential, and non-data part of gradient		potential = 0;		for (size_t layer_idx = 0; layer_idx < count_weights_layers(); layer_idx++) {			//potential -= 0.5 * (hyperparameters[layer_idx * 2] * weights(parameters, layer_idx).squaredNorm() + hyperparameters[layer_idx * 2 + 1] * biases(parameters, layer_idx).squaredNorm());			potential -= 0.5 * (hyperparameters[0] * weights(parameters, layer_idx).squaredNorm() + hyperparameters[1] * biases(parameters, layer_idx).squaredNorm());			// TODO: Debugging here!			//weights(gradient, layer_idx) = (weights(parameters, layer_idx).array() * -hyperparameters[layer_idx * 2]).matrix();			//biases(gradient, layer_idx) = (biases(parameters, layer_idx).array() * -hyperparameters[layer_idx * 2 + 1]).matrix();			weights(gradient, layer_idx) = (weights(parameters, layer_idx).array() * -hyperparameters[0]).matrix();			biases(gradient, layer_idx) = (biases(parameters, layer_idx).array() * -hyperparameters[1]).matrix();		}		/*if (std::isnan(gradient[0])) {			std::cout << gradient[0] << std::endl;		}*/		// Calculate output part of potential and gradient		for (size_t data_idx = 0; data_idx < view.size(); data_idx++) {			// Get the class label for this observation			size_t class_idx = get_nonzero_idx(view.second(data_idx));			// Calculate the output for this sample, and the gradient of the output with respect to the parameters			// gradient_and_output(size_t variable_idx, const Eigen::VectorXd& inputs, const Eigen::VectorXd& parameters, Eigen::VectorXd& outputs, Eigen::VectorXd& gradient_vector)			/*if (std::isnan(temp_gradient_[0])) {				std::cout << temp_gradient_[0] << std::endl;			}*/			log_gradient_and_output(class_idx, view.first(data_idx), parameters, outputs(), temp_gradient_);			//if (outputs()[class_idx] != 0.)				gradient = gradient + temp_gradient_;						/*if (std::isnan(temp_gradient_[0])) {				std::cout << temp_gradient_[0] << std::endl;			}			if (std::isnan(gradient[0])) {				std::cout << gradient[0] << std::endl;			}*/			//if ()			// NOTE: Does it matter here when -E[theta] = -INF?			//potential += log(outputs()[class_idx]);			potential += outputs()[class_idx];		}		// DEBUG: Check that all entries are finite and not NaN		/*if (!std::isfinite(potential)) {			if (std::isnan(potential))				std::cout << "NaN: Potential" << std::endl;			else if (std::isinf(potential))				std::cout << "INF: Potential" << std::endl;		}		for (size_t idx = 0; idx < static_cast<size_t>(gradient.size()); idx++) {			if (!std::isfinite(gradient[idx])) {				if (std::isnan(gradient[idx]))					std::cout << "NaN: Gradient[" << idx << "]" << std::endl;				else if (std::isinf(gradient[idx]))					std::cout << "NaN: Gradient[" << idx << "]" << std::endl;			}		}*/	}

//.........这里部分代码省略.........    ws->import(*sbres);    ws->import(fitModel);    RooPlot *fmgg = mgg.frame();    data.plotOn(fmgg);    fitModel.plotOn(fmgg);    ws->pdf("b_"+tag+"_ext")->plotOn(fmgg,RooFit::LineColor(kRed),RooFit::Range("Full"),RooFit::NormRange("Full"));    fmgg->SetName(tag+"_frame");    ws->import(*fmgg);    delete fmgg;    RooMinuit(*nll).migrad();    RooPlot *fNs = Nsig.frame(0,25);    fNs->SetName(tag+"_Nsig_pll");    RooAbsReal *pll = nll->createProfile(Nsig);    //nll->plotOn(fNs,RooFit::ShiftToZero(),RooFit::LineColor(kRed));    pll->plotOn(fNs);    ws->import(*fNs);    delete fNs;    RooPlot *fmu = mu.frame(125,132);    fmu->SetName(tag+"_mu_pll");    RooAbsReal *pll_mu = nll->createProfile(mu);    pll_mu->plotOn(fmu);    ws->import(*fmu);    delete fmu;  }  RooArgSet weights("weights");  RooArgSet pdfs_bonly("pdfs_bonly");  RooArgSet pdfs_b("pdfs_b");  RooRealVar minAIC("minAIC","",1E10);  //compute AIC stuff  for(auto t = tags.begin(); t!=tags.end(); t++) {    RooAbsPdf *p_bonly = ws->pdf("bonly_"+*t);    RooAbsPdf *p_b = ws->pdf("b_"+*t);    RooFitResult *sb = (RooFitResult*)ws->obj(*t+"_bonly_fitres");    RooRealVar k(*t+"_b_k","",p_bonly->getParameters(RooArgSet(mgg))->getSize());    RooRealVar nll(*t+"_b_minNll","",sb->minNll());    RooRealVar Npts(*t+"_b_N","",blind_data->sumEntries());    RooFormulaVar AIC(*t+"_b_AIC","2*@0+2*@1+2*@1*(@1+1)/(@[email
C++ weprintf函数代码示例
C++ weight_sscanf函数代码示例