Latin Hypercube Samples (lhs)  1.0
R, C++, and Rcpp code to generate Latin hypercube samples
All Classes Namespaces Files Functions Variables Typedefs Macros Pages
utilityLHS.h
Go to the documentation of this file.
1 
21 #ifndef UTILITYLHS_H
22 #define UTILITYLHS_H
23 
24 #include "LHSCommonDefines.h"
25 
26 namespace lhslib
27 {
33  bool isValidLHS(const bclib::matrix<int> & result);
39  bool isValidLHS(const bclib::matrix<double> & result);
40 
47  template <class T>
48  void rank(const std::vector<T> & toRank, std::vector<int> & ranks)
49  {
50  if (toRank.size() != ranks.size())
51  {
52  ranks.resize(toRank.size(), 0);
53  }
54  typename std::vector<T>::const_iterator toRank_it1;
55  typename std::vector<T>::const_iterator toRank_it2;
56  std::vector<int>::iterator ranks_it;
57  for (toRank_it1 = toRank.begin(), ranks_it = ranks.begin();
58  toRank_it1 != toRank.end() && ranks_it != ranks.end();
59  ++toRank_it1, ++ranks_it)
60  {
61  *ranks_it = 0;
62  for (toRank_it2 = toRank.begin(); toRank_it2 != toRank.end(); ++toRank_it2)
63  {
64  if (*toRank_it1 < *toRank_it2)
65  {
66  (*ranks_it)++;
67  }
68  }
69  }
70  }
71 
76  void initializeAvailableMatrix(bclib::matrix<int> & avail);
77 
83  template <class T>
84  void lhsPrint(const bclib::matrix<T> & A)
85  {
86  PRINT_MACRO << "\n";
87  msize_type cols = A.colsize();
88  msize_type rows = A.rowsize();
89  for (msize_type irow = 0; irow < rows; irow++)
90  {
91  for (msize_type jcol = 0; jcol < cols; jcol++)
92  {
93  PRINT_MACRO << A(irow, jcol) << ", ";
94  }
95  PRINT_MACRO << "\n";
96  }
97  }
98 
104  template <class T>
105  struct squareDifference : public std::binary_function<T, T, T> /*arg1, arg2, return */
106  {
113  T operator()(const T & x, const T & y) const
114  {
115  return (x-y) * (x-y);
116  }
117  };
118 
126  template <class T>
127  T calculateDistanceSquared(const std::vector<T> & A, const std::vector<T> & B)
128  {
129  if (A.size() != B.size())
130  {
131  throw std::runtime_error("Inputs of a different size");
132  }
133  // sum = sum + (a-b)*(a-b)
134  T sum = std::inner_product(A.begin(), A.end(), B.begin(), static_cast<T>(0), std::plus<T>(), squareDifference<T>());
135  return sum;
136  }
137 
151  template <class T, bool ISROWWISE>
152  T calculateDistanceSquared(const typename bclib::matrixConstIter<T,ISROWWISE> Abegin,
153  const typename bclib::matrixConstIter<T,ISROWWISE> Aend,
154  const typename bclib::matrixConstIter<T,ISROWWISE> Bbegin)
155  {
156  // sum = sum + (a-b)*(a-b)
157  T sum = std::inner_product(Abegin, Aend, Bbegin, static_cast<T>(0), std::plus<T>(), squareDifference<T>());
158  return sum;
159  }
160 
167  template <class T>
168  void calculateDistance(const bclib::matrix<T> & mat, bclib::matrix<double> & result)
169  {
170  msize_type m_rows = mat.rowsize();
171  if (result.rowsize() != m_rows || result.colsize() != m_rows)
172  {
173  result = bclib::matrix<double>(m_rows, m_rows);
174  }
175  for (msize_type i = 0; i < m_rows - 1; i++)
176  {
177  for (msize_type j = i+1; j < m_rows; j++)
178  {
179  typename bclib::matrix<T>::const_rowwise_iterator rowi_begin = mat.rowwisebegin(i);
180  typename bclib::matrix<T>::const_rowwise_iterator rowi_end = mat.rowwiseend(i);
181  typename bclib::matrix<T>::const_rowwise_iterator rowj_begin = mat.rowwisebegin(j);
182  T sum = calculateDistanceSquared<T, true>(rowi_begin, rowi_end, rowj_begin);
183  result(i,j) = sqrt(static_cast<double>(sum));
184  }
185  }
186  }
187 
193  template <class T, class W>
194  struct invert : public std::unary_function<T, W> /*arg1, return */
195  {
201  W operator()(const T & x) const
202  {
203  if (x != static_cast<T>(0))
204  {
205  return 1.0 / static_cast<W>(x);
206  }
207  else
208  {
209  return static_cast<W>(x);
210  }
211  }
212  };
213 
220  template <class T>
221  double sumInvDistance(const bclib::matrix<T> & A)
222  {
223  // create a matrix to hold the distances
224  bclib::matrix<double> dist = bclib::matrix<double>(A.rowsize(), A.rowsize());
225  // calculate the distances between the rows of A
226  calculateDistance<T>(A, dist);
227  // invert all the distances
228  std::transform<bclib::matrix<double>::iterator>(dist.begin(), dist.end(),
229  dist.begin(), invert<double, double>());
230  // sum the inverted
231  double totalInvDistance = std::accumulate<bclib::matrix<double>::iterator>(dist.begin(), dist.end(), 0.0);
232  return totalInvDistance;
233  }
234 
241  template <class T>
242  double sumInvDistance_deprecated(const bclib::matrix<T> & A)
243  {
244  msize_type nr = A.rowsize();
245  msize_type nc = A.colsize();
246  T oneDistance;
247  T diff;
248  double totalInvDistance = 0.0;
249  /* iterate the row of the first point from 0 to N-2 */
250  for (msize_type irow = 0; irow < nr - 1; irow++)
251  {
252  /* iterate the row the second point from i+1 to N-1 */
253  for (msize_type jrow = (irow + 1); jrow < nr; jrow++)
254  {
255  oneDistance = static_cast<T>(0);
256  /* iterate through the columns, summing the squared differences */
257  for (msize_type kcol = 0; kcol < nc; kcol++)
258  {
259  /* calculate the square of the difference in one dimension between the
260  * points */
261  diff = A(irow,kcol) - A(jrow,kcol);
262  oneDistance += diff * diff;
263  }
264  /* sum the inverse distances */
265  if (oneDistance != 0)
266  {
267  totalInvDistance += (1.0 / sqrt(static_cast<double>(oneDistance)));
268  }
269  }
270  }
271  return totalInvDistance;
272  }
273 
280  template <class T>
281  void copyMatrix(bclib::matrix<T> & copyTo, const bclib::matrix<T> & copyFrom)
282  {
283  if (copyFrom.rowsize() != copyTo.rowsize() ||
284  copyFrom.colsize() != copyTo.colsize() ||
285  copyFrom.isTransposed() != copyTo.isTransposed())
286  {
287  throw std::runtime_error("Matrices are not compatible for a copy");
288  }
289  std::copy<typename bclib::matrix<T>::const_iterator, typename bclib::matrix<T>::iterator>(copyFrom.begin(), copyFrom.end(), copyTo.begin());
290  }
291 
298  template <class T>
299  double calculateSOptimal(const bclib::matrix<T> & mat)
300  {
301  // B[i] <- 1/sum(1/dist(A[, , i]))
302  double sum = sumInvDistance<T>(mat);
303  return 1.0 / sum;
304  }
305 
312  void runif_std(unsigned int n, std::vector<double> & output, bclib::CRandom<double> & oRandom);
313 
323  template <class T1>
324  void runifint(unsigned int n, T1 min, T1 max, std::vector<T1> & output, bclib::CRandom<double> & oRandom)
325  {
326  if (output.size() != n)
327  {
328  output.resize(n);
329  }
330  std::vector<double> r = std::vector<double>(n);
331  runif_std(n, r, oRandom);
332  typename std::vector<T1>::iterator output_it;
333  std::vector<double>::iterator r_it;
334  double range = static_cast<double>(max) + 1.0 - static_cast<double>(min);
335  for (output_it = output.begin(), r_it = r.begin();
336  output_it != output.end() && r_it != r.end(); ++output_it, ++r_it)
337  {
338  *output_it = min + static_cast<T1>(floor((*r_it) * range));
339  }
340  }
341 
350  template <class T1>
351  void runifint(T1 min, T1 max, T1 * output, bclib::CRandom<double> & oRandom)
352  {
353  double r = oRandom.getNextRandom();
354  double range = static_cast<double>(max) + 1.0 - static_cast<double>(min);
355  *output = min + static_cast<T1>(floor((r * range)));
356  }
357 
358 } // end namespace
359 
360 #endif /* UTILITYLHS_H */
361 
lhslib::isValidLHS
bool isValidLHS(const bclib::matrix< int > &result)
Definition: utilityLHS.cpp:25
PRINT_MACRO
#define PRINT_MACRO
Definition: LHSCommonDefines.h:45
lhslib::squareDifference
Definition: utilityLHS.h:106
lhslib::calculateSOptimal
double calculateSOptimal(const bclib::matrix< T > &mat)
Definition: utilityLHS.h:299
lhslib::copyMatrix
void copyMatrix(bclib::matrix< T > &copyTo, const bclib::matrix< T > &copyFrom)
Definition: utilityLHS.h:281
lhslib::sumInvDistance
double sumInvDistance(const bclib::matrix< T > &A)
Definition: utilityLHS.h:221
lhslib::calculateDistanceSquared
T calculateDistanceSquared(const std::vector< T > &A, const std::vector< T > &B)
Definition: utilityLHS.h:127
LHSCommonDefines.h
lhslib::runif_std
void runif_std(unsigned int n, std::vector< double > &output, bclib::CRandom< double > &oRandom)
Definition: utilityLHS.cpp:82
lhslib::runifint
void runifint(unsigned int n, T1 min, T1 max, std::vector< T1 > &output, bclib::CRandom< double > &oRandom)
Definition: utilityLHS.h:324
lhslib::calculateDistance
void calculateDistance(const bclib::matrix< T > &mat, bclib::matrix< double > &result)
Definition: utilityLHS.h:168
lhslib
Definition: geneticLHS.cpp:25
lhslib::lhsPrint
void lhsPrint(const bclib::matrix< T > &A)
Definition: utilityLHS.h:84
lhslib::squareDifference::operator()
T operator()(const T &x, const T &y) const
Definition: utilityLHS.h:113
bclib::CRandom::getNextRandom
virtual T getNextRandom()=0
lhslib::invert
Definition: utilityLHS.h:195
lhslib::invert::operator()
W operator()(const T &x) const
Definition: utilityLHS.h:201
lhslib::msize_type
bclib::matrix< int >::size_type msize_type
Definition: LHSCommonDefines.h:114
lhslib::rank
void rank(const std::vector< T > &toRank, std::vector< int > &ranks)
Definition: utilityLHS.h:48
lhslib::initializeAvailableMatrix
void initializeAvailableMatrix(bclib::matrix< int > &avail)
Definition: utilityLHS.cpp:70
bclib::CRandom< double >
lhslib::sumInvDistance_deprecated
double sumInvDistance_deprecated(const bclib::matrix< T > &A)
Definition: utilityLHS.h:242