guessmath.h 7.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297
  1. ////////////////////////////////////////////////////////////////////////////////
  2. /// \file guessmath.h
  3. /// \brief Mathematical utilites header file.
  4. ///
  5. /// This header file contains:
  6. /// (1) Definitions of constants and common functions used throughout LPJ-GUESS
  7. ///
  8. /// \author Michael Mischurow
  9. /// $Date: 2018-02-02 18:01:35 +0100 (ven, 02 fév 2018) $
  10. ///
  11. ////////////////////////////////////////////////////////////////////////////////
  12. #ifndef LPJ_GUESS_GUESSMATH_H
  13. #define LPJ_GUESS_GUESSMATH_H
  14. #include <assert.h>
  15. #include "archive.h"
  16. #define _USE_MATH_DEFINES
  17. #include <cmath>
  18. #ifndef M_PI
  19. double const PI = 4 * atan(1.0);
  20. #else
  21. double const PI = M_PI;
  22. #undef M_PI
  23. #endif
  24. const double DEGTORAD = PI / 180.;
  25. inline bool negligible(double dval, int limit = 0) {
  26. // Returns true if |dval| < EPSILON, otherwise false
  27. return limit ? fabs(dval) < pow(10.0, limit) : fabs(dval) < 1.0e-30;
  28. }
  29. inline bool largerthanzero(double dval, int limit = 0) {
  30. // Returns true if |dval| < EPSILON, otherwise false
  31. return limit ? dval > pow(10.0, limit) : dval > 1.0e-30;
  32. }
  33. inline bool equal(double dval1, double dval2) {
  34. // Returns true if |dval1-dval2| < EPSILON, otherwise false
  35. return negligible(dval1 - dval2);
  36. }
  37. inline double mean(double* array, int nitem) {
  38. // Returns arithmetic mean of 'nitem' values in 'array'
  39. double sum=0.0;
  40. for (int i=0; i<nitem; sum += array[i++]);
  41. return sum / (double)nitem;
  42. }
  43. /// Gives the mean of just two values
  44. inline double mean(double x, double y) {
  45. return (x+y)/2.0;
  46. }
  47. /// Calculates variation coefficient of values in an array
  48. inline double variation_coefficient(double* data, int n) {
  49. // 0 and 1 will give division with zero.
  50. if (n <= 1) {
  51. return -1;
  52. }
  53. double avg = mean(data, n);
  54. double dev = 0;
  55. for (int i=0; i<n; i++) {
  56. dev += (data[i] - avg) * (data[i] - avg);
  57. }
  58. double std = sqrt(dev / (n-1));
  59. if (std > 0 && avg > 0) { // check that data appear in the array
  60. return std / avg;
  61. }
  62. return 0;
  63. }
  64. /// A short version of Richards curve where:
  65. /** a is the lower asymptote,
  66. * b is the upper asymptote. If a=0 then b is called the carrying capacity,
  67. * c the growth rate,
  68. * d is the time of maximum growth
  69. * Source: https://en.wikipedia.org/wiki/Generalised_logistic_function, 2013-11-11
  70. */
  71. inline double richards_curve(double a, double b, double c, double d, double x) {
  72. return a + (b - a) / (1 + exp(-c * (x - d)));
  73. }
  74. inline void regress(double* x, double* y, int n, double& a, double& b) {
  75. // Performs a linear regression of array y on array x (n values)
  76. // returning parameters a and b in the fitted model: y=a+bx
  77. // Source: Press et al 1986, Sect 14.2
  78. double sx,sy,sxx,sxy,delta;
  79. sx = sy = sxy = sxx = 0.0;
  80. for (int i=0; i<n; i++) {
  81. sx += x[i];
  82. sy += y[i];
  83. sxx+= x[i]*x[i];
  84. sxy+= x[i]*y[i];
  85. }
  86. delta = (double)n*sxx - sx*sx;
  87. a = (sxx*sy - sx*sxy)/delta;
  88. b = ((double)n*sxy-sx*sy)/delta;
  89. }
  90. // Forward declarations needed for the friendship between Historic and operator&
  91. template<typename T, size_t capacity>
  92. class Historic;
  93. template<typename T, size_t capacity>
  94. ArchiveStream& operator&(ArchiveStream& stream,
  95. Historic<T, capacity>& data);
  96. /// Keeps track of historic values of some variable
  97. /** Useful for calculating running means etc.
  98. *
  99. * The class behaves like a queue with a fixed size,
  100. * when a new value is added, and the queue is full,
  101. * the oldest value is overwritten.
  102. */
  103. template<typename T, size_t capacity>
  104. class Historic {
  105. public:
  106. /// The maximum number of elements stored, given as template parameter
  107. static const size_t CAPACITY = capacity;
  108. Historic()
  109. : current_index(0), full(false) {
  110. }
  111. /// Adds a value, overwriting the oldest if full
  112. void add(double value) {
  113. values[current_index] = value;
  114. current_index = (current_index+1) % capacity;
  115. if (current_index == 0) {
  116. full = true;
  117. }
  118. }
  119. /// Returns the number of values stored (0-CAPACITY)
  120. size_t size() const {
  121. return full ? capacity : current_index;
  122. }
  123. /// Calculates arithmetic mean of the stored values
  124. T mean() const {
  125. const size_t nvalues = size();
  126. assert(nvalues != 0);
  127. return sum()/nvalues;
  128. }
  129. /// Sum of stored values
  130. T sum() const {
  131. T result = 0.0;
  132. const size_t nvalues = size();
  133. for (size_t i = 0; i < nvalues; ++i) {
  134. result += values[i];
  135. }
  136. return result;
  137. }
  138. /// Returnes the latest of the stored values
  139. T lastadd() const { //latest
  140. return (*this)[size()-1];
  141. }
  142. /// Returns the maximum of the stored values
  143. T max() const {
  144. T result = -999999.0;
  145. const size_t nvalues = size();
  146. for (size_t i = 0; i < nvalues; ++i) {
  147. if (values[i]>result) {
  148. result = values[i];
  149. }
  150. }
  151. return result;
  152. }
  153. /// Returns the minimum of the stored values
  154. T min() const {
  155. T result = 999999.0;
  156. const size_t nvalues = size();
  157. for (size_t i = 0; i < nvalues; ++i) {
  158. if (values[i]<result) {
  159. result = values[i];
  160. }
  161. }
  162. return result;
  163. }
  164. /// Calculates arithmetic mean of the stored values for a period from current position and backwards nsteps
  165. T periodicmean(const size_t nsteps) const {
  166. assert(nsteps != 0);
  167. if (nsteps > size()){
  168. return mean();
  169. } else {
  170. return periodicsum(nsteps)/nsteps;
  171. }
  172. }
  173. /// Sum of stored values for a period from current position and backwards nsteps
  174. T periodicsum(const size_t nsteps) const {
  175. T result = 0.0;
  176. if (nsteps >= size()){
  177. return sum();
  178. } else {
  179. for (size_t i = size()-1; i >=size()-nsteps; --i) {
  180. result += (*this)[i];
  181. }
  182. }
  183. return result;
  184. }
  185. /// Returns a single value
  186. /** The values are ordered by age, the oldest value has index 0.
  187. *
  188. * \param pos Index of value to retrieve (must be less than size())
  189. */
  190. T operator[](size_t pos) const {
  191. assert(pos < size());
  192. if (full) {
  193. return values[(current_index+pos)%capacity];
  194. }
  195. else {
  196. return values[pos];
  197. }
  198. }
  199. /// Writes all values to a plain buffer
  200. /** The values will be ordered by age, the oldest value has index 0.
  201. *
  202. * \param buffer Array to write to, must have room for at least size() values
  203. */
  204. void to_array(T* buffer) const {
  205. const size_t first_position = full ? current_index : 0;
  206. const size_t nvalues = size();
  207. for (size_t i = 0; i < nvalues; ++i) {
  208. buffer[i] = values[(first_position+i)%capacity];
  209. }
  210. }
  211. friend ArchiveStream& operator&<T, capacity>(ArchiveStream& stream,
  212. Historic<T, capacity>& data);
  213. private:
  214. /// The stored values
  215. T values[capacity];
  216. /// The next position (in the values array) to write to
  217. size_t current_index;
  218. /// Whether we've stored CAPACITY values yet
  219. bool full;
  220. };
  221. /// Serialization support for Historic
  222. /** We have a friend serialization operator instead of
  223. * implementing Serializable, to avoid overhead of a
  224. * vtable in Historic (since serialize() is virtual).
  225. *
  226. * This could be implemented without friend, by only
  227. * using size(), operator[] and add(). The order of
  228. * the elements in the buffer will however not be
  229. * preserved then, and for instance the sum() function
  230. * will then not give _exactly_ the same results
  231. * (due to limited floating point precision).
  232. */
  233. template<typename T, size_t capacity>
  234. ArchiveStream& operator&(ArchiveStream& stream,
  235. Historic<T, capacity>& data) {
  236. stream & data.values
  237. & data.current_index
  238. & data.full;
  239. return stream;
  240. }
  241. #endif // LPJ_GUESS_GUESSMATH_H