generic.cpp
Go to the documentation of this file.
1 //-----------------------------------------------------------------------------------
2 // d-SEAMS is free software: you can redistribute it and/or modify
3 // it under the terms of the GNU General Public License as published by
4 // the Free Software Foundation, either version 3 of the License, or
5 // (at your option) any later version.
6 //
7 // A copy of the GNU General Public License is available at
8 // http://www.gnu.org/licenses/
9 //-----------------------------------------------------------------------------------
10 
11 #include <generic.hpp>
12 #include <iostream>
13 
23  std::string outFile) {
24  std::ofstream outputFile;
25  // Create a new file in the output directory
26  outputFile.open(outFile);
27 
28  if (outputFile.is_open()) {
29  // First line
30  outputFile << "# Frame\tAtomID\tx\ty\tz\tcij\ticeType\n";
31  // Write out all the information out line by line
32  for (int i = 0; i < yCloud->nop; i++) {
33  outputFile << yCloud->currentFrame << "\t" << yCloud->pts[i].atomID
34  << "\t" << yCloud->pts[i].x << "\t" << yCloud->pts[i].y << "\t"
35  << yCloud->pts[i].z << "\t";
36  // Print out cij
37  // for(int c=0; c<yCloud->pts[i].c_ij.size(); c++){outputFile <<
38  // yCloud->pts[i].c_ij[c]<<"\t";} Print out the classifier
39  outputFile << yCloud->pts[i].iceType << "\n";
40  }
41  }
42  // Close the file
43  outputFile.close();
44  return 0;
45 }
46 
67  molSys::PointCloud<molSys::Point<double>, double> *yCloud, int iatomIndex,
68  int jatomIndex, double *x_i, double *y_i, double *z_i, double *x_j,
69  double *y_j, double *z_j) {
70  //
71  double x_iatom, y_iatom, z_iatom;
72  double x_jatom, y_jatom, z_jatom;
73  double x_ij, y_ij, z_ij; // Relative distance
74  std::vector<double> box = yCloud->box;
75  double xPBC, yPBC, zPBC; // Actual unwrapped distance
76 
77  // ----------------------------------------------------------------------
78  // INIT
79  // iatom
80  x_iatom = yCloud->pts[iatomIndex].x;
81  y_iatom = yCloud->pts[iatomIndex].y;
82  z_iatom = yCloud->pts[iatomIndex].z;
83  // jatom
84  x_jatom = yCloud->pts[jatomIndex].x;
85  y_jatom = yCloud->pts[jatomIndex].y;
86  z_jatom = yCloud->pts[jatomIndex].z;
87  // ----------------------------------------------------------------------
88  // GET RELATIVE DISTANCE
89  x_ij = x_iatom - x_jatom;
90  y_ij = y_iatom - y_jatom;
91  z_ij = z_iatom - z_jatom;
92  // ----------------------------------------------------------------------
93  // SHIFT COORDINATES IF REQUIRED
94  // Shift x
95  if (fabs(x_ij) > 0.5 * box[0]) {
96  // Get the actual distance
97  xPBC = box[0] - fabs(x_ij);
98  if (x_ij < 0) {
99  x_jatom = x_iatom - xPBC;
100  } // To the -x side of currentIndex
101  else {
102  x_jatom = x_iatom + xPBC;
103  } // Add to the + side
104  } // Shift nextElement
105  //
106  // Shift y
107  if (fabs(y_ij) > 0.5 * box[1]) {
108  // Get the actual distance
109  yPBC = box[1] - fabs(y_ij);
110  if (y_ij < 0) {
111  y_jatom = y_iatom - yPBC;
112  } // To the -y side of currentIndex
113  else {
114  y_jatom = y_iatom + yPBC;
115  } // Add to the + side
116  } // Shift nextElement
117  //
118  // Shift z
119  if (fabs(z_ij) > 0.5 * box[2]) {
120  // Get the actual distance
121  zPBC = box[2] - fabs(z_ij);
122  if (z_ij < 0) {
123  z_jatom = z_iatom - zPBC;
124  } // To the -z side of currentIndex
125  else {
126  z_jatom = z_iatom + zPBC;
127  } // Add to the + side
128  } // Shift nextElement
129  // ----------------------------------------------------------------------
130  // Assign values
131  *x_i = x_iatom;
132  *y_i = y_iatom;
133  *z_i = z_iatom;
134  *x_j = x_jatom;
135  *y_j = y_jatom;
136  *z_j = z_jatom;
137 
138  return 0;
139 }
140 
151  Eigen::Vector3d eigOO = Eigen::Map<Eigen::Vector3d>(OO.data(), OO.size());
152  Eigen::Vector3d eigOH = Eigen::Map<Eigen::Vector3d>(OH.data(), OH.size());
153  double angle;
154  angle = acos(eigOO.dot(eigOH) / (eigOO.norm() * eigOH.norm()));
155  return angle;
156 }
157 
165  //
166  double avgVal = 0.0; // Average value, excluding the outliers
167  double median; // Median value
168  int n = inpVec.size(); // Number of values
169  std::vector<double> lowerRange,
170  upperRange; // n/2 smallest and n/2 largest numbers
171  double firstQuartile, thirdQuartile; // First and third quartiles
172  double iqr; // Interquartile range
173  double outlierLimLow, outlierLimHigh; // Outliers limit
174  int numOfObservations = 0; // Number of observations used for the average
175  // ----------------------
176  // Calculate the median (the vector is sorted inside the function)
177  median = calcMedian(&inpVec);
178  // ----------------------
179  // Get the n/2 smallest and largest numbers
180  //
181  if (n % 2 == 0) {
182  for (int i = 0; i < n / 2; i++) {
183  // n/2 smallest numbers
184  lowerRange.push_back(inpVec[i]);
185  // n/2 largest numbers
186  upperRange.push_back(inpVec[n / 2 + i]);
187  } // end of loop to fill up the n/2 smallest and n/2 largest
188  } // even
189  else {
190  //
191  int halfN = (n + 1) / 2;
192  // Exclude the median
193  for (int i = 0; i < halfN; i++) {
194  // (n+1)/2 smallest numbers
195  lowerRange.push_back(inpVec[i]);
196  // (n+1)/2 largest numbers
197  upperRange.push_back(inpVec[halfN + i]);
198  } // end of filling up the smallest and largest half-ranges
199  } // for odd numbers
200  // ----------------------
201  // Calculate the first and third quartiles, and interquartile range
202  //
203  // First quartile
204  firstQuartile = calcMedian(&lowerRange);
205  // Third quartile
206  thirdQuartile = calcMedian(&upperRange);
207  // Interquartile range
208  iqr = thirdQuartile - firstQuartile;
209  // ----------------------
210  // Calculate the average without outliers
211  // Outliers are defined as values which
212  // are less than Q1 - 1.5IQR
213  // or greater than Q3 + 1.5IQR
214  //
215  // Get the limits for the outliers
216  outlierLimLow = firstQuartile - 1.5 * iqr;
217  outlierLimHigh = thirdQuartile + 1.5 * iqr;
218  //
219  // Loop through the values in inpVec to get the average, excluding outliers
220  for (int i = 0; i < n; i++) {
221  //
222  if (inpVec[i] < outlierLimLow) {
223  continue;
224  } // lower limit outlier
225  else if (inpVec[i] > outlierLimHigh) {
226  continue;
227  } // higher limit outlier
228  else {
229  // Number of observations added
230  numOfObservations++;
231  // Add to the average
232  avgVal += inpVec[i];
233  } // take the average
234  } // end of loop for getting the average
235  //
236  // Divide by the number of observations used
237  avgVal /= numOfObservations;
238  // ----------------------
239  // This fails if there are not enough observations (ring size = 3)
240  if (numOfObservations == 0) {
241  double sumVal = 0.0;
242  // Loop through all the values and sum
243  for (int i = 0; i <= n; i++) {
244  sumVal += inpVec[i];
245  } // end of sum
246  // Normalize
247  avgVal = sumVal / n;
248  return avgVal;
249  } // for triangular prism blocks
250  // ----------------------
251 
252  return avgVal;
253 }
254 
263  std::vector<double> quat2) {
264  //
265  double prod; // Product of quat1 and conjugate of quat2
266  // The angular distance is
267  // angularDistance = 2*cosInverse(quat1*conj(quat2))
268  prod = quat1[0] * quat2[0] - quat1[1] * quat2[1] - quat1[2] * quat2[2] -
269  quat1[3] * quat2[3];
270  // The angular distance is:
271  double angDist = 2 * acos(prod) * 180.0 / (gen::pi);
272  // Return the angular distance
273  return angDist;
274 }
T close(T... args)
T data(T... args)
File for containing generic or common functions.
double eigenVecAngle(std::vector< double > OO, std::vector< double > OH)
Eigen function for getting the angle (in radians) between the O–O and O-H vectors.
Definition: generic.cpp:150
double calcMedian(std::vector< double > *input)
Inline generic function for calculating the median given a vector of the values.
Definition: generic.hpp:73
const double pi
Definition: generic.hpp:50
double angDistDegQuaternions(std::vector< double > quat1, std::vector< double > quat2)
Definition: generic.cpp:262
int unwrappedCoordShift(molSys::PointCloud< molSys::Point< double >, double > *yCloud, int iatomIndex, int jatomIndex, double *x_i, double *y_i, double *z_i, double *x_j, double *y_j, double *z_j)
Shift particles (unwrapped coordinates)
Definition: generic.cpp:66
int prettyPrintYoda(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::string outFile)
Generic function for printing all the struct information.
Definition: generic.cpp:21
double getAverageWithoutOutliers(std::vector< double > inpVec)
Get the average, after excluding the outliers, using quartiles.
Definition: generic.cpp:164
T is_open(T... args)
T open(T... args)
T push_back(T... args)
T size(T... args)
This contains a collection of points; contains information for a particular frame.
Definition: mol_sys.hpp:166
This contains per-particle information.
Definition: mol_sys.hpp:145