diff --git a/README.md b/README.md index a9a6b0f99b..185975b7a8 100644 --- a/README.md +++ b/README.md @@ -251,3 +251,5 @@ External code/libraries bundled with CPPTRAJ * CPPTRAJ uses code for the [permuted congruent pseudo-random number generator](https://www.pcg-random.org/index.html) PCG32 by Melissa O'Neill and the [Xoshiro 128++ pseudo-random number generator](http://prng.di.unimi.it) by David Blackman and Sebastino Vigna. * The code for quaternion RMSD calculation was adapted from code in [qcprot.c](https://theobald.brandeis.edu/qcp/qcprot.c) originally written by Douglas L. Theobald and Pu Lio (Brandeis University). + +* The code for reading numpy arrays in `src/libnpy` is from [libnpy](https://github.com/llohse/libnpy) written by Leon Merten Lohse et al. (Universität Göttingen). diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx index 6c06584d01..43128a7261 100644 --- a/doc/cpptraj.lyx +++ b/doc/cpptraj.lyx @@ -6253,7 +6253,7 @@ The following COORDS data set commands are available: \begin_layout Standard \align center \begin_inset Tabular - + @@ -6365,6 +6365,26 @@ Write a COORDS set to a file. \begin_inset Text +\begin_layout Plain Layout +crdtransform +\end_layout + +\end_inset + + +\begin_inset Text + +\begin_layout Plain Layout +Transform a COORDS set in one of several ways. +\end_layout + +\end_inset + + + + +\begin_inset Text + \begin_layout Plain Layout createcrd \end_layout @@ -6405,6 +6425,27 @@ Run simple energy minimization on a frame of a COORDS set. \begin_inset Text +\begin_layout Plain Layout +extendedcomp +\end_layout + +\end_inset + + +\begin_inset Text + +\begin_layout Plain Layout +Calculate extended comparison similarity values for each frame in COORDS + set. +\end_layout + +\end_inset + + + + +\begin_inset Text + \begin_layout Plain Layout graft \end_layout @@ -6726,6 +6767,147 @@ crd1 crdout crd1 crd1.pdb multi crdframes 1,10 \end_layout +\begin_layout Subsection +crdtransform +\end_layout + +\begin_layout LyX-Code +crdtransform [name ] +\end_layout + +\begin_layout LyX-Code + { rmsrefine [mask ] [mass] [rmstol ] | +\end_layout + +\begin_layout LyX-Code + normcoords | +\end_layout + +\begin_layout LyX-Code + trim [metric ] [{ntrimmed <#>|cutoff }] +\end_layout + +\begin_layout LyX-Code + [criterion {comp|medoid}]] +\end_layout + +\begin_layout LyX-Code + } +\end_layout + +\begin_deeper +\begin_layout Description + COORDS set to transform. +\end_layout + +\begin_layout Description +[name +\begin_inset space ~ +\end_inset + +] COORDS set to create; if not specified will be modified. +\end_layout + +\begin_layout Description +rmsrefine Do iterative RMS refinement. +\end_layout + +\begin_deeper +\begin_layout Description +[mask +\begin_inset space ~ +\end_inset + +] Mask of atoms to fit during refinement. +\end_layout + +\begin_layout Description +[mass] Mass-weight the refinement. +\end_layout + +\begin_layout Description +[rmstol +\begin_inset space ~ +\end_inset + +] Tolerance (in Ang.) below which RMS-refinement will stop. +\end_layout + +\end_deeper +\begin_layout Description +normcoords Normalize coordinates between 0.0 and 1.0 using the minimum and + maximum coordinate values. +\end_layout + +\begin_layout Description +trim Remove trajectory frames using extended similarity metrics. +\end_layout + +\begin_deeper +\begin_layout Description +[metric +\begin_inset space ~ +\end_inset + +] Metric to use; default MSD. +\end_layout + +\begin_layout Description +[{ntrimmed +\begin_inset space ~ +\end_inset + +<#>|cutoff +\begin_inset space ~ +\end_inset + +} # of frames or fraction of trajectory to trim. +\end_layout + +\begin_layout Description +[criterion +\begin_inset space ~ +\end_inset + +{comp|medoid}] Trim frames by comparitive similarity (i.e. + trim most dissimilar frames) or comparison to medoid (i.e. + trim most dissimilar to medoid frame). +\end_layout + +\end_deeper +\end_deeper +\begin_layout Standard +Transform a COORDS set in one of several ways. + Does not yet work with TRAJ data sets. + The iterative RMS refinement is similar to the procedure outlined by Klem + et al. + (J. + Chem. + Theory Comput. + 2022, 18, 3218−3230). + The extended similarity metrics are those defined by Racz et al. + (J. + Comp.-Aid. + Mol. + Design, 2022, 36, 157-173). +\end_layout + \begin_layout Subsection createcrd \end_layout @@ -6947,6 +7129,359 @@ Perform steepest descent minimization on a frame in a COORDS set using a A simple nonbonded term can be added as well if desired. \end_layout +\begin_layout Subsection +extendedcomp +\begin_inset CommandInset label +LatexCommand label +name "subsec:cpptraj-extendedcomp" + +\end_inset + + +\end_layout + +\begin_layout LyX-Code +extendedcomp [name ] +\end_layout + +\begin_layout LyX-Code + [metric ] [out ] +\end_layout + +\begin_layout LyX-Code + = msd bub fai gle ja jt rt rr sm ss1 ss2 +\end_layout + +\begin_deeper +\begin_layout Description + Input COORDS set. +\end_layout + +\begin_layout Description +[name +\begin_inset space ~ +\end_inset + +] Output data set name containing values. +\end_layout + +\begin_layout Description +[metric +\begin_inset space ~ +\end_inset + +] Metric to use. +\end_layout + +\begin_layout Description +[out +\begin_inset space ~ +\end_inset + +,file>] File to write values to. +\end_layout + +\begin_layout Standard +DataSets generated: +\end_layout + +\begin_layout Description + Set containing similarity values for each COORDS frame. +\end_layout + +\end_deeper +\begin_layout Standard +Calculate extended comparison similarity values for each frame in a COORDS + set. + The extended similarity metrics are those defined by Racz et al. + (J. + Comp.-Aid. + Mol. + Design, 2022, 36, 157-173). + The metrics are as follows: +\end_layout + +\begin_layout Standard +\begin_inset Tabular + + + + + + +\begin_inset Text + +\begin_layout Plain Layout + +\series bold +Keyword +\end_layout + +\end_inset + + +\begin_inset Text + +\begin_layout Plain Layout + +\series bold +Metric +\end_layout + +\end_inset + + + + +\begin_inset Text + +\begin_layout Plain Layout +msd +\end_layout + +\end_inset + + +\begin_inset Text + +\begin_layout Plain Layout +Mean-squared deviation. +\end_layout + +\end_inset + + + + +\begin_inset Text + +\begin_layout Plain Layout +bub +\end_layout + +\end_inset + + +\begin_inset Text + +\begin_layout Plain Layout +Bhattacharyya's U coefficient +\end_layout + +\end_inset + + + + +\begin_inset Text + +\begin_layout Plain Layout +fai +\end_layout + +\end_inset + + +\begin_inset Text + +\begin_layout Plain Layout +Faiman's coefficient +\end_layout + +\end_inset + + + + +\begin_inset Text + +\begin_layout Plain Layout +gle +\end_layout + +\end_inset + + +\begin_inset Text + +\begin_layout Plain Layout +Gleason's coefficient +\end_layout + +\end_inset + + + + +\begin_inset Text + +\begin_layout Plain Layout +ja +\end_layout + +\end_inset + + +\begin_inset Text + +\begin_layout Plain Layout +Jaccard's coefficient +\end_layout + +\end_inset + + + + +\begin_inset Text + +\begin_layout Plain Layout +jt +\end_layout + +\end_inset + + +\begin_inset Text + +\begin_layout Plain Layout +Jaccard-Tanimoto coefficient +\end_layout + +\end_inset + + + + +\begin_inset Text + +\begin_layout Plain Layout +rt +\end_layout + +\end_inset + + +\begin_inset Text + +\begin_layout Plain Layout +Rogers-Tanimoto coefficient +\end_layout + +\end_inset + + + + +\begin_inset Text + +\begin_layout Plain Layout +rr +\end_layout + +\end_inset + + +\begin_inset Text + +\begin_layout Plain Layout +Russell-Rao coefficient +\end_layout + +\end_inset + + + + +\begin_inset Text + +\begin_layout Plain Layout +sm +\end_layout + +\end_inset + + +\begin_inset Text + +\begin_layout Plain Layout +Simpson's coefficient +\end_layout + +\end_inset + + + + +\begin_inset Text + +\begin_layout Plain Layout +ss1 +\end_layout + +\end_inset + + +\begin_inset Text + +\begin_layout Plain Layout +Sokal-Sneath 1 coefficient +\end_layout + +\end_inset + + + + +\begin_inset Text + +\begin_layout Plain Layout +ss2 +\end_layout + +\end_inset + + +\begin_inset Text + +\begin_layout Plain Layout +Sokal-Sneath 2 coefficient +\end_layout + +\end_inset + + + + +\end_inset + + +\end_layout + \begin_layout Subsection graft \end_layout diff --git a/src/Action_Keep.cpp b/src/Action_Keep.cpp index a70016acbf..1a4bba2e88 100644 --- a/src/Action_Keep.cpp +++ b/src/Action_Keep.cpp @@ -257,6 +257,7 @@ Action::RetType Action_Keep::DoAction(int frameNum, ActionFrame& frm) if (err == Action::OK) { keepFrame_.SetFrame(frm.Frm(), atomsToKeep_); frm.SetFrame( &keepFrame_ ); + err = Action::MODIFY_COORDS; } return err; diff --git a/src/Command.cpp b/src/Command.cpp index 158e2d5fb2..2311a4d316 100644 --- a/src/Command.cpp +++ b/src/Command.cpp @@ -54,6 +54,8 @@ #include "Exec_RotateDihedral.h" #include "Exec_SplitCoords.h" #include "Exec_CatCrd.h" +#include "Exec_CrdTransform.h" +#include "Exec_ExtendedComparison.h" // ----- TRAJECTORY ------------------------------------------------------------ #include "Exec_Traj.h" // ----- TOPOLOGY -------------------------------------------------------------- @@ -268,7 +270,9 @@ void Command::Init() { Command::AddCmd( new Exec_CombineCoords(), Cmd::EXE, 1, "combinecrd" ); Command::AddCmd( new Exec_CrdAction(), Cmd::EXE, 1, "crdaction" ); Command::AddCmd( new Exec_CrdOut(), Cmd::EXE, 1, "crdout" ); - Command::AddCmd( new Exec_Emin(), Cmd::EXE, 1, "emin"); // hidden + Command::AddCmd( new Exec_CrdTransform(), Cmd::EXE, 1, "crdtransform" ); + Command::AddCmd( new Exec_Emin(), Cmd::EXE, 1, "emin"); + Command::AddCmd( new Exec_ExtendedComparison(),Cmd::EXE,1, "extendedcomp" ); Command::AddCmd( new Exec_Graft(), Cmd::EXE, 1, "graft"); Command::AddCmd( new Exec_LoadCrd(), Cmd::EXE, 1, "loadcrd" ); Command::AddCmd( new Exec_LoadTraj(), Cmd::EXE, 1, "loadtraj" ); diff --git a/src/DataFile.cpp b/src/DataFile.cpp index 18d599d975..32c6c44ee2 100644 --- a/src/DataFile.cpp +++ b/src/DataFile.cpp @@ -31,6 +31,7 @@ #include "DataIO_Peaks.h" #include "DataIO_NetCDF.h" #include "DataIO_AmberEne.h" +#include "DataIO_Numpy.h" // CONSTRUCTOR DataFile::DataFile() : @@ -83,6 +84,7 @@ const FileTypes::AllocToken DataFile::DF_AllocArray[] = { { "NetCDF data", 0, 0, 0 }, # endif { "Amber Energy File", 0, 0, DataIO_AmberEne::Alloc}, + { "Numpy array", 0, 0, DataIO_Numpy::Alloc}, { "Unknown Data file", 0, 0, 0 } }; @@ -108,6 +110,7 @@ const FileTypes::KeyToken DataFile::DF_KeyArray[] = { { CMATRIX_NETCDF,"nccmatrix", ".nccmatrix" }, { PEAKS, "peaks", ".peaks" }, { AMBERENE, "amberene", ".ene" }, + { NUMPY, "numpy", ".npy" }, { UNKNOWN_DATA, 0, 0 } }; diff --git a/src/DataFile.h b/src/DataFile.h index 4b72424e61..4571381644 100644 --- a/src/DataFile.h +++ b/src/DataFile.h @@ -18,7 +18,7 @@ class DataFile { DATAFILE=0, XMGRACE, GNUPLOT, XPLOR, OPENDX, REMLOG, MDOUT, EVECS, VECTRAJ, XVG, CCP4, CHARMMREPD, CHARMMFASTREP, CHARMMOUT, CPOUT, CHARMMRTFPRM, CMATRIX_BINARY, CMATRIX_NETCDF, PEAKS, - NETCDFDATA, AMBERENE, + NETCDFDATA, AMBERENE, NUMPY, UNKNOWN_DATA }; DataFile(); diff --git a/src/DataIO_Numpy.cpp b/src/DataIO_Numpy.cpp new file mode 100644 index 0000000000..2c6f04cf82 --- /dev/null +++ b/src/DataIO_Numpy.cpp @@ -0,0 +1,157 @@ +#include "DataIO_Numpy.h" +#include "CpptrajStdio.h" + +#include "DataSet_Coords.h" + +#include "libnpy/npy.hpp" + +/// CONSTRUCTOR +DataIO_Numpy::DataIO_Numpy() +{ + +} + +// DataIO_Numpy::ID_DataFormat() +bool DataIO_Numpy::ID_DataFormat(CpptrajFile& infile) +{ + if (infile.OpenFile()) return false; + // Read first 6 bytes + unsigned char magic[6]; + magic[0] = 0; + magic[1] = 0; + magic[2] = 0; + magic[3] = 0; + magic[4] = 0; + magic[5] = 0; + infile.Read(magic, 6); + infile.CloseFile(); + if (magic[0] == 93 && magic[1] == 'N' && magic[2] == 'U' && + magic[3] == 'M' && magic[4] == 'P' && magic[5] == 'Y') + return true; + return false; +} + +// DataIO_Numpy::ReadHelp() +void DataIO_Numpy::ReadHelp() +{ + +} + +// DataIO_Numpy::processReadArgs() +int DataIO_Numpy::processReadArgs(ArgList& argIn) +{ + + return 0; +} + +/** Read 2d numpy array as COORDS set */ +int DataIO_Numpy::read_data_as_coords(std::string const& dsname, DataSetList& dsl, + std::vector const& data, + unsigned long nframes, + unsigned long ncoords) +const +{ + if (ncoords % 3 != 0) { + mprinterr("Internal Error: DataIO_Numpy::read_data_as_coords(): # coords %lu not divisible by 3.\n", ncoords); + return 1; + } + unsigned long natoms = ncoords / 3; + + DataSet* ds = dsl.AddSet( DataSet::COORDS, MetaData(dsname) ); + if (ds == 0) { + mprinterr("Error: Could not allocate COORDS set with name '%s'\n", dsname.c_str()); + return 1; + } + DataSet_Coords* coords = static_cast( ds ); + + // Create a pseudo topology + Topology top; + for (unsigned long iat = 0; iat != natoms; iat++) + top.AddTopAtom( Atom("CA","C"), Residue("XXX", 1, ' ', ' ') ); + top.SetSingleMolecule(); + top.CommonSetup(false, false); + top.Summary(); + + // Set up COORDS set + if (coords->CoordsSetup( top, CoordinateInfo(Box(), false, false, false) ) ) { + mprinterr("Error: Could not set up COORDS set '%s'\n", coords->legend()); + return 1; + } + if (coords->Allocate(DataSet::SizeArray(1, nframes))) { + mprinterr("Error: Could not allocate COORDS set '%s'\n", coords->legend()); + return 1; + } + Frame frm = coords->AllocateFrame(); + + std::vector::const_iterator it = data.begin(); + for (unsigned long ifrm = 0; ifrm != nframes; ifrm++) { + frm.ClearAtoms(); + for (unsigned long iat = 0; iat != natoms; iat++) { + double xyz[3]; + xyz[0] = *(it++); + xyz[1] = *(it++); + xyz[2] = *(it++); + frm.AddXYZ( xyz ); + } + coords->AddFrame( frm ); + } + + + return 0; +} + + +// DataIO_Numpy::ReadData() +int DataIO_Numpy::ReadData(FileName const& fname, DataSetList& dsl, std::string const& dsname) +{ + mprintf("\tReading numpy array from '%s'\n", fname.full()); + npy::npy_data dataIn = npy::read_npy( fname.Full() ); + + //std::vector data = dataIn.data; + //std::vector shape = dataIn.shape; + //bool fortran_order = dataIn.fortran_order; + + mprintf("\tData array size = %zu\n", dataIn.data.size()); + mprintf("\tShape array size = %zu\n", dataIn.shape.size()); + for (std::vector::const_iterator it = dataIn.shape.begin(); + it != dataIn.shape.end(); ++it) + mprintf("\t\t%lu\n", *it); + mprintf("\tFortran order = %i\n", (int)dataIn.fortran_order); + + // Convert to coordinates. Expect 2 dimensions, frame and ncoords + if (dataIn.shape.size() != 2) { + mprinterr("Error: Shape of numpy array is not 2 (%zu)\n", dataIn.shape.size()); + return 1; + } + if (dataIn.fortran_order) { + mprinterr("Error: Cannot yet process numpy array in fortran order.\n"); + return 1; + } + int err = read_data_as_coords(dsname, dsl, dataIn.data, dataIn.shape[0], dataIn.shape[1]); + if (err != 0) { + mprinterr("Error: Could not convert numpy array to COORDS.\n"); + return 1; + } + + return 0; +} + +// DataIO_Numpy::WriteHelp() +void DataIO_Numpy::WriteHelp() +{ + +} + +// DataIO_Numpy::processWriteArgs() +int DataIO_Numpy::processWriteArgs(ArgList& argIn) +{ + + return 0; +} + +// DataIO_Numpy::WriteData() +int DataIO_Numpy::WriteData(FileName const& fname, DataSetList const& dsl) +{ + mprintf("\tWriting to numpy array '%s'\n", fname.full()); + return 1; +} diff --git a/src/DataIO_Numpy.h b/src/DataIO_Numpy.h new file mode 100644 index 0000000000..f2ed55eb34 --- /dev/null +++ b/src/DataIO_Numpy.h @@ -0,0 +1,22 @@ +#ifndef INC_DATAIO_NUMPY_H +#define INC_DATAIO_NUMPY_H +#include "DataIO.h" +/// +class DataIO_Numpy : public DataIO { + public: + DataIO_Numpy(); + static void ReadHelp(); + static void WriteHelp(); + static BaseIOtype* Alloc() { return (BaseIOtype*)new DataIO_Numpy(); } + int processReadArgs(ArgList&); + int ReadData(FileName const&, DataSetList&, std::string const&); + int processWriteArgs(ArgList&); + int WriteData(FileName const&, DataSetList const&); + bool ID_DataFormat(CpptrajFile&); + private: + int read_data_as_coords(std::string const&, DataSetList&, + std::vector const&, + unsigned long, unsigned long) const; + +}; +#endif diff --git a/src/Exec_CrdTransform.cpp b/src/Exec_CrdTransform.cpp new file mode 100644 index 0000000000..475d98a1e1 --- /dev/null +++ b/src/Exec_CrdTransform.cpp @@ -0,0 +1,422 @@ +#include "Exec_CrdTransform.h" +#include "CpptrajStdio.h" +#include "Constants.h" +#include // std::max,min,sort +#include // floor + +/// Get the minimum and maximum coordinates in a given frame, store in min and max +static inline void getMaxMin(Frame const& frmIn, Vec3& max, Vec3& min) { + for (int at = 0; at != frmIn.Natom(); at++) { + const double* xyz = frmIn.XYZ(at); + for (int i = 0; i != 3; i++) { + max[i] = std::max(max[i], xyz[i]); + min[i] = std::min(min[i], xyz[i]); + } + } +} + +/** Normalize coordinates between 0 and 1. */ +int Exec_CrdTransform::normalizeCoords(DataSet_Coords* crdIn, + DataSet_Coords* crdOut) +const +{ + mprintf("\tNormalize coordinates between 0 and 1.\n"); + mprintf("\tInput coords: %s\n", crdIn->legend()); + mprintf("\tOutput coords: %s\n", crdOut->legend()); + // Get max and min X Y and Z + Frame frmIn = crdIn->AllocateFrame(); + crdIn->GetFrame(0, frmIn); + Vec3 xyzmax( frmIn.XYZ(0) ); + Vec3 xyzmin( frmIn.XYZ(0) ); + getMaxMin(frmIn, xyzmax, xyzmin); + for (unsigned int idx = 1; idx < crdIn->Size(); idx++) { + crdIn->GetFrame(idx, frmIn); + getMaxMin(frmIn, xyzmax, xyzmin); + } + mprintf("\tMax: %g %g %g\n", xyzmax[0], xyzmax[1], xyzmax[2]); + mprintf("\tMin: %g %g %g\n", xyzmin[0], xyzmin[0], xyzmin[0]); + // Choose the overall max and min + double max = xyzmax[0]; + max = std::max(max, xyzmax[1]); + max = std::max(max, xyzmax[2]); + double min = xyzmin[0]; + min = std::min(min, xyzmin[1]); + min = std::min(min, xyzmin[2]); + + Vec3 norm( max - min ); + // Protect against bad values + bool hasBadValues = false; + static const char dirStr[] = { 'X', 'Y', 'Z' }; + for (int ii = 0; ii != 3; ii++) { + if (norm[ii] < 0) { + mprinterr("Error: Min value > max value for %c coordinate.\n", dirStr[ii]); + hasBadValues = true; + } + if (norm[ii] < Constants::SMALL) { + mprinterr("Error: Min value == max value for %c coordinate.\n", dirStr[ii]); + hasBadValues = true; + } + } + if (hasBadValues) return 1; + // Transform coords between 0 and 1 + for (unsigned int idx = 0; idx < crdIn->Size(); idx++) { + crdIn->GetFrame(idx, frmIn); + for (int crdidx = 0; crdidx < frmIn.size(); crdidx+=3) { + frmIn[crdidx ] = (frmIn[crdidx ] - min) / norm[0]; + frmIn[crdidx+1] = (frmIn[crdidx+1] - min) / norm[1]; + frmIn[crdidx+2] = (frmIn[crdidx+2] - min) / norm[2]; + } + crdOut->SetCRD(idx, frmIn); + } + + return 0; +} + +/** Transform coordinates by RMS-fitting to an average structure, calculating + * a new average, then RMS-fitting to that average and so on until a + * tolerance is reached. Essentially the procedure described by + * Klem et al. J. Chem. Theory Comput. 2022, 18, 3218−3230. + */ +int Exec_CrdTransform::iterativeRmsRefinement(AtomMask const& maskIn, + bool useMass, + double tolIn, + DataSet_Coords* crdIn, + DataSet_Coords* crdOut) +const +{ + mprintf("\tRMS iterative refinement.\n"); + mprintf("\tInput coords: %s\n", crdIn->legend()); + mprintf("\tOutput coords: %s\n", crdOut->legend()); + mprintf("\tAtom mask: %s\n", maskIn.MaskString()); + mprintf("\tRMS Tolerance: %g Ang.\n", tolIn); + if (useMass) + mprintf("\tMass-weighting on.\n"); + else + mprintf("\tMass-weighting off.\n"); + // Do the initial fit to the first frame. + Frame frmIn = crdIn->AllocateFrame(); + crdIn->GetFrame(0, frmIn); + Frame selectedRef; + selectedRef.SetupFrameFromMask( maskIn, crdIn->Top().Atoms() ); + selectedRef.SetCoordinates( frmIn, maskIn ); + // Ensure reference is centered on the origin + Vec3 refTrans = selectedRef.CenterOnOrigin( useMass ); + // Set up frame for selected incoming atoms + Frame selectedTgt = selectedRef; + // Set up frame to hold average + Frame avgFrm = selectedTgt; + + double currentTol = tolIn + 9999.0; + + unsigned int iteration = 0; + while (currentTol > tolIn) { + avgFrm.ZeroCoords(); + Vec3 tgtTrans(0.0); + Matrix_3x3 rot(0.0); + for (unsigned int idx = 0; idx != crdIn->Size(); idx++) { + crdIn->GetFrame(idx, frmIn); + selectedTgt.SetCoordinates( frmIn, maskIn ); + selectedTgt.RMSD_CenteredRef( selectedRef, rot, tgtTrans, useMass ); + frmIn.Trans_Rot_Trans(tgtTrans, rot, refTrans); + crdOut->SetCRD(idx, frmIn); + avgFrm.AddByMask( frmIn, maskIn ); + } + avgFrm.Divide( (double)crdIn->Size() ); + // Calc RMS of current average to current reference + currentTol = avgFrm.RMSD_CenteredRef( selectedRef, rot, tgtTrans, useMass ); + mprintf("\t%8u %12.4f\n", iteration+1, currentTol); + // Fit the current average TODO is this necessary? + avgFrm.Trans_Rot_Trans(tgtTrans, rot, refTrans); + // Set current average to be new reference + selectedRef = avgFrm; + iteration++; + } + + return 0; +} + +/** Strings corresponding to CriterionType */ +const char* Exec_CrdTransform::CriterionStr_[] = { + "comp_sim", "sim_to_medoid", "No criterion" +}; + + +/** Trim a desired percentage of outliers (most dissimilar) from the COORDS + * data set by calculating the largest complement similarity. + */ +int Exec_CrdTransform::trimOutliers(int n_trimmed, double cutoffIn, + ExtendedSimilarity::MetricType metric, + CriterionType criterion, + DataSet_Coords* crdIn, + DataSet_Coords* crdOut) +const +{ + mprintf("\tTrimming outliers.\n"); + mprintf("\tInput coords: %s\n", crdIn->legend()); + mprintf("\tOutput coords: %s\n", crdOut->legend()); + mprintf("\tUsing metric: %s\n", ExtendedSimilarity::metricStr(metric)); + mprintf("\tCriterion: %s\n", CriterionStr_[criterion]); + unsigned int Ncoords = crdIn->Top().Natom() * 3; + unsigned int Nframes = crdIn->Size(); + Frame frmIn = crdIn->AllocateFrame(); + mprintf("\t'%s' has %u coordinates, %u frames.\n", crdIn->legend(), Ncoords, Nframes); + if (Nframes < 2) { + mprintf("Warning: Less than 2 frames, nothing to trim.\n"); // TODO something with crdOut? + return 0; + } + // Specify n_trimmed or cutoff, but not both. + if (n_trimmed < 0 && cutoffIn < 0) { + mprinterr("Internal Error: Must specify either number to trim or cutoff.\n"); + return 1; + } + if (n_trimmed >= 0 && cutoffIn > 0) { + mprinterr("Error: Must specify either number to trim or cutoff, but not both.\n"); + return 1; + } + // Set cutoff value + unsigned int cutoff; + if (n_trimmed >= 0) { + cutoff = (unsigned int)n_trimmed; + mprintf("\t# to trim: %i\n", n_trimmed); + } else { + cutoff = (unsigned int)(floor(Nframes * cutoffIn)); + mprintf("\tFraction of outliers to remove: %f\n", cutoffIn); + } + mprintf("\tUsing cutoff value: %u\n", cutoff); + // Do extended similarity calculation for each frame + ExtendedSimilarity ExtSim; + if (ExtSim.SetOpts( metric, crdIn->Size(), Ncoords )) { + mprinterr("Error: Extended similarity setup for trim failed.\n"); + return 1; + } + ExtendedSimilarity::Darray csimvals = ExtSim.CalculateCompSim( *crdIn ); + if (csimvals.empty()) { + mprinterr("Error: No comparitive similarity values for trim calculated.\n"); + return 1; + } + // DEBUG +/* CpptrajFile dbg; + dbg.OpenWrite("test.cpptraj.out"); + for (ExtendedSimilarity::Darray::const_iterator it = csimvals.begin(); it != csimvals.end(); ++it) + dbg.Printf("%8li %16.8f\n", it - csimvals.begin(), *it); + dbg.CloseFile();*/ + // Array type for sorting by sim. value while preserving indices + typedef std::pair IdxValPairType; + std::vector comp_sims; + comp_sims.reserve( crdIn->Size() ); + // For sorting IdxValPairType by values + struct IdxValPairCmp { + inline bool operator()(IdxValPairType const& first, IdxValPairType const& second) { + return (first.second < second.second); + } + }; + struct ReverseIdxValPairCmp { + inline bool operator()(IdxValPairType const& first, IdxValPairType const& second) { + return (second.second < first.second); + } + }; + + if (criterion == COMP_SIM) { + // ----- Comp sim ------------------ + // Place values in sortable array and sort lowest to highest. + for (unsigned int idx = 0; idx < crdIn->Size(); idx++) + comp_sims.push_back( IdxValPairType(idx, csimvals[idx]) ); + std::sort(comp_sims.begin(), comp_sims.end(), IdxValPairCmp()); + } else if (criterion == SIM_TO_MEDOID) { + // ----- SIM TO MEDOID ------------- + mprintf("\tMedoid index= %li\n", ExtSim.MedoidIndex() + 1); + if (ExtSim.MedoidIndex() < 0) { + mprinterr("Error: Medoid index is < 0\n"); + return 1; + } + // Get comp sim values to medoid (excluding medoid) + unsigned int medoid_index = (unsigned int)ExtSim.MedoidIndex(); + Frame medoid = crdIn->AllocateFrame(); + crdIn->GetFrame(medoid_index, medoid); + //mprintf("["); + for (unsigned int idx = 0; idx < crdIn->Size(); idx++) { + if (idx != medoid_index) { + //mprintf("DEBUG\n"); + crdIn->GetFrame(idx, frmIn); + comp_sims.push_back( IdxValPairType(idx, ExtSim.CalculateCompSim(frmIn, medoid)) ); + //mprintf(" %10.8g", comp_sims.back().second); + } + } + //mprintf("]\n"); + // Sort highest to lowest + std::sort(comp_sims.begin(), comp_sims.end(), ReverseIdxValPairCmp()); + } + // Remove frames up to the cutoff + std::vector keepFrame(comp_sims.size(), true); + if (debug_ > 0) mprintf("["); + unsigned int nToRemove = 0; + for (unsigned int idx = 0; idx < cutoff; idx++) { + if (debug_ > 0) mprintf(" %u", comp_sims[idx].first+1); + keepFrame[comp_sims[idx].first] = false; + nToRemove++; + } + if (debug_ > 0) mprintf("]\n"); + mprintf("\tRemoving %u frames.\n", nToRemove); + + // Populate the output trajectory + for (unsigned int idx = 0; idx < crdIn->Size(); idx++) { + if (keepFrame[idx]) { + crdIn->GetFrame(idx, frmIn); + crdOut->AddFrame( frmIn ); + } + } + return 0; +} + + +// Exec_CrdTransform::Help() +void Exec_CrdTransform::Help() const +{ + mprintf("\t [name ]\n" + "\t{ rmsrefine [mask ] [mass] [rmstol ] |\n" + "\t normcoords |\n" + "\t trim [metric ] [{ntrimmed <#>|cutoff }]\n" + "\t [criterion {comp|medoid}]]\n" + "\t}\n"); + mprintf(" = %s\n", ExtendedSimilarity::MetricKeys().c_str()); + mprintf(" Transform a trajectory in one of several ways:\n" + " - rmsrefine : Do an iterative RMS refinement of all frames.\n" + " - normcoords : Normalize coordinates between 0.0 and 1.0.\n" + " - trim : Remove trajectory frames using extended similarity metrics.\n"); +} + +// Exec_CrdTransform::Execute() +Exec::RetType Exec_CrdTransform::Execute(CpptrajState& State, ArgList& argIn) +{ + debug_ = State.Debug(); + AtomMask mask; + bool useMass = false; + double rmsTol = -1.0; + + int n_trimmed = -1; + double cutoff = -1.0; + ExtendedSimilarity::MetricType metric = ExtendedSimilarity::NO_METRIC; + CriterionType criterion = NO_CRITERION; + // Determine mode + bool lengthWillBeModified = false; + enum ModeType { RMSREFINE = 0, NORMCOORDS, TRIM, UNSPECIFIED }; + ModeType mode = UNSPECIFIED; + if (argIn.hasKey("rmsrefine")) { + mode = RMSREFINE; + mask.SetMaskString( argIn.GetStringKey("mask") ); + useMass = argIn.hasKey("mass"); + rmsTol = argIn.getKeyDouble("rmstol", 0.0001); + } else if (argIn.hasKey("normcoords")) { + mode = NORMCOORDS; + } else if (argIn.hasKey("trim")) { + lengthWillBeModified = true; + mode = TRIM; + n_trimmed = argIn.getKeyInt("ntrimmed", -1); + cutoff = argIn.getKeyDouble("cutoff", -1.0); + std::string mstr = argIn.GetStringKey("metric"); + if (!mstr.empty()) { + metric = ExtendedSimilarity::TypeFromKeyword( mstr ); + if (metric == ExtendedSimilarity::NO_METRIC) { + mprinterr("Error: Metric '%s' not recognized.\n", mstr.c_str()); + return CpptrajState::ERR; + } + } else { + metric = ExtendedSimilarity::MSD; + } + std::string cstr = argIn.GetStringKey("criterion"); + if (!cstr.empty()) { + if (cstr == "comp") + criterion = COMP_SIM; + else if (cstr == "medoid") + criterion = SIM_TO_MEDOID; + else { + mprinterr("Error: Unrecognized criterion: %s\n", cstr.c_str()); + return CpptrajState::ERR; + } + } else { + criterion = COMP_SIM; + } + } else { + mprinterr("Error: Expected 'trim', 'rmsrefine', or 'normcoords'\n"); + return CpptrajState::ERR; + } + + // Get COORDS set + std::string outname = argIn.GetStringKey("name"); + std::string setname = argIn.GetStringNext(); + if (setname.empty()) { + mprinterr("Error: %s: Specify input COORDS dataset name.\n", argIn.Command()); + return CpptrajState::ERR; + } + DataSet_Coords* CRD = (DataSet_Coords*)State.DSL().FindSetOfGroup( setname, DataSet::COORDINATES ); + if (CRD == 0) { + mprinterr("Error: %s: No COORDS set with name %s found.\n", argIn.Command(), setname.c_str()); + return CpptrajState::ERR; + } + mprintf("\tUsing set '%s'\n", CRD->legend()); + if (CRD->Size() < 1) { + mprinterr("Error: Set '%s' has no frames.\n", CRD->legend()); + return CpptrajState::ERR; + } + if (CRD->Type() == DataSet::TRAJ) { + mprinterr("Error: TRAJ sets not yet supported.\n"); // FIXME + return CpptrajState::ERR; + } + // ----- No more input args processed below here --------- + + // Set up output set if needed TODO FRAMES set + DataSet_Coords* OUT = 0; + if (!outname.empty()) { + OUT = (DataSet_Coords*)State.DSL().AddSet( DataSet::COORDS, MetaData(outname) ); + if (OUT == 0) { + mprinterr("Error: Could not create output coords set %s\n", outname.c_str()); + return CpptrajState::ERR; + } + if (OUT->CoordsSetup( CRD->Top(), CRD->CoordsInfo() )) { + mprinterr("Error: Could not set up output coords set %s\n", OUT->legend()); + return CpptrajState::ERR; + } + if (!lengthWillBeModified) OUT->Allocate( DataSet::SizeArray(1, CRD->Size()) ); + } + bool needToDeleteCRD = false; + if (OUT == 0) { + if (!lengthWillBeModified) + OUT = CRD; + else { + // Need to replace CRD with a new set. Remove CRD from master DSL. + needToDeleteCRD = true; + State.DSL().PopSet( CRD ); + OUT = (DataSet_Coords*)State.DSL().AddSet( DataSet::COORDS, CRD->Meta() ); + if (OUT == 0) { + mprinterr("Error: Could not replace coords set %s\n", CRD->legend()); + return CpptrajState::ERR; + } + if (OUT->CoordsSetup( CRD->Top(), CRD->CoordsInfo() )) { + mprinterr("Error: Could not set up replacement coords set %s\n", OUT->legend()); + return CpptrajState::ERR; + } + } + } + + // Set up mask + if (mask.MaskStringSet()) { + if (CRD->Top().SetupIntegerMask( mask )) { + mprinterr("Error: Could not set up mask.\n"); + return CpptrajState::ERR; + } + mask.MaskInfo(); + } + + int err = 0; + switch (mode) { + case RMSREFINE : err = iterativeRmsRefinement(mask, useMass, rmsTol, CRD, OUT); break; + case NORMCOORDS : err = normalizeCoords(CRD, OUT); break; + // TODO pass in criterion + case TRIM : err = trimOutliers(n_trimmed, cutoff, metric, criterion, CRD, OUT); break; + default : err = 1; break; + } + if (needToDeleteCRD) delete CRD; + if (err != 0) return CpptrajState::ERR; + + return CpptrajState::OK; +} diff --git a/src/Exec_CrdTransform.h b/src/Exec_CrdTransform.h new file mode 100644 index 0000000000..096b60bd8b --- /dev/null +++ b/src/Exec_CrdTransform.h @@ -0,0 +1,27 @@ +#ifndef INC_EXEC_CRDTRANSFORM_H +#define INC_EXEC_CRDTRANSFORM_H +#include "Exec.h" +#include "ExtendedSimilarity.h" +/// Used to transform/condition coordinates +class Exec_CrdTransform : public Exec { + public: + Exec_CrdTransform() : Exec(COORDS), debug_(0) {} + void Help() const; + DispatchObject* Alloc() const { return (DispatchObject*)new Exec_CrdTransform(); } + RetType Execute(CpptrajState&, ArgList&); + private: + enum CriterionType { COMP_SIM = 0, ///< Remove most dissimilar objects based on complement similarity + SIM_TO_MEDOID, ///< Remove most dissimilar objects based on similarity to medoid. + NO_CRITERION }; + + static const char* CriterionStr_[]; + + int iterativeRmsRefinement(AtomMask const&, bool, double, + DataSet_Coords*, DataSet_Coords*) const; + int normalizeCoords(DataSet_Coords*, DataSet_Coords*) const; + int trimOutliers(int, double, ExtendedSimilarity::MetricType, CriterionType, + DataSet_Coords*, DataSet_Coords*) const; + + int debug_; +}; +#endif diff --git a/src/Exec_ExtendedComparison.cpp b/src/Exec_ExtendedComparison.cpp new file mode 100644 index 0000000000..a8b15438d3 --- /dev/null +++ b/src/Exec_ExtendedComparison.cpp @@ -0,0 +1,83 @@ +#include "Exec_ExtendedComparison.h" +#include "CpptrajStdio.h" +#include "ExtendedSimilarity.h" + +// Exec_ExtendedComparison::Help() +void Exec_ExtendedComparison::Help() const +{ + mprintf("\t [name ]\n" + "\t[metric ] [out ]\n"); + mprintf(" = %s\n", ExtendedSimilarity::MetricKeys().c_str()); + mprintf(" Calculate extended comparison similarity values for each trajectory frame.\n"); +} + +// Exec_ExtendedComparison::Execute() +Exec::RetType Exec_ExtendedComparison::Execute(CpptrajState& State, ArgList& argIn) +{ + // Metric + std::string mstr = argIn.GetStringKey("metric"); + ExtendedSimilarity::MetricType metric = ExtendedSimilarity::NO_METRIC; + if (!mstr.empty()) { + metric = ExtendedSimilarity::TypeFromKeyword( mstr ); + if (metric == ExtendedSimilarity::NO_METRIC) { + mprinterr("Error: Metric '%s' not recognized.\n", mstr.c_str()); + return CpptrajState::ERR; + } + } else { + metric = ExtendedSimilarity::MSD; + } + // Get output set + std::string outname = argIn.GetStringKey("name"); + DataSet* out = State.DSL().AddSet( DataSet::DOUBLE, outname, "EXTCOMP" ); + if (out == 0) { + mprinterr("Error: Could not set up output data set for extended comparison.\n"); + return CpptrajState::ERR; + } + out->SetDim(Dimension::X, Dimension(1.0, 1.0, "Frame") ); // TODO make time an option? + // Output file + DataFile* df = State.DFL().AddDataFile( argIn.GetStringKey("out"), argIn ); + if (df != 0) + df->AddDataSet( out ); + // Get COORDS set + std::string setname = argIn.GetStringNext(); + if (setname.empty()) { + mprinterr("Error: %s: Specify input COORDS dataset name.\n", argIn.Command()); + return CpptrajState::ERR; + } + DataSet_Coords* CRD = (DataSet_Coords*)State.DSL().FindSetOfGroup( setname, DataSet::COORDINATES ); + if (CRD == 0) { + mprinterr("Error: %s: No COORDS set with name %s found.\n", argIn.Command(), setname.c_str()); + return CpptrajState::ERR; + } + mprintf("\tUsing set '%s'\n", CRD->legend()); + if (CRD->Size() < 1) { + mprinterr("Error: Set '%s' has no frames.\n", CRD->legend()); + return CpptrajState::ERR; + } + out->Allocate( DataSet::SizeArray(1, CRD->Size()) ); + + mprintf("\tCalculating extended comparison similarity values.\n"); + mprintf("\tInput coords: %s\n", CRD->legend()); + mprintf("\tOutput set: %s\n", out->legend()); + if (df != 0) + mprintf("\tWriting set to file: %s\n", df->DataFilename().full()); + mprintf("\tUsing metric: %s\n", ExtendedSimilarity::metricStr(metric)); + + // Do extended similarity calculation for each frame + // TODO have ExtendedSimilarity return a DataSet_double? + ExtendedSimilarity ExtSim; + if (ExtSim.SetOpts( metric, CRD->Size(), CRD->Top().Natom()*3 )) { + mprinterr("Error: Extended similarity setup failed.\n"); + return CpptrajState::ERR; + } + ExtendedSimilarity::Darray csimvals = ExtSim.CalculateCompSim( *CRD ); + if (csimvals.empty()) { + mprinterr("Error: No comparitive similarity values calculated.\n"); + return CpptrajState::ERR; + } + const double* dptr = &csimvals[0]; + for (unsigned int idx = 0; idx != csimvals.size(); idx++, dptr++) + out->Add(idx, dptr); + + return CpptrajState::OK; +} diff --git a/src/Exec_ExtendedComparison.h b/src/Exec_ExtendedComparison.h new file mode 100644 index 0000000000..01b6ea1021 --- /dev/null +++ b/src/Exec_ExtendedComparison.h @@ -0,0 +1,12 @@ +#ifndef INC_EXEC_EXTENDEDCOMPARISON_H +#define INC_EXEC_EXTENDEDCOMPARISON_H +#include "Exec.h" +/// Calculate extended comparison similarity values for a COORDS set +class Exec_ExtendedComparison : public Exec { + public: + Exec_ExtendedComparison() : Exec(COORDS) {} + void Help() const; + DispatchObject* Alloc() const { return (DispatchObject*)new Exec_ExtendedComparison(); } + RetType Execute(CpptrajState&, ArgList&); +}; +#endif diff --git a/src/ExtendedSimilarity.cpp b/src/ExtendedSimilarity.cpp new file mode 100644 index 0000000000..b7824b7e8a --- /dev/null +++ b/src/ExtendedSimilarity.cpp @@ -0,0 +1,512 @@ +#include "ExtendedSimilarity.h" +#include "DataSet_Coords.h" +#include "CpptrajStdio.h" +#include // ceil, pow, sqrt + +/** CONSTRUCTOR */ +ExtendedSimilarity::ExtendedSimilarity() : + metric_(NO_METRIC), + cthreshType_(NO_THRESHOLD), + c_threshold_(0), + wfactorType_(FRACTION), + power_(0) +{} + +/** Set options - metric, c. threshold type, c. threshold value, weight type, weight value, + * nframes, ncoords. +  */ +int ExtendedSimilarity::SetOpts(MetricType mt, CoincidenceThresholdType ct, double cv, + WeightFactorType wt, double wv, unsigned int Nframes, + unsigned int Ncoords) +{ + metric_ = mt; + cthreshType_ = ct; + c_threshold_ = cv; + wfactorType_ = wt; + power_ = wv; + c_sum_.assign( Ncoords, 0.0 ); + return isValid(Nframes); +} + +/** Set options - metric, nframes, ncoords only; defaults for the rest. */ +int ExtendedSimilarity::SetOpts(MetricType mt, unsigned int Nframes, unsigned int Ncoords) +{ + metric_ = mt; + cthreshType_ = NO_THRESHOLD; + c_threshold_ = 0; + wfactorType_ = FRACTION; + power_ = 0; + c_sum_.assign( Ncoords, 0.0 ); + return isValid(Nframes); +} + +/** Check options and do any common setup. + * \return 0 if options are valid. + */ +int ExtendedSimilarity::isValid(unsigned int n_objects) { + if (metric_ == MSD) { + if (c_sum_.size() % 3 != 0) { + mprinterr("Error: # of coords (%zu) is not divisible by 3; required for MSD.\n", c_sum_.size()); + return 1; + } + natoms_ = c_sum_.size() / 3; + if (natoms_ < 1) { + mprinterr("Error: Similarity options are set up for MSD metric but # atoms < 1.\n"); + return 1; + } + sq_sum_.assign( c_sum_.size(), 0.0 ); + } else if (metric_ == NO_METRIC) { + mprinterr("Error: Similarity options metric not set.\n"); + return 1; + } else { + if (cthreshType_ == N_OBJECTS) { + if ((unsigned int)c_threshold_ >= n_objects) { + mprinterr("Error: c_threshold cannot be equal or greater to n_objects.\n"); + return 1; + } + } else if (cthreshType_ == FRAC_OBJECTS) { + bool in_range = (c_threshold_ > 0 && c_threshold_ < 1); + if (!in_range) { + mprinterr("Error: c_threshold fraction must be between 0 and 1.\n"); + return 1; + } + } + } + return 0; +} + +/** Strings corresponding to MetricType */ +const char* ExtendedSimilarity::MetricStr_[] = { + "Mean-squared deviation", + "Bhattacharyya's U coefficient", + "Faiman's coefficient", + "Gleason's coefficient", + "Jaccard's coefficient", + "Jaccard-Tanimoto coefficient", + "Rogers-Tanimoto coefficient", + "Russell-Rao coefficient", + "Simpson's coefficient", + "Sokal-Sneath 1 coefficient", + "Sokal-Sneath 2 coefficient", + "No metric" +}; + +/** Keywords corresponding to MetricType. */ +const char* ExtendedSimilarity::MetricKeys_[] = { + "msd", + "bub", + "fai", + "gle", + "ja", + "jt", + "rt", + "rr", + "sm", + "ss1", + "ss2", + 0 +}; + +/** \return Character string corresponding to given metric type. */ +const char* ExtendedSimilarity::metricStr(MetricType m) { + return MetricStr_[m]; +} + +/** \return MetricType corresponding to keyword. */ +ExtendedSimilarity::MetricType ExtendedSimilarity::TypeFromKeyword(std::string const& key) { + for (int i = 0; i < (int)NO_METRIC; i++) { + if (key == std::string(MetricKeys_[i])) { + return (MetricType)i; + } + } + return NO_METRIC; +} + +/** \return string containing all metric keywords. */ +std::string ExtendedSimilarity::MetricKeys() { + std::string keys; + for (int i = 0; i < (int)NO_METRIC; i++) { + if (i == 0) + keys.assign( MetricKeys_[i] ); + else + keys.append( " " + std::string(MetricKeys_[i]) ); + } + return keys; +} + +/// For debug, print double array. +static inline void printDarray(std::vector const& arr) { + int col = 0; + mprintf("["); + for (std::vector::const_iterator it = arr.begin(); it != arr.end(); ++it) { + mprintf(" %10.8g", *it); + col++; + if (col == 6) { + mprintf("\n"); + col = 0; + } + } + mprintf("]\n"); +} + +/// For debug, print Frame +static inline void printFrame(Frame const& arr) { + int col = 0; + mprintf("["); + for (int idx = 0; idx != arr.size(); idx++) { + mprintf(" %10.8g", arr[idx]); + col++; + if (col == 6) { + mprintf("\n"); + col = 0; + } + } + mprintf("]\n"); +} +/** Calculate the comparitive similarity value between two frames. */ +double ExtendedSimilarity::CalculateCompSim(Frame const& f1, Frame const& f2) +{ + // Sanity check + if (metric_ == NO_METRIC || c_sum_.empty()) { + mprinterr("Internal Error: ExtendedSimilarity::Comparison() called before SetOpts().\n"); + return 0; + } + //mprintf("["); + //printFrame(f1); + //printFrame(f2); + //mprintf("\n"); + double val = 0; + unsigned int Ncoords = c_sum_.size(); + //c_sum_.assign( Ncoords, 0.0 ); + for (unsigned int icrd = 0; icrd < Ncoords; icrd++) + c_sum_[icrd] = f1[icrd] + f2[icrd]; + if (metric_ == MSD) { + for (unsigned int icrd = 0; icrd < Ncoords; icrd++) + sq_sum_[icrd] = (f1[icrd] * f1[icrd]) + (f2[icrd] * f2[icrd]); + val = msd_condensed(c_sum_, sq_sum_, 2, natoms_); + } else { + val = Comparison(c_sum_, 2); + } + //mprintf("DEBUG val=%10.8g\n",val); + + return val; +} + +/** Calculate the comparitive similarity values for COORDS set. */ +ExtendedSimilarity::Darray ExtendedSimilarity::CalculateCompSim(DataSet_Coords& crdIn) +{ + // Sanity check + if (metric_ == NO_METRIC || c_sum_.empty()) { + mprinterr("Internal Error: ExtendedSimilarity::Comparison() called before SetOpts().\n"); + return Darray(); + } + + unsigned int Nframes = crdIn.Size(); + unsigned int Ncoords = c_sum_.size(); + Frame frmIn = crdIn.AllocateFrame(); + c_sum_.assign( Ncoords, 0.0 ); + if (metric_ == MSD) { + sq_sum_.assign( Ncoords, 0.0 ); + // Get sum and sum squares for each coordinate + for (unsigned int idx = 0; idx < Nframes; idx++) { + crdIn.GetFrame(idx, frmIn); + for (unsigned int icrd = 0; icrd < Ncoords; icrd++) { + c_sum_[icrd] += frmIn[icrd]; + sq_sum_[icrd] += frmIn[icrd] * frmIn[icrd]; + } + } + } else { + // Get sum for each coordinate + for (unsigned int idx = 0; idx < Nframes; idx++) { + crdIn.GetFrame(idx, frmIn); + for (unsigned int icrd = 0; icrd < Ncoords; icrd++) + c_sum_[icrd] += frmIn[icrd]; + } + } + + // For each frame, get the comp. similarity. + Darray comp_sims; + comp_sims.reserve( Nframes ); + Darray c_arr(c_sum_.size(), 0.0); + if (metric_ == MSD) { + Darray sq_arr(sq_sum_.size(), 0.0); + for (unsigned int idx = 0; idx < Nframes; idx++) { + crdIn.GetFrame(idx, frmIn); + for (unsigned int icrd = 0; icrd < Ncoords; icrd++) { + c_arr[icrd] = c_sum_[icrd] - frmIn[icrd]; + sq_arr[icrd] = sq_sum_[icrd] - (frmIn[icrd]*frmIn[icrd]); + } + //printDarray(sq_arr); + //printDarray(c_sum); + double val = msd_condensed(c_arr, sq_arr, Nframes-1, natoms_); + //mprintf("%8u %16.8f\n", idx, val); + comp_sims.push_back( val ); + } + } else { + for (unsigned int idx = 0; idx < Nframes; idx++) { + crdIn.GetFrame(idx, frmIn); + for (unsigned int icrd = 0; icrd < Ncoords; icrd++) + c_arr[icrd] = c_sum_[icrd] - frmIn[icrd]; + //printDarray(c_sum); + double val = Comparison(c_arr, Nframes-1); + //mprintf("%8u %16.8f\n", idx, val); + comp_sims.push_back( val ); + } + } + // Record the medioid (max dissimilarity) + max_dissim_val_ = -1; + max_dissim_idx_ = -1; + for (Darray::const_iterator it = comp_sims.begin(); it != comp_sims.end(); ++it) { + if (*it > max_dissim_val_) { + max_dissim_val_ = *it; + max_dissim_idx_ = it - comp_sims.begin(); + } + } + //mprintf("DEBUG: Max dissim. val %g at idx %li\n", max_dissim_val_, max_dissim_idx_); + + return comp_sims; +} + +/** \return Extended comparison value. Should NOT be called for MSD. + * \param c_sum Column sum of the data + * \param Nframes Number of samples (frames) + */ +double ExtendedSimilarity::Comparison(Darray const& c_sum, unsigned int Nframes) +const +{ + double val = 0; + Counters count; + count = calculate_counters(c_sum, Nframes); + //mprintf("%10.8g %10.8g %10.8g %10.8g %10.8g\n", count.w_a_, count.w_d_, count.a_, count.d_, count.total_dis_); + switch (metric_) { + case MSD : + mprinterr("Internal Error: ExtendedSimilarity::Comparison() called for MSD metric.\n"); + break; + case BUB : + val = (sqrt(count.w_a_ * count.w_d_) + count.w_a_) / + (sqrt(count.a_ * count.d_) + count.a_ + count.total_dis_); + break; + case FAI : + val = (count.w_a_ + 0.5 * count.w_d_) / (count.p_); break; + case GLE : + val = (2 * count.w_a_) / (2 * count.a_ + count.total_dis_); break; + case JA : + val = (3 * count.w_a_) / (3 * count.a_ + count.total_dis_); break; + case JT : + val = (count.w_a_) / (count.a_ + count.total_dis_); break; + case RT : + val = (count.total_w_sim_) / (count.p_ + count.total_dis_); break; + case RR : + val = (count.w_a_) / (count.p_); break; + case SM : + val = (count.total_w_sim_) / (count.p_); break; + case SS1 : + val = (count.w_a_) / (count.a_ + 2 * count.total_dis_); break; + case SS2 : + val = (2 * count.total_w_sim_) / (count.p_ + count.total_sim_); break; + default: + mprinterr("Internal Error: ExtendedSimilarity::Comparison(): Metric '%s' is unhandled.\n", + MetricStr_[metric_]); + } + // NOTE 1 - val is invalid for MSD, but we should not be here for MSD + return 1 - val; +} + +/** Mean-squared deviation + * \param c_sum (Natoms*3) Column sum of the data + * \param sq_sum (Natoms*3) Column sum of the squared data + * \param Nframes Number of samples (frames) + * \param Natoms Number of atoms in the system + */ +double ExtendedSimilarity::msd_condensed(Darray const& c_sum, Darray const& sq_sum, unsigned int Nframes, unsigned int Natoms) +const +{ + if (c_sum.size() != sq_sum.size()) { + mprinterr("Internal Error: ExtendedSimilarity::msd_condensed(): Array sizes are not equal.\n"); + return 0; + } + //msd = np.sum(2 * (N * sq_sum - c_sum ** 2)) / (N * (N - 1)) + double msd = 0; + for (unsigned int idx = 0; idx != c_sum.size(); idx++) + msd += (2 * (Nframes * sq_sum[idx] - (c_sum[idx] * c_sum[idx]))); + msd /= ((double)Nframes * (double(Nframes-1))); + //norm_msd = msd / N_atoms + msd /= Natoms; + return msd; +} + +// ------------------------------------- +// Power weight functions +ExtendedSimilarity::Darray ExtendedSimilarity::f_s_power(Darray const& in, unsigned int n_objects, double p) { + Darray out; + out.reserve(in.size()); + for (Darray::const_iterator d = in.begin(); d != in.end(); ++d) + out.push_back( pow(p, (double)(-(n_objects - *d))) ); + return out; +} + +ExtendedSimilarity::Darray ExtendedSimilarity::f_d_power(Darray const& in, unsigned int n_objects, double p) { + Darray out; + out.reserve(in.size()); + for (Darray::const_iterator d = in.begin(); d != in.end(); ++d) + out.push_back( pow(p, (double)(-(*d - (n_objects % 2)))) ); + return out; +} + +ExtendedSimilarity::Darray ExtendedSimilarity::f_s_frac(Darray const& in, unsigned int n_objects, double p) { + Darray out; + out.reserve(in.size()); + for (Darray::const_iterator d = in.begin(); d != in.end(); ++d) + out.push_back( *d / n_objects ); + return out; +} + +ExtendedSimilarity::Darray ExtendedSimilarity::f_d_frac(Darray const& in, unsigned int n_objects, double p) { + Darray out; + out.reserve(in.size()); + for (Darray::const_iterator d = in.begin(); d != in.end(); ++d) + out.push_back( (*d - (n_objects % 2)) / n_objects ); + return out; +} + +ExtendedSimilarity::Darray ExtendedSimilarity::f_one(Darray const& in, unsigned int n_objects, double p) { + Darray out(in.size(), 1.0); + return out; +} + +// ------------------------------------- +/// For debug, print boolean array +static inline void printBarray(std::vector const& arr) { + int col = 0; + mprintf("["); + for (std::vector::const_iterator it = arr.begin(); it != arr.end(); ++it) { + if (*it) + mprintf(" True"); + else + mprintf(" False"); + col++; + if (col == 12) { + mprintf("\n"); + col = 0; + } + } + mprintf("]\n"); +} + +/// \return count of true elements in boolean array +static inline unsigned int Bsum(std::vector const& arr) { + unsigned int count = 0; + for (std::vector::const_iterator it = arr.begin(); it != arr.end(); ++it) + if (*it) count++; + return count; +} + +/// \return sum over double array +static inline double Dsum(std::vector const& arr) { + double sum = 0; + for (std::vector::const_iterator it = arr.begin(); it != arr.end(); ++it) + sum += *it; + return sum; +} + +/** \return Array containing elements of d for which b is true, times 2 minus n_objects */ +ExtendedSimilarity::Darray ExtendedSimilarity::subArray(Darray const& d, Barray const& b, unsigned int n_objects) +{ + Darray out; + out.reserve(d.size()); + for (unsigned int idx = 0; idx != d.size(); idx++) + if (b[idx]) out.push_back(2 * d[idx] - n_objects); + return out; +} + +/** \return Array containing absolute value of elements of d for which b is true, times 2 minus n_objects */ +ExtendedSimilarity::Darray ExtendedSimilarity::absSubArray(Darray const& d, Barray const& b, unsigned int n_objects) +{ + Darray out; + out.reserve(d.size()); + for (unsigned int idx = 0; idx != d.size(); idx++) + if (b[idx]) out.push_back( fabs(2 * d[idx] - n_objects) ); + return out; +} + +/** Calculate 1-similarity, 0-similarity, and dissimilarity counters. + * \param c_total Column sum of the data (c_sum) + * \param n_objects Number of samples (frames) + */ +ExtendedSimilarity::Counters + ExtendedSimilarity::calculate_counters(Darray const& c_total, unsigned int n_objects) +const +{ + // Assign c_threshold + unsigned int c_threshold; + switch (cthreshType_) { + case NO_THRESHOLD : c_threshold = n_objects % 2; break; + case DISSIMILAR : c_threshold = ceil(n_objects / 2); break; + case N_OBJECTS : c_threshold = (unsigned int)c_threshold_; break; + case FRAC_OBJECTS : c_threshold = (unsigned int)(c_threshold_ * n_objects); break; + } + // Set w_factor + typedef Darray (*WgtFxnType)(Darray const&, unsigned int, double); + WgtFxnType f_s; // Similarity function + WgtFxnType f_d; // Dissimilarity function + + switch(wfactorType_) { + case POWER: + f_s = f_s_power; + f_d = f_d_power; + break; + case FRACTION: + f_s = f_s_frac; + f_d = f_d_frac; + break; + case OTHER: + f_s = f_one; + f_d = f_one; + break; + } + // A indices + Barray a_indices; + a_indices.reserve(c_total.size()); + for (Darray::const_iterator it = c_total.begin(); it != c_total.end(); ++it) + a_indices.push_back( 2 * *it - n_objects > c_threshold ); + //printBarray( a_indices ); + // D indices + Barray d_indices; + d_indices.reserve(c_total.size()); + for (Darray::const_iterator it = c_total.begin(); it != c_total.end(); ++it) + d_indices.push_back( n_objects - 2 * *it > c_threshold ); + //printBarray( d_indices ); + // Dissimilarity indices + Barray dis_indices; + dis_indices.reserve(c_total.size()); + for (Darray::const_iterator it = c_total.begin(); it != c_total.end(); ++it) + dis_indices.push_back( fabs( 2 * *it - n_objects) <= c_threshold ); + //printBarray( dis_indices ); + + Counters count; + count.a_ = Bsum(a_indices); + count.d_ = Bsum(d_indices); + count.total_dis_ = Bsum(dis_indices); + + //mprintf("%g %g %g\n", count.a_, count.d_, count.total_dis_); + + Darray a_w_array = f_s( subArray(c_total, a_indices, n_objects), n_objects, power_ ); + //printDarray( a_w_array ); + Darray d_w_array = f_s( absSubArray(c_total, d_indices, n_objects), n_objects, power_ ); + //printDarray( d_w_array ); + Darray total_w_dis_array = f_d( absSubArray(c_total, dis_indices, n_objects), n_objects, power_ ); + //printDarray( total_w_dis_array ); + + count.w_a_ = Dsum( a_w_array ); + count.w_d_ = Dsum( d_w_array ); + count.total_w_dis_ = Dsum( total_w_dis_array ); + //mprintf("%10.8f %10.8f %10.8f\n", count.w_a_, count.w_d_, count.total_w_dis_); + + count.total_sim_ = count.a_ + count.d_; + count.total_w_sim_ = count.w_a_ + count.w_d_; + count.p_ = count.total_sim_ + count.total_dis_; + count.w_p_ = count.total_w_sim_ + count.total_w_dis_; + //mprintf("%8g %10.8f %8g %10.8f\n", count.total_sim_, count.total_w_sim_, count.p_, count.w_p_); + + return count; +} diff --git a/src/ExtendedSimilarity.h b/src/ExtendedSimilarity.h new file mode 100644 index 0000000000..41f7e5c31f --- /dev/null +++ b/src/ExtendedSimilarity.h @@ -0,0 +1,116 @@ +#ifndef INC_EXTENDED_SIMILARITY_H +#define INC_EXTENDED_SIMILARITY_H +#include +#include +// Fwd declares +class DataSet_Coords; +class Frame; +/// Implements extended similarity comparisons +class ExtendedSimilarity { + public: + typedef std::vector Darray; + /// Metric types. Sync with MetricStr_ and MetricKeys_ + enum MetricType { + MSD = 0, ///< Mean-squared deviation + BUB, ///< Bhattacharyya's U coefficient + FAI, ///< Faiman's coefficient + GLE, ///< Gleason's coefficient + JA, ///< Jaccard's coefficient + JT, ///< Jaccard-Tanimoto coefficient + RT, ///< Rogers-Tanimoto coefficient + RR, ///< Russell-Rao coefficient + SM, ///< Simpson's coefficient + SS1, ///< Sokal-Sneath 1 coefficient + SS2, ///< Sokal-Sneath 2 coefficient + NO_METRIC + }; + /// Coincidence threshold types + enum CoincidenceThresholdType { + NO_THRESHOLD = 0, ///< Default, c_threshold = n_objects % 2 + DISSIMILAR, ///< Dissimilar, c_threshold = ceil(n_objects/2) + N_OBJECTS, ///< Target number of objects (< total number of objects) + FRAC_OBJECTS ///< Fraction of total number of objects + }; + /// Weight factor types + enum WeightFactorType { + FRACTION = 0, ///< similarity = d[k]/n_objects, dissimilarity = 1 - (d[k] - n_objects % 2)/n_objects + POWER, ///< similarity = n^-(n_objects - d[k]), dissimilarity = n^-(d[k] - n_objects % 2) + OTHER ///< similarity = dissimilarity = 1 + }; + + /// CONSTRUCTOR + ExtendedSimilarity(); + /// Set options - metric, c. threshold type, c. threshold value, weight type, weight value, #frames, #coords + int SetOpts(MetricType, CoincidenceThresholdType, double, WeightFactorType, double, unsigned int, unsigned int); + /// Set options - metric, #frames, #coords + int SetOpts(MetricType, unsigned int, unsigned int); + + /// \return Char string corresponding to given MetricType + static const char* metricStr(MetricType); + /// \return Type corresponding to given keyword + static MetricType TypeFromKeyword(std::string const&); + /// \return string containing all metric keywords + static std::string MetricKeys(); + + /// Calculate complimentary similarity over given COORDS set + Darray CalculateCompSim(DataSet_Coords&); + /// Calculate complimentary similarity between two frames + double CalculateCompSim(Frame const&, Frame const&); + + /// \return Medoid frame index from last call to Darray CalculateCompSim + long int MedoidIndex() const { return max_dissim_idx_; } + private: + typedef std::vector Barray; + /// Hold counters from calculate_counters + class Counters { + public: + Counters() : a_(0), w_a_(0), d_(0), w_d_(0), total_sim_(0), total_w_sim_(0), + total_dis_(0), total_w_dis_(0), p_(0), w_p_(0) {} + + double a_; + double w_a_; + double d_; + double w_d_; + double total_sim_; + double total_w_sim_; + double total_dis_; + double total_w_dis_; + double p_; + double w_p_; + }; + /// Descriptive strings corresponding to MetricType + static const char* MetricStr_[]; + /// Keywords corresponding to MetricType + static const char* MetricKeys_[]; + + /// \return 1 if set up is invalid + int isValid(unsigned int); + /// \return Extended comparison value for given arrays + double Comparison(Darray const&, unsigned int) const; + /// Calculate MSD from sum and squared sum arrays + double msd_condensed(Darray const&, Darray const&, unsigned int, unsigned int) const; + /// \return Sub-array based on values of given boolean array + static inline Darray subArray(Darray const&, Barray const&, unsigned int); + /// \return Absolute value sub-array based on values of given boolean array + static inline Darray absSubArray(Darray const&, Barray const&, unsigned int); + /// Calculate 1-similarity, 0-similarity, and dissimilarity counters from sum array + Counters calculate_counters(Darray const&, unsigned int) const; + + static Darray f_s_power(Darray const&, unsigned int, double); + static Darray f_d_power(Darray const&, unsigned int, double); + static Darray f_s_frac(Darray const&, unsigned int, double); + static Darray f_d_frac(Darray const&, unsigned int, double); + static Darray f_one(Darray const&, unsigned int, double); + + MetricType metric_; ///< Desired metric + Darray c_sum_; ///< Sum array + Darray sq_sum_; ///< Sum of squares array (MSD) + unsigned int natoms_; ///< Number of atoms (MSD) + CoincidenceThresholdType cthreshType_; ///< Coincidence threshold type + double c_threshold_; ///< Coincidence threshold value + WeightFactorType wfactorType_; ///< Weight factor type + double power_; ///< Power for weight factor power type + double max_dissim_val_; ///< Max dissimilarity value (medioid) + long int max_dissim_idx_; ///< Max dissimilarity index (medioid) +}; +#endif diff --git a/src/RPNcalc.cpp b/src/RPNcalc.cpp index 64f5cabef2..6af5bf28f7 100644 --- a/src/RPNcalc.cpp +++ b/src/RPNcalc.cpp @@ -849,6 +849,11 @@ int RPNcalc::TokenLoop(DataSetList& DSL) const { // DataSet OP Value DataSet* ds1 = Dval[1].DS(); double d2 = Dval[0].Value(); + if (ds1 == 0) { + mprinterr("Error: Expected data set '%s' value, but data set is null.\n", + T->Description()); + return 1; + } if (debug_ > 0) mprintf("DEBUG: '%s' [%s] '%f' => 'TEMP:%li'\n", ds1->legend(), T->Description(), d2, T-tokens_.begin()); diff --git a/src/Topology.cpp b/src/Topology.cpp index 39a5ccc6f4..703e866f3f 100644 --- a/src/Topology.cpp +++ b/src/Topology.cpp @@ -1175,6 +1175,16 @@ int Topology::DetermineMolecules() { return 0; } +/** Put all atoms in a single molecule. Mostly intended for cases + * where you want a pseudo-topology and do not really care about + * molecule info. + */ +int Topology::SetSingleMolecule() { + molecules_.clear(); + molecules_.push_back( Molecule(0, Natom()) ); + return 0; +} + // ----------------------------------------------------------------------------- // Topology::DetermineNumExtraPoints() void Topology::DetermineNumExtraPoints() { diff --git a/src/Topology.h b/src/Topology.h index b6946d4d69..7f0e67ca6e 100644 --- a/src/Topology.h +++ b/src/Topology.h @@ -98,6 +98,8 @@ class Topology { int NresInMol(int) const; /// Determine molecules based on bond information int DetermineMolecules(); + /// Designate all atoms as part of a single molecule. + int SetSingleMolecule(); // ----- Bond-specific routines -------------- size_t Nbonds() const { return bonds_.size()+bondsh_.size(); } BondArray const& Bonds() const { return bonds_; } diff --git a/src/Version.h b/src/Version.h index da68c718d0..4c7559c4c1 100644 --- a/src/Version.h +++ b/src/Version.h @@ -12,7 +12,7 @@ * Whenever a number that precedes is incremented, all subsequent * numbers should be reset to 0. */ -#define CPPTRAJ_INTERNAL_VERSION "V6.20.4" +#define CPPTRAJ_INTERNAL_VERSION "V6.20.5" /// PYTRAJ relies on this #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION #endif diff --git a/src/cpptrajdepend b/src/cpptrajdepend index 322b7e0562..e2c8c578a7 100644 --- a/src/cpptrajdepend +++ b/src/cpptrajdepend @@ -186,7 +186,7 @@ ClusterMap.o : ClusterMap.cpp AssociatedData.h ClusterMap.h Constants.h CpptrajF Cmd.o : Cmd.cpp Cmd.h DispatchObject.h CmdInput.o : CmdInput.cpp CmdInput.h StringRoutines.h CmdList.o : CmdList.cpp Cmd.h CmdList.h DispatchObject.h -Command.o : Command.cpp Action.h ActionFrameCounter.h ActionList.h ActionState.h ActionTopWriter.h Action_Align.h Action_Angle.h Action_AreaPerMol.h Action_AtomMap.h Action_AtomicCorr.h Action_AtomicFluct.h Action_AutoImage.h Action_Average.h Action_AvgBox.h Action_Bounds.h Action_Box.h Action_Center.h Action_Channel.h Action_CheckChirality.h Action_CheckStructure.h Action_Closest.h Action_ClusterDihedral.h Action_Contacts.h Action_CreateCrd.h Action_CreateReservoir.h Action_DNAionTracker.h Action_DSSP.h Action_Density.h Action_Diffusion.h Action_Dihedral.h Action_DihedralRMS.h Action_Dipole.h Action_DistRmsd.h Action_Distance.h Action_Energy.h Action_Esander.h Action_FilterByData.h Action_FixAtomOrder.h Action_FixImagedBonds.h Action_GIST.h Action_Grid.h Action_GridFreeEnergy.h Action_HydrogenBond.h Action_Image.h Action_InfraredSpectrum.h Action_Jcoupling.h Action_Keep.h Action_LESsplit.h Action_LIE.h Action_LipidOrder.h Action_MakeStructure.h Action_Mask.h Action_Matrix.h Action_MinImage.h Action_Molsurf.h Action_MultiDihedral.h Action_MultiPucker.h Action_MultiVector.h Action_NAstruct.h Action_NMRrst.h Action_NativeContacts.h Action_OrderParameter.h Action_Outtraj.h Action_PairDist.h Action_Pairwise.h Action_Principal.h Action_Projection.h Action_Pucker.h Action_Radgyr.h Action_Radial.h Action_RandomizeIons.h Action_Remap.h Action_ReplicateCell.h Action_Rmsd.h Action_Rotate.h Action_RunningAvg.h Action_STFC_Diffusion.h Action_Scale.h Action_SetVelocity.h Action_Spam.h Action_Strip.h Action_Surf.h Action_SymmetricRmsd.h Action_Temperature.h Action_Time.h Action_ToroidalDiffusion.h Action_Translate.h Action_Unstrip.h Action_Unwrap.h Action_Vector.h Action_VelocityAutoCorr.h Action_Volmap.h Action_Volume.h Action_Watershell.h Action_XtalSymm.h Analysis.h AnalysisList.h AnalysisState.h Analysis_AmdBias.h Analysis_AutoCorr.h Analysis_Average.h Analysis_CalcDiffusion.h Analysis_Clustering.h Analysis_ConstantPHStats.h Analysis_Corr.h Analysis_CrankShaft.h Analysis_CrdFluct.h Analysis_CrossCorr.h Analysis_CurveFit.h Analysis_Divergence.h Analysis_EvalPlateau.h Analysis_FFT.h Analysis_HausdorffDistance.h Analysis_Hist.h Analysis_IRED.h Analysis_Integrate.h Analysis_KDE.h Analysis_Lifetime.h Analysis_LowestCurve.h Analysis_Matrix.h Analysis_MeltCurve.h Analysis_Modes.h Analysis_MultiHist.h Analysis_Multicurve.h Analysis_Overlap.h Analysis_PhiPsi.h Analysis_Regression.h Analysis_RemLog.h Analysis_Rms2d.h Analysis_RmsAvgCorr.h Analysis_Rotdif.h Analysis_RunningAvg.h Analysis_Slope.h Analysis_Spline.h Analysis_State.h Analysis_Statistics.h Analysis_TI.h Analysis_TICA.h Analysis_Timecorr.h Analysis_VectorMath.h Analysis_Wavelet.h ArgList.h Array1D.h ArrayIterator.h AssociatedData.h Atom.h AtomMap.h AtomMask.h AtomType.h AxisType.h BaseIOtype.h Box.h BoxArgs.h BufferedLine.h CharMask.h Cluster/Algorithm.h Cluster/BestReps.h Cluster/CentroidArray.h Cluster/Cframes.h Cluster/Control.h Cluster/DrawGraph.h Cluster/List.h Cluster/Metric.h Cluster/MetricArray.h Cluster/Node.h Cluster/Sieve.h Cluster/Silhouette.h ClusterMap.h Cmd.h CmdInput.h CmdList.h Command.h CompactFrameArray.h ComplexArray.h Constants.h Constraints.h ControlBlock.h ControlBlock_For.h CoordinateInfo.h Corr.h Cph.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataFilter.h DataSet.h DataSetList.h DataSet_1D.h DataSet_2D.h DataSet_3D.h DataSet_Coords.h DataSet_Coords_CRD.h DataSet_Coords_REF.h DataSet_GridFlt.h DataSet_MatrixDbl.h DataSet_MatrixFlt.h DataSet_Mesh.h DataSet_Modes.h DataSet_RemLog.h DataSet_Vector.h DataSet_double.h DataSet_float.h DataSet_integer.h DataSet_integer_mem.h DataSet_pH.h DataSet_string.h Deprecated.h DiffusionResults.h DihedralSearch.h Dimension.h DispatchObject.h Energy.h Energy_Sander.h EnsembleIn.h EnsembleOutList.h Ewald.h EwaldOptions.h Ewald_ParticleMesh.h ExclusionArray.h Exec.h Exec_AddMissingRes.h Exec_Analyze.h Exec_Calc.h Exec_CatCrd.h Exec_Change.h Exec_ClusterMap.h Exec_CombineCoords.h Exec_Commands.h Exec_CompareClusters.h Exec_CompareEnergy.h Exec_CompareTop.h Exec_CrdAction.h Exec_CrdOut.h Exec_CreateSet.h Exec_DataFile.h Exec_DataFilter.h Exec_DataSetCmd.h Exec_Emin.h Exec_Flatten.h Exec_GenerateAmberRst.h Exec_Graft.h Exec_Help.h Exec_HmassRepartition.h Exec_LoadCrd.h Exec_LoadTraj.h Exec_ParallelAnalysis.h Exec_ParmBox.h Exec_ParmSolvent.h Exec_ParmStrip.h Exec_ParmWrite.h Exec_ParseTiming.h Exec_PermuteDihedrals.h Exec_Precision.h Exec_PrepareForLeap.h Exec_PrintData.h Exec_Random.h Exec_ReadData.h Exec_ReadEnsembleData.h Exec_ReadInput.h Exec_RotateDihedral.h Exec_RunAnalysis.h Exec_ScaleDihedralK.h Exec_SequenceAlign.h Exec_Set.h Exec_Show.h Exec_SortEnsembleData.h Exec_SplitCoords.h Exec_System.h Exec_Top.h Exec_Traj.h Exec_UpdateParameters.h Exec_ViewRst.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h GIST_PME.h Grid.h GridAction.h GridBin.h HistBin.h Hungarian.h ImageOption.h ImageTypes.h InputTrajCommon.h MapAtom.h MaskArray.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NC_Routines.h NameType.h NetcdfFile.h OnlineVarT.h OutputTrajCommon.h PDBfile.h PairList.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h PubFFT.h Pucker.h Pucker_PuckerMask.h Pucker_PuckerSearch.h Pucker_PuckerToken.h RPNcalc.h Random.h Range.h ReferenceAction.h ReferenceFrame.h RemdReservoirNC.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h Spline.h SplineFxnTable.h StructureCheck.h SymbolExporting.h SymmetricRmsdCalc.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajectoryFile.h Trajin.h TrajinList.h TrajoutList.h Trajout_Single.h TypeNameHolder.h Unit.h Vec3.h cuda_kernels/GistCudaSetup.cuh helpme_standalone.h molsurf.h +Command.o : Command.cpp Action.h ActionFrameCounter.h ActionList.h ActionState.h ActionTopWriter.h Action_Align.h Action_Angle.h Action_AreaPerMol.h Action_AtomMap.h Action_AtomicCorr.h Action_AtomicFluct.h Action_AutoImage.h Action_Average.h Action_AvgBox.h Action_Bounds.h Action_Box.h Action_Center.h Action_Channel.h Action_CheckChirality.h Action_CheckStructure.h Action_Closest.h Action_ClusterDihedral.h Action_Contacts.h Action_CreateCrd.h Action_CreateReservoir.h Action_DNAionTracker.h Action_DSSP.h Action_Density.h Action_Diffusion.h Action_Dihedral.h Action_DihedralRMS.h Action_Dipole.h Action_DistRmsd.h Action_Distance.h Action_Energy.h Action_Esander.h Action_FilterByData.h Action_FixAtomOrder.h Action_FixImagedBonds.h Action_GIST.h Action_Grid.h Action_GridFreeEnergy.h Action_HydrogenBond.h Action_Image.h Action_InfraredSpectrum.h Action_Jcoupling.h Action_Keep.h Action_LESsplit.h Action_LIE.h Action_LipidOrder.h Action_MakeStructure.h Action_Mask.h Action_Matrix.h Action_MinImage.h Action_Molsurf.h Action_MultiDihedral.h Action_MultiPucker.h Action_MultiVector.h Action_NAstruct.h Action_NMRrst.h Action_NativeContacts.h Action_OrderParameter.h Action_Outtraj.h Action_PairDist.h Action_Pairwise.h Action_Principal.h Action_Projection.h Action_Pucker.h Action_Radgyr.h Action_Radial.h Action_RandomizeIons.h Action_Remap.h Action_ReplicateCell.h Action_Rmsd.h Action_Rotate.h Action_RunningAvg.h Action_STFC_Diffusion.h Action_Scale.h Action_SetVelocity.h Action_Spam.h Action_Strip.h Action_Surf.h Action_SymmetricRmsd.h Action_Temperature.h Action_Time.h Action_ToroidalDiffusion.h Action_Translate.h Action_Unstrip.h Action_Unwrap.h Action_Vector.h Action_VelocityAutoCorr.h Action_Volmap.h Action_Volume.h Action_Watershell.h Action_XtalSymm.h Analysis.h AnalysisList.h AnalysisState.h Analysis_AmdBias.h Analysis_AutoCorr.h Analysis_Average.h Analysis_CalcDiffusion.h Analysis_Clustering.h Analysis_ConstantPHStats.h Analysis_Corr.h Analysis_CrankShaft.h Analysis_CrdFluct.h Analysis_CrossCorr.h Analysis_CurveFit.h Analysis_Divergence.h Analysis_EvalPlateau.h Analysis_FFT.h Analysis_HausdorffDistance.h Analysis_Hist.h Analysis_IRED.h Analysis_Integrate.h Analysis_KDE.h Analysis_Lifetime.h Analysis_LowestCurve.h Analysis_Matrix.h Analysis_MeltCurve.h Analysis_Modes.h Analysis_MultiHist.h Analysis_Multicurve.h Analysis_Overlap.h Analysis_PhiPsi.h Analysis_Regression.h Analysis_RemLog.h Analysis_Rms2d.h Analysis_RmsAvgCorr.h Analysis_Rotdif.h Analysis_RunningAvg.h Analysis_Slope.h Analysis_Spline.h Analysis_State.h Analysis_Statistics.h Analysis_TI.h Analysis_TICA.h Analysis_Timecorr.h Analysis_VectorMath.h Analysis_Wavelet.h ArgList.h Array1D.h ArrayIterator.h AssociatedData.h Atom.h AtomMap.h AtomMask.h AtomType.h AxisType.h BaseIOtype.h Box.h BoxArgs.h BufferedLine.h CharMask.h Cluster/Algorithm.h Cluster/BestReps.h Cluster/CentroidArray.h Cluster/Cframes.h Cluster/Control.h Cluster/DrawGraph.h Cluster/List.h Cluster/Metric.h Cluster/MetricArray.h Cluster/Node.h Cluster/Sieve.h Cluster/Silhouette.h ClusterMap.h Cmd.h CmdInput.h CmdList.h Command.h CompactFrameArray.h ComplexArray.h Constants.h Constraints.h ControlBlock.h ControlBlock_For.h CoordinateInfo.h Corr.h Cph.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataFilter.h DataSet.h DataSetList.h DataSet_1D.h DataSet_2D.h DataSet_3D.h DataSet_Coords.h DataSet_Coords_CRD.h DataSet_Coords_REF.h DataSet_GridFlt.h DataSet_MatrixDbl.h DataSet_MatrixFlt.h DataSet_Mesh.h DataSet_Modes.h DataSet_RemLog.h DataSet_Vector.h DataSet_double.h DataSet_float.h DataSet_integer.h DataSet_integer_mem.h DataSet_pH.h DataSet_string.h Deprecated.h DiffusionResults.h DihedralSearch.h Dimension.h DispatchObject.h Energy.h Energy_Sander.h EnsembleIn.h EnsembleOutList.h Ewald.h EwaldOptions.h Ewald_ParticleMesh.h ExclusionArray.h Exec.h Exec_AddMissingRes.h Exec_Analyze.h Exec_Calc.h Exec_CatCrd.h Exec_Change.h Exec_ClusterMap.h Exec_CombineCoords.h Exec_Commands.h Exec_CompareClusters.h Exec_CompareEnergy.h Exec_CompareTop.h Exec_CrdAction.h Exec_CrdOut.h Exec_CrdTransform.h Exec_CreateSet.h Exec_DataFile.h Exec_DataFilter.h Exec_DataSetCmd.h Exec_Emin.h Exec_ExtendedComparison.h Exec_Flatten.h Exec_GenerateAmberRst.h Exec_Graft.h Exec_Help.h Exec_HmassRepartition.h Exec_LoadCrd.h Exec_LoadTraj.h Exec_ParallelAnalysis.h Exec_ParmBox.h Exec_ParmSolvent.h Exec_ParmStrip.h Exec_ParmWrite.h Exec_ParseTiming.h Exec_PermuteDihedrals.h Exec_Precision.h Exec_PrepareForLeap.h Exec_PrintData.h Exec_Random.h Exec_ReadData.h Exec_ReadEnsembleData.h Exec_ReadInput.h Exec_RotateDihedral.h Exec_RunAnalysis.h Exec_ScaleDihedralK.h Exec_SequenceAlign.h Exec_Set.h Exec_Show.h Exec_SortEnsembleData.h Exec_SplitCoords.h Exec_System.h Exec_Top.h Exec_Traj.h Exec_UpdateParameters.h Exec_ViewRst.h ExtendedSimilarity.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h GIST_PME.h Grid.h GridAction.h GridBin.h HistBin.h Hungarian.h ImageOption.h ImageTypes.h InputTrajCommon.h MapAtom.h MaskArray.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NC_Routines.h NameType.h NetcdfFile.h OnlineVarT.h OutputTrajCommon.h PDBfile.h PairList.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h PubFFT.h Pucker.h Pucker_PuckerMask.h Pucker_PuckerSearch.h Pucker_PuckerToken.h RPNcalc.h Random.h Range.h ReferenceAction.h ReferenceFrame.h RemdReservoirNC.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h Spline.h SplineFxnTable.h StructureCheck.h SymbolExporting.h SymmetricRmsdCalc.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajectoryFile.h Trajin.h TrajinList.h TrajoutList.h Trajout_Single.h TypeNameHolder.h Unit.h Vec3.h cuda_kernels/GistCudaSetup.cuh helpme_standalone.h molsurf.h CompactFrameArray.o : CompactFrameArray.cpp Box.h CompactFrameArray.h CoordinateInfo.h CpptrajStdio.h Matrix_3x3.h Parallel.h ReplicaDimArray.h Vec3.h ComplexArray.o : ComplexArray.cpp ArrayIterator.h ComplexArray.h Constraints.o : Constraints.cpp ArgList.h Atom.h AtomMask.h AtomType.h Box.h CharMask.h Constants.h Constraints.h CoordinateInfo.h CpptrajStdio.h FileName.h Frame.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h Topology.h TypeNameHolder.h Unit.h Vec3.h @@ -199,7 +199,7 @@ CpptrajFile.o : CpptrajFile.cpp CpptrajFile.h CpptrajStdio.h FileIO.h FileIO_Bzi CpptrajState.o : CpptrajState.cpp Action.h ActionList.h ActionState.h Action_CreateCrd.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h CompactFrameArray.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_CRD.h DataSet_Coords_REF.h DataSet_Coords_TRJ.h DataSet_Topology.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleNavigator.h EnsembleOutList.h FileIO.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h ProgressBar.h Random.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajFrameIndex.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h CpptrajStdio.o : CpptrajStdio.cpp Parallel.h CurveFit.o : CurveFit.cpp CurveFit.h -DataFile.o : DataFile.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h BufferedLine.h Cluster/Cframes.h Cluster/Cmatrix_NC.h Constants.h CoordinateInfo.h Cph.h CpptrajFile.h CpptrajStdio.h DataFile.h DataIO.h DataIO_AmberEne.h DataIO_CCP4.h DataIO_CharmmFastRep.h DataIO_CharmmOutput.h DataIO_CharmmRepLog.h DataIO_CharmmRtfPrm.h DataIO_Cmatrix_Binary.h DataIO_Cmatrix_NC.h DataIO_Cpout.h DataIO_Evecs.h DataIO_Gnuplot.h DataIO_Grace.h DataIO_Mdout.h DataIO_NetCDF.h DataIO_OpenDx.h DataIO_Peaks.h DataIO_RemLog.h DataIO_Std.h DataIO_VecTraj.h DataIO_XVG.h DataIO_Xplor.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_PairwiseCache.h DataSet_PairwiseCache_MEM.h DataSet_RemLog.h Dimension.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajectoryFile.h TypeNameHolder.h Unit.h Vec3.h +DataFile.o : DataFile.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h BufferedLine.h Cluster/Cframes.h Cluster/Cmatrix_NC.h Constants.h CoordinateInfo.h Cph.h CpptrajFile.h CpptrajStdio.h DataFile.h DataIO.h DataIO_AmberEne.h DataIO_CCP4.h DataIO_CharmmFastRep.h DataIO_CharmmOutput.h DataIO_CharmmRepLog.h DataIO_CharmmRtfPrm.h DataIO_Cmatrix_Binary.h DataIO_Cmatrix_NC.h DataIO_Cpout.h DataIO_Evecs.h DataIO_Gnuplot.h DataIO_Grace.h DataIO_Mdout.h DataIO_NetCDF.h DataIO_Numpy.h DataIO_OpenDx.h DataIO_Peaks.h DataIO_RemLog.h DataIO_Std.h DataIO_VecTraj.h DataIO_XVG.h DataIO_Xplor.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_PairwiseCache.h DataSet_PairwiseCache_MEM.h DataSet_RemLog.h Dimension.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajectoryFile.h TypeNameHolder.h Unit.h Vec3.h DataFileList.o : DataFileList.cpp ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h PDBfile.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h DataFilter.o : DataFilter.cpp ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataFilter.h DataSet.h DataSetList.h DataSet_1D.h DataSet_2D.h DataSet_3D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_integer.h Dimension.h FileIO.h FileName.h FileTypes.h Frame.h GridBin.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h DataIO.o : DataIO.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataIO.h DataSet.h DataSetList.h DataSet_2D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_MatrixDbl.h Dimension.h FileIO.h FileName.h Frame.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h @@ -217,6 +217,7 @@ DataIO_Gnuplot.o : DataIO_Gnuplot.cpp ArgList.h Array1D.h AssociatedData.h Atom. DataIO_Grace.o : DataIO_Grace.cpp ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h BufferedLine.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataIO.h DataIO_Grace.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_double.h DataSet_string.h Dimension.h FileIO.h FileName.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h DataIO_Mdout.o : DataIO_Mdout.cpp ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h BufferedLine.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataIO.h DataIO_Mdout.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_double.h Dimension.h FileIO.h FileName.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h DataIO_NetCDF.o : DataIO_NetCDF.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Cluster/Cframes.h ComplexArray.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataIO.h DataIO_NetCDF.h DataSet.h DataSetList.h DataSet_1D.h DataSet_2D.h DataSet_3D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_Mat3x3.h DataSet_MatrixDbl.h DataSet_Mesh.h DataSet_Modes.h DataSet_PairwiseCache.h DataSet_PairwiseCache_MEM.h DataSet_Vector.h DataSet_Vector_Scalar.h DataSet_string.h DataSet_unsignedInt.h Dimension.h FileIO.h FileName.h Frame.h GridBin.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NC_Routines.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h Spline.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h Version.h +DataIO_Numpy.o : DataIO_Numpy.cpp ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataIO.h DataIO_Numpy.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h FileIO.h FileName.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h libnpy/npy.hpp DataIO_OpenDx.o : DataIO_OpenDx.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h BufferedLine.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataIO.h DataIO_OpenDx.h DataSet.h DataSetList.h DataSet_3D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_GridDbl.h DataSet_GridFlt.h Dimension.h FileIO.h FileName.h Frame.h Grid.h GridBin.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h ProgressBar.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h DataIO_Peaks.o : DataIO_Peaks.cpp ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataIO.h DataIO_Peaks.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_Vector_Scalar.h Dimension.h FileIO.h FileName.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h DataIO_RemLog.o : DataIO_RemLog.cpp ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h BufferedLine.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataIO.h DataIO_RemLog.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_RemLog.h Dimension.h FileIO.h FileName.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h ProgressBar.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h @@ -292,11 +293,13 @@ Exec_CompareEnergy.o : Exec_CompareEnergy.cpp Action.h ActionList.h ActionState. Exec_CompareTop.o : Exec_CompareTop.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_CompareTop.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h Exec_CrdAction.o : Exec_CrdAction.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Cmd.h CmdList.h Command.h CompactFrameArray.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_CRD.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_CrdAction.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h ProgressBar.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h Exec_CrdOut.o : Exec_CrdOut.cpp Action.h ActionFrameCounter.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_CrdOut.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h OutputTrajCommon.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h ProgressBar.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajectoryFile.h Trajin.h TrajinList.h TrajoutList.h Trajout_Single.h TypeNameHolder.h Unit.h Vec3.h +Exec_CrdTransform.o : Exec_CrdTransform.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_CrdTransform.h ExtendedSimilarity.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h Exec_CreateSet.o : Exec_CreateSet.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_CreateSet.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h RPNcalc.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h Exec_DataFile.o : Exec_DataFile.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_DataFile.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h Exec_DataFilter.o : Exec_DataFilter.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataFilter.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_DataFilter.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h ProgressBar.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h Exec_DataSetCmd.o : Exec_DataSetCmd.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h ComplexArray.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_1D.h DataSet_2D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_Mat3x3.h DataSet_MatrixDbl.h DataSet_Mesh.h DataSet_Vector.h DataSet_string.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_DataSetCmd.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h Spline.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h Exec_Emin.o : Exec_Emin.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h CharMask.h Constants.h Constraints.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnergyArray.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_Emin.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MdOpts.h MetaData.h Minimize_SteepestDescent.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h PotentialFunction.h PotentialTerm.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h +Exec_ExtendedComparison.o : Exec_ExtendedComparison.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_ExtendedComparison.h ExtendedSimilarity.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h Exec_Flatten.o : Exec_Flatten.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_1D.h DataSet_2D.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_Flatten.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h OnlineVarT.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h Exec_GenerateAmberRst.o : Exec_GenerateAmberRst.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h DistRoutines.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_GenerateAmberRst.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h ImageOption.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TorsionRoutines.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h Exec_Graft.o : Exec_Graft.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_Graft.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h @@ -331,6 +334,7 @@ Exec_Top.o : Exec_Top.cpp Action.h ActionList.h ActionState.h Analysis.h Analysi Exec_Traj.o : Exec_Traj.cpp Action.h ActionFrameCounter.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_Traj.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h Exec_UpdateParameters.o : Exec_UpdateParameters.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_Parameters.h DataSet_Topology.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_UpdateParameters.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h Exec_ViewRst.o : Exec_ViewRst.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h BufferedLine.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_ViewRst.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h ViewRst.h +ExtendedSimilarity.o : ExtendedSimilarity.cpp AssociatedData.h Atom.h AtomMask.h AtomType.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSet_Coords.h Dimension.h ExtendedSimilarity.h FileIO.h FileName.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h TextFormat.h Topology.h TypeNameHolder.h Unit.h Vec3.h ExternalFxn.o : ExternalFxn.cpp Random.h FileIO_Bzip2.o : FileIO_Bzip2.cpp CpptrajStdio.h FileIO.h FileIO_Bzip2.h FileIO_Gzip.o : FileIO_Gzip.cpp CpptrajStdio.h FileIO.h FileIO_Gzip.h diff --git a/src/cpptrajfiles b/src/cpptrajfiles index cca042ad42..224bdaae47 100644 --- a/src/cpptrajfiles +++ b/src/cpptrajfiles @@ -190,6 +190,7 @@ COMMON_SOURCES= \ DataIO_Peaks.cpp \ DataIO_Mdout.cpp \ DataIO_NetCDF.cpp \ + DataIO_Numpy.cpp \ DataIO_RemLog.cpp \ DataIO_Std.cpp \ DataIO_VecTraj.cpp \ @@ -262,11 +263,13 @@ COMMON_SOURCES= \ Exec_CompareTop.cpp \ Exec_CrdAction.cpp \ Exec_CrdOut.cpp \ + Exec_CrdTransform.cpp \ Exec_CreateSet.cpp \ Exec_DataFile.cpp \ Exec_DataFilter.cpp \ Exec_DataSetCmd.cpp \ Exec_Emin.cpp \ + Exec_ExtendedComparison.cpp \ Exec_Flatten.cpp \ Exec_GenerateAmberRst.cpp \ Exec_Graft.cpp \ @@ -301,6 +304,7 @@ COMMON_SOURCES= \ Exec_Traj.cpp \ Exec_UpdateParameters.cpp \ Exec_ViewRst.cpp \ + ExtendedSimilarity.cpp \ File_TempName.cpp \ FileIO_Bzip2.cpp \ FileIO_Gzip.cpp \ diff --git a/src/libnpy/LICENSE.sdx b/src/libnpy/LICENSE.sdx new file mode 100644 index 0000000000..c620e5f0c5 --- /dev/null +++ b/src/libnpy/LICENSE.sdx @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2021 Leon Merten Lohse + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/src/libnpy/README.md b/src/libnpy/README.md new file mode 100644 index 0000000000..ec31a167f3 --- /dev/null +++ b/src/libnpy/README.md @@ -0,0 +1,117 @@ +# libnpy + +libnpy is a simple C++ library for reading and writing of numpy's [.npy files](https://docs.scipy.org/doc/numpy/neps/npy-format.html). + +Refer to [format.py](https://github.com/numpy/numpy/blob/master/numpy/lib/format.py) for a detailed description of the .npy format. + +This libraries primary purpose is *writing* numerical data easily and efficiently into the .npy format. +It also allows *reading* .npy files, although only a very limited subset of data types are supported. + +## Features + - Writing C++ vectors (std::vector) to .npy files + - Reading (some) simple .npy files into C++ vectors + +## Supported data types and mapping + Only *scalar* *numeric* data types are supported. There is no natural way to represent more complex objects and parsing the header becomes tremendously more complex. + Supported types: + - unsigned integer + - signed integer + - floating point + - complex floating point (std::complex, ...) + +## Usage +libnpy is a header only library. You only need to download `npy.hpp` into your include path. Neither special compiler flags nor a specific build system are required. + +Optional: If you use meson, you can use the provided `meson.build` file to declare the dependency on libnpy. + +The API has changed in the last release. The old C-style API is still available, but might get removed in the future. + +### Reading data: +```c++ +#include "npy.hpp" +#include +#include + +int main() { + + const std::string path {"data.npy"}; + npy::npy_data d = npy::read_npy(path); + + std::vector data = d.data; + std::vector shape = d.shape; + bool fortran_order = d.fortran_order; +} + +``` + +### Writing data: +```c++ +#include "npy.hpp" +#include +#include + +int main() { + const std::vector data{1, 2, 3, 4, 5, 6}; + + npy::npy_data d; + d.data = data; + d.shape = {2, 3}; + d.fortran_order = false; // optional + + const std::string path{"out.npy"}; + write_npy(path, d); +} + +``` + +This will involve an additional copy of the data, which might be undesireable for larger data. The copy can be avoided by using `npy::npy_data_ptr` as follows. + +```c++ +#include "npy.hpp" +#include +#include + +int main() { + const std::vector data{1, 2, 3, 4, 5, 6}; + + npy::npy_data_ptr d; + d.data_ptr = data.data(); + d.shape = {2, 3}; + d.fortran_order = false; // optional + + const std::string path{"out.npy"}; + write_npy(path, d); +} + +``` + +See `test/` for further examples. +C++14 is required. + +## Tests +The tests can be build with `meson>=0.55` and depend on catch2. +``` +cd tests +meson setup builddir +meson test -Cbuilddir +``` + +## Known limitations +1. Only a few data types are supported. + +2. The numpy header is a literal Python dictionary and the Python syntax is very permissive. libnpy's parser was only tested with numpy's implemenation of the .npy format. + +## Contributing +Feel free to send me a pull request, open an issue, or contact me directly. + +The code is formatted with clang-format. +Please test your changes by running the tests and static analysis. +Meson automatically builds a target for clang-tidy: +``` +cd tests +meson setup builddir +ninja -C builddir clang-tidy +``` + +## License +The project is licensed under the [MIT](LICENSE) license diff --git a/src/libnpy/npy.hpp b/src/libnpy/npy.hpp new file mode 100644 index 0000000000..5f9d0fac31 --- /dev/null +++ b/src/libnpy/npy.hpp @@ -0,0 +1,599 @@ +/* + Copyright 2017-2023 Leon Merten Lohse + + Permission is hereby granted, free of charge, to any person obtaining a copy + of this software and associated documentation files (the "Software"), to deal + in the Software without restriction, including without limitation the rights + to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + copies of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be included in + all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + SOFTWARE. +*/ + +#ifndef NPY_HPP_ +#define NPY_HPP_ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace npy { + +/* Compile-time test for byte order. + If your compiler does not define these per default, you may want to define + one of these constants manually. + Defaults to little endian order. */ +#if defined(__BYTE_ORDER) && __BYTE_ORDER == __BIG_ENDIAN || defined(__BIG_ENDIAN__) || defined(__ARMEB__) || \ + defined(__THUMBEB__) || defined(__AARCH64EB__) || defined(_MIBSEB) || defined(__MIBSEB) || defined(__MIBSEB__) +const bool big_endian = true; +#else +const bool big_endian = false; +#endif + +const size_t magic_string_length = 6; +const std::array magic_string = {'\x93', 'N', 'U', 'M', 'P', 'Y'}; + +const char little_endian_char = '<'; +const char big_endian_char = '>'; +const char no_endian_char = '|'; + +constexpr std::array endian_chars = {little_endian_char, big_endian_char, no_endian_char}; +constexpr std::array numtype_chars = {'f', 'i', 'u', 'c'}; + +constexpr char host_endian_char = (big_endian ? big_endian_char : little_endian_char); + +/* npy array length */ +using ndarray_len_t = unsigned long int; +using shape_t = std::vector; + +using version_t = std::pair; + +struct dtype_t { + char byteorder; + char kind; + unsigned int itemsize; + + inline std::string str() const { + std::stringstream ss; + ss << byteorder << kind << itemsize; + return ss.str(); + } + + inline std::tuple tie() const { + return std::tie(byteorder, kind, itemsize); + } +}; + +struct header_t { + dtype_t dtype; + bool fortran_order; + shape_t shape; +}; + +inline void write_magic(std::ostream &ostream, version_t version) { + ostream.write(magic_string.data(), magic_string_length); + ostream.put(version.first); + ostream.put(version.second); +} + +inline version_t read_magic(std::istream &istream) { + std::array buf{}; + istream.read(buf.data(), sizeof(buf)); + + if (!istream) { + throw std::runtime_error("io error: failed reading file"); + } + + if (!std::equal(magic_string.begin(), magic_string.end(), buf.begin())) + throw std::runtime_error("this file does not have a valid npy format."); + + version_t version; + version.first = buf[magic_string_length]; + version.second = buf[magic_string_length + 1]; + + return version; +} + +const std::unordered_map dtype_map = { + {std::type_index(typeid(float)), {host_endian_char, 'f', sizeof(float)}}, + {std::type_index(typeid(double)), {host_endian_char, 'f', sizeof(double)}}, + {std::type_index(typeid(long double)), {host_endian_char, 'f', sizeof(long double)}}, + {std::type_index(typeid(char)), {no_endian_char, 'i', sizeof(char)}}, + {std::type_index(typeid(signed char)), {no_endian_char, 'i', sizeof(signed char)}}, + {std::type_index(typeid(short)), {host_endian_char, 'i', sizeof(short)}}, + {std::type_index(typeid(int)), {host_endian_char, 'i', sizeof(int)}}, + {std::type_index(typeid(long)), {host_endian_char, 'i', sizeof(long)}}, + {std::type_index(typeid(long long)), {host_endian_char, 'i', sizeof(long long)}}, + {std::type_index(typeid(unsigned char)), {no_endian_char, 'u', sizeof(unsigned char)}}, + {std::type_index(typeid(unsigned short)), {host_endian_char, 'u', sizeof(unsigned short)}}, + {std::type_index(typeid(unsigned int)), {host_endian_char, 'u', sizeof(unsigned int)}}, + {std::type_index(typeid(unsigned long)), {host_endian_char, 'u', sizeof(unsigned long)}}, + {std::type_index(typeid(unsigned long long)), {host_endian_char, 'u', sizeof(unsigned long long)}}, + {std::type_index(typeid(std::complex)), {host_endian_char, 'c', sizeof(std::complex)}}, + {std::type_index(typeid(std::complex)), {host_endian_char, 'c', sizeof(std::complex)}}, + {std::type_index(typeid(std::complex)), {host_endian_char, 'c', sizeof(std::complex)}}}; + +// helpers +inline bool is_digits(const std::string &str) { return std::all_of(str.begin(), str.end(), ::isdigit); } + +template +inline bool in_array(T val, const std::array &arr) { + return std::find(std::begin(arr), std::end(arr), val) != std::end(arr); +} + +inline dtype_t parse_descr(std::string typestring) { + if (typestring.length() < 3) { + throw std::runtime_error("invalid typestring (length)"); + } + + char byteorder_c = typestring.at(0); + char kind_c = typestring.at(1); + std::string itemsize_s = typestring.substr(2); + + if (!in_array(byteorder_c, endian_chars)) { + throw std::runtime_error("invalid typestring (byteorder)"); + } + + if (!in_array(kind_c, numtype_chars)) { + throw std::runtime_error("invalid typestring (kind)"); + } + + if (!is_digits(itemsize_s)) { + throw std::runtime_error("invalid typestring (itemsize)"); + } + unsigned int itemsize = std::stoul(itemsize_s); + + return {byteorder_c, kind_c, itemsize}; +} + +namespace pyparse { + +/** + Removes leading and trailing whitespaces + */ +inline std::string trim(const std::string &str) { + const std::string whitespace = " \t"; + auto begin = str.find_first_not_of(whitespace); + + if (begin == std::string::npos) return ""; + + auto end = str.find_last_not_of(whitespace); + + return str.substr(begin, end - begin + 1); +} + +inline std::string get_value_from_map(const std::string &mapstr) { + size_t sep_pos = mapstr.find_first_of(":"); + if (sep_pos == std::string::npos) return ""; + + std::string tmp = mapstr.substr(sep_pos + 1); + return trim(tmp); +} + +/** + Parses the string representation of a Python dict + + The keys need to be known and may not appear anywhere else in the data. + */ +inline std::unordered_map parse_dict(std::string in, const std::vector &keys) { + std::unordered_map map; + + if (keys.size() == 0) return map; + + in = trim(in); + + // unwrap dictionary + if ((in.front() == '{') && (in.back() == '}')) + in = in.substr(1, in.length() - 2); + else + throw std::runtime_error("Not a Python dictionary."); + + std::vector> positions; + + for (auto const &value : keys) { + size_t pos = in.find("'" + value + "'"); + + if (pos == std::string::npos) throw std::runtime_error("Missing '" + value + "' key."); + + std::pair position_pair{pos, value}; + positions.push_back(position_pair); + } + + // sort by position in dict + std::sort(positions.begin(), positions.end()); + + for (size_t i = 0; i < positions.size(); ++i) { + std::string raw_value; + size_t begin{positions[i].first}; + size_t end{std::string::npos}; + + std::string key = positions[i].second; + + if (i + 1 < positions.size()) end = positions[i + 1].first; + + raw_value = in.substr(begin, end - begin); + + raw_value = trim(raw_value); + + if (raw_value.back() == ',') raw_value.pop_back(); + + map[key] = get_value_from_map(raw_value); + } + + return map; +} + +/** + Parses the string representation of a Python boolean + */ +inline bool parse_bool(const std::string &in) { + if (in == "True") return true; + if (in == "False") return false; + + throw std::runtime_error("Invalid python boolan."); +} + +/** + Parses the string representation of a Python str + */ +inline std::string parse_str(const std::string &in) { + if ((in.front() == '\'') && (in.back() == '\'')) return in.substr(1, in.length() - 2); + + throw std::runtime_error("Invalid python string."); +} + +/** + Parses the string represenatation of a Python tuple into a vector of its items + */ +inline std::vector parse_tuple(std::string in) { + std::vector v; + const char seperator = ','; + + in = trim(in); + + if ((in.front() == '(') && (in.back() == ')')) + in = in.substr(1, in.length() - 2); + else + throw std::runtime_error("Invalid Python tuple."); + + std::istringstream iss(in); + + for (std::string token; std::getline(iss, token, seperator);) { + v.push_back(token); + } + + return v; +} + +template +inline std::string write_tuple(const std::vector &v) { + if (v.size() == 0) return "()"; + + std::ostringstream ss; + + if (v.size() == 1) { + ss << "(" << v.front() << ",)"; + } else { + const std::string delimiter = ", "; + // v.size() > 1 + ss << "("; + std::copy(v.begin(), v.end() - 1, std::ostream_iterator(ss, delimiter.c_str())); + ss << v.back(); + ss << ")"; + } + + return ss.str(); +} + +inline std::string write_boolean(bool b) { + if (b) + return "True"; + else + return "False"; +} + +} // namespace pyparse + +inline header_t parse_header(std::string header) { + /* + The first 6 bytes are a magic string: exactly "x93NUMPY". + The next 1 byte is an unsigned byte: the major version number of the file + format, e.g. x01. The next 1 byte is an unsigned byte: the minor version + number of the file format, e.g. x00. Note: the version of the file format + is not tied to the version of the numpy package. The next 2 bytes form a + little-endian unsigned short int: the length of the header data HEADER_LEN. + The next HEADER_LEN bytes form the header data describing the array's + format. It is an ASCII string which contains a Python literal expression of + a dictionary. It is terminated by a newline ('n') and padded with spaces + ('x20') to make the total length of the magic string + 4 + HEADER_LEN be + evenly divisible by 16 for alignment purposes. The dictionary contains + three keys: + + "descr" : dtype.descr + An object that can be passed as an argument to the numpy.dtype() + constructor to create the array's dtype. "fortran_order" : bool Whether the + array data is Fortran-contiguous or not. Since Fortran-contiguous arrays + are a common form of non-C-contiguity, we allow them to be written directly + to disk for efficiency. "shape" : tuple of int The shape of the array. For + repeatability and readability, this dictionary is formatted using + pprint.pformat() so the keys are in alphabetic order. + */ + + // remove trailing newline + if (header.back() != '\n') throw std::runtime_error("invalid header"); + header.pop_back(); + + // parse the dictionary + std::vector keys{"descr", "fortran_order", "shape"}; + auto dict_map = npy::pyparse::parse_dict(header, keys); + + if (dict_map.size() == 0) throw std::runtime_error("invalid dictionary in header"); + + std::string descr_s = dict_map["descr"]; + std::string fortran_s = dict_map["fortran_order"]; + std::string shape_s = dict_map["shape"]; + + std::string descr = npy::pyparse::parse_str(descr_s); + dtype_t dtype = parse_descr(descr); + + // convert literal Python bool to C++ bool + bool fortran_order = npy::pyparse::parse_bool(fortran_s); + + // parse the shape tuple + auto shape_v = npy::pyparse::parse_tuple(shape_s); + + shape_t shape; + for (auto item : shape_v) { + auto dim = static_cast(std::stoul(item)); + shape.push_back(dim); + } + + return {dtype, fortran_order, shape}; +} + +inline std::string write_header_dict(const std::string &descr, bool fortran_order, const shape_t &shape) { + std::string s_fortran_order = npy::pyparse::write_boolean(fortran_order); + std::string shape_s = npy::pyparse::write_tuple(shape); + + return "{'descr': '" + descr + "', 'fortran_order': " + s_fortran_order + ", 'shape': " + shape_s + ", }"; +} + +inline void write_header(std::ostream &out, const header_t &header) { + std::string header_dict = write_header_dict(header.dtype.str(), header.fortran_order, header.shape); + + size_t length = magic_string_length + 2 + 2 + header_dict.length() + 1; + + version_t version{1, 0}; + if (length >= 255 * 255) { + length = magic_string_length + 2 + 4 + header_dict.length() + 1; + version = {2, 0}; + } + size_t padding_len = 16 - length % 16; + std::string padding(padding_len, ' '); + + // write magic + write_magic(out, version); + + // write header length + if (version == version_t{1, 0}) { + auto header_len = static_cast(header_dict.length() + padding.length() + 1); + + std::array header_len_le16{static_cast((header_len >> 0) & 0xff), + static_cast((header_len >> 8) & 0xff)}; + out.write(reinterpret_cast(header_len_le16.data()), 2); + } else { + auto header_len = static_cast(header_dict.length() + padding.length() + 1); + + std::array header_len_le32{ + static_cast((header_len >> 0) & 0xff), static_cast((header_len >> 8) & 0xff), + static_cast((header_len >> 16) & 0xff), static_cast((header_len >> 24) & 0xff)}; + out.write(reinterpret_cast(header_len_le32.data()), 4); + } + + out << header_dict << padding << '\n'; +} + +inline std::string read_header(std::istream &istream) { + // check magic bytes an version number + version_t version = read_magic(istream); + + uint32_t header_length = 0; + if (version == version_t{1, 0}) { + std::array header_len_le16{}; + istream.read(reinterpret_cast(header_len_le16.data()), 2); + header_length = (header_len_le16[0] << 0) | (header_len_le16[1] << 8); + + if ((magic_string_length + 2 + 2 + header_length) % 16 != 0) { + // TODO(llohse): display warning + } + } else if (version == version_t{2, 0}) { + std::array header_len_le32{}; + istream.read(reinterpret_cast(header_len_le32.data()), 4); + + header_length = + (header_len_le32[0] << 0) | (header_len_le32[1] << 8) | (header_len_le32[2] << 16) | (header_len_le32[3] << 24); + + if ((magic_string_length + 2 + 4 + header_length) % 16 != 0) { + // TODO(llohse): display warning + } + } else { + throw std::runtime_error("unsupported file format version"); + } + + auto buf_v = std::vector(header_length); + istream.read(buf_v.data(), header_length); + std::string header(buf_v.data(), header_length); + + return header; +} + +inline ndarray_len_t comp_size(const shape_t &shape) { + ndarray_len_t size = 1; + for (ndarray_len_t i : shape) size *= i; + + return size; +} + +template +struct npy_data { + std::vector data = {}; + shape_t shape = {}; + bool fortran_order = false; +}; + +template +struct npy_data_ptr { + const Scalar *data_ptr = nullptr; + shape_t shape = {}; + bool fortran_order = false; +}; + +template +inline npy_data read_npy(std::istream &in) { + std::string header_s = read_header(in); + + // parse header + header_t header = parse_header(header_s); + + // check if the typestring matches the given one + const dtype_t dtype = dtype_map.at(std::type_index(typeid(Scalar))); + + if (header.dtype.tie() != dtype.tie()) { + throw std::runtime_error("formatting error: typestrings not matching"); + } + + // compute the data size based on the shape + auto size = static_cast(comp_size(header.shape)); + + npy_data data; + + data.shape = header.shape; + data.fortran_order = header.fortran_order; + + data.data.resize(size); + + // read the data + in.read(reinterpret_cast(data.data.data()), sizeof(Scalar) * size); + + return data; +} + +template +inline npy_data read_npy(const std::string &filename) { + std::ifstream stream(filename, std::ifstream::binary); + if (!stream) { + throw std::runtime_error("io error: failed to open a file."); + } + + return read_npy(stream); +} + +template +inline void write_npy(std::ostream &out, const npy_data &data) { + // static_assert(has_typestring::value, "scalar type not + // understood"); + const dtype_t dtype = dtype_map.at(std::type_index(typeid(Scalar))); + + header_t header{dtype, data.fortran_order, data.shape}; + write_header(out, header); + + auto size = static_cast(comp_size(data.shape)); + + out.write(reinterpret_cast(data.data.data()), sizeof(Scalar) * size); +} + +template +inline void write_npy(const std::string &filename, const npy_data &data) { + std::ofstream stream(filename, std::ofstream::binary); + if (!stream) { + throw std::runtime_error("io error: failed to open a file."); + } + + write_npy(stream, data); +} + +template +inline void write_npy(std::ostream &out, const npy_data_ptr &data_ptr) { + const dtype_t dtype = dtype_map.at(std::type_index(typeid(Scalar))); + + header_t header{dtype, data_ptr.fortran_order, data_ptr.shape}; + write_header(out, header); + + auto size = static_cast(comp_size(data_ptr.shape)); + + out.write(reinterpret_cast(data_ptr.data_ptr), sizeof(Scalar) * size); +} + +template +inline void write_npy(const std::string &filename, const npy_data_ptr &data_ptr) { + std::ofstream stream(filename, std::ofstream::binary); + if (!stream) { + throw std::runtime_error("io error: failed to open a file."); + } + + write_npy(stream, data_ptr); +} + +// old interface + +// NOLINTBEGIN(*-avoid-c-arrays) +template +inline void SaveArrayAsNumpy(const std::string &filename, bool fortran_order, unsigned int n_dims, + const unsigned long shape[], const Scalar *data) { + const npy_data_ptr ptr{data, {shape, shape + n_dims}, fortran_order}; + + write_npy(filename, ptr); +} + +template +inline void SaveArrayAsNumpy(const std::string &filename, bool fortran_order, unsigned int n_dims, + const unsigned long shape[], const std::vector &data) { + SaveArrayAsNumpy(filename, fortran_order, n_dims, shape, data.data()); +} + +template +inline void LoadArrayFromNumpy(const std::string &filename, std::vector &shape, bool &fortran_order, + std::vector &data) { + const npy_data n_data = read_npy(filename); + + shape = n_data.shape; + fortran_order = n_data.fortran_order; + + std::copy(n_data.data.begin(), n_data.data.end(), std::back_inserter(data)); +} + +template +inline void LoadArrayFromNumpy(const std::string &filename, std::vector &shape, + std::vector &data) { + bool fortran_order = false; + LoadArrayFromNumpy(filename, shape, fortran_order, data); +} +// NOLINTEND(*-avoid-c-arrays) + +} // namespace npy + +#endif // NPY_HPP_ diff --git a/test/Test_RefineTrajectory/RunTest.sh b/test/Test_RefineTrajectory/RunTest.sh new file mode 100755 index 0000000000..8316c7c6f9 --- /dev/null +++ b/test/Test_RefineTrajectory/RunTest.sh @@ -0,0 +1,111 @@ +#!/bin/bash + +. ../MasterTest.sh + +CleanFiles cpptraj.in rms.dat refinedcoords.crd normcoords.crd normcoords2.crd \ + trimcoords.crd trimcoords2.crd + +INPUT='-i cpptraj.in' + +TESTNAME='Trajectory refinement tests' + +Requires netcdf + +# -------------------------------------- +RefineWithLoop() { + cat > cpptraj.in < cpptraj.in < cpptraj.in < cpptraj.in < cpptraj.in < cpptraj.in <