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

// Action_Outtraj::Init()
Action::RetType Action_Outtraj::Init(ArgList& actionArgs, ActionInit& init, int debugIn)
{
  // Set up output traj
  outtraj_.SetDebug(debugIn);
  std::string trajfilename = actionArgs.GetStringNext();
  if (trajfilename.empty()) {
    mprinterr("Error: No filename given.\nError: Usage: ");
    Help();
    return Action::ERR;
  }
  associatedParm_ = init.DSL().GetTopology(actionArgs);
  if (associatedParm_ == 0) {
    mprinterr("Error: Could not get associated topology for %s\n",trajfilename.c_str());
    return Action::ERR;
  }
  // If maxmin, get the name of the dataset as well as the max and min values.
  double lastmin = 0.0;
  double lastmax = 0.0;
  while ( actionArgs.Contains("maxmin") ) {
    std::string datasetName = actionArgs.GetStringKey("maxmin");
    if (!datasetName.empty()) {
      DataSet* dset = init.DSL().GetDataSet(datasetName);
      if (dset==0) {
        mprintf("Error: maxmin: Could not get dataset %s\n",datasetName.c_str());
        return Action::ERR;
      } else {
        // Currently only allow int, float, or double datasets
        if (dset->Type() != DataSet::INTEGER &&
            dset->Type() != DataSet::FLOAT &&
            dset->Type() != DataSet::DOUBLE) 
        {
          mprinterr("Error: maxmin: Only int, float, or double dataset (%s) supported.\n",
                  datasetName.c_str());
          return Action::ERR;
        }
        Dsets_.push_back( (DataSet_1D*)dset );
        Max_.push_back( actionArgs.getKeyDouble("max",lastmax) );
        Min_.push_back( actionArgs.getKeyDouble("min",lastmin) );
        lastmax = Max_.back();
        lastmin = Min_.back();
      }
    } else {
      mprinterr("Error: maxmin Usage: maxmin <setname> max <max> min <min>\n");
      return Action::ERR;
    }
  }
  // Initialize output trajectory with remaining arguments
  if ( outtraj_.InitEnsembleTrajWrite(trajfilename, actionArgs.RemainingArgs(), 
                                      TrajectoryFile::UNKNOWN_TRAJ, init.DSL().EnsembleNum()) ) 
    return Action::ERR;
  isSetup_ = false;

  mprintf("    OUTTRAJ: Writing frames associated with topology '%s'\n", associatedParm_->c_str());
  for (unsigned int ds = 0; ds < Dsets_.size(); ++ds)
    mprintf("\tmaxmin: Printing trajectory frames based on %g <= %s <= %g\n",
            Min_[ds], Dsets_[ds]->legend(), Max_[ds]);

  return Action::OK;
} 

/// Return true if set is an atomic type (i.e. int, double, float).
static bool GoodCalcType(DataSet const& ds) {
  if (ds.Type() == DataSet::DOUBLE || 
      ds.Type() == DataSet::FLOAT || 
      ds.Type() == DataSet::INT)
    return true;
  mprinterr("Error: DataSet %s is not a valid type for this calc.\n",
            ds.Name().c_str());
  return false;
}

TEST(DataSetCoreTest, DefaultsOk)
{
    DataSet dataset;
    EXPECT_EQ(DataSet::GENERIC, dataset.Type());
    EXPECT_FALSE(dataset.CreatedAt().empty());
    EXPECT_FALSE(dataset.MetaType().empty());
    EXPECT_FALSE(dataset.TimeStampedName().empty());
    EXPECT_FALSE(dataset.UniqueId().empty());

    EXPECT_EQ(0, dataset.TimeStampedName().find("pacbio_dataset_"));

    EXPECT_TRUE(dataset.Format().empty());
    EXPECT_TRUE(dataset.ModifiedAt().empty());
    EXPECT_TRUE(dataset.Name().empty());
    EXPECT_TRUE(dataset.ResourceId().empty());
    EXPECT_TRUE(dataset.Tags().empty());
    EXPECT_EQ(0, dataset.ExternalResources().Size());
    EXPECT_EQ(0, dataset.Filters().Size());
    EXPECT_EQ(0, dataset.SubDataSets().Size());

    EXPECT_EQ(string{"3.0.1"}, dataset.Version());
}

