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

Analysis::RetType Analysis_AutoCorr::Setup(ArgList& analyzeArgs, DataSetList* datasetlist,
                            TopologyList* PFLin, DataFileList* DFLin, int debugIn)
{
  const char* calctype;

  std::string setname = analyzeArgs.GetStringKey("name");
  DataFile* outfile = DFLin->AddDataFile( analyzeArgs.GetStringKey("out"), analyzeArgs );
  lagmax_ = analyzeArgs.getKeyInt("lagmax",-1);
  calc_covar_ = !analyzeArgs.hasKey("nocovar");
  usefft_ = !analyzeArgs.hasKey("direct");
  // Select datasets from remaining args
  ArgList dsetArgs = analyzeArgs.RemainingArgs();
  for (ArgList::const_iterator dsa = dsetArgs.begin(); dsa != dsetArgs.end(); ++dsa)
    dsets_ += datasetlist->GetMultipleSets( *dsa );
  if (dsets_.empty()) {
    mprinterr("Error: autocorr: No data sets selected.\n");
    return Analysis::ERR;
  }
  // If setname is empty generate a default name
  if (setname.empty())
    setname = datasetlist->GenerateDefaultName( "autocorr" );
  // Setup output datasets
  int idx = 0;
  MetaData md( setname );
  for (DataSetList::const_iterator DS = dsets_.begin(); DS != dsets_.end(); ++DS) {
    md.SetIdx( idx++ );
    DataSet* dsout = datasetlist->AddSet( DataSet::DOUBLE, md );
    if (dsout==0) return Analysis::ERR;
    dsout->SetLegend( (*DS)->Meta().Legend() );
    outputData_.push_back( dsout );
    // Add set to output file
    if (outfile != 0) outfile->AddDataSet( outputData_.back() );
  }
 
  if (calc_covar_)
    calctype = "covariance";
  else
    calctype = "correlation";
 
  mprintf("    AUTOCORR: Calculating auto-%s for %i data sets:\n", calctype, dsets_.size());
  dsets_.List();
  if (lagmax_!=-1)
    mprintf("\tLag max= %i\n", lagmax_);
  if ( !setname.empty() )
    mprintf("\tSet name: %s\n", setname.c_str() );
  if ( outfile != 0 )
    mprintf("\tOutfile name: %s\n", outfile->DataFilename().base());
  if (usefft_)
    mprintf("\tUsing FFT to calculate %s.\n", calctype);
  else
    mprintf("\tUsing direct method to calculate %s.\n", calctype);

  return Analysis::OK;
}

// Exec_DataSetCmd::Execute()
Exec::RetType Exec_DataSetCmd::Execute(CpptrajState& State, ArgList& argIn) {
  RetType err = CpptrajState::OK;
  if (argIn.Contains("legend")) {         // Set legend for one data set
    std::string legend = argIn.GetStringKey("legend");
    DataSet* ds = State.DSL().GetDataSet( argIn.GetStringNext() );
    if (ds == 0) return CpptrajState::ERR;
    mprintf("\tChanging legend '%s' to '%s'\n", ds->legend(), legend.c_str());
    ds->SetLegend( legend );
  // ---------------------------------------------
  } else if (argIn.hasKey("outformat")) { // Change double precision set output format
    err = ChangeOutputFormat(State, argIn);
  // ---------------------------------------------
  } else if (argIn.hasKey("remove")) {    // Remove data sets by various criteria
    err = Remove(State, argIn);
  // ---------------------------------------------
  } else if (argIn.hasKey("makexy")) {    // Combine values from two sets into 1
    err = MakeXY(State, argIn);
  // ---------------------------------------------
  } else if (argIn.hasKey("make2d")) {    // Create 2D matrix from 1D set
    err = Make2D(State, argIn);
  // ---------------------------------------------
  } else if (argIn.hasKey("vectorcoord")) { // Extract vector X/Y/Z coord as new set
    err = VectorCoord(State, argIn);
  // ---------------------------------------------
  } else if (argIn.hasKey("filter")) {    // Filter points in data set to make new data set
    err = Filter(State, argIn);
  // ---------------------------------------------
  } else if (argIn.hasKey("cat")) {       // Concatenate two or more data sets
    err = Concatenate(State, argIn);
  // ---------------------------------------------
  } else if (argIn.hasKey("droppoints")) { // Drop points from set
    err = ModifyPoints(State, argIn, true);
  // ---------------------------------------------
  } else if (argIn.hasKey("keeppoints")) { // Keep points in set
    err = ModifyPoints(State, argIn, false);
  // ---------------------------------------------
  } else if (argIn.hasKey("dim")) {        // Modify dimension of set(s)
    err = ChangeDim(State, argIn);
  // ---------------------------------------------
  } else if (argIn.hasKey("invert")) {     // Invert set(s) X/Y, create new sets
    err = InvertSets(State, argIn);
  // ---------------------------------------------
  } else {                                // Default: change mode/type for one or more sets.
    err = ChangeModeType(State, argIn);
  }
  return err;
}

