C++ MCWalkerConfiguration::end方法代码示例

 //NOTE: weights are not handled nicely.
 //weight should be done carefully, not valid for DMC
 //will add a function to MCWalkerConfiguration to track total weight
 void EstimatorManager::accumulate(MCWalkerConfiguration& W)
 {
   BlockWeight += W.getActiveWalkers();
   RealType norm=1.0/W.getGlobalNumWalkers();
   for(int i=0; i< Estimators.size(); i++) 
     Estimators[i]->accumulate(W,W.begin(),W.end(),norm);
 }

  int WalkerControlBase::branch(int iter, MCWalkerConfiguration& W, RealType trigger) {

    int prefill_numwalkers = sortWalkers(W);

    measureProperties(iter);
    W.EnsembleProperty=EnsembleProperty;

    //un-biased variance but we use the saimple one
    //W.EnsembleProperty.Variance=(e2sum*wsum-esum*esum)/(wsum*wsum-w2sum);

    ////add to the accumData for block average: REMOVE THIS
    //accumData[ENERGY_INDEX]     += curData[ENERGY_INDEX]*wgtInv;
    //accumData[ENERGY_SQ_INDEX]  += curData[ENERGY_SQ_INDEX]*wgtInv;
    //accumData[WALKERSIZE_INDEX] += curData[WALKERSIZE_INDEX];
    //accumData[WEIGHT_INDEX]     += curData[WEIGHT_INDEX];

    int nw_tot = copyWalkers(W);

    //set Weight and Multiplicity to default values
    MCWalkerConfiguration::iterator it(W.begin()),it_end(W.end());
    while(it != it_end) {
      (*it)->Weight= 1.0;
      (*it)->Multiplicity=1.0;
      ++it;
    }

    //set the global number of walkers
    W.setGlobalNumWalkers(nw_tot);

    return prefill_numwalkers;
  }

/** accumulate Local energies and collectables
 * @param W ensemble
 */
void EstimatorManager::accumulate(MCWalkerConfiguration& W)
{
  BlockWeight += W.getActiveWalkers();
  RealType norm=1.0/W.getGlobalNumWalkers();
  for(int i=0; i< Estimators.size(); i++)
    Estimators[i]->accumulate(W,W.begin(),W.end(),norm);
  if(Collectables)//collectables are normalized by QMC drivers
    Collectables->accumulate_all(W.Collectables,1.0);
}

/**  Write the set of walker configurations to the HDF5 file.  
 *@param W set of walker configurations
 *@param ic the number of frames
 *
 * \if ic==-1
 *  use only the last frame for a restart
 * \else if ic>=0
 *  use ic frames from the file for opitimizations
 */
bool  
HDFWalkerInput0::put(MCWalkerConfiguration& W, int ic){

  if(Counter<0) return false;

  int selected = ic;
  if(ic<0) {
    XMLReport("Will use the last set from " << NumSets << " of configurations.")
    selected = NumSets-1;
  }

  typedef MCWalkerConfiguration::PosType PosType;
  typedef MCWalkerConfiguration::PropertyContainer_t ProtertyContainer_t;

  typedef Matrix<PosType>  PosContainer_t;

  int nwt = 0;
  int npt = 0;
  //2D array of PosTypes (x,y,z) indexed by (walker,particle)
  PosContainer_t Pos_temp;

  //open the group
  char GrpName[128];
  sprintf(GrpName,"config%04d",selected);
  hid_t group_id = H5Gopen(h_config,GrpName);
    
  HDFAttribIO<PosContainer_t> Pos_in(Pos_temp);
  //read the dataset
  Pos_in.read(group_id,"coord");
  //close the group
  H5Gclose(group_id);

  /*check to see if the number of walkers and particles is  consistent with W */
  int nptcl = Pos_temp.cols();
  nwt = Pos_temp.rows();

  int curWalker = W.getActiveWalkers();
  if(curWalker) {
    LOGMSG("Adding " << nwt << " walkers to " << curWalker)
    W.createWalkers(nwt);
  } else {
    W.resize(nwt,nptcl); 
  }

  //assign configurations to W
  int iw=0;
  MCWalkerConfiguration::iterator it = W.begin()+curWalker; 
  MCWalkerConfiguration::iterator it_end = W.end(); 
  while(it != it_end) {
    std::copy(Pos_temp[iw],Pos_temp[iw+1], (*it)->R.begin());
    ++it;++iw;
  }

  return true;
}