/** Replace all variables in given ArgList with their values. */
ArgList VariableArray::ReplaceVariables(ArgList const& argIn, DataSetList const& DSL, int debug)
{
  if (debug > 0) mprintf("DEBUG: Before variable replacement:  [%s]\n", argIn.ArgLine());
  ArgList modCmd = argIn;
  for (int n = 0; n < modCmd.Nargs(); n++) {
    size_t pos = modCmd[n].find("$");
    while (pos != std::string::npos) {
      // Argument is/contains a variable. Find first non-alphanumeric char
      size_t len = 1;
      for (size_t pos1 = pos+1; pos1 < modCmd[n].size(); pos1++, len++)
        if (!isalnum(modCmd[n][pos1])) break;
      std::string var_in_arg = modCmd[n].substr(pos, len);
      // See if variable occurs in CurrentVars_
      Varray::const_iterator vp = CurrentVars_.begin();
      for (; vp != CurrentVars_.end(); ++vp)
        if (vp->first == var_in_arg) break;
      // If found replace with value from CurrentVars_
      if (vp != CurrentVars_.end()) {
        if (debug > 0)
          mprintf("DEBUG: Replaced variable '%s' with value '%s'\n",
                  var_in_arg.c_str(), vp->second.c_str());
        std::string arg = modCmd[n];
        arg.replace(pos, vp->first.size(), vp->second);
        modCmd.ChangeArg(n, arg);
      } else {
        // Not found in CurrentVars_; see if this is a DataSet.
        for (size_t pos1 = pos+len; pos1 < modCmd[n].size(); pos1++, len++)
          if (!isalnum(modCmd[n][pos1]) &&
              modCmd[n][pos1] != '[' &&
              modCmd[n][pos1] != ':' &&
              modCmd[n][pos1] != ']' &&
              modCmd[n][pos1] != '_' &&
              modCmd[n][pos1] != '-' &&
              modCmd[n][pos1] != '%')
            break;
        var_in_arg = modCmd[n].substr(pos+1, len-1);
        DataSet* ds = DSL.GetDataSet( var_in_arg );
        if (ds == 0) {
          mprinterr("Error: Unrecognized variable in command: %s\n", var_in_arg.c_str());
          return ArgList();
        } else {
          if (ds->Type() != DataSet::STRING && ds->Group() != DataSet::SCALAR_1D) {
            mprinterr("Error: Only 1D data sets supported.\n");
            return ArgList();
          }
          if (ds->Size() < 1) {
            mprinterr("Error: Set is empty.\n");
            return ArgList();
          }
          if (ds->Size() > 1)
            mprintf("Warning: Only using first value.\n");
          std::string value;
          if (ds->Type() == DataSet::STRING)
            value = (*((DataSet_string*)ds))[0];
          else
            value = doubleToString(((DataSet_1D*)ds)->Dval(0));
          if (debug > 0)
            mprintf("DEBUG: Replaced variable '$%s' with value '%s' from DataSet '%s'\n",
                    var_in_arg.c_str(), value.c_str(), ds->legend());
          std::string arg = modCmd[n];
          arg.replace(pos, var_in_arg.size()+1, value);
          modCmd.ChangeArg(n, arg);
        }
      }
      pos = modCmd[n].find("$");
    } // END loop over this argument
  }
  return modCmd;
}

/** Given an array of already set up data sets (from e.g. input data files)
  * and optional X values, add the sets to this list if they do not exist
  * or append any existing sets.
  */
