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

// Analysis_CrdFluct::CalcBfactors()
void Analysis_CrdFluct::CalcBfactors( Frame SumCoords, Frame SumCoords2, double Nsets,
                                      DataSet& outset )
{
  SumCoords.Divide(Nsets);
  SumCoords2.Divide(Nsets);
  //SumCoords2 = SumCoords2 - (SumCoords * SumCoords);
  SumCoords *= SumCoords;
  SumCoords2 -= SumCoords;
  AtomMask::const_iterator maskat = mask_.begin();
  if (bfactor_) {
    // Set up b factor normalization
    // B-factors are (8/3)*PI*PI * <r>**2 hence we do not sqrt the fluctuations
    double bfac = (8.0/3.0)*Constants::PI*Constants::PI;
    for (int i = 0; i < SumCoords2.size(); i+=3) {
      double fluct = (SumCoords2[i] + SumCoords2[i+1] + SumCoords2[i+2]) * bfac;
      outset.Add( *maskat, &fluct );
      ++maskat; 
    }
  } else {
    // Atomic fluctuations
    for (int i = 0; i < SumCoords2.size(); i+=3) {
      double fluct = SumCoords2[i] + SumCoords2[i+1] + SumCoords2[i+2];
      if (fluct > 0)
        outset.Add( *maskat, &fluct );
      ++maskat;
    }
  }
}

DataSet Model::Run(int totalTime, int totalEnergy)
{
  DataSet data;
  
  for(int timeStep = 0; timeStep < totalTime; timeStep++)
    {
      int populationCount = 0;
      int deathCount = 0;
      int reproductionCount = 0;
      
      // initialise random seed
      srand(time(NULL));

      // distribute energy
      for(int energyCount = 0; energyCount < totalEnergy; energyCount++)
	{
	  int index = rand() % population.size();
	  population[index].Consume(1);
	}

      // check for expiration and reproduction
      // and age all organisms 
      for(vector<Organism>::iterator it = population.begin();
	  it != population.end(); ++it)
	{
	  if(!it->Expire())
	    {
	      populationCount++;
	      if(it->Reproduce())
		{
		  reproductionCount++;
		}
	      it->Age();
	    }else
	    {
	      // Erase the expired entry
	      it = population.erase(it);
	      deathCount++;
	    }
	}

      // add new organisms as determined by 
      // reporduction count
      for(int newOrganisms = 0;
	  newOrganisms < reproductionCount;
	  newOrganisms++)
	{
	  Organism o;
	  population.push_back(o);
	}
      data.Add( populationCount,reproductionCount,deathCount);
    }
  return data;
}

/** For per-residue RMSD only. Sync the per-residue RMSD data set since
  * it is not part of the master DataSetList in Cpptraj. Setup output
  * file options. Calculate averages if requested.
  */
void Action_Rmsd::Print() {
  if (!perres_ || PerResRMSD_.empty()) return;
  // Per-residue output
  if (perresout_ != 0) {
    // Add data sets to perresout
    for (std::vector<DataSet*>::iterator set = PerResRMSD_.begin();
                                         set != PerResRMSD_.end(); ++set)
      perresout_->AddSet(*set);
    // Set output file to be inverted if requested
    if (perresinvert_) 
      perresout_->ProcessArgs("invert");
    mprintf("    RMSD: Per-residue: Writing data for %zu residues to %s\n",
            PerResRMSD_.size(), perresout_->DataFilename().base());
  }

  // Average
  if (perresavg_ != 0) {
    int Nperres = (int)PerResRMSD_.size();
    // Use the per residue rmsd dataset list to add one more for averaging
    DataSet* PerResAvg = masterDSL_->AddSetAspect(DataSet::DOUBLE, rmsd_->Name(), "Avg");
    // another for stdev
    DataSet* PerResStdev = masterDSL_->AddSetAspect(DataSet::DOUBLE, rmsd_->Name(), "Stdev");
    // Add the average and stdev datasets to the master datafile list
    perresavg_->AddSet(PerResAvg);
    perresavg_->AddSet(PerResStdev);
    perresavg_->ProcessArgs("xlabel Residue");
    // For each residue, get the average rmsd
    double stdev = 0;
    double avg = 0;
    for (int pridx = 0; pridx < Nperres; pridx++) {
      avg = DS_Math::Avg(*PerResRMSD_[pridx], &stdev );
      int dsidx = PerResRMSD_[pridx]->Idx() - 1;
      PerResAvg->Add(dsidx, &avg);
      PerResStdev->Add(dsidx,&stdev);
    }
  }
}

// Exec_DataSetCmd::VectorCoord()
Exec::RetType Exec_DataSetCmd::VectorCoord(CpptrajState& State, ArgList& argIn) {
  // Keywords
  std::string name = argIn.GetStringKey("name");
  int idx;
  if (argIn.hasKey("X"))
    idx = 0;
  else if (argIn.hasKey("Y"))
    idx = 1;
  else if (argIn.hasKey("Z"))
    idx = 2;
  else {
    mprinterr("Error: 'vectorcoord' requires specifying X, Y, or Z.\n");
    return CpptrajState::ERR;
  }
  // Data set
  DataSet* ds1 = State.DSL().GetDataSet( argIn.GetStringNext() );
  if (ds1 == 0) return CpptrajState::ERR;
  if (ds1->Type() != DataSet::VECTOR) {
    mprinterr("Error: 'vectorcoord' only works with vector data sets.\n");
    return CpptrajState::ERR;
  }
  if (ds1->Size() < 1) {
    mprinterr("Error: '%s' is empty.\n", ds1->legend());
    return CpptrajState::ERR;
  }
  // Create output set.
  static const char* XYZ[3] = { "X", "Y", "Z" };
  DataSet* out = State.DSL().AddSet( DataSet::DOUBLE, name, "COORD");
  if (out == 0) return CpptrajState::ERR;
  // Extract data
  mprintf("\tExtracting %s coordinate from vector %s to %s\n",
          XYZ[idx], ds1->legend(), out->Meta().PrintName().c_str());
  DataSet_Vector const& vec = static_cast<DataSet_Vector const&>( *ds1 );
  for (unsigned int n = 0; n != vec.Size(); n++) {
    double d = vec.VXYZ(n)[idx];
    out->Add( n, &d );
  }
  return CpptrajState::OK;
}

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