diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx index 224e24ffec..5a0c25aae6 100644 --- a/doc/cpptraj.lyx +++ b/doc/cpptraj.lyx @@ -2911,7 +2911,7 @@ status open \begin_layout Plain Layout \align center \begin_inset Tabular - + @@ -3822,7 +3822,7 @@ Used for cluster analysis. - + \begin_inset Text \begin_layout Plain Layout @@ -3831,7 +3831,7 @@ Pairwise Cache (NetCDF) \end_inset - + \begin_inset Text \begin_layout Plain Layout @@ -3840,7 +3840,7 @@ Pairwise Cache (NetCDF) \end_inset - + \begin_inset Text \begin_layout Plain Layout @@ -3849,7 +3849,7 @@ nccmatrix \end_inset - + \begin_inset Text \begin_layout Plain Layout @@ -3858,13 +3858,60 @@ pairwise distances \end_inset - + \begin_inset Text \begin_layout Plain Layout Used for cluster analysis. \end_layout +\end_inset + + + + +\begin_inset Text + +\begin_layout Plain Layout +NetCDF Data +\end_layout + +\end_inset + + +\begin_inset Text + +\begin_layout Plain Layout +.nc +\end_layout + +\end_inset + + +\begin_inset Text + +\begin_layout Plain Layout +netcdf +\end_layout + +\end_inset + + +\begin_inset Text + +\begin_layout Plain Layout +All data +\end_layout + +\end_inset + + +\begin_inset Text + +\begin_layout Plain Layout +Only state info saved for pH data. +\end_layout + \end_inset diff --git a/src/Action_Bounds.cpp b/src/Action_Bounds.cpp index 691d556772..4d1ada65c1 100644 --- a/src/Action_Bounds.cpp +++ b/src/Action_Bounds.cpp @@ -8,6 +8,7 @@ Action_Bounds::Action_Bounds() : outfile_(0), offset_(1), + gridIsSetup_(false), grid_(0), ds_xmin_(0), ds_ymin_(0), @@ -44,6 +45,7 @@ Action::RetType Action_Bounds::Init(ArgList& actionArgs, ActionInit& init, int d if (mask_.SetMaskString( actionArgs.GetMaskNext() )) return Action::ERR; std::string dsname = actionArgs.GetStringKey("name"); offset_ = actionArgs.getKeyInt("offset", 1); + gridIsSetup_ = false; if (dxyz_[0] > -1.0) { // Set up grid if (dsname.empty()) { @@ -127,14 +129,42 @@ int Action_Bounds::SyncAction() { trajComm_.ReduceMaster( &buf, max_, 3, MPI_DOUBLE, MPI_MAX ); if (trajComm_.Master()) std::copy( buf, buf+3, max_ ); + // Set up the grid here so all processes have it. + trajComm_.MasterBcast(min_, 3, MPI_DOUBLE); + trajComm_.MasterBcast(max_, 3, MPI_DOUBLE); + calc_center_and_bins(); return 0; } #endif +/** Calculate bounds center and grid bins based on min/max/spacing. */ +int Action_Bounds::calc_center_and_bins() { + if (gridIsSetup_) return 0; + for (int i = 0; i < 3; i++) { + if (dxyz_[i] > 0.0) { + center_[i] = (max_[i] + min_[i]) / 2.0; + long int nbins = (long int)ceil( (max_[i] - min_[i]) / dxyz_[i] ) + (long int)offset_; + if (nbins < 0) { + mprinterr("Error: Overflow in bounds dim %i: max=%f min=%f spacing=%f\n", + i, max_[i], min_[i], dxyz_[i]); + return 1; + } + nxyz_[i] = (size_t)nbins; + } + } + if (grid_ != 0) { + DataSet_3D& grid3d = static_cast( *grid_ ); + if (grid3d.Allocate_N_C_D( nxyz_[0], nxyz_[1], nxyz_[2], center_, dxyz_ )) { + mprinterr("Error: Could not allocate grid %s\n", grid_->legend()); + return 1; + } + } + gridIsSetup_ = true; + return 0; +} + void Action_Bounds::Print() { static const char cXYZ[3] = {'X', 'Y', 'Z'}; - Vec3 center; - size_t nxyz[3]; mprintf(" BOUNDS: Output to %s\n", outfile_->Filename().full()); ds_xmin_->Add(0, min_ ); @@ -143,19 +173,14 @@ void Action_Bounds::Print() { ds_xmax_->Add(0, max_ ); ds_ymax_->Add(0, max_+1); ds_zmax_->Add(0, max_+2); + + calc_center_and_bins(); + for (int i = 0; i < 3; i++) { outfile_->Printf("%f < %c < %f", min_[i], cXYZ[i], max_[i]); - if (dxyz_[i] > 0.0) { - center[i] = (max_[i] + min_[i]) / 2.0; - long int nbins = (long int)ceil( (max_[i] - min_[i]) / dxyz_[i] ) + (long int)offset_; - nxyz[i] = (size_t)nbins; - outfile_->Printf("\tCenter= %f Bins=%zu", center[i], nxyz[i]); - } + if (dxyz_[i] > 0.0) + outfile_->Printf("\tCenter= %f Bins=%zu", center_[i], nxyz_[i]); outfile_->Printf("\n"); } - if (grid_ != 0) { - DataSet_3D& grid3d = static_cast( *grid_ ); - if (grid3d.Allocate_N_C_D( nxyz[0], nxyz[1], nxyz[2], center, dxyz_ )) - mprinterr("Error: Could not allocate grid %s\n", grid_->legend()); - } + } diff --git a/src/Action_Bounds.h b/src/Action_Bounds.h index daced4a07e..1e17a4500a 100644 --- a/src/Action_Bounds.h +++ b/src/Action_Bounds.h @@ -1,6 +1,7 @@ #ifndef INC_ACTION_BOUNDS_H #define INC_ACTION_BOUNDS_H #include "Action.h" +#include // size_t /// Report the min/max XYZ values for atoms in mask. class Action_Bounds : public Action { public: @@ -16,13 +17,19 @@ class Action_Bounds : public Action { Parallel::Comm trajComm_; # endif void Print(); - AtomMask mask_; - CpptrajFile* outfile_; - double max_[3]; - double min_[3]; - Vec3 dxyz_; - int offset_; - DataSet* grid_; + + int calc_center_and_bins(); + + AtomMask mask_; ///< Mask of atoms to obtain the bounds of. + CpptrajFile* outfile_; ///< File to print bounds to. + double max_[3]; ///< Minimum extent of any atom in mask. + double min_[3]; ///< Maximum extent of any atom in mask. + Vec3 dxyz_; ///< Grid spacing (if creating a grid from bounds). + size_t nxyz_[3]; ///< # grid bins in each dimension based on spacing. + Vec3 center_; ///< Grid center based on min/max. + int offset_; ///< # bins offset for grid. + bool gridIsSetup_; ///< True if grid has already been set up for bounds. + DataSet* grid_; ///< Hold grid created from bounds. DataSet* ds_xmin_; DataSet* ds_ymin_; DataSet* ds_zmin_; diff --git a/src/Cluster/Cframes.h b/src/Cluster/Cframes.h index f69ba13e1e..ff6e2e4280 100644 --- a/src/Cluster/Cframes.h +++ b/src/Cluster/Cframes.h @@ -28,6 +28,7 @@ class Cframes { void clear() { frames_.clear(); } void reserve(std::size_t n) { frames_.reserve(n); } Iarray const& Data() const { return frames_; } + int* Ptr() { return &frames_[0]; } bool empty() const { return frames_.empty(); } bool HasFrame(int) const; diff --git a/src/Cluster/Cmatrix_NC.cpp b/src/Cluster/Cmatrix_NC.cpp index 3a08929cc4..b74d1034aa 100644 --- a/src/Cluster/Cmatrix_NC.cpp +++ b/src/Cluster/Cmatrix_NC.cpp @@ -29,7 +29,7 @@ Cpptraj::Cluster::Cmatrix_NC::~Cmatrix_NC() { /** \return true if NetCDF file has CPPTRAJ_CMATRIX conventions. */ bool Cpptraj::Cluster::Cmatrix_NC::IsCpptrajCmatrix(int NCID) { - return (NC::GetAttrText(NCID, "Conventions") == "CPPTRAJ_CMATRIX"); + return (NC::GetConventions(NCID) == NC::NC_CPPTRAJCMATRIX); } #else /* BINTRAJ */ @@ -251,7 +251,7 @@ int Cpptraj::Cluster::Cmatrix_NC::CreateCmatrix(FileName const& fname, unsigned actualFrames_VID_ = -1; // Attributes - if (NC::CheckErr(nc_put_att_text(ncid_, NC_GLOBAL, "Conventions", 15, "CPPTRAJ_CMATRIX"))) + if (NC::PutConventions(ncid_, NC::NC_CPPTRAJCMATRIX)) return 1; if (NC::CheckErr(nc_put_att_text(ncid_, NC_GLOBAL, "Version", 3, "1.0"))) return 1; diff --git a/src/DataFile.cpp b/src/DataFile.cpp index 768560fbd7..18d599d975 100644 --- a/src/DataFile.cpp +++ b/src/DataFile.cpp @@ -29,6 +29,7 @@ #include "DataIO_Cmatrix_Binary.h" #include "DataIO_Cmatrix_NC.h" #include "DataIO_Peaks.h" +#include "DataIO_NetCDF.h" #include "DataIO_AmberEne.h" // CONSTRUCTOR @@ -76,10 +77,16 @@ const FileTypes::AllocToken DataFile::DF_AllocArray[] = { { "Pairwise Cache (NetCDF)", 0, 0, 0 }, # endif { "Peaks", 0, 0, DataIO_Peaks::Alloc }, +# ifdef BINTRAJ + { "NetCDF data", 0, 0, DataIO_NetCDF::Alloc }, +# else + { "NetCDF data", 0, 0, 0 }, +# endif { "Amber Energy File", 0, 0, DataIO_AmberEne::Alloc}, { "Unknown Data file", 0, 0, 0 } }; +/** Types that support reads. */ const FileTypes::KeyToken DataFile::DF_KeyArray[] = { { DATAFILE, "dat", ".dat" }, { XMGRACE, "grace", ".agr" }, @@ -96,6 +103,7 @@ const FileTypes::KeyToken DataFile::DF_KeyArray[] = { { CHARMMREPD, "charmmrepd", ".exch" }, { CHARMMOUT, "charmmout", ".charmmout"}, { CHARMMRTFPRM, "charmmrtfprm", ".rtfprm"}, + { NETCDFDATA, "netcdf", ".nc" }, { CMATRIX_BINARY,"cmatrix", ".cmatrix" }, { CMATRIX_NETCDF,"nccmatrix", ".nccmatrix" }, { PEAKS, "peaks", ".peaks" }, @@ -120,6 +128,7 @@ const FileTypes::KeyToken DataFile::DF_WriteKeyArray[] = { { CMATRIX_NETCDF,"nccmatrix", ".nccmatrix" }, { CHARMMRTFPRM, "charmmrtfprm", ".prm" }, { PEAKS, "peaks", ".peaks" }, + { NETCDFDATA, "netcdf", ".nc" }, { UNKNOWN_DATA, 0, 0 } }; @@ -333,13 +342,15 @@ int DataFile::AddDataSet(DataSet* dataIn) { dataIn->legend()); } } else { - if ((int)dataIn->Ndim() != dimension_) { - mprinterr("Error: DataSets in DataFile %s have dimension %i\n" - "Error: Attempting to add set %s of dimension %zu\n", - filename_.base(), dimension_, - dataIn->legend(), dataIn->Ndim()); - return Error("Error: Adding DataSets with different dimensions to same file" - " is currently unsupported.\n"); + if (Type() != NETCDFDATA) { + if ((int)dataIn->Ndim() != dimension_) { + mprinterr("Error: DataSets in DataFile %s have dimension %i\n" + "Error: Attempting to add set %s of dimension %zu\n", + filename_.base(), dimension_, + dataIn->legend(), dataIn->Ndim()); + return Error("Error: Adding DataSets with different dimensions to same file" + " is currently unsupported.\n"); + } } if (!dataio_->CheckValidFor(*dataIn)) { mprinterr("Error: DataSet '%s' is not valid for DataFile '%s' format.\n", @@ -501,7 +512,7 @@ int DataFile::WriteSetsToFile(FileName const& fname, DataSetList& setsToWrite) dftimer.Total()); # endif if (err > 0) - mprinterr("Error writing %iD Data to %s\n", dimension_, fname.base()); + mprinterr("Error: Could not write data sets to %s\n", fname.base()); } return err; } diff --git a/src/DataFile.h b/src/DataFile.h index 2f901aedcf..4b72424e61 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, - AMBERENE, + NETCDFDATA, AMBERENE, UNKNOWN_DATA }; DataFile(); diff --git a/src/DataIO_NetCDF.cpp b/src/DataIO_NetCDF.cpp new file mode 100644 index 0000000000..55b779fff3 --- /dev/null +++ b/src/DataIO_NetCDF.cpp @@ -0,0 +1,2037 @@ +#include "DataIO_NetCDF.h" +#include "CpptrajStdio.h" +#ifdef BINTRAJ +# include +# include +# include +# include "DataSet_1D.h" +# include "NC_Routines.h" +# include "StringRoutines.h" +# include "Version.h" +// DataSets +# include "DataSet_Mesh.h" +# include "DataSet_MatrixDbl.h" +# include "DataSet_Modes.h" +# include "DataSet_3D.h" +# include "DataSet_string.h" +# include "DataSet_Vector.h" +# include "DataSet_Vector_Scalar.h" +# include "DataSet_PairwiseCache_MEM.h" +# include "DataSet_unsignedInt.h" +#endif + +/// CONSTRUCTOR +DataIO_NetCDF::DataIO_NetCDF() : + DataIO(true, true, true), // Valid for 1D, 2D, 3D + ncid_(-1) +{ + SetValid( DataSet::MODES ); + SetValid( DataSet::VECTOR_SCALAR ); + SetValid( DataSet::PMATRIX_MEM ); +} + +// DataIO_NetCDF::ID_DataFormat() +bool DataIO_NetCDF::ID_DataFormat(CpptrajFile& infile) +{ + if (infile.OpenFile()) return false; + // Read first 3 bytes + unsigned char magic[3]; + magic[0] = 0; + magic[1] = 0; + magic[2] = 0; + infile.Read(magic, 3); + infile.CloseFile(); + if (magic[0] == 'C' && magic[1] == 'D' && magic[2] == 'F') { +# ifdef BINTRAJ + // Check that we have cpptraj conventions + return (NC::GetConventions(infile.Filename().Full()) == NC::NC_CPPTRAJDATA); +# else + mprintf("Warning: '%s' is a NetCDF file but CPPTRAJ was compiled without NetCDF support.\n", + infile.Filename().full()); +# endif + } + return false; +} + +// ----------------------------------------------------------------------------- +/** Hold info for a netcdf variable. */ +class DataIO_NetCDF::NcVar { + public: + /// CONSTRUCTOR - blank + NcVar() : vid_(-999), vtype_(0), hasBeenRead_(false) {} + /// CONSTRUCTOR + NcVar(int vidIn, nc_type vtypeIn, const char* vnameIn, int ndims, const int* dimids) : + vid_(vidIn), vtype_(vtypeIn), vname_(vnameIn), hasBeenRead_(false), dimIds_(ndims) + { + for (int i = 0; i < ndims; i++) dimIds_[i] = dimids[i]; + } + /// COPY CONSTRUCTOR + NcVar(NcVar const& rhs) : + vid_(rhs.vid_), vtype_(rhs.vtype_), vname_(rhs.vname_), hasBeenRead_(rhs.hasBeenRead_), + dimIds_(rhs.dimIds_) {} + /// ASSIGNMENT + NcVar& operator=(NcVar const& rhs) { + if (this == &rhs) return *this; + vid_ = rhs.vid_; + vtype_ = rhs.vtype_; + vname_ = rhs.vname_; + hasBeenRead_ = rhs.hasBeenRead_; + dimIds_ = rhs.dimIds_; + return *this; + } + + int VID() const { return vid_; } + nc_type Vtype() const { return vtype_; } + std::string const& Vname() const { return vname_; } + const char* vname() const { return vname_.c_str(); } + bool HasBeenRead() const { return hasBeenRead_; } + unsigned int Ndims() const { return dimIds_.size(); } + int DimId(unsigned int idx) const { return dimIds_[idx]; } + bool Empty() const { return (vid_ == -999); } + void Print() const { + mprintf("DEBUG:\tVariable %i - '%s' type %i dimIds:", vid_, vname_.c_str(), (int)vtype_); + for (std::vector::const_iterator it = dimIds_.begin(); it != dimIds_.end(); ++it) + mprintf(" %i", *it); + mprintf("\n"); + } + + void MarkRead() { hasBeenRead_ = true; } + private: + int vid_; ///< Netcdf variable id + nc_type vtype_; ///< Netcdf variable type + std::string vname_; ///< Variable name + bool hasBeenRead_; ///< True if the variable has been read in. + std::vector dimIds_; ///< Netcdf dimension ids for this var. +}; +// ----------------------------------------------------------------------------- +/** Hold info for a netcdf dimension. */ +class DataIO_NetCDF::NcDim { + public: + /// CONSTRUCTOR - blank + NcDim() : did_(-999), size_(0) {} + /// CONSTRUCTOR + NcDim(int did, std::string const& lbl, unsigned int sze) : + did_(did), label_(lbl), size_(sze) {} + /// COPY CONSTRUCTOR + NcDim(NcDim const& rhs) : did_(rhs.did_), label_(rhs.label_), size_(rhs.size_) {} + /// ASSIGNMENT + NcDim& operator=(NcDim const& rhs) { + if (this == &rhs) return *this; + did_ = rhs.did_; + label_ = rhs.label_; + size_ = rhs.size_; + return *this; + } + + int DID() const { return did_; } + std::string const& Label() const { return label_; } + unsigned int Size() const { return size_; } + bool Empty() const { return (did_ == -999); } + void Print() const { + mprintf("DEBUG:\tDimension %i - '%s' (%u)\n", did_, label_.c_str(), size_); + } + private: + int did_; ///< Netcdf dimension ID + std::string label_; ///< Dimension label + unsigned int size_; ///< Dimension size +}; + +// ----------------------------------------------------------------------------- + +// DataIO_NetCDF::ReadHelp() +void DataIO_NetCDF::ReadHelp() +{ + +} + +// DataIO_NetCDF::processReadArgs() +int DataIO_NetCDF::processReadArgs(ArgList& argIn) +{ + + return 0; +} + +/// Get integer attribute from variable +/** \return 1 if an error occurred, 0 if found, and -1 if not present. */ +static inline int GetVarIntAtt(int& ival, const char* desc, int ncid, int varid) +{ + ival = -1; + int ncerr = nc_get_att_int(ncid, varid, desc, &ival); + if (ncerr != NC_NOERR) { + if (ncerr != NC_ENOTATT) { + NC::CheckErr(ncerr); + mprinterr("Error: Could not get '%s' attribute.\n", desc); + return 1; + } + return -1; + } + return 0; +} + +/// Get double attribute from variable +static inline int GetVarDblAtt(double& dval, const char* desc, int ncid, int varid) +{ + dval = 0; + int ncerr = nc_get_att_double(ncid, varid, desc, &dval); + if (ncerr != NC_NOERR) { + if (ncerr != NC_ENOTATT) { + NC::CheckErr(ncerr); + mprinterr("Error: Could not get '%s' attribute.\n", desc); + return 1; + } + return -1; + } + return 0; +} + +/// Get double attribute array from variable +static inline int GetVarDblArrayAtt(double* darray, const char* desc, int ncid, int varid) +{ + int ncerr = nc_get_att_double(ncid, varid, desc, darray); + if (ncerr != NC_NOERR) { + if (ncerr != NC_ENOTATT) { + NC::CheckErr(ncerr); + mprinterr("Error: Could not get '%s' attribute array.\n", desc); + return 1; + } + return -1; + } + return 0; +} + +/// Get DataSet metadata from a variable +static inline MetaData GetVarMetaData(int& errStat, int ncid, int varid) +{ + MetaData meta; + errStat = 0; + // Filename + std::string att = NC::GetAttrText(ncid, varid, "filename"); + if (!att.empty()) + meta.SetFileName( att ); + // Name + att = NC::GetAttrText(ncid, varid, "name"); + if (att.empty()) { + mprinterr("Error: 'name' attribute missing for VID %i\n", varid); + errStat = 1; + return meta; + } + meta.SetName( att ); + // Aspect + att = NC::GetAttrText(ncid, varid, "aspect"); + if (!att.empty()) + meta.SetAspect( att ); + // Legend + att = NC::GetAttrText(ncid, varid, "legend"); + if (!att.empty()) + meta.SetLegend( att ); + // Index + int ival; + int ret = GetVarIntAtt(ival, "index", ncid, varid); + if (ret == 1) { + errStat = 1; + return meta; + } else if (ret == 0) { + meta.SetIdx( ival ); + } + // Ensemble number + ret = GetVarIntAtt(ival, "ensemblenum", ncid, varid); + if (ret == 1) { + errStat = 1; + return meta; + } else if (ret == 0) { + meta.SetEnsembleNum( ival ); + } + // TODO TimeSeries? + // Scalar mode + MetaData::scalarMode smode = MetaData::UNKNOWN_MODE; + att = NC::GetAttrText(ncid, varid, "scalarmode"); + if (!att.empty()) + smode = MetaData::ModeFromKeyword( att ); + meta.SetScalarMode( smode ); + // Scalar type + MetaData::scalarType stype = MetaData::UNDEFINED; + att = NC::GetAttrText(ncid, varid, "scalartype"); + if (!att.empty()) + stype = MetaData::TypeFromKeyword( att, smode ); + meta.SetScalarType( stype ); + + return meta; +} + +/// \return Variable index information from attributes +std::vector DataIO_NetCDF::getVarIndexInfo(int& errStat, int ncid, int varid) +const +{ +// idxVarIds.clear(); + errStat = 0; + std::vector Dims; + // Get # of index dimensions + int nindexvar = 0; + int ret = GetVarIntAtt(nindexvar, "nindexvar", ncid, varid); + if (ret != 0) { + mprinterr("Error: Missing 'nindexvar' attribute.\n"); + errStat = 1; + return Dims; + } + if (debug_ > 0) + mprintf("DEBUG: nindexvar= %i\n", nindexvar); + if (nindexvar < 1) return Dims; + + Dims.resize( nindexvar ); + for (int i = 0; i < nindexvar; i++) { + Dimension& dim = Dims[i]; + // Expect either label and indexid, or + // label, min, and step + std::string suffix( integerToString(i) ); + std::string label = "label" + suffix; + std::string att = NC::GetAttrText(ncid, varid, label.c_str()); + if (att.empty()) { + mprintf("Warning: Index variable %i missing 'label' attribute.\n", i); + } else { + dim.SetLabel( att ); + } + // Check for index + std::string str = "indexid" + suffix; + int indexid = -1; + int ret = GetVarIntAtt(indexid, str.c_str(), ncid, varid); + if (ret == 1) { + mprinterr("Error: Getting '%s' attribute.\n", str.c_str()); + errStat = 1; + return Dims; + } else if (ret == 0) { +// // index present +// idxVarIds[i] = indexid; + // TODO setting dim for backwards compatibility. This should + // not be necessary eventually when index vars are better handled. + dim.ChangeMin(1.0); + dim.ChangeStep(1.0); + } else if (ret == -1) { + // No index, check for min and step + std::string min = "min" + suffix; + std::string step = "step" + suffix; + double dval; + ret = GetVarDblAtt(dval, min.c_str(), ncid, varid); + if (ret != 0) { + mprinterr("Error: '%s' attribute is missing for varid %i.\n", min.c_str(), varid); + errStat = 1; + return Dims; + } + dim.ChangeMin( dval ); + ret = GetVarDblAtt(dval, step.c_str(), ncid, varid); + if (ret != 0) { + mprinterr("Error: '%s' attribute is missing for varid %i.\n", step.c_str(), varid); + errStat = 1; + return Dims; + } + dim.ChangeStep( dval ); + } + } + return Dims; +} + +/** \return the length of the dimension of a 1D variable. */ +int DataIO_NetCDF::get_1D_var_dimlen(size_t& len, int varid, const char* desc) const { + len = 0; + // Get the first dimension dim id + int lengthDim[1]; + if (NC::CheckErr(nc_inq_vardimid(ncid_, varid, lengthDim))) { + mprinterr("Error: Could not get dimension id for %s var.\n", desc); + return 1; + } + // Get lengths dim size + if (NC::CheckErr(nc_inq_dimlen(ncid_, lengthDim[0], &len))) { + mprinterr("Error: Could not get size of dimension for %s var.\n", desc); + return 1; + } + return 0; +} + +/// Shortcut for getting dimension length from Dimensions_ array +unsigned int DataIO_NetCDF::dimLen(int did) const { + return Dimensions_[did].Size(); +} + +// ----------------------------------------------------------------------------- +/** Read CPPTRAJ XY mesh set with CPPTRAJ conventions. */ +int DataIO_NetCDF::readData_1D_xy(DataSet* ds, NcVar const& yVar, VarArray& Vars) const { + // ----- XY Mesh --------------- + size_t start[1]; + size_t count[1]; + start[0] = 0; + count[0] = dimLen(yVar.DimId(0)); + // Get the Y var id + int xvarid; + int ret = GetVarIntAtt(xvarid, "indexid0", ncid_, yVar.VID()); + if (ret != 0) { + mprinterr("Error: No 'Yid' attribute for XY set '%s'.\n", yVar.vname()); + return 1; + } + DataSet_Mesh& set = static_cast( *ds ); + set.Resize(count[0]); + DataSet_Mesh::Darray& Xvals = set.SetMeshX(); + DataSet_Mesh::Darray& Yvals = set.SetMeshY(); + // Get X + if (NC::CheckErr(nc_get_vara(ncid_, xvarid, start, count, (void*)(&Xvals[0])))) { + mprinterr("Error: Could not get X values for XY set.\n"); + return 1; + } + Vars[xvarid].MarkRead(); + // Get Y + if (NC::CheckErr(nc_get_vara(ncid_, yVar.VID(), start, count, (void*)(&Yvals[0])))) { + mprinterr("Error: Could not get Y values for XY set.\n"); + return 1; + } + Vars[yVar.VID()].MarkRead(); + return 0; +} + +/** Read unsigned integer set with CPPTRAJ conventions. Separate routine from + * readData_1D needed because netcdf classic has no unsigned int type. + */ +int DataIO_NetCDF::readData_1D_unsignedInt(DataSet* ds, NcVar const& yVar, VarArray& Vars) const { + size_t start[1], count[1]; + count[0] = 1; + + unsigned int nelts = dimLen(yVar.DimId(0)); + + DataSet_unsignedInt& uintSet = static_cast( *ds ); + uintSet.Allocate( DataSet::SizeArray(1, nelts) ); + + for (unsigned int idx = 0; idx < nelts; idx++) { + start[0] = idx; + int ival; + if (NC::CheckErr(nc_get_vara(ncid_, yVar.VID(), start, count, &ival))) { + mprinterr("Error: Coult not get element %u for unsigned int set '%s'\n", idx, uintSet.legend()); + return 1; + } + uintSet.AddElement( ival ); + } + Vars[yVar.VID()].MarkRead(); + return 0; +} + +/** Read vector set with CPPTRAJ conventions. */ +int DataIO_NetCDF::readData_1D_vector(DataSet* ds, NcVar const& yVar, VarArray& Vars) const { + size_t start[2]; + start[1] = 0; + size_t count[2]; + count[0] = 1; + count[1] = 3; + + unsigned int nvecs = dimLen(yVar.DimId(0)); + + DataSet_Vector& vecSet = static_cast( *ds ); + //vecSet.Allocate(DataSet::SizeArray(1, nvecs)); + + // Get the xyz values + vecSet.Resize( nvecs ); + for (unsigned int idx = 0; idx < nvecs; idx++) { + start[0] = idx; + if (NC::CheckErr(nc_get_vara(ncid_, yVar.VID(), start, count, vecSet[idx].Dptr()))) { + mprinterr("Error: Could not get vector xyz values for %u '%s'\n", idx, vecSet.legend()); + return 1; + } + } + Vars[yVar.VID()].MarkRead(); + // Get the origin values + int originVarid; + int ret = GetVarIntAtt(originVarid, "originsid", ncid_, yVar.VID()); + if (ret == 1) { + mprinterr("Error: Problem getting 'originsid' attribute for '%s'\n", vecSet.legend()); + return 1; + } else if (ret == 0) { + vecSet.ResizeOrigins( nvecs ); + for (unsigned int idx = 0; idx < nvecs; idx++) { + start[0] = idx; + if (NC::CheckErr(nc_get_vara(ncid_, originVarid, start, count, vecSet.ModifyOxyz(idx).Dptr()))) { + mprinterr("Error: Could not get vector origin values for %u '%s'\n", idx, vecSet.legend()); + return 1; + } + } + Vars[originVarid].MarkRead(); + } + return 0; +} + +/** Read vector/scalar set with CPPTRAJ conventions. */ +int DataIO_NetCDF::readData_1D_vector_scalar(DataSet* ds, NcVar const& yVar, VarArray& Vars) const { + size_t start[2]; + start[1] = 0; + size_t count[2]; + count[0] = 1; + count[1] = 3; + + unsigned int nvecs = dimLen(yVar.DimId(0)); + + DataSet_Vector_Scalar& vecSet = static_cast( *ds ); + //vecSet.Allocate(DataSet::SizeArray(1, nvecs)); + + // Get the xyz values + vecSet.Resize( nvecs ); + for (unsigned int idx = 0; idx < nvecs; idx++) { + start[0] = idx; + if (NC::CheckErr(nc_get_vara(ncid_, yVar.VID(), start, count, vecSet.ModifyVec(idx).Dptr()))) { + mprinterr("Error: Could not get vector xyz values for %u '%s'\n", idx, vecSet.legend()); + return 1; + } + } + Vars[yVar.VID()].MarkRead(); + + // Get the scalar values + int scalarVid; + int ret = GetVarIntAtt(scalarVid, "scalarid", ncid_, yVar.VID()); + if (ret != 0) { + mprinterr("Error: Could not get 'scalarid' attribute for '%s'\n", vecSet.legend()); + return 1; + } + + start[0] = 0; + count[0] = nvecs; + if (NC::CheckErr(nc_get_vara(ncid_, scalarVid, start, count, vecSet.ValPtr()))) { + mprinterr("Error: Could not get scalar values for '%s'\n", vecSet.legend()); + return 1; + } + Vars[scalarVid].MarkRead(); + + return 0; +} + +/** Read string set with CPPTRAJ conventions. */ +int DataIO_NetCDF::readData_1D_string(DataSet* ds, NcVar const& yVar, VarArray& Vars) const { + // Get the lengths var id + int lengthsVarid; + int ret = GetVarIntAtt(lengthsVarid, "strlengthsid", ncid_, yVar.VID()); + if (ret != 0) { + mprinterr("Error: No 'strlengthsid' attribute for string set '%s'.\n", yVar.vname()); + return 1; + } + // Get the string lengths array dimension size + size_t nstrings; + if (get_1D_var_dimlen(nstrings, lengthsVarid, "string lengths")) return 1; + + DataSet_string& strSet = static_cast( *ds ); + strSet.Resize( nstrings ); + + size_t cstart[1], nstart[1]; + size_t ccount[1], ncount[1]; + cstart[0] = 0; + nstart[0] = 0; + ncount[0] = 1; + std::vector buffer; + for (size_t idx = 0; idx < nstrings; idx++) { + //if (debug_ > 1) mprintf("DEBUG: read string %zu", idx); + // Get the string length + int len; + if (NC::CheckErr(nc_get_vara_int(ncid_, lengthsVarid, nstart, ncount, &len))) { + mprinterr("Error: Could not get length of string %zu\n", idx); + return 1; + } + //if (debug_ > 1) mprintf(" len=%i", len); + if ((size_t)len+1 > buffer.size()) + buffer.resize( len+1 ); + nstart[0]++; + // Get the string + ccount[0] = len; + if (NC::CheckErr(nc_get_vara_text(ncid_, yVar.VID(), cstart, ccount, &buffer[0]))) { + mprinterr("Error: Could not get string %zu\n", idx); + return 1; + } + buffer[len] = '\0'; + //if (debug_ > 1) mprintf(" str='%s'\n", &buffer[0]); + cstart[0] += (size_t)len; + // Assign to DataSet + strSet.SetElement( idx, &buffer[0] ); + } + // Mark vars as read + Vars[yVar.VID()].MarkRead(); + Vars[lengthsVarid].MarkRead(); + + return 0; +} + +/** Read 1D array with CPPTRAJ conventions. */ +int DataIO_NetCDF::readData_1D(DataSet* ds, NcVar const& yVar, VarArray& Vars) const { + // ----- 1D Scalar ------------- + size_t start[1]; + size_t count[1]; + start[0] = 0; + count[0] = dimLen(yVar.DimId(0)); + + DataSet_1D& set = static_cast( *ds ); + set.Resize( count[0] ); + if (NC::CheckErr(nc_get_vara(ncid_, yVar.VID(), start, count, (void*)(set.Yptr())))) { + mprinterr("Error: Could not get values for set.\n"); + return 1; + } + Vars[yVar.VID()].MarkRead(); + return 0; +} + +/** Read cluster pairwise matrix with CPPTRAJ conventions. */ +int DataIO_NetCDF::readData_cluster_pwmatrix(DataSet* ds, NcVar const& matrixVar, VarArray& Vars) const { + // Get n_original_frames attribute + int n_original_frames; + int ret = GetVarIntAtt(n_original_frames, "n_original_frames", ncid_, matrixVar.VID()); + if (ret != 0) { + mprinterr("Error: Could not get 'n_original_frames' attribute.\n"); + return 1; + } + // Get metric description + std::string metric_description = NC::GetAttrText(ncid_, matrixVar.VID(), "MetricDescription"); + // Get the sieve attribute + int sieve; + ret = GetVarIntAtt(sieve, "sieve", ncid_, matrixVar.VID()); + if (ret != 0) { + mprinterr("Error: Could not get 'sieve' attribute.\n"); + return 1; + } + // Get the frames being clustered + Cpptraj::Cluster::Cframes frames_to_cluster; + int actual_framesid = -1; + if (sieve != 1) { + ret = GetVarIntAtt(actual_framesid, "indexid0", ncid_, matrixVar.VID()); + if (ret != 0) { + mprinterr("Error: Could not get 'indexid0' attribute (for clustered frames variable ID) but sieve != 1.\n"); + return 1; + } + // TODO check that indexid0 == indexid1? + // Get the actual_frames dimension length (i.e. # rows) TODO check against nrows attribute? + size_t n_actual_frames; + if (get_1D_var_dimlen(n_actual_frames, actual_framesid, "actual # frames")) return 1; + // Read actual frames + frames_to_cluster.assign(n_actual_frames, -1); + size_t start[1], count[1]; + start[0] = 0; + count[0] = n_actual_frames; + if (NC::CheckErr(nc_get_vara(ncid_, actual_framesid, start, count, frames_to_cluster.Ptr()))) { + mprinterr("Error: Could not read actual frames array.\n"); + return 1; + } + Vars[actual_framesid].MarkRead(); + } else { + // Sieve 1 (i.e. no sieveing) + frames_to_cluster.reserve(n_original_frames); + for (int frm = 0; frm < n_original_frames; frm++) + frames_to_cluster.push_back( frm ); + } + if (debug_ > 1) { + mprintf("DEBUG: Frames in pairwise matrix:"); + for (Cpptraj::Cluster::Cframes::const_iterator it = frames_to_cluster.begin(); + it != frames_to_cluster.end(); ++it) + mprintf(" %i", *it); + mprintf("\n"); + } + // Allocate matrix + DataSet_PairwiseCache_MEM& pwmatrix = static_cast( *ds ); + if (pwmatrix.SetupCache( n_original_frames, frames_to_cluster, sieve, metric_description )) + { + mprinterr("Error: Failed to set up cluster pairwise matrix.\n"); + return 1; + } + // Read matrix + size_t start[1], count[1]; + start[0] = 0; + count[0] = dimLen(matrixVar.DimId(0)); + if (NC::CheckErr(nc_get_vara(ncid_, matrixVar.VID(), start, count, pwmatrix.Ptr()))) { + mprinterr("Error: Could not read pairwise matrix variable.\n"); + return 1; + } + + return 0; +} + +/** Read 2D matrix with CPPTRAJ conventions. */ +int DataIO_NetCDF::readData_2D(DataSet* ds, NcVar const& matVar, VarArray& Vars) const { + // ----- 2D Matrix ------------- + // Get nrows/ncols + int ncols, nrows; + int ret = GetVarIntAtt(ncols, "ncols", ncid_, matVar.VID()); + if (ret != 0) { + mprinterr("Error: Could not get 'ncols'.\n"); + return 1; + } + ret = GetVarIntAtt(nrows, "nrows", ncid_, matVar.VID()); + if (ret != 0) { + mprinterr("Error: Could not get 'nrows'.\n"); + return 1; + } + // Get matrix kind + std::string mkind = NC::GetAttrText(ncid_, matVar.VID(), "matrixkind"); + // Allocate + DataSet_2D& mat = static_cast( *ds ); + int allocErr = 0; + if (mkind == "full") + allocErr = mat.Allocate2D(ncols, nrows); + else if (mkind == "half") + allocErr = mat.AllocateHalf(ncols); + else if (mkind == "tri") + allocErr = mat.AllocateTriangle(ncols); + else { + mprinterr("Error: Urecognized matrix kind: %s\n", mkind.c_str()); + return 1; + } + if (allocErr != 0) { + mprinterr("Error: Could not allocate matrix.\n"); + return 1; + } + // Read values + size_t start[1]; + size_t count[1]; + start[0] = 0; + count[0] = dimLen(matVar.DimId(0)); + if (NC::CheckErr(nc_get_vara(ncid_, matVar.VID(), start, count, (void*)(mat.MatrixPtr())))) { + mprinterr("Error: Could not get values for matrix.\n"); + return 1; + } + Vars[matVar.VID()].MarkRead(); + // Check for nsnapshots + int nsnapshots = 0; + ret = GetVarIntAtt(nsnapshots, "nsnapshots", ncid_, matVar.VID()); + if (ret == 1) + return 1; + else if (ret == 0) { + DataSet_MatrixDbl& dmat = static_cast( mat ); + dmat.SetNsnapshots( nsnapshots ); + } + // Check for vectid + int vectVarId = -1; + ret = GetVarIntAtt(vectVarId, "vectid", ncid_, matVar.VID()); + if (ret == 1) return 1; + if (ret == 0) { + mprintf("\tMatrix has diagonal vector data.\n"); + if (mat.Type() != DataSet::MATRIX_DBL) { + mprinterr("Error: Variable has vect data but set is not double matrix.\n"); + return 1; + } + unsigned int vectLength = dimLen(Vars[vectVarId].DimId(0)); + DataSet_MatrixDbl& dmat = static_cast( mat ); + dmat.AllocateVector( vectLength ); + count[0] = vectLength; + if (NC::CheckErr(nc_get_vara(ncid_, vectVarId, start, count, (void*)(&dmat.V1()[0])))) { + mprinterr("Error: Could not get vect for matrix.\n"); + return 1; + } + Vars[vectVarId].MarkRead(); + } + // Check for massid + int massVarId = -1; + ret = GetVarIntAtt(massVarId, "massid", ncid_, matVar.VID()); + if (ret == 1) return 1; + if (ret == 0) { + mprintf("\tMatrix has mass data.\n"); + if (mat.Type() != DataSet::MATRIX_DBL) { + mprinterr("Error: Variable has mass data but set is not double matrix.\n"); + return 1; + } + unsigned int massLength = dimLen(Vars[massVarId].DimId(0)); + DataSet_MatrixDbl& dmat = static_cast( mat ); + dmat.AllocateMass( massLength ); + count[0] = massLength; + if (NC::CheckErr(nc_get_vara(ncid_, massVarId, start, count, (void*)(&dmat.M1()[0])))) { + mprinterr("Error: Could not get mass for matrix.\n"); + return 1; + } + Vars[massVarId].MarkRead(); + } + return 0; +} + +/** Read 3D grid with CPPTRAJ conventions. */ +int DataIO_NetCDF::readData_3D(DataSet* ds, NcVar const& gridVar, VarArray& Vars) const { + // ----- 3D Grid --------------------- + // Get nx/ny/nz + int nx, ny, nz; + int ret = GetVarIntAtt(nx, "nx", ncid_, gridVar.VID()); + if (ret != 0) { + mprinterr("Error: Could not get 'nx'.\n"); + return 1; + } + ret = GetVarIntAtt(ny, "ny", ncid_, gridVar.VID()); + if (ret != 0) { + mprinterr("Error: Could not get 'ny'.\n"); + return 1; + } + ret = GetVarIntAtt(nz, "nz", ncid_, gridVar.VID()); + if (ret != 0) { + mprinterr("Error: Could not get 'nz'.\n"); + return 1; + } + // Get origin + double oxyz[3]; + ret = GetVarDblArrayAtt(oxyz, "origin", ncid_, gridVar.VID()); + if (ret == 1) { + mprinterr("Error: Could not get attribute 'origin'.\n"); + return 1; + } else if (ret == -1) { + mprinterr("Error: Missing attribute 'origin'.\n"); + return 1; + } + // Get ucell + double ucell[9]; + ret = GetVarDblArrayAtt(ucell, "ucell", ncid_, gridVar.VID()); + if (ret == 1) { + mprinterr("Error: Could not get attribute 'ucell'.\n"); + return 1; + } else if (ret == -1) { + mprinterr("Error: Missing attribute 'ucell'.\n"); + return 1; + } + Box gridBox; + if (gridBox.SetupFromUcell( ucell )) { + mprinterr("Error: Could not set up grid unit cell.\n"); + return 1; + } + + // Allocate + DataSet_3D& grid = static_cast( *ds ); + int allocErr = grid.Allocate_N_O_Box( nx, ny, nz, Vec3(oxyz), gridBox ); + if (allocErr != 0) { + mprinterr("Error: Could not allocate grid.\n"); + return 1; + } + // Read values + size_t start[1]; + size_t count[1]; + start[0] = 0; + count[0] = dimLen(gridVar.DimId(0)); + if (NC::CheckErr(nc_get_vara(ncid_, gridVar.VID(), start, count, (void*)(grid.GridPtr())))) { + mprinterr("Error: Could not get values for grid.\n"); + return 1; + } + Vars[gridVar.VID()].MarkRead(); + // DEBUG + if (debug_ > 0) + grid.GridInfo(); + + return 0; +} + +/** Read modes data with CPPTRAJ conventions. */ +int DataIO_NetCDF::readData_modes(DataSet* ds, NcVar const& modesVar, VarArray& Vars) const { + // ----- Modes ----------------- + unsigned int n_eigenvalues = dimLen(modesVar.DimId(0)); + // Get the eigenvectors variable ID + int vectorsVarId = -1; + int ret = GetVarIntAtt(vectorsVarId, "eigenvectorsid", ncid_, modesVar.VID()); + if (ret != 0) { + mprinterr("Error: Modes data missing 'eigenvectorsid'.\n"); + return 1; + } + unsigned int evectorLength = dimLen(Vars[vectorsVarId].DimId(0)); + // Get the avg. coords variable ID + int coordsVarId = -1; + ret = GetVarIntAtt(coordsVarId, "avgcoordsid", ncid_, modesVar.VID()); + if (ret != 0) { + mprinterr("Error: Modes data missing 'avgcoordsid'.\n"); + return 1; + } + unsigned int avgCoordsLength = dimLen(Vars[coordsVarId].DimId(0)); + // Get the mass variable ID if present + int massVarId = -1; + unsigned int massLength = 0; + ret = GetVarIntAtt(massVarId, "massid", ncid_, modesVar.VID()); + if (ret == 1) { + mprinterr("Error: Could not get 'massid'.\n"); + return 1; + } else if (ret == 0) { + massLength = dimLen(Vars[massVarId].DimId(0)); + } + mprintf("\tModes: # values= %u, evector size= %u, avg coords size= %u, mass size = %u\n", + n_eigenvalues, evectorLength, avgCoordsLength, massLength); + // Allocate the modes set + DataSet_Modes& modes = static_cast( *ds ); + if (modes.AllocateModes(n_eigenvalues, evectorLength, avgCoordsLength, massLength)) { + mprinterr("Error: Could not allocate memory for modes set.\n"); + return 1; + } + // Read modes data + size_t start[1]; + size_t count[1]; + start[0] = 0; + count[0] = dimLen(modesVar.DimId(0)); + if (NC::CheckErr(nc_get_vara(ncid_, modesVar.VID(), start, count, modes.EvalPtr()))) { + mprinterr("Error: Could not read eigenvalues.\n"); + return 1; + } + count[0] = evectorLength; + if (NC::CheckErr(nc_get_vara(ncid_, vectorsVarId, start, count, modes.EvectPtr()))) { + mprinterr("Error: Could not read eigenvectors.\n"); + return 1; + } + count[0] = avgCoordsLength; + if (NC::CheckErr(nc_get_vara(ncid_, coordsVarId, start, count, modes.AvgFramePtr()))) { + mprinterr("Error: Could not read avg. coords.\n"); + return 1; + } + if (massVarId != -1 && massLength > 0) { + count[0] = massLength; + if (NC::CheckErr(nc_get_vara(ncid_, massVarId, start, count, modes.MassPtr()))) { + mprinterr("Error: Could not read masses.\n"); + return 1; + } + } + // Mark variables as read + Vars[modesVar.VID()].MarkRead(); + Vars[vectorsVarId].MarkRead(); + Vars[coordsVarId].MarkRead(); + if (massVarId != -1) + Vars[massVarId].MarkRead(); + return 0; +} + +/** Read variable with CPPTRAJ conventions. */ +int DataIO_NetCDF::read_cpptraj_vars(DataSetList& dsl, std::string const& dsname, VarArray& Vars) +const +{ + for (VarArray::const_iterator var = Vars.begin(); var != Vars.end(); ++var) + { + if (var->HasBeenRead()) continue; + + // Get the description + std::string desc = NC::GetAttrText(ncid_, var->VID(), "description"); + // Get the type from the description + DataSet::DataType dtype = DataSet::TypeFromDescription( desc ); + if (debug_ > 0) + mprintf("\t%s Description: %s (%i)\n", var->vname(), desc.c_str(), (int)dtype); + // Get metaData + int errStat = 0; + MetaData meta = GetVarMetaData( errStat, ncid_, var->VID() ); + if (errStat != 0) { + mprinterr("Error: Could not set up meta data for variable '%s'\n", var->vname()); + return 1; + } + // Replace name if user has specified a name + if (!user_specified_name_.empty()) { + // First, check if overwriting the name will clash with an existing set. + MetaData checkMeta = meta; + checkMeta.SetName( user_specified_name_ ); + DataSet* checkSet = dsl.CheckForSet( checkMeta ); + if (checkSet != 0) { + mprinterr("Error: Using specified name '%s' clashes with existing set '%s'\n", + checkMeta.PrintName().c_str(), checkSet->Meta().PrintName().c_str()); + return 1; + } + meta.SetName( user_specified_name_ ); + } + // Get DataSet dimensions + errStat = 0; + std::vector Dims = getVarIndexInfo(errStat, ncid_, var->VID()); + if (errStat != 0) return 1; + if (debug_ > 0) + for (std::vector::const_iterator dim = Dims.begin(); dim != Dims.end(); ++dim) + mprintf("DEBUG:\t Var %s dim %s min %f step %f\n", var->vname(), dim->label(), dim->Min(), dim->Step()); + // Add the set + DataSet* ds = dsl.AddSet( dtype, meta ); + if (ds == 0) { + mprinterr("Error: Could not allocate set '%s'\n", meta.PrintName().c_str()); + return 1; + } + // Set DataSet dimensions + unsigned int idx = 0; + for (std::vector::const_iterator dim = Dims.begin(); dim != Dims.end(); ++dim) + ds->SetDim(idx++, *dim); + // Check netcdf variable dimensions + if (debug_ > 0) + mprintf("DEBUG: %s dim length %u\n", var->vname(), dimLen(var->DimId(0)) ); + if (dtype == DataSet::XYMESH) { + if (readData_1D_xy(ds, *var, Vars)) + return 1; + } else if (dtype == DataSet::VECTOR) { + if (readData_1D_vector(ds, *var, Vars)) + return 1; + } else if (dtype == DataSet::VECTOR_SCALAR) { + if (readData_1D_vector_scalar(ds, *var, Vars)) + return 1; + } else if (dtype == DataSet::STRING) { + if (readData_1D_string(ds, *var, Vars)) + return 1; + } else if (ds->Group() == DataSet::SCALAR_1D) { + if (readData_1D(ds, *var, Vars)) + return 1; + } else if (dtype == DataSet::PMATRIX_MEM) { + if (readData_cluster_pwmatrix(ds, *var, Vars)) + return 1; + } else if (ds->Group() == DataSet::MATRIX_2D) { + if (readData_2D(ds, *var, Vars)) + return 1; + } else if (ds->Group() == DataSet::GRID_3D) { + if (readData_3D(ds, *var, Vars)) + return 1; + } else if ( dtype == DataSet::MODES ) { + if (readData_modes(ds, *var, Vars)) + return 1; + } else { + mprinterr("Error: Cannot read type '%s' yet.\n", desc.c_str()); + return 1; + } + } // END loop over Vars + return 0; +} + +// DataIO_NetCDF::ReadData() +int DataIO_NetCDF::ReadData(FileName const& fname, DataSetList& dsl, std::string const& dsname) +{ +# ifdef BINTRAJ + // Check if the user specified a data set name. The default is to use the + // file base name, optionally plus '_' + if ( dsname.find(fname.Base()) == std::string::npos ) { + mprintf("\tUser has specified a data set name.\n" + "\tThis will replace any existing data set names in NetCDF file.\n"); + user_specified_name_ = dsname; + } else { + mprintf("\tUser has not specified a data set name. Using data set names in NetCDF file.\n"); + user_specified_name_.clear(); + } + + if (NC::CheckErr( nc_open( fname.full(), NC_NOWRITE, &ncid_ ) != NC_NOERR )) { + mprinterr("Error: Could not open NetCDF data file '%s'\n", fname.full()); + return 1; + } + + // Check if we have CPPTRAJ conventions + bool hasCpptrajConventions = (NC::GetConventions(ncid_) == NC::NC_CPPTRAJDATA); + if (hasCpptrajConventions) + mprintf("\tNetCDF data file has CPPTRAJ conventions.\n"); + else { + mprinterr("Error: NetCDF data file has unknown conventions.\n"); + return 1; + } + + // Get the number of dimensions, variables, attributes, and ID of the + // unlimited dimension (if any). + int ndimsp, nvarsp, ngattsp,unlimdimidp; + if (NC::CheckErr( nc_inq(ncid_, &ndimsp, &nvarsp, &ngattsp, &unlimdimidp) )) { + mprinterr("Error: Could not get NetCDF data file information.\n"); + return 1; + } + if (debug_ > 0) + mprintf("DEBUG: '%s' : ndimsp=%i nvarsp=%i ngattsp=%i unlimdimidp=%i\n", + fname.full(), ndimsp, nvarsp, ngattsp, unlimdimidp); + char varName[NC_MAX_NAME+1]; + // Get the length of all dimensions + Dimensions_.clear(); + Dimensions_.reserve(ndimsp); + for (int idim = 0; idim < ndimsp; idim++) { + size_t diml; + if (NC::CheckErr(nc_inq_dim(ncid_, idim, varName, &diml))) { + mprinterr("Error: Could not get length of NetCDF data dimension %i\n", idim); + return 1; + } + Dimensions_.push_back( NcDim( idim, varName, diml ) ); + if (debug_ > 0) + Dimensions_.back().Print(); + } + + VarArray AllVars; + AllVars.reserve( nvarsp ); + + // Loop over all variables in the NetCDF file. + nc_type varType = 0; // Variable type + int nVarDims = -1; // # variable dimensions + int varDimIds[NC_MAX_VAR_DIMS]; // variable dim ids + int nVarAttributes = -1; // number of variable attributes + for (int ivar = 0; ivar < nvarsp; ivar++) { + if (NC::CheckErr(nc_inq_var(ncid_, ivar, varName, &varType, &nVarDims, varDimIds, &nVarAttributes))) { + mprinterr("Error: Could not get NetCDF data variable name %i\n", ivar); + return 1; + } + if (debug_ > 0) + mprintf("DEBUG:\tVariable %i - '%s', %i dims, %i attributes\n", ivar, varName, nVarDims, nVarAttributes); + AllVars.push_back( NcVar(ivar, varType, varName, nVarDims, varDimIds) ); + } + + for (VarArray::const_iterator it = AllVars.begin(); it != AllVars.end(); ++it) + mprintf("\t %i (%s)\n", it->VID(), it->vname()); + + if (read_cpptraj_vars(dsl, dsname, AllVars)) return 1; + + nc_close( ncid_ ); + return 0; +# else + return 1; +# endif +} + +// ============================================================================= +// DataIO_NetCDF::WriteHelp() +void DataIO_NetCDF::WriteHelp() +{ + +} + +// DataIO_NetCDF::processWriteArgs() +int DataIO_NetCDF::processWriteArgs(ArgList& argIn) +{ + + return 0; +} + +// ----------------------------------------------------------------------------- +/** Hold a pool of pointers to DataSets in the list. They will be marked off + * as they are used. + */ +class DataIO_NetCDF::SetPool { + public: + /// CONSTRUCTOR - place sets from DataSetList in this pool + SetPool(DataSetList const& dsl) : nUsed_(0) { + sets_.reserve( dsl.size() ); + isUsed_.assign( dsl.size(), false ); + for (DataSetList::const_iterator it = dsl.begin(); it != dsl.end(); ++it) + sets_.push_back( *it ); + } + /// \return set at idx + DataSet const* Set(unsigned int idx) const { return sets_[idx]; } + /// \return Number of sets + unsigned int Nsets() const { return sets_.size(); } + /// \return true if set at idx has been used + bool IsUsed(unsigned int idx) const { return isUsed_[idx]; } + /// \return true if all sets have been marked as used + bool AllUsed() const { return (nUsed_ == isUsed_.size()); } + /// Print unused sets to stdout + void PrintUnused() const { + for (unsigned int idx = 0; idx < sets_.size(); idx++) + if (!isUsed_[idx]) mprintf("\tUnused: %s\n", sets_[idx]->legend()); + } + /// Mark set at idx as used + void MarkUsed(unsigned int idx) { + isUsed_[idx] = true; + nUsed_++; + } + private: + std::vector sets_; + std::vector isUsed_; + unsigned int nUsed_; +}; +// ----------------------------------------------------------------------------- + +/** Hold a pointer to DataSet and its original index. Used to refer back to + * original DataSetList from sets in a SetPool. + */ +class DataIO_NetCDF::Set { + public: + //Set(DataSet const* ds, int oidx) : ds_(ds), oidx_(oidx) {} + Set(DataSet const* ds) : ds_(ds) {} + + DataSet const* DS() const { return ds_; } + //int OriginalIdx() const { return oidx_; } + private: + DataSet const* ds_; + //int oidx_; +}; +// ----------------------------------------------------------------------------- + +#ifdef BINTRAJ +/// Tell netcdf file to end define mode +static inline int EndDefineMode(int ncid) { + // End netcdf definitions + if (NC::CheckErr(nc_enddef(ncid))) { + mprinterr("NetCDF data error on ending definitions."); + return 1; + } + return 0; +} + +/// Tell netcdf file to enter define mode +static inline int EnterDefineMode(int ncid) { + int err = nc_redef(ncid); + if (err != NC_NOERR && err != NC_EINDEFINE) { + NC::CheckErr( err ); + return 1; + } + return 0; +} + +/// Add a string as text attribute to a variable +static inline int AddDataSetStringAtt(std::string const& str, const char* desc, int ncid, int varid) +{ + if (!str.empty()) { + if (NC::CheckErr(nc_put_att_text(ncid, varid, desc, str.size(), str.c_str()) )) + { + mprinterr("Error: Writing attribute %s.\n", desc); + return 1; + } + } + return 0; +} + +/// Add an integer as integer attribute +static inline int AddDataSetIntAtt(int ival, const char* desc, int ncid, int varid) +{ + // NOTE: -1 means not set TODO define that in MetaData + if (ival != -1) { + if (NC::CheckErr(nc_put_att_int(ncid, varid, desc, NC_INT, 1, &ival))) + { + mprinterr("Error: Writing attribute %s.\n", desc); + return 1; + } + } + return 0; +} + +/// Add a double as double attribute +static inline int AddDataSetDblAtt(double dval, const char* desc, int ncid, int varid) +{ + if (NC::CheckErr(nc_put_att_double(ncid, varid, desc, NC_DOUBLE, 1, &dval))) + { + mprinterr("Error: Writing attribute %s.\n", desc); + return 1; + } + return 0; +} + +/// Add DataSet metadata to a variable +static inline int AddDataSetMetaData(MetaData const& meta, int ncid, int varid) +{ + // Filename + if (AddDataSetStringAtt(meta.Fname().Full(), "filename", ncid, varid)) return 1; + // Name + if (AddDataSetStringAtt(meta.Name(), "name", ncid, varid)) return 1; + // Aspect + if (AddDataSetStringAtt(meta.Aspect(), "aspect", ncid, varid)) return 1; + // Legend + if (AddDataSetStringAtt(meta.Legend(), "legend", ncid, varid)) return 1; + // Index + if (AddDataSetIntAtt(meta.Idx(), "index", ncid, varid)) return 1; + // Ensemble number + if (AddDataSetIntAtt(meta.EnsembleNum(), "ensemblenum", ncid, varid)) return 1; + // TODO TimeSeries? + // Scalar type + if (meta.ScalarType() != MetaData::UNDEFINED) { + if (AddDataSetStringAtt(meta.TypeString(), "scalartype", ncid, varid)) return 1; + } + // Scalar mode + if (meta.ScalarMode() != MetaData::UNKNOWN_MODE) { + if (AddDataSetStringAtt(meta.ModeString(), "scalarmode", ncid, varid)) return 1; + } + + return 0; +} + +/// Add DataSet index information to a target variable +/** \param ds Input DataSet + * \param indexVarIds Optional index variable IDs, one for each dimension. + * \param ncid File ncid + * \param varid Target variable id + */ +static inline int AddDataSetIndexInfo(DataSet const* ds, std::vector const& indexVarIds, + int ncid, int varid) +{ + // Add number of index variables + if (AddDataSetIntAtt( ds->Ndim(), "nindexvar", ncid, varid )) return 1; + // Loop over index variables + for (int idx = 0 ; idx < (int)ds->Ndim(); idx++) + { + Dimension const& dim = ds->Dim(idx); + std::string suffix( integerToString(idx) ); + if (!indexVarIds.empty() && indexVarIds[idx] > -1) { + // Index variable ID present for this dimension + std::string indexid( "indexid" + suffix ); + if (AddDataSetIntAtt(indexVarIds[idx], indexid.c_str(), ncid, varid)) return 1; + } else { + // No index variable ID for this dimension + std::string min( "min" + suffix); + std::string step( "step" + suffix); + if (AddDataSetDblAtt(dim.Min(), min.c_str(), ncid, varid)) return 1; + if (AddDataSetDblAtt(dim.Step(), step.c_str(), ncid, varid)) return 1; + } + std::string label("label" + suffix); + if (AddDataSetStringAtt(dim.Label(), label.c_str(), ncid, varid)) return 1; + } + return 0; +} + +/// Add DataSet MetaData, index information (with optional index var IDs), and description to target variable +/** \param ds Set with MetaData to add. + * \param indexVarIds Optional index variable IDs, one for each dimension. + * \param ncid NetCDF id to add to. + * \param varid NetCDF variable id to add to. + */ +static inline int AddDataSetInfo(DataSet const* ds, std::vector const& indexVarIds, + int ncid, int varid) +{ + // Add DataSet metadata as attributes + if (AddDataSetMetaData(ds->Meta(), ncid, varid)) return 1; + // Add index info + if (AddDataSetIndexInfo(ds, indexVarIds, ncid, varid)) return 1; + // Store the description + if (AddDataSetStringAtt(ds->description(), "description", ncid, varid)) return 1; + return 0; +} + +/// Add DataSet MetaData, index information (no index var IDs), and description to target variable +/** \param ds Set with MetaData to add. + * \param ncid NetCDF id to add to. + * \param varid NetCDF variable id to add to. + */ +static inline int AddDataSetInfo(DataSet const* ds, int ncid, int varid) { + return AddDataSetInfo(ds, std::vector(), ncid, varid); +} + +/** Define dimension. Ensure name is unique by appending an index. + * Dimension label will be '