Add ability for 'datafilter' command to modify data sets.

doc/cpptraj.lyx

@@ -7589,13 +7589,17 @@ name "subsec:cpptraj-datafilter"
 \end_layout
 
 \begin_layout LyX-Code
-datafilter <dataset arg> min <min> max <max> [out <file> [name <setname>]]
+datafilter {<dataset arg> min <min> max <max> ...} [out <file>] [name <setname>]
 \begin_inset Separator latexpar
 \end_inset
 
 
 \end_layout
 
+\begin_layout LyX-Code
+           [{multi | filterset <set> [newset <newname>]}]
+\end_layout
+
 \begin_deeper
 \begin_layout Description
 <dataset
@@ -7618,7 +7622,7 @@ max
 \begin_inset space ~
 \end_inset
 
-<max> Data set and min/max cutoffs to use; can specify more than once.
+<max> Data set name and min/max cutoffs to use; can specify more than one.
 \end_layout
 
 \begin_layout Description
@@ -7634,9 +7638,37 @@ max
 \begin_inset space ~
 \end_inset
 
-<setname>] Name of filter data set.
+<setname>] Name of filter data set containing 1 when cutoffs satisfied,
+ 0 otherwise.
+\end_layout
+
+\begin_layout Description
+[multi] Filter each set separately instead of all together (creates filter
+ set for each input set).
+ Cannot be used with 'filterset'.
+\end_layout
+
+\begin_layout Description
+[filterset
+\begin_inset space ~
+\end_inset
+
+<set>] If specified, <set> will be filtered to only contain data that satisfies
+ cutoffs.
+ Cannot be used with 'multi'.
+\end_layout
+
+\begin_deeper
+\begin_layout Description
+[newset
+\begin_inset space ~
+\end_inset
+
+<newname>] If specified a new set will be created from 'filterset' instead
+ of replacing 'filterset'.
 \end_layout
 
+\end_deeper
 \end_deeper
 \begin_layout Standard
 Create a data set (optionally named <setname>) containing 1 for data within
@@ -7658,9 +7690,38 @@ Create a data set (optionally named <setname>) containing 1 for data within
 <max>
 \series default
  argument, and can be as many as there are specified data sets.
- For example, to read in data from two separate files (d1.dat and a1.dat)
- and generate a filter data set named FILTER having 1 when d1 is between
- 0.0 and 3.0 and a1 is between 135.0 and 180.0:
+ If 
+\series bold
+'multi'
+\series default
+ is specified then only filter data sets will be created for each data set
+ instead.
+ If 
+\series bold
+'filterset'
+\series default
+ is specified, the specified 
+\series bold
+<set>
+\series default
+ will be modified to only contain '1' frames; cannot be used with 
+\series bold
+'multi'
+\series default
+.
+ If 
+\series bold
+'newset'
+\series default
+ is also specified, a new set will be created containing the '1' frames
+ instead.
+ The 'filterset' functionality only works for 1D scalar sets.
+\end_layout
+
+\begin_layout Standard
+For example, to read in data from two separate files (d1.dat and a1.dat) and
+ generate a filter data set named FILTER having 1 when d1 is between 0.0
+ and 3.0 and a1 is between 135.0 and 180.0:
 \end_layout
 
 \begin_layout LyX-Code

src/Action_FilterByData.h

