C++ DataFile::SetupDatafile方法代码示例

/** Create new DataFile, or return existing DataFile. */
DataFile* DataFileList::AddDataFile(FileName const& nameIn, ArgList& argIn,
                                    DataFile::DataFormatType typeIn)
{
  // If no filename, no output desired
  if (nameIn.empty()) return 0;
  FileName fname( nameIn );
  // Append ensemble number if set.
  //rprintf("DEBUG: Setting up data file '%s' with ensembleNum %i\n", nameIn.base(), ensembleNum_);
  if (ensembleNum_ != -1)
    fname.Append( "." + integerToString(ensembleNum_) );
  // Check if filename in use by CpptrajFile.
  CpptrajFile* cf = GetCpptrajFile(fname);
  if (cf != 0) {
    mprinterr("Error: Data file name '%s' already in use by text output file '%s'.\n",
              fname.full(), cf->Filename().full());
    return 0;
  }
  // Check if this filename already in use
  DataFile* Current = GetDataFile(fname);
  // If no DataFile associated with name, create new DataFile
  if (Current==0) {
    Current = new DataFile();
    if (Current->SetupDatafile(fname, argIn, typeIn, debug_)) {
      mprinterr("Error: Setting up data file %s\n", fname.full());
      delete Current;
      return 0;
    }
    fileList_.push_back(Current);
  } else {
    // Set debug level
    Current->SetDebug(debug_);
    // If a type was specified, make sure it matches.
    if (typeIn != DataFile::UNKNOWN_DATA && typeIn != Current->Type()) {
      mprinterr("Error: '%s' is type %s but has been requested as type %s.\n",
                Current->DataFilename().full(), Current->FormatString(),
                DataFile::FormatString( typeIn ));
      return 0;
    }
    // Check for keywords that do not match file type
    DataFile::DataFormatType kType = DataFile::GetFormatFromArg( argIn );
    if (kType != DataFile::UNKNOWN_DATA && kType != Current->Type())
      mprintf("Warning: %s is type %s but type %s keyword specified; ignoring keyword.\n",
              Current->DataFilename().full(), Current->FormatString(),
              DataFile::FormatString( kType ));
    // Process Arguments
    if (!argIn.empty())
      Current->ProcessArgs( argIn );
  }
  return Current;
}

void Action_Spam::Print() {
  // Print the spam info file if we didn't do pure water
  if (!purewater_) {
    // Warn about any overflows
    if (overflow_)
      mprinterr("Warning: SPAM: Some frames had a box too small for the cutoff.\n");

    // Print information about each missing peak
    infofile_->Printf("# There are %d density peaks and %d frames\n\n",
                (int)peaks_.size(), Nframes_);
    // Loop over every Data set
    for (unsigned int i = 0; i < peakFrameData_.size(); i++) {
      // Skip peaks with 0 unoccupied sites
      if (peakFrameData_[i].size() == 0) continue;
      // Find out how many double-occupied frames there are
      int ndouble = 0;
      for (unsigned int j = 0; j < peakFrameData_[i].size(); j++)
        if (peakFrameData_[i][j] < 0) ndouble++;
      infofile_->Printf("# Peak %u has %d omitted frames (%d double-occupied)\n",
                  i, peakFrameData_[i].size(), ndouble);
      for (unsigned int j = 0; j < peakFrameData_[i].size(); j++) {
        if (j > 0 && j % 10 == 0) infofile_->Printf("\n");
        infofile_->Printf("%7d ", peakFrameData_[i][j]);
      }
      infofile_->Printf("\n\n");
    }
  }

  // Print the summary file with the calculated SPAM energies
  if (!summaryfile_.empty()) {
    // Not enabled yet -- just print out the data files.
    mprinterr("Warning: SPAM: SPAM calculation not yet enabled.\n");
    if (datafile_.empty()) datafile_ = summaryfile_;
  }
  // Now print the energy data
  if (!datafile_.empty()) {
    // Now write the data file with all of the SPAM energies
    DataFile dfl;
    ArgList dummy;
    dfl.SetupDatafile(datafile_, dummy, 0);
    for (int i = 0; i < (int)myDSL_.size(); i++) {
      dfl.AddDataSet(myDSL_[i]);
    }
    dfl.WriteDataOut();
  }
}