// Analysis_FFT::Setup()
Analysis::RetType Analysis_FFT::Setup(ArgList& analyzeArgs, DataSetList* datasetlist, DataFileList* DFLin, int debugIn)
{
  std::string setname = analyzeArgs.GetStringKey("name");
  DataFile* outfile = DFLin->AddDataFile(analyzeArgs.GetStringKey("out"), analyzeArgs);
  dt_ = analyzeArgs.getKeyDouble("dt",1.0);
  // Select datasets from remaining args
  if (input_dsets_.AddSetsFromArgs( analyzeArgs.RemainingArgs(), *datasetlist )) {
    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 setname is empty generate a default name
  if (setname.empty())
    setname = datasetlist->GenerateDefaultName( "FFT" );
  // Setup output datasets.
  int idx = 0;
  if ( input_dsets_.size() == 1 )
    idx = -1; // Only one input set, no need to refer to it by index
  for ( Array1D::const_iterator DS = input_dsets_.begin(); 
                                DS != input_dsets_.end(); ++DS) 
  {
    DataSet* dsout = datasetlist->AddSet( DataSet::DOUBLE, MetaData(setname, idx++) );
    if (dsout==0) return Analysis::ERR;
    dsout->SetLegend( (*DS)->Meta().Legend() );
    output_dsets_.push_back( (DataSet_1D*)dsout );
    if (outfile != 0) outfile->AddDataSet( dsout );
  }

  mprintf("    FFT: Calculating FFT for %u data sets.\n", input_dsets_.size());
  mprintf("\tTime step: %f\n", dt_);
  if ( !setname.empty() )
    mprintf("\tSet name: %s\n", setname.c_str() );
  if ( outfile != 0 )
    mprintf("\tOutfile name: %s\n", outfile->DataFilename().base());

  return Analysis::OK;
}

Analysis::RetType Analysis_Integrate::Setup(ArgList& analyzeArgs, AnalysisSetup& setup, int debugIn)
{
  std::string setname = analyzeArgs.GetStringKey("name");
  outfile_ = setup.DFL().AddDataFile(analyzeArgs.GetStringKey("out"), 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;
  }

  // Set up output datasets
  if (outfile_ != 0) {
    for (Array1D::const_iterator dsIn = input_dsets_.begin();
                                 dsIn != input_dsets_.end(); ++dsIn)
    {
      DataSet* ds = setup.DSL().AddSet(DataSet::XYMESH, setname, "Int");
      if (ds == 0) return Analysis::ERR;
      ds->SetLegend( "Int(" + (*dsIn)->Meta().Legend() + ")" );
      outfile_->AddDataSet( ds );
      output_dsets_.push_back( (DataSet_Mesh*)ds );
    }
  }
  
  mprintf("    INTEGRATE: Calculating integral of %i data sets.\n",
          input_dsets_.size());
  if (outfile_ != 0) {
    if (!setname.empty())
      mprintf("\tOutput set name: %s\n", setname.c_str());
    mprintf("\tOutfile name: %s\n", outfile_->DataFilename().base());
  }
  //for (Array1D::const_iterator set = input_dsets_.begin(); set != input_dsets_.end(); ++set)
  //  mprintf("\t%s\n", (*set)->legend());
  return Analysis::OK;
}

/** Perform setup required for per residue rmsd calculation.
  * Need to set up a target mask, reference mask, and dataset for each
  * residue specified in ResRange.
  * NOTE: Residues in the range arguments from user start at 1, internal
  *       res nums start from 0.
  */
int Action_Rmsd::perResSetup(Topology* currentParm, Topology* RefParm) {
  Range tgt_range;
  Range ref_range;

  NumResidues_ = currentParm->FinalSoluteRes();

  // If no target range previously specified do all solute residues
  if (ResRange_.Empty()) 
    tgt_range.SetRange(1, NumResidues_+1);
  else
    tgt_range = ResRange_;

  // If the reference range is empty, set it to match the target range
  if (RefRange_.Empty()) 
    ref_range = tgt_range;
  else
    ref_range = RefRange_;

  // Check that the number of reference residues matches number of target residues
  if (tgt_range.Size() != ref_range.Size()) {
    mprintf("Warning: RMSD: perres: Number of residues %i does not match\n",tgt_range.Size());
    mprintf("Warning:       number of reference residues %i.\n",ref_range.Size());
    return 1;
  }

  // Setup a dataset, target mask, and reference mask, for each residue.
  // Since we will only calculate per res rmsd for residues that can be
  // successfully set up, keep track of that as well.
  //mprinterr("DEBUG: Setting up %i masks and data for %s\n",nres,currentParm->parmName);
  resizeResMasks();
  resIsActive_.reserve(NumResidues_);
  resIsActive_.assign(NumResidues_, false);
  int N = -1; // Set to -1 since increment is at top of loop
  Range::const_iterator ref_it = ref_range.begin();
  for (Range::const_iterator tgt_it = tgt_range.begin();
                             tgt_it != tgt_range.end(); ++tgt_it)
  {
    int tgtRes = *tgt_it;
    int refRes = *ref_it;
    ++ref_it;
    // Check if either the residue num or the reference residue num out of range.
    if ( tgtRes < 1 || tgtRes > NumResidues_) {
      mprintf("Warning: Rmsd: perres: Specified residue # %i is out of range.\n",
              tgtRes);
      continue;
    }
    if ( refRes < 1 || refRes > NumResidues_ ) {
      mprintf("Warning: Rmsd: perres: Specified reference residue # %i is out of range.\n",
              refRes);
      continue;
    }
    ++N;
    // Create dataset for res - if already present this returns null 
    DataSet* prDataSet = masterDSL_->AddSetIdxAspect( DataSet::DOUBLE, rmsd_->Name(), tgtRes, "res");
    prDataSet->SetLegend( currentParm->TruncResNameNum(tgtRes-1) );
    PerResRMSD_.push_back( prDataSet );

    // Setup mask strings. Note that masks are based off user residue nums
    std::string tgtArg = ":" + integerToString(tgtRes) + perresmask_;
    tgtResMask_[N].SetMaskString(tgtArg.c_str());
    std::string refArg = ":" + integerToString(refRes) + perresmask_;
    refResMask_[N].SetMaskString(refArg.c_str());
    //mprintf("DEBUG: RMSD: PerRes: Mask %s RefMask %s\n",tgtArg,refArg);

    // Setup the reference mask
    if (RefParm->SetupIntegerMask(refResMask_[N])) {
      mprintf("      perres: Could not setup reference mask for residue %i\n",refRes);
      continue;
    }
    if (refResMask_[N].None()) {
      mprintf("      perres: No atoms selected for reference residue %i\n",refRes);
      continue;
    }

    // Setup the target mask
    if (currentParm->SetupIntegerMask(tgtResMask_[N])) {
      mprintf("      perres: Could not setup target mask for residue %i\n",tgtRes);
      continue;
    }
    if (tgtResMask_[N].None()) {
      mprintf("      perres: No atoms selected for target residue %i\n",tgtRes);
      continue;
    }

    // Check that # atoms in target and reference masks match
    if (tgtResMask_[N].Nselected() != refResMask_[N].Nselected()) {
      mprintf("      perres: Res %i: # atoms in Tgt [%i] != # atoms in Ref [%i]\n",
              tgtRes,tgtResMask_[N].Nselected(),refResMask_[N].Nselected());
      continue;
    }

    // Indicate that these masks were properly set up
    resIsActive_[N]=true;
  }   
//.........这里部分代码省略.........

//.........这里部分代码省略.........
                 *ba < setup.Top().Res(*res).LastAtom() )
            {
              if ( !selected[ *ba - res_first_atom ] &&
                   setup.Top()[ *ba ].Element() == Atom::HYDROGEN )
                heavyAtomGroup.push_back( *ba );
            }
          if (!heavyAtomGroup.empty())
            potentialSites.push_back( Site(*res, heavyAtomGroup) );
        }
      }
    }
    mprintf("\t%zu potential NOE sites:\n", potentialSites.size());
    for (SiteArray::const_iterator site = potentialSites.begin();
                                   site != potentialSites.end(); ++site)
    {
      mprintf("  %u\tRes %i:", site - potentialSites.begin(), site->ResNum()+1);
      for (unsigned int idx = 0; idx != site->Nindices(); ++idx)
        mprintf(" %s", setup.Top().TruncAtomNameNum( site->Idx(idx) ).c_str());
      mprintf("\n");
    }
    if (noeArray_.empty()) {
      size_t siteArraySize = 0;
      // Set up all potential NOE pairs. Keep track of size.
      for (SiteArray::const_iterator site1 = potentialSites.begin();
                                     site1 != potentialSites.end(); ++site1)
      {
        for (SiteArray::const_iterator site2 = site1 + 1;
                                       site2 != potentialSites.end(); ++site2)
        {
          if (site1->ResNum() != site2->ResNum()) {
            std::string legend = site1->SiteLegend(setup.Top()) + "--" +
                                 site2->SiteLegend(setup.Top());
            DataSet* ds = 0;
            if (series_) {
              ds = masterDSL_->AddSet(DataSet::FLOAT,
                                      MetaData(setname_, "foundNOE", noeArray_.size()));
              if (ds == 0) return Action::ERR;
              // Construct a data set name.
              ds->SetLegend(legend);
            }
            noeArray_.push_back( NOEtype(*site1, *site2, ds, legend) );
            siteArraySize += (2 * sizeof(int) * site1->Nindices()) +
                             (2 * sizeof(int) * site2->Nindices());
          }
        }
      }
      numNoePairs_ = noeArray_.size();
      size_t siteSize = sizeof(int) + (2 * sizeof(Iarray)) + sizeof(Site);
      size_t noeSize = (2 * siteSize) + sizeof(DataSet*) + sizeof(double) 
                       + sizeof(NOEtype);
      if (series_) noeSize += sizeof(std::vector<float>);
      size_t noeArraySize = (noeSize * numNoePairs_) + siteArraySize;
      if (series_)
        noeArraySize += (setup.Nframes() * numNoePairs_ * sizeof(float));
      mprintf("\t%zu potential NOE pairs. Estimated memory usage is %s\n",
              numNoePairs_, ByteString(noeArraySize, BYTE_DECIMAL).c_str());
    } else if (numNoePairs_ != potentialSites.size()) {
      mprinterr("Warning: Found NOE matrix has already been set up for %zu potential\n"
                "Warning:   NOEs, but %zu NOEs currently found.\n", numNoePairs_,
                potentialSites.size());
      return Action::SKIP;
    }
  }
  // ---------------------------------------------
  // Set up NOEs specified on the command line
  if (!Pairs_.empty()) {
    if (!specifiedNOEs_.empty()) {
      mprintf("Warning: Specifying NOEs currently only works with first topology used.\n");
      return Action::SKIP;
    }
    for (MaskPairArray::iterator mp = Pairs_.begin(); mp != Pairs_.end(); mp++) {
      if (setup.Top().SetupIntegerMask( mp->first )) return Action::ERR;
      int res1 = CheckSameResidue(setup.Top(), mp->first);
      if (res1 < 0) continue;
      if (setup.Top().SetupIntegerMask( mp->second )) return Action::ERR;
      int res2 = CheckSameResidue(setup.Top(), mp->second);
      if (res2 < 0) continue;
      Site site1( res1, mp->first.Selected() );
      Site site2( res2, mp->second.Selected() );
      std::string legend = site1.SiteLegend(setup.Top()) + "--" +
                           site2.SiteLegend(setup.Top());
      DataSet* ds = 0;
      if (series_) {
        ds = masterDSL_->AddSet(DataSet::FLOAT,
                                MetaData(setname_, "specNOE", specifiedNOEs_.size()));
        if (ds == 0) return Action::ERR;
        ds->SetLegend(legend);
      }
      specifiedNOEs_.push_back( NOEtype(site1, site2, ds, legend) );
    }
  } 
  // Set up imaging info for this parm
  Image_.SetupImaging( setup.CoordInfo().TrajBox().Type() );
  if (Image_.ImagingEnabled())
    mprintf("\tImaged.\n");
  else
    mprintf("\tImaging off.\n");

  return Action::OK;  
}