int 
WalkerControlMPI::branch(int iter, MCWalkerConfiguration& W, RealType trigger) {

  std::fill(curData.begin(),curData.end(),0);
  //std::fill(NumPerNode.begin(),NumPerNode.end(),0);
  sortWalkers(W);

  //update the number of walkers for this node
  curData[LE_MAX+MyContext]=NumWalkers;

  int nw = copyWalkers(W);

  myComm->allreduce(curData);

  RealType wgtInv(1.0/curData[WEIGHT_INDEX]);
  accumData[ENERGY_INDEX]     += curData[ENERGY_INDEX]*wgtInv;
  accumData[ENERGY_SQ_INDEX]  += curData[ENERGY_SQ_INDEX]*wgtInv;
  accumData[WALKERSIZE_INDEX] += curData[WALKERSIZE_INDEX];
  accumData[WEIGHT_INDEX]     += curData[WEIGHT_INDEX];

  Cur_pop=0;
  for(int i=0, j=LE_MAX; i<NumContexts; i++,j++) {
    Cur_pop+= NumPerNode[i]=static_cast<int>(curData[j]);
  }

  swapWalkersSimple(W); 

  //Do not need to use a trigger.
  //Cur_min=Nmax; 
  //Cur_max=0; 
  //Cur_pop=0;
  //for(int i=0, j=LE_MAX; i<NumContexts; i++,j++) {
  //  Cur_pop+= NumPerNode[i]=static_cast<int>(curData[j]);
  //  Cur_min = std::min(Cur_min,NumPerNode[i]);
  //  Cur_max = std::max(Cur_max,NumPerNode[i]);
  //}
  //int max_diff = std::max(Cur_max*NumContexts-Cur_pop,Cur_pop-Cur_min*NumContexts);
  //double diff_pop = static_cast<double>(max_diff)/static_cast<double>(Cur_pop);
  //if(diff_pop > trigger) { 
  //  swapWalkersSimple(W); 
  //  //swapWalkersMap(W); 
  //}

  //set Weight and Multiplicity to default values
  MCWalkerConfiguration::iterator it(W.begin()),it_end(W.end());
  while(it != it_end) {
    (*it)->Weight= 1.0;
    (*it)->Multiplicity=1.0;
    ++it;
  }

  //update the global number of walkers
  W.setGlobalNumWalkers(Cur_pop);
  return Cur_pop;
}

/**
   @brief accumulate data for all the estimators
*/
void 
ScalarEstimatorManager::accumulate(const MCWalkerConfiguration& W) {

  for(MCWalkerConfiguration::const_iterator it = W.begin(); 
      it != W.end(); it++){    
    RealType wgt = (*it)->Properties(WEIGHT);
    WeightSum += wgt;
    for(int i=0; i< Estimators.size(); i++) 
      Estimators[i]->accumulate(**it,wgt);
  }
}

 void WalkerControlBase::setWalkerID(MCWalkerConfiguration& walkers)
 {
   start(); //do the normal start
   MCWalkerConfiguration::iterator wit(walkers.begin());
   MCWalkerConfiguration::iterator wit_end(walkers.end());
   for(; wit != wit_end; ++wit)
   {
     if((*wit)->ID==0) 
     {
       (*wit)->ID=(++NumWalkersCreated)*NumContexts+MyContext;
       (*wit)->ParentID=(*wit)->ID;
     }
   }
 }

/** accumulate data for all the estimators
 */
