C++ DataSet::ModifyDim方法代码示例

// Analysis_TI::Setup()
Analysis::RetType Analysis_TI::Setup(ArgList& analyzeArgs, AnalysisSetup& setup, int debugIn)
{
  debug_ = debugIn;
  int nq = analyzeArgs.getKeyInt("nq", 0);
  ArgList nskipArg(analyzeArgs.GetStringKey("nskip"), ","); // Comma-separated
  avg_increment_ = analyzeArgs.getKeyInt("avgincrement", -1);
  avg_max_ = analyzeArgs.getKeyInt("avgmax", -1);
  avg_skip_ = analyzeArgs.getKeyInt("avgskip", 0);
  n_bootstrap_pts_ = analyzeArgs.getKeyInt("bs_pts", -1);
  n_bootstrap_samples_ = analyzeArgs.getKeyInt("bs_samples", 0);
  bootstrap_seed_ = analyzeArgs.getKeyInt("bs_seed", -1);
  bootstrap_fac_ = analyzeArgs.getKeyDouble("bs_fac", 0.75);
  if (!nskipArg.empty()) {
    avgType_ = SKIP;
    // Specified numbers of points to skip
    nskip_.clear();
    for (int i = 0; i != nskipArg.Nargs(); i++) {
      nskip_.push_back( nskipArg.getNextInteger(0) );
      if (nskip_.back() < 0) nskip_.back() = 0;
    }
  } else if (avg_increment_ > 0)
    avgType_ = INCREMENT;
  else if (n_bootstrap_samples_ > 0)
    avgType_ = BOOTSTRAP;
  else
    avgType_ = AVG;
  masterDSL_ = setup.DslPtr();
  // Get lambda values
  ArgList xArgs(analyzeArgs.GetStringKey("xvals"), ","); // Also comma-separated
  if (!xArgs.empty()) {
    xval_.clear();
    for (int i = 0; i != xArgs.Nargs(); i++)
      xval_.push_back( xArgs.getNextDouble(0.0) );
  }
  std::string setname = analyzeArgs.GetStringKey("name");
  DataFile* outfile = setup.DFL().AddDataFile(analyzeArgs.GetStringKey("out"), analyzeArgs);
  curveout_ = setup.DFL().AddDataFile(analyzeArgs.GetStringKey("curveout"), analyzeArgs);
  // Select datasets from remaining args
  if (input_dsets_.AddSetsFromArgs( analyzeArgs.RemainingArgs(), setup.DSL() )) {
    mprinterr("Error: Could not add data sets.\n");
    return Analysis::ERR;
  }
  if (input_dsets_.empty()) {
    mprinterr("Error: No input data sets.\n");
    return Analysis::ERR;
  }
  if (SetQuadAndWeights(nq)) return Analysis::ERR;
  // Determine integration mode
  if (nq > 0)
    mode_ = GAUSSIAN_QUAD;
  else
    mode_ = TRAPEZOID;
  // Check that # abscissas matches # data sets
  if (xval_.size() != input_dsets_.size()) {
     mprinterr("Error: Expected %zu data sets for integration, got %zu\n",
               input_dsets_.size(), xval_.size());
    return Analysis::ERR;
  }
  // Set up output data sets
  DataSet::DataType dtype = DataSet::DOUBLE;
  if (avgType_ == SKIP || avgType_ == INCREMENT)
    dtype = DataSet::XYMESH;
  dAout_ = setup.DSL().AddSet(dtype, setname, "TI");
  if (dAout_ == 0) return Analysis::ERR;
  if (outfile != 0) outfile->AddDataSet( dAout_ );
  MetaData md(dAout_->Meta().Name(), "TIcurve");
  if (avgType_ == AVG) {
    // Single curve
    curve_.push_back( setup.DSL().AddSet(DataSet::XYMESH, md) );
    if (curve_.back() == 0) return Analysis::ERR;
    curve_.back()->ModifyDim(Dimension::X).SetLabel("Lambda");
    if (curveout_ != 0) curveout_->AddDataSet( curve_.back() );
    if (outfile != 0) outfile->ProcessArgs("noxcol");
  } else if (avgType_ == SKIP) {
    // As many curves as skip values
    for (Iarray::const_iterator it = nskip_.begin(); it != nskip_.end(); ++it) {
      md.SetIdx( *it );
      DataSet* ds = setup.DSL().AddSet(DataSet::XYMESH, md);
      if (ds == 0) return Analysis::ERR;
      ds->ModifyDim(Dimension::X).SetLabel("Lambda");
      ds->SetLegend( md.Name() + "_Skip" + integerToString(*it) );
      if (curveout_ != 0) curveout_->AddDataSet( ds );
      curve_.push_back( ds );
    }
  } else if (avgType_ == BOOTSTRAP) {
    // As many curves as resamples
    for (int nsample = 0; nsample != n_bootstrap_samples_; nsample++) {
      md.SetIdx(nsample);
      DataSet* ds = setup.DSL().AddSet(DataSet::XYMESH, md);
      if (ds == 0) return Analysis::ERR;
      ds->ModifyDim(Dimension::X).SetLabel("Lambda");
      ds->SetLegend( md.Name() + "_Sample" + integerToString(nsample) );
      if (curveout_ != 0) curveout_->AddDataSet( ds );
      curve_.push_back( ds );
    }
    // Standard devation of avg free energy over samples
    dA_SD_ = setup.DSL().AddSet(DataSet::DOUBLE, MetaData(md.Name(), "SD"));
    if (dA_SD_ == 0) return Analysis::ERR;
    if (outfile != 0) {
//.........这里部分代码省略.........

// Analysis_TI::Calc_Increment()
int Analysis_TI::Calc_Increment() {
  // Determine max points if not given.
  int maxpts = avg_max_;
  if (maxpts == -1) {
    for (unsigned int idx = 0; idx != input_dsets_.size(); idx++) {
      DataSet_1D const& ds = static_cast<DataSet_1D const&>( *(input_dsets_[idx]) );
      if (maxpts == -1)
        maxpts = (int)ds.Size();
      else if (maxpts != (int)ds.Size()) {
        mprintf("Warning: # points in '%s' (%zu) is different than %i.\n",
                ds.legend(), ds.Size(), maxpts);
        maxpts = std::min( maxpts, (int)ds.Size() );
        mprintf("Warning:   Will only use %i points.\n", maxpts);
      }
    }
  }
  if (maxpts < 1) {
    mprinterr("Error: Max points to use is < 1.\n");
    return 1;
  }
  if (avg_skip_ >= maxpts) {
    mprinterr("Error: 'avgskip' (%i) > max (%i).\n", avg_skip_, maxpts);
    return 1;
  }
  // sum: Hold the results of integration for each curve (increment)
  Darray sum;
  // points: Hold point values at which each avg is being calculated
  Iarray points;
  // Loop over input data sets. 
  for (unsigned int idx = 0; idx != input_dsets_.size(); idx++) {
    DataSet_1D const& ds = static_cast<DataSet_1D const&>( *(input_dsets_[idx]) );
    if (CheckSet(ds)) return 1; 
    // Calculate averages for each increment
    Darray avg;
    Iarray increments;
    int count = 0;
    int endpt = maxpts -1;
    double currentSum = 0.0;
    if (debug_ > 0) mprintf("DEBUG: Lambda %g\n", xval_[idx]);
    for (int pt = avg_skip_; pt != maxpts; pt++)
    {
      currentSum += ds.Dval(pt);
      count++;
      if (count == avg_increment_ || pt == endpt) {
        avg.push_back( currentSum / ((double)(pt - avg_skip_ + 1)) );
        increments.push_back(pt+1);
        if (debug_ > 0)
          mprintf("DEBUG:\t\tAvg from %i to %i: %g\n", avg_skip_+1, pt+1, avg.back());
        count = 0;
      }
    }
    if (sum.empty()) {
      sum.resize(avg.size());
      points = increments;
    } else if (sum.size() != avg.size()) {
      mprinterr("Error: Different # of increments for set '%s'; got %zu, expected %zu.\n",
                ds.legend(), avg.size(), sum.size());
      return 1;
    }
    // Create increment curve data sets
    if (curve_.empty()) {
      MetaData md(dAout_->Meta().Name(), "TIcurve");
      for (unsigned int j = 0; j != avg.size(); j++) {
        md.SetIdx( increments[j] );
        DataSet* ds = masterDSL_->AddSet(DataSet::XYMESH, md);
        if (ds == 0) return Analysis::ERR;
        ds->ModifyDim(Dimension::X).SetLabel("Lambda");
        ds->SetLegend( md.Name() + "_Skip" + integerToString(increments[j]) );
        if (curveout_ != 0) curveout_->AddDataSet( ds );
        curve_.push_back( ds );
      }
    }
    for (unsigned int j = 0; j != avg.size(); j++) {
      DataSet_Mesh& CR = static_cast<DataSet_Mesh&>( *(curve_[j]) );
      CR.AddXY(xval_[idx], avg[j]);
      if (mode_ == GAUSSIAN_QUAD)
        sum[j] += (wgt_[idx] * avg[j]);
    }
  } // END loop over data sets
  if (mode_ == TRAPEZOID) Integrate_Trapezoid(sum);
  // Store final integration values
  DataSet_Mesh& DA = static_cast<DataSet_Mesh&>( *dAout_ );
  DA.ModifyDim(Dimension::X).SetLabel("Point");
  for (unsigned int j = 0; j != points.size(); j++)
    DA.AddXY(points[j], sum[j]);

  return 0;
}