Analysis::RetType Analysis_AutoCorr::Setup(ArgList& analyzeArgs, AnalysisSetup& setup, int debugIn)
{
  const char* calctype;

  std::string setname = analyzeArgs.GetStringKey("name");
  DataFile* outfile = setup.DFL().AddDataFile( analyzeArgs.GetStringKey("out"), analyzeArgs );
  lagmax_ = analyzeArgs.getKeyInt("lagmax",-1);
  calc_covar_ = !analyzeArgs.hasKey("nocovar");
  usefft_ = !analyzeArgs.hasKey("direct");
  // Select datasets from remaining args
  dsets_.clear();
  ArgList dsetArgs = analyzeArgs.RemainingArgs();
  for (ArgList::const_iterator dsa = dsetArgs.begin(); dsa != dsetArgs.end(); ++dsa) {
    DataSetList setsIn = setup.DSL().GetMultipleSets( *dsa );
    for (DataSetList::const_iterator ds = setsIn.begin(); ds != setsIn.end(); ++ds) {
      if ( (*ds)->Group() != DataSet::SCALAR_1D && (*ds)->Type() != DataSet::VECTOR )
        mprintf("Warning: Set '%s' type not supported in AUTOCORR - skipping.\n",
                (*ds)->legend());
      else
        dsets_.push_back( *ds );
    }
  }
  if (dsets_.empty()) {
    mprinterr("Error: No data sets selected.\n");
    return Analysis::ERR;
  }
  // If setname is empty generate a default name
  if (setname.empty())
    setname = setup.DSL().GenerateDefaultName( "autocorr" );
  // Setup output datasets
  MetaData md( setname );
  for (unsigned int idx = 0; idx != dsets_.size(); idx++) {
    md.SetIdx( idx );
    DataSet* dsout = setup.DSL().AddSet( DataSet::DOUBLE, md );
    if (dsout==0) return Analysis::ERR;
    dsout->SetLegend( dsets_[idx]->Meta().Legend() );
    outputData_.push_back( dsout );
    // Add set to output file
    if (outfile != 0) outfile->AddDataSet( outputData_.back() );
  }
 
  if (calc_covar_)
    calctype = "covariance";
  else
    calctype = "correlation";
 
  mprintf("    AUTOCORR: Calculating auto-%s for %i data sets:\n\t", calctype, dsets_.size());
  for (unsigned int idx = 0; idx != dsets_.size(); ++idx)
    mprintf(" %s", dsets_[idx]->legend());
  mprintf("\n");
  if (lagmax_!=-1)
    mprintf("\tLag max= %i\n", lagmax_);
  if ( !setname.empty() )
    mprintf("\tSet name: %s\n", setname.c_str() );
  if ( outfile != 0 )
    mprintf("\tOutfile name: %s\n", outfile->DataFilename().base());
  if (usefft_)
    mprintf("\tUsing FFT to calculate %s.\n", calctype);
  else
    mprintf("\tUsing direct method to calculate %s.\n", calctype);

  return Analysis::OK;
}