void 
ScalarEstimatorManager::accumulate(const MCWalkerConfiguration& W) {

  RealType wgt_sum=0;
  MCWalkerConfiguration::const_iterator it(W.begin()),it_end(W.end());
  while(it != it_end) {
    RealType wgt = (*it)->Weight;
    wgt_sum+= wgt;
    for(int i=0; i< Estimators.size(); i++) Estimators[i]->accumulate(**it,wgt);
    ++it;
  }

  MyData[WEIGHT_INDEX]+=wgt_sum;
}

 void EstimatorManager::getCurrentStatistics(MCWalkerConfiguration& W
     , RealType& eavg, RealType& var)
 {
   LocalEnergyOnlyEstimator energynow;
   energynow.clear();
   energynow.accumulate(W,W.begin(),W.end(),1.0);
   vector<RealType> tmp(3);
   tmp[0]= energynow.scalars[0].result();
   tmp[1]= energynow.scalars[0].result2();
   tmp[2]= energynow.scalars[0].count();
   myComm->allreduce(tmp);
   eavg=tmp[0]/tmp[2];
   var=tmp[1]/tmp[2]-eavg*eavg;
 }

 void Write2XYZ(MCWalkerConfiguration& W)
 {
   ofstream fout("bad.xyz");
   MCWalkerConfiguration::iterator it(W.begin());
   MCWalkerConfiguration::iterator it_end(W.end());
   int nptcls(W.getTotalNum());
   while(it != it_end) {
     fout << nptcls << endl 
       << "# E = " << (*it)->Properties(LOCALENERGY) 
       << " Wgt= " << (*it)->Weight << endl;
     for(int i=0; i<nptcls; i++)
       fout << "H " << (*it)->R[i] << endl;
     ++it;
   }
 }

int WalkerControlBase::doNotBranch(int iter, MCWalkerConfiguration& W)
{
    MCWalkerConfiguration::iterator it(W.begin());
    MCWalkerConfiguration::iterator it_end(W.end());
    RealType esum=0.0,e2sum=0.0,wsum=0.0,ecum=0.0, w2sum=0.0;
    RealType r2_accepted=0.0,r2_proposed=0.0;
    for(; it!=it_end; ++it)
    {
        r2_accepted+=(*it)->Properties(R2ACCEPTED);
        r2_proposed+=(*it)->Properties(R2PROPOSED);
        RealType e((*it)->Properties(LOCALENERGY));
        int nc= std::min(static_cast<int>((*it)->Multiplicity),MaxCopy);
        RealType wgt((*it)->Weight);
        esum += wgt*e;
        e2sum += wgt*e*e;
        wsum += wgt;
        w2sum += wgt*wgt;
        ecum += e;
    }

    //temp is an array to perform reduction operations
    std::fill(curData.begin(),curData.end(),0);

    curData[ENERGY_INDEX]=esum;
    curData[ENERGY_SQ_INDEX]=e2sum;
    curData[WALKERSIZE_INDEX]=W.getActiveWalkers();
    curData[WEIGHT_INDEX]=wsum;
    curData[EREF_INDEX]=ecum;
    curData[R2ACCEPTED_INDEX]=r2_accepted;
    curData[R2PROPOSED_INDEX]=r2_proposed;

    myComm->allreduce(curData);

    measureProperties(iter);
    trialEnergy=EnsembleProperty.Energy;
    W.EnsembleProperty=EnsembleProperty;

    //return the current data
    return W.getGlobalNumWalkers();
}

int 
WalkerReconfigurationMPI::branch(int iter, MCWalkerConfiguration& W, RealType trigger) {

  int nwkept = swapWalkers(W);

  RealType wgtInv(1.0/curData[WEIGHT_INDEX]);
  accumData[ENERGY_INDEX]     += curData[ENERGY_INDEX]*wgtInv;
  accumData[ENERGY_SQ_INDEX]  += curData[ENERGY_SQ_INDEX]*wgtInv;
  accumData[WALKERSIZE_INDEX] += nwkept;
  //accumData[WALKERSIZE_INDEX] += curData[WALKERSIZE_INDEX];
  accumData[WEIGHT_INDEX]     += curData[WEIGHT_INDEX];

  //set Weight and Multiplicity to default values
  MCWalkerConfiguration::iterator it(W.begin()),it_end(W.end());
  while(it != it_end) {
    (*it)->Weight= 1.0;
    (*it)->Multiplicity=1.0;
    ++it;
  }

  return nwkept;
}