int DataSetList::AddOrAppendSets(std::string const& XlabelIn, Darray const& Xvals,
                                 DataListType const& Sets)
{
  if (debug_ > 0)
    mprintf("DEBUG: Calling AddOrAppendSets for %zu sets, %zu X values, Xlabel= %s.\n",
            Sets.size(), Xvals.size(), XlabelIn.c_str());
  if (Sets.empty()) return 0; // No error for now.
  // If no X label assume 'Frame' for backwards compat.
  std::string Xlabel;
  if (XlabelIn.empty())
    Xlabel.assign("Frame");
  else
    Xlabel = XlabelIn;
  Dimension Xdim;
  // First determine if X values increase monotonically with a regular step
  bool isMonotonic = true;
  double xstep = 1.0;
  if (Xvals.size() > 1) {
    xstep = (Xvals.back() - Xvals.front()) / (double)(Xvals.size() - 1);
    for (Darray::const_iterator X = Xvals.begin()+2; X != Xvals.end(); ++X)
      if ((*X - *(X-1)) - xstep > Constants::SMALL) {
        isMonotonic = false;
        break;
      }
    // Set dim even for non-monotonic sets so Label is correct. FIXME is this ok?
    Xdim = Dimension( Xvals.front(), xstep, Xlabel );
  } else
    // No X values. set generic X dim.
    Xdim = Dimension(1.0, 1.0, Xlabel);
  if (debug_ > 0) {
    mprintf("DEBUG: xstep %g xmin %g\n", Xdim.Step(), Xdim.Min());
    if (isMonotonic) mprintf("DEBUG: Xdim is monotonic.\n");
  }
  for (const_iterator ds = Sets.begin(); ds != Sets.end(); ++ds) {
    if (*ds == 0) continue; // TODO: Warn about blanks?
    if (debug_ > 0) mprintf("DEBUG: AddOrAppend set '%s'", (*ds)->legend());
    if (isMonotonic) (*ds)->SetDim(Dimension::X, Xdim);
    DataSet* existingSet = CheckForSet( (*ds)->Meta() );
    if (existingSet == 0) {
      // New set. If set is scalar 1D but X values are not monotonic,
      // convert to XY mesh if necessary.
      if ( !isMonotonic &&
           (*ds)->Group() == DataSet::SCALAR_1D &&
           (*ds)->Type() != DataSet::XYMESH )
      {
        DataSet* xyptr = AddSet( DataSet::XYMESH, (*ds)->Meta() );
        if (xyptr == 0) {
          mprinterr("Error: Could not convert set %s to XY mesh.\n", (*ds)->legend());
          continue; // TODO exit instead?
        }
        if ( (*ds)->Size() != Xvals.size() ) { // Sanity check
          mprinterr("Error: # of X values does not match set %s size.\n", (*ds)->legend());
          continue;
        }
        DataSet_1D const& set = static_cast<DataSet_1D const&>( *(*ds) );
        DataSet_Mesh& xy = static_cast<DataSet_Mesh&>( *xyptr );
        for (unsigned int i = 0; i != set.Size(); i++)
          xy.AddXY( Xvals[i], set.Dval(i) );
        xy.SetDim(Dimension::X, Xdim);
        if (debug_ > 0) mprintf(", New set, converted to XY-MESH\n");
        // Free memory since set has now been converted.
        delete *ds;
      } else { // No conversion. Just add.
        (*ds)->SetDim(Dimension::X, Xdim);
        AddSet( *ds );
        if (debug_ > 0) mprintf(", New set\n");
      }
    } else {
      if (debug_ > 0) mprintf(", appending to existing set\n");
      // Set exists. Try to append this set to existing set.
      bool canAppend = true;
      if ( (*ds)->Group() == DataSet::GENERIC ) {
        // For GENERIC sets, base Type must match. // TODO: Should this logic be in DataSets?
        if ( (*ds)->Type() != existingSet->Type() )
          canAppend = false;
      } else {
        // For non-GENERIC sets, Group must match
        if ( (*ds)->Group() != existingSet->Group() )
          canAppend = false;
      }
      if (!canAppend)
        mprinterr("Error: Cannot append set of type %s to set of type %s\n",
                  DataArray[(*ds)->Type()].Description,
                  DataArray[existingSet->Type()].Description);
      // If cannot append or attempt to append fails, rename and add as new set.
      if (!canAppend || existingSet->Append( *ds )) { // TODO Dimension check?
        if (canAppend)
          mprintf("Warning: Append currently not supported for type %s\n",
                  DataArray[existingSet->Type()].Description);
        MetaData md = (*ds)->Meta();
        md.SetName( GenerateDefaultName("X") );
        mprintf("Warning: Renaming %s to %s\n", (*ds)->Meta().PrintName().c_str(),
                md.PrintName().c_str());
        (*ds)->SetMeta( md );
       
        AddSet( *ds );
//.........这里部分代码省略.........

// Exec_UpdateParameters::Execute()
Exec::RetType Exec_UpdateParameters::Execute(CpptrajState& State, ArgList& argIn)
{
  std::string dsname = argIn.GetStringKey("setname");
  if (dsname.empty()) {
    mprinterr("Error: Specify parameter set.\n");
    return CpptrajState::ERR;
  }
  DataSet* ds = State.DSL().GetDataSet( dsname );
  if (ds == 0) {
    mprinterr("Error: Parameter data set '%s' not found.\n", dsname.c_str());
    return CpptrajState::ERR;
  }
  if (ds->Type() != DataSet::PARAMETERS) {
    mprinterr("Error: Set '%s' is not a parameter data set.\n", ds->legend());
    return CpptrajState::ERR;
  }
  DataSet_Parameters const& prm = static_cast<DataSet_Parameters const&>( *ds );
  Topology* dstop = State.DSL().GetTopology( argIn );
  if (dstop == 0) {
    mprinterr("Error: No topology specified.\n");
    return CpptrajState::ERR;
  }
  Topology& top = static_cast<Topology&>( *dstop );

  mprintf("\tUpdating parameters in topology '%s' using those in set '%s'\n",
          top.c_str(), prm.legend());
  // Get list of existing parameters.
  // We assume a parameter in topology is never repeated, so as soon
  // as it is found we can move to the next one. 
  ParmHolder<int> ParmIndices;
  // Bond parameters
  BondTypes(ParmIndices, top, top.Bonds());
  BondTypes(ParmIndices, top, top.BondsH());
  for (ParmHolder<BondParmType>::const_iterator it1 = prm.BP().begin();
                                                it1 != prm.BP().end(); ++it1)
  {
    for (ParmHolder<int>::const_iterator it0 = ParmIndices.begin();
                                         it0 != ParmIndices.end(); ++it0)
    {
      if (it1->first == it0->first)
      {
        BondParmType& bp = top.SetBondParm(it0->second);
        mprintf("\tUpdating bond parameter %s - %s from %f %f to %f %f\n",
                *(it0->first[0]), *(it0->first[1]), bp.Rk(), bp.Req(),
                it1->second.Rk(), it1->second.Req());
        bp = it1->second;
        break;
      }
    }
  }
  ParmIndices.clear();
  // Angle parameters
  AngleTypes(ParmIndices, top, top.Angles());
  AngleTypes(ParmIndices, top, top.AnglesH());
  for (ParmHolder<AngleParmType>::const_iterator it1 = prm.AP().begin();
                                                 it1 != prm.AP().end(); ++it1)
  {
    for (ParmHolder<int>::const_iterator it0 = ParmIndices.begin();
                                         it0 != ParmIndices.end(); ++it0)
    {
      if (it1->first == it0->first)
      {
        AngleParmType& bp = top.SetAngleParm(it0->second);
        mprintf("\tUpdating angle parameter %s - %s - %s from %f %f to %f %f\n",
                *(it0->first[0]), *(it0->first[1]), *(it0->first[2]), bp.Tk(), bp.Teq(),
                it1->second.Tk(), it1->second.Teq());
        bp = it1->second;
        break;
      }
    }
  }
  ParmIndices.clear();
  // Dihedral parameters
  DihedralTypes(ParmIndices, top, top.Dihedrals());
  DihedralTypes(ParmIndices, top, top.DihedralsH());
  for (ParmHolder<DihedralParmType>::const_iterator it1 = prm.DP().begin();
                                                    it1 != prm.DP().end(); ++it1)
  {
    for (ParmHolder<int>::const_iterator it0 = ParmIndices.begin();
                                         it0 != ParmIndices.end(); ++it0)
    {
      if (it1->first == it0->first)
      {
        DihedralParmType& bp = top.SetDihedralParm(it0->second);
        mprintf("\tUpdating dihedral parameter %s - %s - %s  %s from %f %f %f to %f %f %f\n",
                *(it0->first[0]), *(it0->first[1]), *(it0->first[2]), *(it0->first[3]),
                bp.Pk(), bp.Pn(), bp.Phase(),
                it1->second.Pk(), it1->second.Pn(), it1->second.Phase());
        bp = it1->second;
        break;
      }
    }
  }
  ParmIndices.clear();
  if (top.Chamber().HasChamber()) {
    ChamberParmType& chm = top.SetChamber();
    // UB parameters
    BondTypes(ParmIndices, top, chm.UB());
    for (ParmHolder<BondParmType>::const_iterator it1 = prm.UB().begin();
//.........这里部分代码省略.........

/** Calculate time correlation between two DataSets.
  * \D1 DataSet to calculate correlation for.
  * \D2 DataSet to calculate correlation to.
  * \Ct DataSet to store time correlation fn, must be DOUBLE.
  * \lagmaxIn Max lag to calculate corr. -1 means use size of dataset.
  * \calccovar If true calculate covariance (devation from avg).
  * \return 0 on success, 1 on error.
  */
int DS_Math::CrossCorr( DataSet& D1, DataSet& D2, DataSet& Ct, int lagmaxIn, 
                        bool calccovar, bool usefft )
{
  int lagmax;
  double ct;
  // Check if D1 and D2 are valid types
  if ( !GoodCalcType(D1) ) return 1;
  if ( !GoodCalcType(D2) ) return 1;
  // Check that D1 and D2 have same # data points.
  int Nelements = D1.Size();
  if (Nelements != D2.Size()) {
    mprinterr("Error: CrossCorr: # elements in dataset %s (%i) not equal to\n", 
              D1.Legend().c_str(), Nelements);
    mprinterr("Error:            # elements in dataset %s (%i)\n", 
              D2.Legend().c_str(), D2.Size());
    return 1;
  }
  if (Nelements < 2) {
    mprinterr("Error: CrossCorr: # elements is less than 2 (%i)\n", Nelements);
    return 1;
  }
  // Check return dataset type
  if ( Ct.Type() != DataSet::DOUBLE ) {
    mprinterr("Internal Error: CrossCorr: Ct must be of type DataSet::DOUBLE.\n");
    return 1;
  }
  // Check if lagmaxIn makes sense. Set default lag to be Nelements 
  // if not specified.
  if (lagmaxIn == -1)
    lagmax = Nelements;
  else if (lagmaxIn > Nelements) {
    mprintf("Warning: CrossCorr [%s][%s]: max lag (%i) > Nelements (%i), setting to Nelements.\n",
            D1.Legend().c_str(), D2.Legend().c_str(), lagmaxIn, Nelements);
    lagmax = Nelements;
  } else
    lagmax = lagmaxIn;
  // If calculating covariance calculate averages
  double avg1 = 0;
  double avg2 = 0;
  if ( calccovar ) {
    avg1 = Avg(D1);
    avg2 = Avg(D2);
  }
  // Calculate correlation
  double norm = 1.0;
  if ( usefft ) {
    // Calc using FFT
    CorrF_FFT pubfft1(Nelements);
    ComplexArray data1 = pubfft1.Array();
    data1.PadWithZero(Nelements);
    for (int i = 0; i < Nelements; ++i)
      data1[i*2] = D1.Dval(i) - avg1;
    if (&D2 == &D1)
      pubfft1.AutoCorr(data1);
    else {
      // Populate second dataset if different
      ComplexArray data2 = pubfft1.Array();
      data2.PadWithZero(Nelements);
      for (int i = 0; i < Nelements; ++i)
        data2[i*2] = D2.Dval(i) - avg2;
      pubfft1.CrossCorr(data1, data2);
    }
    // Put real components of data1 in output DataSet
    norm = 1.0 / fabs( data1[0] );
    for (int i = 0; i < lagmax; ++i) {
      ct = data1[i*2] * norm;
      Ct.Add(i, &ct);
    }
  } else {
    // Direct calc
    for (int lag = 0; lag < lagmax; ++lag) {
      ct = 0;
      int jmax = Nelements - lag;
      for (int j = 0; j < jmax; ++j)
        ct += ((D1.Dval(j) - avg1) * (D2.Dval(j+lag) - avg2));
      if (lag == 0) {
        if (ct != 0)
          norm = fabs( ct );
      }
      ct /= norm;
      Ct.Add(lag, &ct);
    }
  }
  return 0;
}

/** Syntax: dataset invert <set arg0> ... name <new name> */
Exec::RetType Exec_DataSetCmd::InvertSets(CpptrajState& State, ArgList& argIn) {
  DataSetList& DSL = State.DSL();
  // Get keywords
  DataSet* inputNames = 0;
  std::string dsname = argIn.GetStringKey("legendset");
  if (!dsname.empty()) {
    inputNames = DSL.GetDataSet( dsname );
    if (inputNames == 0) {
      mprinterr("Error: Name set '%s' not found.\n", dsname.c_str());
      return CpptrajState::ERR;
    }
    if (inputNames->Type() != DataSet::STRING) {
      mprinterr("Error: Set '%s' does not contain strings.\n", inputNames->legend());
      return CpptrajState::ERR;
    }
    mprintf("\tUsing names from set '%s' as legends for inverted sets.\n", inputNames->legend());
  }
  dsname = argIn.GetStringKey("name");
  if (dsname.empty()) {
    mprinterr("Error: 'invert' requires that 'name <new set name>' be specified.\n");
    return CpptrajState::ERR;
  }
  mprintf("\tNew sets will be named '%s'\n", dsname.c_str());
  DataFile* outfile = State.DFL().AddDataFile( argIn.GetStringKey("out"), argIn );
  if (outfile != 0)
    mprintf("\tNew sets will be output to '%s'\n", outfile->DataFilename().full());
  // TODO determine type some other way
  DataSet::DataType outtype = DataSet::DOUBLE;
  // Get input DataSets
  std::vector<DataSet_1D*> input_sets; 
  std::string dsarg = argIn.GetStringNext();
  while (!dsarg.empty()) {
    DataSetList sets = DSL.GetMultipleSets( dsarg );
    for (DataSetList::const_iterator ds = sets.begin(); ds != sets.end(); ++ds)
    {
      if ( (*ds)->Group() != DataSet::SCALAR_1D ) {
        mprintf("Warning: '%s': Inversion only supported for 1D scalar data sets.\n",
                (*ds)->legend());
      } else {
        if (!input_sets.empty()) {
          if ( (*ds)->Size() != input_sets.back()->Size() ) {
            mprinterr("Error: Set '%s' has different size (%zu) than previous set (%zu)\n",
                      (*ds)->legend(), (*ds)->Size(), input_sets.back()->Size());
            return CpptrajState::ERR;
          }
        }
        input_sets.push_back( (DataSet_1D*)*ds );
      }
    }
    dsarg = argIn.GetStringNext();
  }
  if (input_sets.empty()) {
    mprinterr("Error: No sets selected.\n");
    return CpptrajState::ERR;
  }
  if (inputNames != 0 && inputNames->Size() != input_sets.front()->Size()) {
    mprinterr("Error: Name set '%s' size (%zu) differs from # data points (%zu).\n",
              inputNames->legend(), inputNames->Size(), input_sets.front()->Size());
    return CpptrajState::ERR;
  }
  mprintf("\t%zu input sets; creating %zu output sets.\n",
          input_sets.size(), input_sets.front()->Size());
  // Need an output data set for each point in input sets
  std::vector<DataSet*> output_sets;
  int column = 1;
  for (int idx = 0; idx != (int)input_sets[0]->Size(); idx++, column++) {
    DataSet* ds = 0;
    ds = DSL.AddSet(outtype, MetaData(dsname, column));
    if (ds == 0) return CpptrajState::ERR;
    if (inputNames != 0)
      ds->SetLegend( (*((DataSet_string*)inputNames))[idx] );
    output_sets.push_back( ds );
    if (outfile != 0) outfile->AddDataSet( ds );
  }
  // Create a data set containing names of each input data set
  DataSet* nameset = DSL.AddSet(DataSet::STRING, MetaData(dsname, column));
  if (nameset == 0) return CpptrajState::ERR;
  if (inputNames != 0)
    nameset->SetLegend("Names");
  if (outfile != 0) outfile->AddDataSet( nameset );
  // Populate output data sets
  for (int jdx = 0; jdx != (int)input_sets.size(); jdx++)
  {
    DataSet_1D const& INP = static_cast<DataSet_1D const&>( *(input_sets[jdx]) );
    nameset->Add( jdx, INP.legend() );
    for (unsigned int idx = 0; idx != INP.Size(); idx++)
    {
      double dval = INP.Dval( idx );
      output_sets[idx]->Add( jdx, &dval );
    }
  }

  return CpptrajState::OK;
}

// Action_CreateReservoir::Init()
Action::RetType Action_CreateReservoir::Init(ArgList& actionArgs, ActionInit& init, int debugIn)
{
# ifndef BINTRAJ
  mprinterr("Error: NetCDF reservoir requires NetCDF support. Recompile with -DBINTRAJ.\n");
  return Action::ERR;
# endif
# ifdef MPI
  if (init.TrajComm().Size() > 1) {
    mprinterr("Error: 'createreservoir' action does not work with > 1 process (%i processes currently).\n", init.TrajComm().Size());
    return Action::ERR;
  }
# endif
  // Get keywords
  filename_.SetFileName( actionArgs.GetStringNext() );
  if (filename_.empty()) {
    mprinterr("Error: createreservoir: No filename specified.\n");
    return Action::ERR;
  }
  reservoirT_ = actionArgs.getKeyDouble("temp0", -1.0);
  if (reservoirT_ < 0.0) {
    mprinterr("Error: Reservoir temperature must be specified and cannot be < 0.0\n");
    return Action::ERR;
  }
  iseed_ = actionArgs.getKeyInt("iseed", 0);
  if (iseed_ < 1) {
    mprinterr("Error: Reservoir random seed must be specified and > 0\n");
    return Action::ERR;
  }
  useVelocity_ = !actionArgs.hasKey("novelocity");
  useForce_ = !actionArgs.hasKey("noforce");
  // Get parm for reservoir traj
  original_trajparm_ = init.DSL().GetTopology( actionArgs );
  if (original_trajparm_ == 0) {
    mprinterr("Error: createreservoir: no topology.\n");
    return Action::ERR;
  }
  // Get energy data set
  std::string eneDsname = actionArgs.GetStringKey("ene");
  DataSet* dstmp = init.DSL().GetDataSet( eneDsname );
  if (dstmp == 0) {
    mprinterr("Error: could not get energy data set %s\n", eneDsname.c_str());
    return Action::ERR;
  }
  if (dstmp->Type() != DataSet::FLOAT &&
      dstmp->Type() != DataSet::DOUBLE &&
      dstmp->Type() != DataSet::XYMESH)
  {
    mprinterr("Error: energy data set %s must be type FLOAT, DOUBLE, or XYMESH.\n",
              dstmp->legend());
    return Action::ERR;
  }
  if (dstmp->Ndim() != 1) {
    mprinterr("Error: energy data set is not 1D (%zu)\n", dstmp->Ndim());
    return Action::ERR;
  }
  ene_ = static_cast<DataSet_1D*>( dstmp );
  // Get bin data set
  std::string binDSname = actionArgs.GetStringKey("bin");
  if (!binDSname.empty()) {
    dstmp = init.DSL().GetDataSet( binDSname );
    if (dstmp == 0) {
      mprinterr("Error: could not get bin data set %s\n", binDSname.c_str());
      return Action::ERR;
    } else if (dstmp->Ndim() != 1) {
      mprinterr("Error: bin data set must be one dimensional.\n");
      return Action::ERR;
    }
    bin_ = static_cast<DataSet_1D*>( dstmp );
  }
  trajIsOpen_ = false;
  nframes_ = 0;
  // Setup output reservoir file
  reservoir_.SetDebug( debugIn );
  // Set title
  title_ = actionArgs.GetStringKey("title");
  if (title_.empty())
    title_.assign("Cpptraj Generated structure reservoir");

  mprintf("    CREATERESERVOIR: '%s', energy data '%s'", filename_.full(), ene_->legend());
  if (bin_ != 0)
    mprintf(", bin data '%s'", bin_->legend());
  mprintf("\n\tTitle: %s\n", title_.c_str());
  mprintf("\tReservoir temperature= %.2f, random seed= %i\n", reservoirT_, iseed_);
  if (useVelocity_)
    mprintf("\tVelocities will be written to reservoir if present.\n");
  else
    mprintf("\tVelocities will not be written to reservoir.\n");
  if (useForce_)
    mprintf("\tForces will be written to reservoir if present.\n");
  else
    mprintf("\tForces will not be written to reservoir.\n");
  mprintf("\tTopology: %s\n", original_trajparm_->c_str());
  return Action::OK;
}

// DataIO_Std::WriteDataNormal()
int DataIO_Std::WriteDataNormal(CpptrajFile& file, DataSetList const& Sets) {
  // Assume all 1D data sets.
  if (Sets.empty() || CheckAllDims(Sets, 1)) return 1;
  // For this output to work the X-dimension of all sets needs to match.
  // The most important things for output are min and step so just check that.
  // Use X dimension of set 0 for all set dimensions.
  CheckXDimension( Sets );
  // TODO: Check for empty dim.
  DataSet* Xdata = Sets[0];
  Dimension const& Xdim = static_cast<Dimension const&>( Xdata->Dim(0) );
  int xcol_width = Xdim.Label().size() + 1; // Only used if hasXcolumn_
  if (xcol_width < 8) xcol_width = 8;
  int xcol_precision = 3;

  // Determine size of largest DataSet.
  size_t maxFrames = DetermineMax( Sets );

  // Set up X column.
  TextFormat x_col_format(XcolFmt());
  if (hasXcolumn_) {
    if (XcolPrecSet()) {
      xcol_width = XcolWidth();
      x_col_format = TextFormat(XcolFmt(), XcolWidth(), XcolPrec());
    } else {
      // Create format string for X column based on dimension in first data set.
      // Adjust X col precision as follows: if the step is set and has a 
      // fractional component set the X col width/precision to either the data
      // width/precision or the current width/precision, whichever is larger. If
      // the set is XYMESH but step has not been set (so we do not know spacing 
      // between X values) use default precision. Otherwise the step has no
      // fractional component so make the precision zero.
      double step_i;
      double step_f = modf( Xdim.Step(), &step_i );
      double min_f  = modf( Xdim.Min(),  &step_i );
      if (Xdim.Step() > 0.0 && (step_f > 0.0 || min_f > 0.0)) {
        xcol_precision = std::max(xcol_precision, Xdata->Format().Precision());
        xcol_width = std::max(xcol_width, Xdata->Format().Width());
      } else if (Xdata->Type() != DataSet::XYMESH)
        xcol_precision = 0;
      x_col_format.SetCoordFormat( maxFrames, Xdim.Min(), Xdim.Step(), xcol_width, xcol_precision );
    }
  } else {
    // If not writing an X-column, no leading space for the first dataset.
    Xdata->SetupFormat().SetFormatAlign( TextFormat::RIGHT );
  }

  // Write header to buffer
  std::vector<int> Original_Column_Widths;
  if (writeHeader_) {
    // If x-column present, write x-label
    if (hasXcolumn_)
      WriteNameToBuffer( file, Xdim.Label(), xcol_width, true );
    // To prevent truncation of DataSet legends, adjust the width of each
    // DataSet if necessary.
    for (DataSetList::const_iterator ds = Sets.begin(); ds != Sets.end(); ++ds) {
      // Record original column widths in case they are changed.
      Original_Column_Widths.push_back( (*ds)->Format().Width() );
      int colLabelSize;
      if (ds == Sets.begin() && !hasXcolumn_)
        colLabelSize = (int)(*ds)->Meta().Legend().size() + 1;
      else
        colLabelSize = (int)(*ds)->Meta().Legend().size();
      //mprintf("DEBUG: Set '%s', fmt width= %i, colWidth= %i, colLabelSize= %i\n",
      //        (*ds)->legend(), (*ds)->Format().Width(), (*ds)->Format().ColumnWidth(),
      //        colLabelSize);
      if (colLabelSize >= (*ds)->Format().ColumnWidth())
        (*ds)->SetupFormat().SetFormatWidth( colLabelSize );
    }
    // Write dataset names to header, left-aligning first set if no X-column
    DataSetList::const_iterator set = Sets.begin();
    if (!hasXcolumn_)
      WriteNameToBuffer( file, (*set)->Meta().Legend(), (*set)->Format().ColumnWidth(), true  );
    else
      WriteNameToBuffer( file, (*set)->Meta().Legend(), (*set)->Format().ColumnWidth(), false );
    ++set;
    for (; set != Sets.end(); ++set) 
      WriteNameToBuffer( file, (*set)->Meta().Legend(), (*set)->Format().ColumnWidth(), false );
    file.Printf("\n"); 
  }

  // Write Data
  DataSet::SizeArray positions(1);
  for (positions[0] = 0; positions[0] < maxFrames; positions[0]++) {
    // Output Frame for each set
    if (hasXcolumn_)
      file.Printf( x_col_format.fmt(), Xdata->Coord(0, positions[0]) );
    for (DataSetList::const_iterator set = Sets.begin(); set != Sets.end(); ++set) 
      (*set)->WriteBuffer(file, positions);
    file.Printf("\n"); 
  }
  // Restore original column widths if necessary
  if (!Original_Column_Widths.empty())
    for (unsigned int i = 0; i != Original_Column_Widths.size(); i++)
      Sets[i]->SetupFormat().SetFormatWidth( Original_Column_Widths[i] );
  return 0;
}