// Cluster_DBSCAN::ComputeKdistMap()
void Cluster_DBSCAN::ComputeKdistMap( Range const& Kvals, 
                                      std::vector<int> const& FramesToCluster ) const
{
  int pt1_idx, pt2_idx, d_idx, point;
  mprintf("\tCalculating Kdist map for %s\n", Kvals.RangeArg());
  double* kdist_array; // Store distance of pt1 to every other point.
  int nframes = (int)FramesToCluster.size();
  // Ensure all Kdist points are within proper range
  Range::const_iterator kval;
  for (kval = Kvals.begin(); kval != Kvals.end(); ++kval)
    if (*kval < 1 || *kval >= nframes) {
      mprinterr("Error: Kdist value %i is out of range (1 <= Kdist < %i)\n",
                 *kval, nframes);
      return;
    }
  int nvals = (int)Kvals.Size();
  double** KMAP; // KMAP[i] has the ith nearest point for each point.
  KMAP = new double*[ nvals ];
  for (int i = 0; i != nvals; i++)
    KMAP[i] = new double[ nframes ];
  ParallelProgress progress( nframes );
# ifdef _OPENMP
# pragma omp parallel private(pt1_idx, pt2_idx, d_idx, kval, point, kdist_array) firstprivate(progress)
  {
  progress.SetThread( omp_get_thread_num() );
#endif
  kdist_array = new double[ nframes ];
# ifdef _OPENMP
# pragma omp for
# endif
  for (pt1_idx = 0; pt1_idx < nframes; pt1_idx++) // X
  {
    progress.Update( pt1_idx );
    point = FramesToCluster[pt1_idx];
    d_idx = 0;
    // Store distances from pt1 to pt2
    for (pt2_idx = 0; pt2_idx != nframes; pt2_idx++)
      kdist_array[d_idx++] = FrameDistances_.GetFdist(point, FramesToCluster[pt2_idx]);
    // Sort distances; will be smallest to largest
    std::sort( kdist_array, kdist_array + nframes );
    // Save the distance of specified nearest neighbors to this point.
    d_idx = 0;
    for (kval = Kvals.begin(); kval != Kvals.end(); ++kval) // Y
      KMAP[d_idx++][pt1_idx] = kdist_array[ *kval ];
  }
  delete[] kdist_array;
# ifdef _OPENMP
  } // END omp parallel
# endif
  progress.Finish();
  // Sort all of the individual kdist plots, smallest to largest.
  for (int i = 0; i != nvals; i++)
    std::sort(KMAP[i], KMAP[i] + nframes);
  // Save in matrix, largest to smallest.
  DataSet_MatrixDbl kmatrix;
  kmatrix.Allocate2D( FramesToCluster.size(), Kvals.Size() );
  for (int y = 0; y != nvals; y++) {
    for (int x = nframes - 1; x != -1; x--)
      kmatrix.AddElement( KMAP[y][x] );
    delete[] KMAP[y];
  }
  delete[] KMAP;
  // Write matrix to file
  DataFile outfile;
  ArgList outargs("usemap");
  outfile.SetupDatafile(k_prefix_ + "Kmatrix.gnu", outargs, debug_);
  outfile.AddDataSet( (DataSet*)&kmatrix );
  outfile.WriteDataOut();
  // Write out the largest and smallest values for each K.
  // This means for each value of K the point with the furthest Kth-nearest
  // neighbor etc.
  CpptrajFile maxfile;
  if (maxfile.OpenWrite(k_prefix_ + "Kmatrix.max.dat")) return;
  maxfile.Printf("%-12s %12s %12s\n", "#Kval", "MaxD", "MinD");
  d_idx = 0;
  for (kval = Kvals.begin(); kval != Kvals.end(); ++kval, d_idx++)
    maxfile.Printf("%12i %12g %12g\n", *kval, kmatrix.GetElement(0, d_idx),
                   kmatrix.GetElement(nframes-1, d_idx));
  maxfile.CloseFile();
}