utils.cpp
Go to the documentation of this file.
1 
9 #include <target/utils.hpp>
10 #include <algorithm>
11 
12 namespace target {
13 
14  // Complex step derivative
15  arma::mat deriv(cx_func f, arma::vec theta) {
16  arma::cx_vec thetac = arma::conv_to<arma::cx_vec>::from(theta);
17  arma::cx_mat val0 = f(thetac);
18  unsigned n = val0.n_elem;
19  unsigned p = theta.n_elem;
20  arma::mat res(n, p);
21  double h = DBL_MIN;
22  cx_dbl h0 = cx_dbl(0, h);
23  for (unsigned i=0; i < p; i++) {
24  arma::cx_vec theta0 = thetac;
25  theta0[i] += h0;
26  arma::mat val = imag(f(theta0))/h;
27  for (unsigned j=0; j < n; j++)
28  res(j, i) = val[j];
29  }
30  return(res);
31  }
32 
33 
34  // Find index of clusters/ids in sorted integer vector
35  arma::umat clusterid(const arma::uvec &id) {
36  unsigned ncl = 1;
37  unsigned n = id.size();
38  unsigned id0 = id(0);
39  for (unsigned i=0; i < n; i++) {
40  if (id(i) != id0) {
41  ncl++;
42  id0 = id(i);
43  }
44  }
45  arma::umat res(ncl, 2); // column 1: index, column 2: size
46  res.fill(0);
47  unsigned cl = 0;
48  id0 = id(0);
49  for (unsigned i=0; i < n; i++) {
50  if (id(i) != id0) {
51  cl++;
52  res(cl, 0) = i; // index
53  id0 = id(i);
54  }
55  res(cl, 1)++; // size
56  }
57  return res;
58  }
59 
60 
61  arma::mat groupsum(const arma::mat &x,
62  const arma::uvec &cluster,
63  bool reduce) {
64  unsigned ncl = cluster.n_elem;
65  unsigned n = ncl;
66  if (!reduce) {
67  n = x.n_rows;
68  }
69  unsigned stop;
70  arma::mat res(n, x.n_cols);
71  arma::rowvec tmp(x.n_cols);
72  for (unsigned i=0; i < ncl; i++) { // Iterate over individuals
73  unsigned start = cluster(i);
74  if (i == (ncl-1)) {
75  stop = x.n_rows;
76  } else {
77  stop = cluster(i+1);
78  }
79  tmp.fill(0);
80  for (unsigned j=start; j < stop; j++) {
81  tmp += x.row(j);
82  }
83  if (reduce) {
84  res.row(i) = tmp;
85  } else {
86  for (unsigned j=start; j < stop; j++) {
87  res.row(j) = tmp;
88  }
89  }
90  }
91  return(res);
92  }
93 
94  arma::mat interpolate(const arma::mat &input, double tau, bool locf) {
95  arma::vec time = input.col(0);
96  unsigned n = time.n_elem;
97  double t0 = time(0);
98  double tn = time(n-1);
99  unsigned N = std::ceil((tn-t0)/tau)+1;
100  arma::mat input2(N, input.n_cols);
101  unsigned cur = 0;
102  input2.row(0) = input.row(0);
103  double curtime = t0;
104  arma::rowvec slope(input.n_cols);
105  if (locf) {
106  slope.fill(0); slope(0) = 1;
107  } else {
108  slope = (input.row(cur+1)-input.row(cur))/(time(cur+1)-time(cur));
109  }
110  for (unsigned i=0; i < N-1; i++) {
111  while (time(cur+1) < curtime) {
112  cur++;
113  if (cur == (n-1)) break;
114  if (!locf)
115  slope = (input.row(cur+1)-input.row(cur))/(time(cur+1)-time(cur));
116  }
117  double delta = curtime-time(cur);
118  input2.row(i) = input.row(cur) + slope*delta;
119  curtime += tau;
120  }
121  tau = tn-input2(N-2, 0);
122  input2.row(N-1) = input.row(input.n_rows-1);
123  return( input2 );
124  }
125 
126 
127  void fastpattern(const arma::umat &y,
128  arma::umat &pattern,
129  arma::uvec &group,
130  unsigned categories) {
131  unsigned n = y.n_rows;
132  unsigned k = y.n_cols;
133 
134  arma::uvec mygroup(n);
135  double lognpattern = k*std::log(static_cast<double>(categories));
136  unsigned npattern = n;
137  if (lognpattern<std::log(static_cast<double>(npattern))) {
138  npattern = (unsigned) pow(static_cast<double>(categories),
139  static_cast<double>(k));
140  }
141  arma::umat mypattern(npattern, k);
142  mypattern.fill(1);
143  unsigned K = 0;
144 
145  for (unsigned i=0; i < n; i++) {
146  arma::urowvec ri = y.row(i);
147  bool found = false;
148  for (unsigned j=0; j < K; j++) {
149  if (sum(ri != mypattern.row(j)) == 0) {
150  found = true;
151  mygroup(i) = j;
152  break;
153  }
154  }
155  if (!found) {
156  K++;
157  mypattern.row(K-1) = ri;
158  mygroup(i) = K-1;
159  }
160  }
161  pattern = mypattern.rows(0, K-1);
162  group = mygroup;
163  }
164 
165  arma::umat fastapprox(arma::vec time, // sorted times
166  const arma::vec &newtime,
167  bool equal,
168  // type: (0: nearedst, 1: right, 2: left)
169  unsigned type) {
170  arma::uvec idx(newtime.n_elem);
171  arma::uvec eq(newtime.n_elem);
172 
173  double vmax = time(time.n_elem-1);
174  arma::mat::iterator it;
175  double upper = 0.0;
176  int pos;
177  for (unsigned i=0; i < newtime.n_elem; i++) {
178  eq[i] = 0;
179  if (newtime(i) > vmax) {
180  pos = time.n_elem-1;
181  } else {
182  it = std::lower_bound(time.begin(), time.end(), newtime(i));
183  upper = *it;
184  if (it == time.begin()) {
185  pos = 0;
186  if (equal && (newtime(i) == upper)) { eq(i) = 1; }
187  } else {
188  pos = static_cast<int>(it-time.begin());
189  if (type == 0 && std::fabs(newtime(i)-time(pos-1)) <
190  std::fabs(newtime(i)-time(pos)))
191  pos -= 1;
192  if (equal && (newtime(i) == upper)) { eq[i] = pos+1; }
193  }
194  }
195  if (type == 2 && newtime(i) < upper) pos--;
196  idx(i) = pos;
197  }
198  if (equal) {
199  return( arma::join_horiz(idx, eq) );
200  }
201  return( std::move(idx) );
202  }
203 
204 
205  double SupTest(const arma::vec &D) {
206  return arma::max(arma::abs(D));
207  }
208 
209  double L2Test(const arma::vec &D, const arma::vec &t) {
210  arma::vec delta(t.n_elem);
211  for (unsigned i=0; i < t.n_elem-1; i++) delta(i) = t[i+1]-t[i];
212  delta(delta.n_elem-1) = 0;
213  return std::sqrt(sum(delta % D % D));
214  }
215 
216  // Comparison of ECDF of (x1,...,xn) with null CDF
217  // evaluated in G = (G(x1),...,G(xn))
218  double CramerVonMises(const arma::vec &x, arma::vec G) {
219  // back to original order of input data:
220  arma::uvec ord = arma::stable_sort_index(x);
221  G = G.elem(ord);
222  unsigned n = G.n_elem;
223  double res = 1/(double)(12*n);
224  for (unsigned i=0; i < G.n_elem; i++) {
225  double val = (2*i-1)/static_cast<double>((2*n)-G(i));
226  res += val*val;
227  }
228  return res;
229  }
230 
231 
232  arma::mat const EmptyMat = arma::mat();
233  arma::vec const EmptyVec = arma::vec();
234 
235 
236  // Foreground colors are in form of 3x, bacground are 4x
237  const char* COL_RESET = "\x1b[0m";
238  const char* COL_DEF = "\x1b[39m";
239  const char* BLACK = "\x1b[30m";
240  const char* RED = "\x1b[31m";
241  const char* MAGENTA = "\x1b[35m";
242  const char* YELLOW = "\x1b[33m";
243  const char* GREEN = "\x1b[32m";
244  const char* BLUE = "\x1b[34m";
245  const char* CYAN = "\x1b[36m";
246  const char* WHITE = "\x1b[37m";
247  const char* GRAY = "\x1b[90m";
248  const char* LRED = "\x1b[91m";
249  const char* LGREEN = "\x1b[92m";
250  const char* LYELLOW = "\x1b[93m";
251  const char* LBLUE = "\x1b[94m";
252  const char* LMAGENTA = "\x1b[95m";
253  const char* LCYAN = "\x1b[96m";
254  const char* LWHITE = "\x1b[97m";
255 
256 } // namespace target
Various utility functions and constants.