C++ DataSetList::AddSet方法代码示例

int Action_Esander::AddSet(Energy_Sander::Etype typeIn, DataSetList& DslIn, DataFile* outfile,
                           std::string const& setname)
{
  Esets_[typeIn] = DslIn.AddSet(DataSet::DOUBLE, MetaData(setname, Energy_Sander::Easpect(typeIn)));
  if (Esets_[typeIn] == 0) return 1;
  if (outfile != 0) outfile->AddDataSet( Esets_[typeIn] );
  return 0;
}

int Action_Energy::AddSet(Etype typeIn, DataSetList& DslIn, DataFile* outfile,
                          std::string const& setname)
{
  Energy_[typeIn] = DslIn.AddSet(DataSet::DOUBLE, MetaData(setname, Estring[typeIn]));
  if (Energy_[typeIn] == 0) return 1;
  if (outfile != 0) outfile->AddDataSet( Energy_[typeIn] );
  return 0;
}

// DataIO_OpenDx::ReadData()
int DataIO_OpenDx::ReadData(FileName const& fname, 
                            DataSetList& datasetlist, std::string const& dsname)
{
  // Add grid data set. Default to float for now.
  DataSet* ds = datasetlist.AddSet( DataSet::GRID_FLT, dsname, "GRID" );
  if (ds==0) return 1;
  if (LoadGrid(fname.full(), *ds)) {
    // Load failed. Erase grid data set.
    datasetlist.RemoveSet( ds );
    return 1;
  }
  return 0;
}

// Analysis_Hist::Setup()
Analysis::RetType Analysis_Hist::ExternalSetup(DataSet_1D* dsIn, std::string const& histname,
                                       int setidx, std::string const& outfilenameIn,
                                       bool minArgSetIn, double minIn,
                                       bool maxArgSetIn, double maxIn,
                                       double stepIn, int binsIn, double tempIn,
                                       NormMode normIn,
                                       DataSetList& datasetlist, DataFileList& DFLin)
{
  debug_ = 0;
  if (dsIn == 0) return Analysis::ERR;
  outfilename_ = outfilenameIn;
  outfile_ = DFLin.AddDataFile(outfilename_);
  Temp_ = tempIn; 
  if (Temp_ != -1.0)
    calcFreeE_ = true;
  else
    calcFreeE_ = false;
  gnuplot_ = false;
  normalize_ = normIn;
  circular_ = false;
  nativeOut_ = false;
  minArgSet_ = minArgSetIn;
  if (minArgSet_)
    default_min_ = minIn;
  maxArgSet_ = maxArgSetIn;
  if (maxArgSet_)
    default_max_ = maxIn;
  default_step_ = stepIn;
  default_bins_ = binsIn;
  calcAMD_ = false;
  amddata_ = 0;

  dimensionArgs_.push_back( ArgList(dsIn->Meta().Legend()) ); // Needed for dim label
  histdata_.push_back( dsIn );
  N_dimensions_ = 1;
  std::string setname = histname;
  std::string htype;
  if (calcFreeE_)
    htype = "FreeE_";
  else
    htype = "Hist_";
  if (setname.empty())
    setname = datasetlist.GenerateDefaultName(htype + dsIn->Meta().Name());
  hist_ = datasetlist.AddSet( DataSet::DOUBLE, MetaData(setname, dsIn->Meta().Aspect(), setidx) );
  if (hist_ == 0) return Analysis::ERR;
  hist_->SetLegend(htype + dsIn->Meta().Legend());
  if (outfile_ != 0) outfile_->AddDataSet( hist_ );
  return Analysis::OK;
}