int 
GlobalWalkerControl::branch(int iter, MCWalkerConfiguration& W, RealType trigger) {

  std::fill(NumPerNode.begin(),NumPerNode.end(),0);

  sortWalkers(W);

  NumPerNode[MyContext] = NumWalkers;

  int nw = copyWalkers(W);

  //wait until everynode comes here
  OHMMS::Controller->barrier();
  gsum(NumPerNode,0);

  Cur_min=Nmax; Cur_max=0; Cur_pop=0;
  for(int i=0; i<NumContexts; i++) {
    Cur_pop+= NumPerNode[i];
    Cur_min = std::min(Cur_min,NumPerNode[i]);
    Cur_max = std::max(Cur_max,NumPerNode[i]);
  }

  int max_diff = std::max(Cur_max*NumContexts-Cur_pop,Cur_pop-Cur_min*NumContexts);
  double diff_pop = static_cast<double>(max_diff)/static_cast<double>(Cur_pop);

  if(diff_pop > trigger) { swapWalkersMap(W); }

  //set Weight and Multiplicity to default values
  MCWalkerConfiguration::iterator it(W.begin()),it_end(W.end());
  while(it != it_end) {
    (*it)->Weight= 1.0;
    (*it)->Multiplicity=1.0;
    ++it;
  }

  return Cur_pop;
}

  /** evaluate curData and mark the bad/good walkers
   */
  int WalkerControlBase::sortWalkers(MCWalkerConfiguration& W) {

    MCWalkerConfiguration::iterator it(W.begin());

    vector<Walker_t*> bad,good_rn;
    vector<int> ncopy_rn;
    NumWalkers=0;
    MCWalkerConfiguration::iterator it_end(W.end());
    RealType esum=0.0,e2sum=0.0,wsum=0.0,ecum=0.0, w2sum=0.0, besum=0.0, bwgtsum=0.0;
    RealType r2_accepted=0.0,r2_proposed=0.0;
    int nrn(0),ncr(0);
    while(it != it_end) 
    {
      bool inFN=(((*it)->ReleasedNodeAge)==0);
      if ((*it)->ReleasedNodeAge==1) ncr+=1;
      int nc= std::min(static_cast<int>((*it)->Multiplicity),MaxCopy);
      r2_accepted+=(*it)->Properties(R2ACCEPTED);
      r2_proposed+=(*it)->Properties(R2PROPOSED);
      RealType e((*it)->Properties(LOCALENERGY));
      RealType bfe((*it)->Properties(ALTERNATEENERGY));
      RealType rnwgt(0.0);
      if (inFN)
        rnwgt=((*it)->Properties(SIGN));
      
//       RealType wgt((*it)->Weight);
      RealType wgt(0.0);
      if (inFN)
        wgt=((*it)->Weight); 
      
      esum += wgt*e;
      e2sum += wgt*e*e;
      wsum += wgt;
      w2sum += wgt*wgt;
      ecum += e;
      besum += bfe*rnwgt*wgt;
      bwgtsum += rnwgt*wgt;

      if((nc) && (inFN))
      {
        NumWalkers += nc;
        good_w.push_back(*it);
        ncopy_w.push_back(nc-1);
      }
      else if (nc)
      {
        NumWalkers += nc;
        nrn+=nc;
        good_rn.push_back(*it);
        ncopy_rn.push_back(nc-1);
      }
      else
      {
        bad.push_back(*it);
      }
      ++it;
    }

    //temp is an array to perform reduction operations
    std::fill(curData.begin(),curData.end(),0);

    //update curData
    curData[ENERGY_INDEX]=esum;
    curData[ENERGY_SQ_INDEX]=e2sum;
    curData[WALKERSIZE_INDEX]=W.getActiveWalkers();
    curData[WEIGHT_INDEX]=wsum;
    curData[EREF_INDEX]=ecum;
    curData[R2ACCEPTED_INDEX]=r2_accepted;
    curData[R2PROPOSED_INDEX]=r2_proposed;
    curData[FNSIZE_INDEX]=static_cast<RealType>(good_w.size());
    curData[RNONESIZE_INDEX]=static_cast<RealType>(ncr);
    curData[RNSIZE_INDEX]=nrn;
    curData[B_ENERGY_INDEX]=besum;
    curData[B_WGT_INDEX]=bwgtsum;
    
    ////this should be move
    //W.EnsembleProperty.NumSamples=curData[WALKERSIZE_INDEX];
    //W.EnsembleProperty.Weight=curData[WEIGHT_INDEX];
    //W.EnsembleProperty.Energy=(esum/=wsum);
    //W.EnsembleProperty.Variance=(e2sum/wsum-esum*esum);
    //W.EnsembleProperty.Variance=(e2sum*wsum-esum*esum)/(wsum*wsum-w2sum);

    //remove bad walkers empty the container
    for(int i=0; i<bad.size(); i++) delete bad[i];
    if (!WriteRN)
    {
     if(good_w.empty()) {
      app_error() << "All the walkers have died. Abort. " << endl;
      APP_ABORT("WalkerControlBase::sortWalkers");
     }

     int sizeofgood = good_w.size();
     //check if the projected number of walkers is too small or too large
     if(NumWalkers>Nmax) {
      int nsub=0;
      int nsub_target=(NumWalkers-nrn)-static_cast<int>(0.9*Nmax);
      int i=0;
      while(i< sizeofgood && nsub<nsub_target) {
        if(ncopy_w[i]) {ncopy_w[i]--; nsub++;}
//.........这里部分代码省略.........

  int WalkerControlBase::doNotBranch(int iter, MCWalkerConfiguration& W)
  {
    MCWalkerConfiguration::iterator it(W.begin());
    MCWalkerConfiguration::iterator it_end(W.end());
    RealType esum=0.0,e2sum=0.0,wsum=0.0,ecum=0.0, w2sum=0.0, besum=0.0, bwgtsum=0.0;
    RealType r2_accepted=0.0,r2_proposed=0.0;
    int nrn(0),ncr(0),nfn(0),ngoodfn(0);
    for(; it!=it_end;++it)
    {
      bool inFN=(((*it)->ReleasedNodeAge)==0);
      int nc= std::min(static_cast<int>((*it)->Multiplicity),MaxCopy);
      
      if ((*it)->ReleasedNodeAge==1) ncr+=1;
      else if ((*it)->ReleasedNodeAge==0) 
      {
        nfn+=1;
        ngoodfn+=nc;
      }
      
      r2_accepted+=(*it)->Properties(R2ACCEPTED);
      r2_proposed+=(*it)->Properties(R2PROPOSED);
      RealType e((*it)->Properties(LOCALENERGY));
      RealType bfe((*it)->Properties(ALTERNATEENERGY));
      RealType rnwgt(0.0);
      if (inFN)
        rnwgt=((*it)->Properties(SIGN));
      else
        nrn+=nc;
      //       RealType wgt((*it)->Weight);
      RealType wgt(0.0);
      if (inFN)
        wgt=((*it)->Weight); 
      
      esum += wgt*e;
      e2sum += wgt*e*e;
      wsum += wgt;
      w2sum += wgt*wgt;
      ecum += e;
      besum += bfe*rnwgt*wgt;
      bwgtsum += rnwgt*wgt;
    }

    //temp is an array to perform reduction operations
    std::fill(curData.begin(),curData.end(),0);

    curData[ENERGY_INDEX]=esum;
    curData[ENERGY_SQ_INDEX]=e2sum;
    curData[WALKERSIZE_INDEX]=W.getActiveWalkers();
    curData[WEIGHT_INDEX]=wsum;
    curData[EREF_INDEX]=ecum;
    curData[R2ACCEPTED_INDEX]=r2_accepted;
    curData[R2PROPOSED_INDEX]=r2_proposed;
    curData[FNSIZE_INDEX]=static_cast<RealType>(nfn);
    curData[RNONESIZE_INDEX]=static_cast<RealType>(ncr);
    curData[RNSIZE_INDEX]=nrn;
    curData[B_ENERGY_INDEX]=besum;
    curData[B_WGT_INDEX]=bwgtsum;

    myComm->allreduce(curData);

    measureProperties(iter);
    trialEnergy=EnsembleProperty.Energy;
    W.EnsembleProperty=EnsembleProperty;

    return W.getActiveWalkers();
  }