@@ -12,6 +12,8 @@ class Action_FilterByData : public Action {
     size_t DetermineFrames() const;
     Action::RetType Init(ArgList&, ActionInit&, int);
     Action::RetType DoAction(int, ActionFrame&);
+    /// \return the DataSet containing 1 for OK, 0 for filtered out.
+    DataSet* FilterSet() const { return maxmin_; }
   private:
     Action::RetType Setup(ActionSetup&) { return Action::OK; }
     void Print();

src/Exec_DataFilter.cpp

@@ -1,17 +1,39 @@
 #include "Exec_DataFilter.h"
-#include "CpptrajStdio.h"
 #include "Action_FilterByData.h"
+#include "CpptrajStdio.h"
+#include "DataSet_integer.h"
 #include "ProgressBar.h"
 
 void Exec_DataFilter::Help() const {
-  mprintf("\t<dataset arg> min <min> max <max> [out <file> [name <setname>]]\n"
+  mprintf("\t{<dataset arg> min <min> max <max> ...} [out <file>] [name <setname>]\n"
+          "\t[{multi | filterset <set> [newset <newname>]}]\n"
           "  Create a data set (optionally named <setname>) containing 1 for\n"
           "  data within given <min> and <max> criteria for each specified\n"
           "  data set. There must be at least one <min> and <max> argument,\n"
-          "  and can be as many as there are specified data sets.\n");
+          "  and can be as many as there are specified data sets.\n"
+          "  If 'multi' is specified then only filter data sets will be created for each\n"
+          "  data set instead.\n"
+          "  If 'filterset' is specified, the specified <set> will be modified\n"
+          "  to only contain '1' frames; cannot be used with 'multi'. If 'newset'\n"
+          "  is also specified, a new set will be created containing the '1' frames instead.\n"
+          "  The 'filterset' functionality only works for 1D scalar sets.\n");
 }
 
+// Exec_DataFilter::Execute()
 Exec::RetType Exec_DataFilter::Execute(CpptrajState& State, ArgList& argIn) {
+
+  // Get args specific to this exec.
+  std::string filterSetName = argIn.GetStringKey("filterset");
+  std::string newSetName;
+  if (!filterSetName.empty()) {
+    // Not intended to work with 'multi'
+    if (argIn.hasKey("multi")) {
+      mprinterr("Error: 'filterset' can not be used with 'multi' keyword.\n");
+      return CpptrajState::ERR;
+    }
+    newSetName = argIn.GetStringKey("newset");
+  }
+  // Init and set up Action_FilterByData
   Action_FilterByData filterAction;
   ActionInit state(State.DSL(), State.DFL());
   if (filterAction.Init(argIn, state, State.Debug()) != Action::OK)
@@ -21,13 +43,68 @@ Exec::RetType Exec_DataFilter::Execute(CpptrajState& State, ArgList& argIn) {
     mprinterr("Error: No data to filter. All sets must contain some data.\n");
     return CpptrajState::ERR;
   }
+  // Get the set to be filtered if needed.
+  DataSet_1D* SetToBeFiltered = 0;
+  DataSet_1D* FilteredSet = 0;
+  DataSet_integer* FilterSet = 0;
+  if (!filterSetName.empty()) {
+    DataSet* ds = filterAction.FilterSet();
+    if (ds == 0 || ds->Type() != DataSet::INTEGER) {
+      mprinterr("Internal Error: FilterSet was not created.\n");
+      return CpptrajState::ERR;
+    }
+    FilterSet = (DataSet_integer*)ds;
+    ds = State.DSL().GetDataSet( filterSetName );
+    if (ds == 0) {
+      mprinterr("Error: Set to be filtered '%s' not found.\n", filterSetName.c_str());
+      return CpptrajState::ERR;
+    }
+    if (ds->Group() != DataSet::SCALAR_1D)
+    {
+      mprinterr("Error: '%s' is not 1D scalar.\n", ds->legend());
+      return CpptrajState::ERR;
+    }
+    if (ds->Size() < nframes) {
+      mprinterr("Error: Set to be filtered size (%zu) is too small (%zu)\n",
+                ds->Size(), nframes);
+      return CpptrajState::ERR;
+    }
+    SetToBeFiltered = (DataSet_1D*)ds;
+    // Create new filter set
+    FilteredSet = (DataSet_1D*)DataSetList::Allocate(SetToBeFiltered->Type());
+    MetaData md(SetToBeFiltered->Meta());
+    if (!newSetName.empty()) {
+      md.SetName(newSetName);
+      ds = State.DSL().CheckForSet( md );
+      if (ds != 0) {
+        mprinterr("Error: New set name '%s' already exists.\n", newSetName.c_str());
+        return CpptrajState::ERR;
+      }
+      mprintf("\tA new filtered set will be created from set '%s'\n", SetToBeFiltered->legend());
+    } else
+      mprintf("\tFiltering set '%s'\n", SetToBeFiltered->legend());
+    FilteredSet->SetMeta( md );
+    // Do not allocate here to save memory.
+  }
+
   ProgressBar progress( nframes );
+  int newidx = 0;
   for (size_t frame = 0; frame != nframes; frame++) {
     progress.Update( frame );
     // Filter does not need frame but does need trajout num.
     ActionFrame frm(0, frame);
     filterAction.DoAction(frame, frm);
+    if (SetToBeFiltered != 0) {
+      if ((*FilterSet)[frame] == 1)
+        FilteredSet->Add(newidx++, SetToBeFiltered->VoidPtr(frame));
+    }
   }
+  // Add/replace the filtered set if necessary
+  if (SetToBeFiltered != 0) {
+    if (newSetName.empty())
+      State.DSL().RemoveSet( SetToBeFiltered );
+    State.DSL().AddSet( FilteredSet );
+  } 
   // Trigger master datafile write just in case
   State.MasterDataFileWrite();
   return CpptrajState::OK;

src/FindDepend.cpp

@@ -36,6 +36,7 @@ enum FileType { SOURCE = 0, HEADER };
 /** \return true if this header should be ignored. */
 bool IgnoreHeader(const char* headername) {
   if (strcmp(headername,"mpi.h")==0) return true;
+  if (strcmp(headername,"mkl.h")==0) return true;
   if (strcmp(headername,"zlib.h")==0) return true;
   if (strcmp(headername,"bzlib.h")==0) return true;
   if (strcmp(headername,"netcdf.h")==0) return true;

src/Version.h

@@ -12,7 +12,7 @@
  * Whenever a number that precedes <revision> is incremented, all subsequent
  * numbers should be reset to 0.
  */
-#define CPPTRAJ_INTERNAL_VERSION "V4.25.11"
+#define CPPTRAJ_INTERNAL_VERSION "V4.25.12"
 /// PYTRAJ relies on this
 #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION
 #endif

src/cpptrajdepend

@@ -251,7 +251,7 @@ Exec_CrdAction.o : Exec_CrdAction.cpp Action.h ActionList.h ActionState.h Analys
 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 AtomExtra.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 SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajectoryFile.h Trajin.h TrajinList.h TrajoutList.h Trajout_Single.h TypeNameHolder.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 AtomExtra.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 SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.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 AtomExtra.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 StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Vec3.h
-Exec_DataFilter.o : Exec_DataFilter.cpp Action.h ActionList.h ActionState.h Action_FilterByData.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h Array1D.h AssociatedData.h Atom.h AtomExtra.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_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 SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Vec3.h
+Exec_DataFilter.o : Exec_DataFilter.cpp Action.h ActionList.h ActionState.h Action_FilterByData.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h Array1D.h AssociatedData.h Atom.h AtomExtra.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_Coords.h DataSet_Coords_REF.h DataSet_integer.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 SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.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 AtomExtra.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_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 Spline.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.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 AtomExtra.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 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 SymbolExporting.h TextFormat.h Timer.h Topology.h TorsionRoutines.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Vec3.h
 Exec_Help.o : Exec_Help.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h BondSearch.h Box.h Cmd.h CmdList.h Command.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_Help.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 ParmFile.h ParmIO.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajectoryFile.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Vec3.h

test/Test_Filter/A2.filtered.dat.save

@@ -0,0 +1,7 @@
+#                 A:2
+       1       0.0000
+       2       0.6360
+       3       0.8260
+       4       0.9030
+       5       0.8820
+       6       1.1830

test/Test_Filter/RunTest.sh

@@ -3,7 +3,7 @@
 . ../MasterTest.sh
 
 # Clean
-CleanFiles filter.in filter.crd filter.dat datafilter.dat
+CleanFiles filter.in filter.crd filter.dat datafilter.dat A2.filtered.dat
 
 INPUT='filter.in'
 TOP='../tz2.truncoct.parm7'
@@ -30,9 +30,12 @@ fi
 cat > filter.in <<EOF
 readdata ../Test_Diffusion/diff_a.xmgr.save index 1 name A as dat
 datafilter A:2 A:3 min 0.0 max 1.2 out datafilter.dat multi name FA
+datafilter A:2 min 0.0 max 1.2 filterset A:2
+writedata A2.filtered.dat A:2
 EOF
 RunCpptraj "Data filter test."
 DoTest datafilter.dat.save datafilter.dat
+DoTest A2.filtered.dat.save A2.filtered.dat
 
 EndTest