Analysis::RetType Analysis_KDE::ExternalSetup(DataSet_1D* dsIn, std::string const& histname,
                                       int setidx, std::string const& outfilenameIn,
                                       bool minArgSetIn, double minIn,
                                       bool maxArgSetIn, double maxIn,
                                       double stepIn, int binsIn, double tempIn,
                                       DataSetList& datasetlist, DataFileList& DFLin)
{
  if (dsIn == 0) return Analysis::ERR;
  data_ = dsIn;
  q_data_ = 0;
  kldiv_ = 0;
  amddata_ = 0;
  bandwidth_ = -1.0;
  minArgSet_ = minArgSetIn;
  if (minArgSet_)
    default_min_ = minIn;
  maxArgSet_ = maxArgSetIn;
  if (maxArgSet_)
    default_max_ = maxIn;
  default_step_ = stepIn;
  default_bins_ = binsIn;
  Temp_ = tempIn;
  if (Temp_ != -1.0)
    calcFreeE_ = true;
  else
    calcFreeE_ = false;
  std::string setname = histname;
  std::string htype;
  if (calcFreeE_)
    htype = "FreeE_";
  else
    htype = "KDE_";
  if (setname.empty())
    setname = datasetlist.GenerateDefaultName(htype + dsIn->Meta().Name());
  DataFile* outfile = DFLin.AddDataFile( outfilenameIn );
  output_ = datasetlist.AddSet(DataSet::DOUBLE, MetaData(setname, dsIn->Meta().Aspect(), setidx));
  if (output_ == 0) return Analysis::ERR;
  output_->SetLegend(htype + dsIn->Meta().Legend());
  if (outfile != 0) outfile->AddDataSet( output_ );
  return Analysis::OK;
}

//  Exec_SortEnsembleData::Sort_pH_Data()
int Exec_SortEnsembleData::Sort_pH_Data(DataSetList const& setsToSort, DataSetList& OutputSets,
                                        unsigned int maxFrames)