// 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;
}

// Analysis_TI::Setup()
Analysis::RetType Analysis_TI::Setup(ArgList& analyzeArgs, DataSetList* datasetlist, DataFileList* DFLin, int debugIn)
{
  int nq = analyzeArgs.getKeyInt("nq", 0);
  ArgList nskipArg(analyzeArgs.GetStringKey("nskip"), ","); // Comma-separated
  if (nskipArg.empty())
    nskip_.resize(1, 0);
  else {
    nskip_.clear();
    for (int i = 0; i != nskipArg.Nargs(); i++) {
      nskip_.push_back( nskipArg.getNextInteger(0) );
      if (nskip_.back() < 0) nskip_.back() = 0;
    }
  }
  std::string setname = analyzeArgs.GetStringKey("name");
  DataFile* outfile = DFLin->AddDataFile(analyzeArgs.GetStringKey("out"), analyzeArgs);
  DataFile* curveout = DFLin->AddDataFile(analyzeArgs.GetStringKey("curveout"), analyzeArgs);
  // Select datasets from remaining args
  if (input_dsets_.AddSetsFromArgs( analyzeArgs.RemainingArgs(), *datasetlist )) {
    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;
  if (quad_.size() != input_dsets_.size()) {
    mprinterr("Error: Expected %zu data sets based on nq, got %zu\n",
              quad_.size(), input_dsets_.size());
    return Analysis::ERR;
  }
  dAout_ = datasetlist->AddSet(DataSet::XYMESH, setname, "TI");
  if (dAout_ == 0) return Analysis::ERR;
  if (outfile != 0) outfile->AddDataSet( dAout_ );
  MetaData md(dAout_->Meta().Name(), "TIcurve");
  for (Iarray::const_iterator it = nskip_.begin(); it != nskip_.end(); ++it) {
    md.SetIdx( *it );
    DataSet* ds = datasetlist->AddSet(DataSet::XYMESH, md);
    if (ds == 0) return Analysis::ERR;
    ds->SetLegend( md.Name() + "_Skip" + integerToString(*it) );
    if (curveout != 0) curveout->AddDataSet( ds );
    curve_.push_back( ds );
  }

  mprintf("    TI: Calculating TI using Gaussian quadrature with %zu points.\n",
          quad_.size());
  mprintf("\t%6s %8s %8s %s\n", "Point", "Abscissa", "Weight", "SetName");
  for (unsigned int i = 0; i != quad_.size(); i++)
    mprintf("\t%6i %8.5f %8.5f %s\n", i, quad_[i], wgt_[i], input_dsets_[i]->legend());
  if (nskip_.front() > 0) {
    mprintf("\tSkipping first");
    for (Iarray::const_iterator it = nskip_.begin(); it != nskip_.end(); ++it)
      mprintf(" %i", *it);
    mprintf(" data points for <DV/DL> calc.\n");
  }
  mprintf("\tResults saved in set '%s'\n", dAout_->legend());
  mprintf("\tTI curve(s) saved in set(s)");
  for (DSarray::const_iterator ds = curve_.begin(); ds != curve_.end(); ++ds)
    mprintf(" '%s'", (*ds)->legend());
  mprintf("\n");
  if (outfile != 0) mprintf("\tResults written to '%s'\n", outfile->DataFilename().full());
  if (curveout!= 0) mprintf("\tTI curve written to '%s'\n", curveout->DataFilename().full());
  return Analysis::OK;
}

// Action_Density::Init()
Action::RetType Action_Density::Init(ArgList& actionArgs, ActionInit& init, int debugIn)
{
# ifdef MPI
  trajComm_ = init.TrajComm();
# endif
  DataFile* outfile = init.DFL().AddDataFile(actionArgs.GetStringKey("out"), actionArgs);

  std::string dsname = actionArgs.GetStringKey("name");

  if (actionArgs.hasKey("x") ) {
    axis_ = DX;
    area_coord_[0] = DY;
    area_coord_[1] = DZ;
  } else if (actionArgs.hasKey("y") ) {
    axis_ = DY;
    area_coord_[0] = DX;
    area_coord_[1] = DZ;
  } else if (actionArgs.hasKey("z") ) {
    axis_ = DZ;
    area_coord_[0] = DX;
    area_coord_[1] = DY;
  }

  property_ = NUMBER;
  if      (actionArgs.hasKey("number") )   property_ = NUMBER;
  else if (actionArgs.hasKey("mass") )     property_ = MASS;
  else if (actionArgs.hasKey("charge") )   property_ = CHARGE;
  else if (actionArgs.hasKey("electron") ) property_ = ELECTRON;

  binType_ = CENTER;
  if      (actionArgs.hasKey("bincenter")) binType_ = CENTER;
  else if (actionArgs.hasKey("binedge")  ) binType_ = EDGE;

  delta_ = actionArgs.getKeyDouble("delta", 0.01);
  if (delta_ <= 0) {
    mprinterr("Error: Delta must be > 0.0\n");
    return Action::ERR;
  }

  // for compatibility with ptraj, ignored because we rely on the atom code to
  // do the right thing, see Atom.{h,cpp}
  if (actionArgs.hasKey("efile"))
    mprintf("Warning: The 'efile' keyword is deprecated.\n");

  // read the rest of the command line as a series of masks
  std::string maskstr;

  unsigned int idx = 1;
  while ( (maskstr = actionArgs.GetMaskNext() ) != emptystring) {
    masks_.push_back( AtomMask(maskstr) );
    if (dsname.empty())
      dsname = init.DSL().GenerateDefaultName("DENSITY");
    MetaData MD(dsname, "avg", idx);
    MD.SetTimeSeries( MetaData::NOT_TS );
    // Hold average density
    DataSet* ads = init.DSL().AddSet( DataSet::DOUBLE, MD );
    if (ads == 0) return Action::ERR;
    ads->SetLegend( NoWhitespace(masks_.back().MaskExpression()) );
    AvSets_.push_back( ads );
    if (outfile != 0) outfile->AddDataSet( ads );
    // Hold SD density
    MD.SetAspect("sd");
    DataSet* sds = init.DSL().AddSet( DataSet::DOUBLE, MD );
    if (sds == 0) return Action::ERR;
    sds->SetLegend( NoWhitespace("sd(" + masks_.back().MaskExpression() + ")") );
    SdSets_.push_back( sds );
    if (outfile != 0) outfile->AddDataSet( sds );
#   ifdef MPI
    ads->SetNeedsSync( false ); // Populated in Print()
    sds->SetNeedsSync( false );
#   endif
    idx++;
  }
  if (masks_.empty()) {
    // If no masks assume we want total system density.
    if (dsname.empty())
      dsname = actionArgs.GetStringNext();
    density_ = init.DSL().AddSet(DataSet::DOUBLE, dsname, "DENSITY");
    if (density_ == 0) return Action::ERR;
    if (outfile != 0) outfile->AddDataSet( density_ );
    image_.InitImaging( true );
    // Hijack delta for storing sum of masses
    delta_ = 0.0;
  } else {
    // Density selected by mask(s) along an axis
    density_ = 0;
    histograms_.resize(masks_.size() );
  }

  mprintf("    DENSITY:");
  if (density_ == 0) {
    const char* binStr[] = {"center", "edge"};
    mprintf(" Determining %s density for %zu masks.\n", PropertyStr_[property_], masks_.size());
    mprintf("\troutine version: %s\n", ROUTINE_VERSION_STRING);
    mprintf("\tDelta is %f\n", delta_);
    mprintf("\tAxis is %s\n", AxisStr_[axis_]);
    mprintf("\tData set name is '%s'\n", dsname.c_str());
    mprintf("\tData set aspect [avg] holds the mean, aspect [sd] holds standard deviation.\n");
    mprintf("\tBin coordinates will be to bin %s.\n", binStr[binType_]);
//.........这里部分代码省略.........

Analysis::RetType Analysis_Lifetime::Setup(ArgList& analyzeArgs, DataSetList* datasetlist,
                            TopologyList* PFLin, DataFileList* DFLin, int debugIn)
{
  // Get Keywords
  DataFile* outfile = DFLin->AddDataFile(analyzeArgs.GetStringKey("out"), analyzeArgs);
  if (outfile != 0) outfile->ProcessArgs("noemptyframes");
  DataFile* maxfile = 0;
  DataFile* avgfile = 0;
  std::string setname = analyzeArgs.GetStringKey("name");
  windowSize_ = analyzeArgs.getKeyInt("window", -1);
  averageonly_ = analyzeArgs.hasKey("averageonly");
  if (!averageonly_ && outfile != 0) {
    maxfile = DFLin->AddDataFile("max." + outfile->DataFilename().Full(), analyzeArgs);
    maxfile->ProcessArgs("noemptyframes");
    avgfile = DFLin->AddDataFile("avg." + outfile->DataFilename().Full(), analyzeArgs);
    avgfile->ProcessArgs("noemptyframes");
  }
  cumulative_ = analyzeArgs.hasKey("cumulative");
  deltaAvg_ = analyzeArgs.hasKey("delta");
  cut_ = analyzeArgs.getKeyDouble("cut", 0.5);
  // Select datasets from remaining args
  ArgList dsetArgs = analyzeArgs.RemainingArgs();
  for (ArgList::const_iterator dsa = dsetArgs.begin(); dsa != dsetArgs.end(); ++dsa)
    inputDsets_ += datasetlist->GetMultipleSets( *dsa );
  if (inputDsets_.empty()) {
    mprinterr("Error: lifetime: No data sets selected.\n");
    return Analysis::ERR;
  }
  // Sort input datasets
  inputDsets_.sort();

  // Create output datasets
  if ( windowSize_ != -1) {
    if (setname.empty()) 
      setname = datasetlist->GenerateDefaultName( "lifetime" );
    int didx = 0;
    for (DataSetList::const_iterator set = inputDsets_.begin(); set != inputDsets_.end(); ++set)
    {
      DataSet* outSet = datasetlist->AddSetIdx( DataSet::FLOAT, setname, didx );
      if (outSet==0) {
        mprinterr("Error: lifetime: Could not allocate output set for %s\n", 
                  (*set)->Legend().c_str());
        return Analysis::ERR;
      }
      outSet->SetLegend( (*set)->Legend() );
      outputDsets_.push_back( outSet );
      if (outfile != 0) outfile->AddSet( outSet );
      if (!averageonly_) {
        // MAX
        // FIXME: CHeck for nullS
        outSet = datasetlist->AddSetIdxAspect( DataSet::INT, setname, didx, "max" );
        outSet->SetLegend( (*set)->Legend() );
        maxDsets_.push_back( outSet );
        if (maxfile != 0) maxfile->AddSet( outSet );
        // AVG
        outSet = datasetlist->AddSetIdxAspect( DataSet::FLOAT, setname, didx, "avg" );
        outSet->SetLegend( (*set)->Legend() );
        avgDsets_.push_back( outSet );
        if (avgfile != 0) avgfile->AddSet( outSet );
      }
      ++didx;
    }
  } else if (outfile != 0) {
    mprinterr("Error: Output file name specified but no window size given ('window <N>')\n");
    return Analysis::ERR;
  }

  if (!averageonly_)
    mprintf("    LIFETIME: Calculating average lifetime using a cutoff of %f", cut_);
  else
    mprintf("    LIFETIME: Calculating only averages");
  mprintf(" of data in %i sets\n", inputDsets_.size());
  if (debugIn > 0)
    inputDsets_.List();
  if (windowSize_ != -1) {
    mprintf("\tAverage of data over windows will be saved to sets named %s\n",
            setname.c_str());
    mprintf("\tWindow size for averaging: %i\n", windowSize_);
    if (cumulative_)
      mprintf("\tCumulative averages will be saved.\n");
    if (deltaAvg_)
      mprintf("\tChange of average from previous average will be saved.\n");
    if (outfile != 0) {
      mprintf("\tOutfile: %s", outfile->DataFilename().base());
      if (!averageonly_)
        mprintf(", %s, %s", maxfile->DataFilename().base(), avgfile->DataFilename().base());
      mprintf("\n");
    }
  }

  return Analysis::OK;
}

Analysis::RetType Analysis_State::Analyze() {
  // Only process as much data as is in the smallest data set.
  size_t nframes = States_.front().DS().Size();
  for (StateArray::const_iterator state = States_.begin(); state != States_.end(); ++state)
  {
    if ( state != States_.begin() && nframes != state->DS().Size() )
      mprintf("Warning: Set '%s' for state '%s' has a different # of frames (%zu)"
              " than previous sets (%zu).\n", state->DS().legend(), state->id(),
              state->DS().Size(), nframes);
    nframes = std::min( nframes, state->DS().Size() );
  }
  mprintf("\tProcessing %zu frames.\n", nframes);
  if (nframes < 1) return Analysis::ERR;

  std::vector<Transition> Status;
  Status.reserve( States_.size() + 1 ); // +1 for state -1, undefined
  for (int i = 0; i != (int)States_.size() + 1; i++) {
    DataSet* ds = masterDSL_->AddSet(DataSet::DOUBLE, 
                                     MetaData(state_data_->Meta().Name(), "sCurve", i));
    if (ds == 0) return Analysis::ERR;
    if (curveOut_ != 0) curveOut_->AddDataSet( ds );
    ds->SetLegend( StateName(i-1) );
    Status.push_back( Transition((DataSet_double*)ds) );
  }

  int last_state = -2;
  size_t last_State_Start = 0;
  size_t final_frame = nframes - 1;
  for (size_t frm = 0; frm != nframes; frm++) {
    // Determine which state we are in.
    int state_num = -1;
    for (StateArray::const_iterator state = States_.begin(); state != States_.end(); ++state)
    {
      double dVal = state->DS().Dval( frm );
      if (dVal >= state->Min() && dVal < state->Max()) { // TODO: Periodic
        if (state_num != -1)
          mprintf("Warning: Frame %zu already defined as state '%s', also matches state '%s'.\n",
                  frm, States_[state_num].id(), state->id());
        else
          state_num = (int)(state - States_.begin());
      }
    }
    state_data_->Add( frm, &state_num );

    // Determine if there has been a transition.
    if (last_state == -2)
      // First frame, no possible transition yet.
      last_state = state_num;
    else if (state_num != last_state || frm == final_frame) {
      int length = (int)(frm - last_State_Start);
      if (state_num != last_state) {
        // There has been a transition from last_state to state_num.
        StatePair sPair(last_state, state_num);
        TransMapType::iterator entry = TransitionMap_.find( sPair );
        if (entry == TransitionMap_.end()) {
          // New transition
          DataSet* ds = masterDSL_->AddSet( DataSet::DOUBLE,
                                            MetaData(state_data_->Meta().Name(), "tCurve",
                                                     TransitionMap_.size()) );
          if (ds == 0) return Analysis::ERR;
          if (curveOut_ != 0) curveOut_->AddDataSet( ds );
          ds->SetLegend( StateName(last_state) + "->" + StateName(state_num) );
          TransitionMap_.insert( TransPair(sPair, Transition(length, (DataSet_double*)ds)) );
        } else
          // Update previous transition
          entry->second.Update(length);
      }
      // Update single state information.
      Status[last_state + 1].Update(length); 
      last_State_Start = frm;
      last_state = state_num;
    }
  }

  // DEBUG: Print single state info.
  if (debug_ > 0) {
    mprintf("  States:\n");
    mprintf("\t%i: %s  max= %i  sum= %i  n= %i  Avg= %g\n", -1, "Undefined",
            Status.front().Max(), Status.front().Sum(),
            Status.front().Nlifetimes(), Status.front().Avg());
    StateArray::const_iterator state = States_.begin();
    for (std::vector<Transition>::const_iterator s = Status.begin() + 1;
                                                 s != Status.end(); ++s, ++state)
    mprintf("\t%u: %s  max= %i  sum= %i  n= %i  Avg= %g\n", state - States_.begin(),
            state->id(), s->Max(), s->Sum(), s->Nlifetimes(), s->Avg());
    // DEBUG: Print transitions.
    mprintf("  Transitions:\n");
    for (TransMapType::const_iterator trans = TransitionMap_.begin();
                                      trans != TransitionMap_.end(); ++trans)
      mprintf("\t%i -> %i: max= %i  sum= %i  n= %i  Avg= %g\n", trans->first.first,
              trans->first.second, trans->second.Max(), trans->second.Sum(),
              trans->second.Nlifetimes(), trans->second.Avg());
  }
  stateOut_->Printf("%-8s %12s %12s %12s %s\n", "#Index", "N", "Average", "Max", "State");
  for (int idx = 0; idx != (int)Status.size(); idx++)
    stateOut_->Printf("%-8i %12i %12.4f %12i %s\n", idx-1, Status[idx].Nlifetimes(),
                      Status[idx].Avg(), Status[idx].Max(), stateName(idx-1));
  transOut_->Printf("%-12s %12s %12s %s\n", "#N", "Average", "Max", "Transition");
  for (TransMapType::const_iterator trans = TransitionMap_.begin();
                                    trans != TransitionMap_.end(); ++trans)
//.........这里部分代码省略.........

/** 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;
}

Analysis::RetType Analysis_Spline::Setup(ArgList& analyzeArgs, DataSetList* datasetlist, DataFileList* DFLin, int debugIn)
{
  std::string setname = analyzeArgs.GetStringKey("name");
  outfile_ = DFLin->AddDataFile(analyzeArgs.GetStringKey("out"), analyzeArgs);
  meshsize_ = analyzeArgs.getKeyInt("meshsize", 0);
  meshfactor_ = -1.0;
  if (meshsize_ < 3) {
    meshfactor_ = analyzeArgs.getKeyDouble("meshfactor", -1.0);
    if (meshfactor_ < Constants::SMALL) {
      mprinterr("Error: Either meshsize must be specified and > 2, or meshfactor must be\n"
                "Error:   specified and > 0.0\n");
      return Analysis::ERR;
    }
  }
  if (analyzeArgs.Contains("meshmin")) {
    meshmin_ = analyzeArgs.getKeyDouble("meshmin", 0.0);
    useDefaultMin_ = true;
  } else
    useDefaultMin_ = false;
  if (analyzeArgs.Contains("meshmax")) {
    meshmax_ = analyzeArgs.getKeyDouble("meshmax", -1.0);
    useDefaultMax_ = true;
  } else
    useDefaultMax_ = false;
  if (useDefaultMin_ && useDefaultMax_ && meshmax_ < meshmin_) {
    mprinterr("Error: meshmax must be > meshmin\n");
    return Analysis::ERR;
  }
  // Select datasets from remaining args
  if (input_dsets_.AddSetsFromArgs( analyzeArgs.RemainingArgs(), *datasetlist )) {
    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;
  }

  // Set up output datasets
  Dimension Xdim( meshmin_, (meshmax_ - meshmin_) / (double)meshsize_ );
  for (Array1D::const_iterator dsIn = input_dsets_.begin();
                               dsIn != input_dsets_.end(); ++dsIn)
  {
    DataSet* ds = datasetlist->AddSet(DataSet::XYMESH, setname, "Spline");
    if (ds == 0) return Analysis::ERR;
    ds->SetLegend( "Spline(" + (*dsIn)->Meta().Legend() + ")" );
    // TODO: Set individually based on input_dsets_
    ds->SetDim(Dimension::X, Xdim);
    if (outfile_ != 0) outfile_->AddDataSet( ds );
    output_dsets_.push_back( (DataSet_Mesh*)ds );
  }

  mprintf("    SPLINE: Applying cubic splining to %u data sets\n", input_dsets_.size());
  if (meshfactor_ < 0)
    mprintf("\tMesh size= %i\n", meshsize_);
  else
    mprintf("\tMesh size will be input set size multiplied by %f\n", meshfactor_);
  if (useDefaultMin_)
    mprintf("\tMesh min= %f,", meshmin_);
  else
    mprintf("\tMesh min will be input set min,");
  if (useDefaultMax_)
    mprintf(" Mesh max= %f\n", meshmax_);
  else
    mprintf(" Mesh max will be input set max.\n");
  if (outfile_ != 0) {
    if (!setname.empty())
      mprintf("\tOutput set name: %s\n", setname.c_str());
    mprintf("\tOutfile name: %s\n", outfile_->DataFilename().base());
  }
  //for (Array1D::const_iterator set = input_dsets_.begin(); set != input_dsets_.end(); ++set)
  //  mprintf("\t%s\n", (*set)->legend());
  return Analysis::OK;
}