const
{
  // Cast sets back to DataSet_PHREMD
  typedef std::vector<DataSet_PHREMD*> Parray;
  Parray PHsets;
  for (DataSetList::const_iterator ds = setsToSort.begin(); ds != setsToSort.end(); ++ds)
    PHsets.push_back( (DataSet_PHREMD*)*ds );

  // Gather initial pH data values, ensure no duplicates
  typedef std::vector<double> Darray;
  Darray pHvalues;
# ifdef MPI
  pHvalues.resize( Parallel::Ensemble_Size() );
  Darray phtmp;
  for (Parray::const_iterator ds = PHsets.begin(); ds != PHsets.end(); ++ds)
    phtmp.push_back( (*ds)->Initial_pH() );
  if (comm_.AllGather(&phtmp[0], phtmp.size(), MPI_DOUBLE, &pHvalues[0])) {
    rprinterr("Error: Gathering pH values.\n");
    return 1;
  }
# else
  for (Parray::const_iterator ds = PHsets.begin(); ds != PHsets.end(); ++ds)
    pHvalues.push_back( (*ds)->Initial_pH() );
# endif
  ReplicaInfo::Map<double> pH_map;
  if (pH_map.CreateMap( pHvalues )) {
    rprinterr("Error: Duplicate pH value detected (%.2f) in ensemble.\n", pH_map.Duplicate());
    return 1;
  }
  Darray sortedPH;
  mprintf("\tInitial pH values:");
  for (ReplicaInfo::Map<double>::const_iterator ph = pH_map.begin(); ph != pH_map.end(); ++ph)
  {
    mprintf(" %6.2f", ph->first);
    sortedPH.push_back( ph->first );
  }
  mprintf("\n");

  // Create sets to hold sorted pH values. Create a set for each pH value
  // and each residue. Final output sets will be PH0R0, PH0R1, PH1R0, ...
  // TODO check that residue info all the same
  DataSet_PHREMD::Rarray const& Residues = PHsets[0]->Residues();
  int defaultState = 0;
# ifdef MPI
  if ( PHsets[0]->Type() == DataSet::PH_IMPL)
    defaultState = -1;
# endif
  if (debug_ > 0)
    rprintf("DEBUG: Sorting %u frames for %zu sets, %zu pH values.\n",
            maxFrames, PHsets.size(), sortedPH.size());
  for (unsigned int idx = 0; idx != sortedPH.size(); idx++) {
    OutputSets.SetEnsembleNum( idx );
    for (unsigned int res = 0; res != Residues.size(); ++res) {
      MetaData md(PHsets[0]->Meta().Name(), Residues[res].Name().Truncated(), Residues[res].Num());
      DataSet_pH* out = (DataSet_pH*)OutputSets.AddSet( DataSet::PH, md );
      if (out==0) return 1;
      //out->SetLegend( "pH " + doubleToString( sortedPH[idx] ) );
      out->Set_Solvent_pH( sortedPH[idx] );
      out->SetResidueInfo( Residues[res] );
      out->SetTimeValues(PHsets[0]->Time());
      out->Resize( maxFrames, defaultState );
    }
  }

  // ---------------------------------------------
  if ( PHsets[0]->Type() == DataSet::PH_EXPL) {
    // Loop over unsorted sets
    for (Parray::const_iterator ds = PHsets.begin(); ds != PHsets.end(); ++ds)
    {
      DataSet_PHREMD_Explicit* in = (DataSet_PHREMD_Explicit*)*ds;
      unsigned int phidx = 0;
      for (unsigned int n = 0; n < maxFrames; n++)
      {
        float phval = in->pH_Values()[n];
        int setidx = pH_map.FindIndex( phval ) * Residues.size();
        //rprintf("DEBUG: %6u Set %10s pH= %6.2f going to %2i\n", n+1, in->legend(), phval, idx);
        //mflush();
        for (unsigned int res = 0; res < in->Residues().size(); res++, setidx++, phidx++)
        {
          DataSet_pH* out = (DataSet_pH*)OutputSets[setidx];
          //if (res == 0 && idx == 0) {
          //  rprintf("DEBUG: Frame %3u res %2u State %2i pH %6.2f\n", 
          //          n, res, in->Res(res).State(n), phval);
          //  mflush();
          //}
          out->SetState(n, in->ResStates()[phidx], in->RecordType(n));
        }
      }
    } // END loop over unsorted sets
#   ifdef MPI
    // Now we need to reduce down each set onto the thread where it belongs.
    if (Parallel::World().Size() > 1) {
      for (int idx = 0; idx != (int)OutputSets.size(); idx++) {
        DataSet_pH* out = (DataSet_pH*)OutputSets[idx];
        int ensembleRank = Parallel::MemberEnsCommRank( out->Meta().EnsembleNum() );
        //rprintf("DEBUG: Reduce set %s to rank %i\n", out->legend(), ensembleRank);
        out->Reduce( comm_, ensembleRank );
//.........这里部分代码省略.........

//.........这里部分代码省略.........
        // Perform endian swapping on header if necessary
        endian_swap(buffer.i, 56);
    }

    // Print DEBUG info
    if (debug_ > 0) {
        mprintf("DEBUG: Columns=%i  Rows=%i  Sections=%i\n", buffer.i[0], buffer.i[1], buffer.i[2]);
        mprintf("DEBUG: Mode=%i\n", buffer.i[3]);
        mprintf("DEBUG: Offsets: C=%i  R=%i  S=%i\n", buffer.i[4], buffer.i[5], buffer.i[6]);
        mprintf("DEBUG: NXYZ={ %i %i %i }\n", buffer.i[7], buffer.i[8], buffer.i[9]);
        mprintf("DEBUG: Box XYZ={ %f %f %f }  ABG={ %f %f %f }\n",
                buffer.f[10], buffer.f[11], buffer.f[12],
                buffer.f[13], buffer.f[14], buffer.f[15]);
        mprintf("DEBUG: Map: ColAxis=%i  RowAxis=%i  SecAxis=%i\n",
                buffer.i[16], buffer.i[17], buffer.i[18]);
        mprintf("DEBUG: SpaceGroup#=%i  SymmOpBytes=%i  SkewFlag=%i\n",
                buffer.i[22], buffer.i[23], buffer.i[24]);
        const int* MSKEW = buffer.i + 25;
        mprintf("DEBUG: Skew matrix: %i %i %i\n"
                "                    %i %i %i\n"
                "                    %i %i %i\n", MSKEW[0], MSKEW[1], MSKEW[2], MSKEW[3],
                MSKEW[4], MSKEW[5], MSKEW[6], MSKEW[7], MSKEW[8]);
        const int* TSKEW = buffer.i + 34;
        mprintf("DEBUG: Skew translation: %i %i %i\n", TSKEW[0], TSKEW[1], TSKEW[2]);
        mprintf("DEBUG: Nlabels=%i\n", buffer.i[55]);
    }
    // Check input data. Only support mode 2 for now.
    if (buffer.i[3] != 2) {
        mprinterr("Error: Mode %i; currently only mode 2 for CCP4 files is supported.\n", buffer.i[3]);
        return 1;
    }
    // Check offsets.
    if (buffer.i[4] != 0 || buffer.i[5] != 0 || buffer.i[6] != 0)
        mprintf("Warning: Non-zero offsets present. This is not yet supported and will be ignored.\n");
    // Check that mapping is col=x row=y section=z
    if (buffer.i[16] != 1 || buffer.i[17] != 2 || buffer.i[18] != 3) {
        mprinterr("Error: Currently only support cols=X, rows=Y, sections=Z\n");
        return 1;
    }
    if (buffer.i[24] != 0) {
        mprintf("Warning: Skew information present but not yet supported and will be ignored.\n");
        return 1;
    }

    // Read 10 80 character text labels
    char Labels[801];
    Labels[800] = '\0';
    infile.Read( Labels, 200*wSize );
    mprintf("\t%s\n", Labels);
    // Symmetry records: operators separated by * and grouped into 'lines' of 80 characters
    int NsymmRecords = buffer.i[23] / 80;
    if (NsymmRecords > 0) {
        char symBuffer[81];
        mprintf("\t%i symmetry records.\n", NsymmRecords);
        for (int ib = 0; ib != NsymmRecords; ib++) {
            infile.Gets( symBuffer, 80 );
            mprintf("\t%s\n", symBuffer);
        }
    }

    // Add grid data set. Default to float for now.
    DataSet* gridDS = datasetlist.AddSet( DataSet::GRID_FLT, dsname, "GRID" );
    if (gridDS == 0) return 1;
    DataSet_GridFlt& grid = static_cast<DataSet_GridFlt&>( *gridDS );
    // Allocate grid from dims and spacing. FIXME OK to assume zero origin?
    if (grid.Allocate_N_O_Box( buffer.i[7], buffer.i[8], buffer.i[9],
                               Vec3(0.0), Box(buffer.f + 10) ) != 0)
    {
        mprinterr("Error: Could not allocate grid.\n");
        return 1;
    }
    // FIXME: Grids are currently indexed so Z is fastest changing.
    //        Should be able to change indexing in grid DataSet.
    size_t mapSize = buffer.i[7] * buffer.i[8] * buffer.i[9];
    mprintf("\tCCP4 map has %zu elements\n", mapSize);
    mprintf("\tDensity: Min=%f  Max=%f  Mean=%f  RMS=%f\n",
            buffer.f[19], buffer.f[20], buffer.f[21], buffer.f[54]);
    std::vector<float> mapbuffer( mapSize );
    int mapBytes = mapSize * wSize;
    int numRead = infile.Read( &mapbuffer[0], mapBytes );
    if (numRead < 1) {
        mprinterr("Error: Could not read CCP4 map data.\n");
        return 1;
    } else if (numRead < mapBytes)
        mprintf("Warning: Expected %i bytes, read only %i bytes\n", mapBytes, numRead);
    if (isBigEndian) endian_swap(&mapbuffer[0], mapSize);

    // FIXME: Place data into grid DataSet with correct ordering.
    int gidx = 0;
    int NXY = buffer.i[7] * buffer.i[8];
    for (int ix = 0; ix != buffer.i[7]; ix++)
        for (int iy = 0; iy != buffer.i[8]; iy++)
            for (int iz = 0; iz != buffer.i[9]; iz++) {
                int midx = (iz * NXY) + (iy * buffer.i[7]) + ix;
                grid[gidx++] = mapbuffer[midx];
            }

    infile.CloseFile();
    return 0;
}

// DataIO_Std::Read_3D()
int DataIO_Std::Read_3D(std::string const& fname, 
                        DataSetList& datasetlist, std::string const& dsname)
{
  BufferedLine buffer;
  if (buffer.OpenFileRead( fname )) return 1;
  mprintf("\tData will be read as 3D grid: X Y Z Value\n");
  if (binCorners_)
    mprintf("\tAssuming X Y Z are bin corners\n");
  else
    mprintf("\tAssuming X Y Z are bin centers\n");
  const char* ptr = buffer.Line();
  // Check if #counts is present
  if (strncmp(ptr,"#counts",7)==0) {
    mprintf("\tReading grid dimensions.\n");
    unsigned int counts[3];
    sscanf(ptr+7,"%u %u %u", counts, counts+1, counts+2);
    for (int i = 0; i < 3; i++) {
      if (dims_[i] == 0)
        dims_[i] = counts[i];
      else if (dims_[i] != (size_t)counts[i])
        mprintf("Warning: Specified size for dim %i (%zu) differs from size in file (%u)\n",
                i, dims_[i], counts[i]);
    }
    ptr = buffer.Line();
  }
  if (dims_[0] == 0 || dims_[1] == 0 || dims_[2] == 0) {
    mprinterr("Error: 'dims' not specified for 'read3d' and no dims in file\n");
    return 1;
  }
  // Check if #origin is present
  if (strncmp(ptr,"#origin",7)==0) {
    mprintf("\tReading grid origin.\n");
    double oxyz[3];
    sscanf(ptr+7,"%lf %lf %lf", oxyz, oxyz+1, oxyz+2);
    for (int i = 0; i < 3; i++) {
      if (!originSpecified_)
        origin_[i] = oxyz[i];
      else if (origin_[i] != oxyz[i])
        mprintf("Warning: Specified origin for dim %i (%g) differs from origin in file (%g)\n",
                i, origin_[i], oxyz[i]);
    }
    ptr = buffer.Line();
  }
  // Check if #delta is present
  bool nonortho = false;
  Box gridBox;
  if (strncmp(ptr,"#delta",6)==0) {
    mprintf("\tReading grid deltas.\n");
    double dvals[9];
    int ndvals = sscanf(ptr+6,"%lf %lf %lf %lf %lf %lf %lf %lf %lf", dvals,
                        dvals+1, dvals+2, dvals+3, dvals+4, dvals+5,
                        dvals+6, dvals+7, dvals+8);
    if (ndvals == 3) {
      for (int i = 0; i < 3; i++) {
        if (!deltaSpecified_)
          delta_[i] = dvals[i];
        else if (delta_[i] != dvals[i])
          mprintf("Warning: Specified delta for dim %i (%g) differs from delta in file (%g)\n",
                  i, delta_[i], dvals[i]);
      }
    } else {
      nonortho = true;
      dvals[0] *= (double)dims_[0]; dvals[1] *= (double)dims_[0]; dvals[2] *= (double)dims_[0];
      dvals[3] *= (double)dims_[1]; dvals[4] *= (double)dims_[1]; dvals[5] *= (double)dims_[1];
      dvals[6] *= (double)dims_[2]; dvals[7] *= (double)dims_[2]; dvals[8] *= (double)dims_[2];
      gridBox = Box(Matrix_3x3(dvals));
    }
    ptr = buffer.Line();
  }
  // Get or allocate data set
  DataSet::DataType dtype;
  if (prec_ == DOUBLE) {
    dtype = DataSet::GRID_DBL;
    mprintf("\tGrid is double precision.\n");
  } else {
    dtype = DataSet::GRID_FLT;
    mprintf("\tGrid is single precision.\n");
  }
  MetaData md( dsname );
  DataSet_3D* ds = 0;
  DataSet* set = datasetlist.CheckForSet( md );
  if (set == 0) {
    ds = (DataSet_3D*)datasetlist.AddSet(dtype, dsname);
    if (ds == 0) return 1;
    int err = 0;
    if (nonortho)
      err = ds->Allocate_N_O_Box(dims_[0], dims_[1], dims_[2], origin_, gridBox);
    else
      err = ds->Allocate_N_O_D(dims_[0], dims_[1], dims_[2], origin_, delta_);
    if (err != 0) return 1;
  } else {
    mprintf("\tAppending to existing set '%s'\n", set->legend());
    if (set->Group() != DataSet::GRID_3D) {
      mprinterr("Error: Set '%s' is not a grid set, cannot append.\n", set->legend());
      return 1;
    }
    ds = (DataSet_3D*)set;
    // Check that dimensions line up. TODO check origin etc too?
    if (dims_[0] != ds->NX() ||
//.........这里部分代码省略.........

/** Read cluster matrix file. Can only get here if file has already been
  * determined to be in the proper format, so do no further error checking.
  * Expected format:
  *   <int> <int> <name>
  */
int DataIO_Std::ReadCmatrix(FileName const& fname,
                            DataSetList& datasetlist, std::string const& dsname)
{
  // Allocate output data set
  DataSet* ds = datasetlist.AddSet( DataSet::CMATRIX, dsname );
  if (ds == 0) return 1;
  DataSet_Cmatrix_MEM& Mat = static_cast<DataSet_Cmatrix_MEM&>( *ds );
  // Buffer file
  BufferedLine buffer;
  if (buffer.OpenFileRead( fname )) return 1;
  // Read past title. See if optional 'nframes' key is there.
  const char* ptr = buffer.Line();
  ArgList header;
  header.SetList(ptr+1, SEPARATORS );
  int nframes = header.getKeyInt("nframes", -1);
  // Need to keep track of frame indices so we can check for sieving.
  std::vector<char> sieveStatus;
  if (nframes > 0)
    sieveStatus.assign(nframes, 'T');
  // Keep track of matrix values.
  std::vector<float> Vals;
  // Read file
  bool checkSieve = true;
  int f1 = -1, f2 = -1, firstf1 = -1;
  float val = 0;
  while ( (ptr = buffer.Line()) != 0 )
  {
    if (checkSieve) {
      sscanf(ptr, "%i %i %f", &f1, &f2, &val);
      if (f2 > (int)sieveStatus.size())
        sieveStatus.resize(f2, 'T');
      if (firstf1 == -1) {
        // First values.
        sieveStatus[f1-1] = 'F';
        sieveStatus[f2-1] = 'F';
        firstf1 = f1;
      } else if (f1 > firstf1) {
          checkSieve = false;
      } else {
        sieveStatus[f2-1] = 'F';
      }
    } else {
      sscanf(ptr, "%*i %*i %f", &val);
    }
    Vals.push_back( val );
  }
  // DEBUG
  //mprintf("Sieved array:\n");
  //for (unsigned int i = 0; i < sieveStatus.size(); i++)
  //  mprintf("\t%6u %c\n", i+1, sieveStatus[i]);
  // Try to determine if sieve is random or not.
  int sieveDelta = 1;
  f1 = -1;
  f2 = -1;
  int actual_nrows = 0;
  for (int i = 0; i < (int)sieveStatus.size(); i++) {
    if (sieveStatus[i] == 'F') {
      actual_nrows++;
      if (sieveDelta != -2) {
        if (f1 == -1) {
          f1 = i;
        } else if (f2 == -1) {
          sieveDelta = i - f1;
          f1 = i;
          f2 = i;
        } else {
          int newDelta = i - f1;
          if (newDelta != sieveDelta) {
            // Random. No need to calculate sieveDelta anymore.
            sieveDelta = -2;
          }
          f1 = i;
        }
      }
    }
  }
  if (sieveDelta == -2) {
    // Random sieve. Try to figure out original sieve value.
    int o_frames = (int)sieveStatus.size();
    int o_sieve_value = o_frames / actual_nrows;
    if ( (o_frames % actual_nrows) != 0 )
      o_sieve_value++;
    sieveDelta = -o_sieve_value;
  }
  if (debug_ > 0)
    mprintf("DEBUG: sieve %i, actual_nrows= %i\n", sieveDelta, actual_nrows);
  if (sieveDelta != 1 && nframes == -1)
    mprintf("Warning: Pairwise distance matrix file contains sieved frames but\n"
            "Warning:   number of original frames is not present in file - this\n"
            "Warning:   may lead to ignored frames in cluster output. Please add\n"
            "Warning:   'nframes <# original frames>' to the pairwise distance\n"
            "Warning:   matrix file header, e.g. '#F1 F2 pw.dat nframes 1000'.\n");
  
  // Save cluster matrix
  if (Mat.Allocate( DataSet::SizeArray(1, actual_nrows) )) return 1;
//.........这里部分代码省略.........