src/sk_hist.c

   1 /* Apache 2.0 INS-AMU 2015 */
   2
   3 #include "stdlib.h"
   4 #include "string.h"
   5 #include "math.h"
   6
   7 #include "sk_util.h"
   8 #include "sk_hist.h"
   9 #include "sk_malloc.h"
  10
  11 struct sk_hist
  12 {
  13         /*           nu+1   nu    nu    nu   nd, maxvi */
  14         int nd, nu, *lim, *len, *pos, *uvi, *vi, *vi2i, maxvi;
  15         /*   sum(len)  nu     nd                */
  16         double *buf, *maxd, *del, dt, t;
  17 };
  18
  19 sk_hist *sk_hist_alloc() {
  20         return sk_malloc(sizeof(sk_hist));
  21 }
  22
  23 int sk_hist_get_maxvi(sk_hist *h) {
  24         return h->maxvi;
  25 }
  26
  27 int sk_hist_get_vi2i(sk_hist *h, int vi) {
  28         return h->vi2i[vi];
  29 }
  30
  31 static void setup_buffer_structure(sk_hist *h)
  32 {
  33         int i, j, ui;
  34         double maxd;
  35
  36         h->maxd = sk_malloc (sizeof(double) * h->nu);
  37         h->lim = sk_malloc (sizeof(double) * (h->nu + 1));
  38         h->len = sk_malloc (sizeof(double) * h->nu);
  39         h->pos = sk_malloc (sizeof(double) * h->nu);
  40
  41         /* vi2i requires max(uvi) then filling in vi2i[ui]=i */
  42         h->maxvi = 0;
  43         for (i=0; i<h->nu; i++)
  44                 if (h->uvi[i] > h->maxvi)
  45                         h->maxvi = h->uvi[i];
  46
  47         h->vi2i = sk_malloc (sizeof(int) * (h->maxvi + 1));
  48
  49         for (i=0; i<h->nu; i++)
  50         {
  51                 ui = h->uvi[i];
  52                 h->vi2i[ui] = i;
  53                 maxd = 0.0;
  54                 for (j=0; j<h->nd; j++)
  55                         if (h->vi[j]==ui && h->del[j]>maxd)
  56                                 maxd = h->del[j];
  57                 h->maxd[i] = maxd;
  58                 h->len[i] = ceil(maxd / h->dt) + 2;
  59                 h->pos[i] = 0;
  60                 if (i==0)
  61                         h->lim[i] = 0;
  62                 else
  63                         h->lim[i] = h->lim[i-1] + h->len[i-1];
  64         }
  65         h->lim[h->nu] = h->lim[h->nu-1] + h->len[h->nu-1];
  66 }
  67
  68 void sk_hist_init(sk_hist * restrict h, int nd, int * restrict vi, double * restrict vd, double t0, double dt)
  69 {
  70         h->dt = dt;
  71         h->t = t0;
  72         h->nd = nd;
  73
  74         /* identify unique delayed variable indices
  75          * NB: this allocates h->uci
  76          */
  77         sk_util_uniqi(nd, vi, &(h->nu), &(h->uvi));
  78
  79         /* alloc & copy the delay values */
  80         h->del = sk_malloc (sizeof(double) * nd);
  81         memcpy(h->del, vd, nd * sizeof(double));
  82
  83         /* and indices */
  84         h->vi = sk_malloc (sizeof(int) * nd);
  85         memcpy(h->vi, vi, nd * sizeof(int));
  86
  87         /* setup maxd, off, len, pos */
  88         setup_buffer_structure(h);
  89
  90         /* alloc buf */
  91         h->buf = sk_malloc (sizeof(double) * h->lim[h->nu]);
  92 }
  93
  94 void sk_hist_free(sk_hist *h)
  95 {
  96         sk_free(h->uvi);
  97         sk_free(h->del);
  98         sk_free(h->vi);
  99         sk_free(h->maxd);
 100         sk_free(h->lim);
 101         sk_free(h->len);
 102         sk_free(h->pos);
 103         sk_free(h->vi2i);
 104         sk_free(h->buf);
 105         sk_free(h);
 106 }
 107
 108 void sk_hist_fill(sk_hist * restrict h, sk_hist_filler filler, void * restrict fill_data)
 109 {
 110         int i, j, ui, o, n, *vi;
 111         double *t;
 112
 113         n = h->lim[h->nu];
 114         t = sk_malloc (sizeof(double) * n);
 115         vi = sk_malloc (sizeof(int) * n);
 116
 117         /* expand indices per buffer element */
 118         for (i=0; i<h->nu; i++)
 119         {
 120                 ui = h->uvi[i];
 121                 for (j=h->lim[i]; j<h->lim[i+1]; j++)
 122                         vi[j] = ui;
 123         }
 124
 125         /* evaluate time per buffer element */
 126         for (i=0; i<h->nu; i++)
 127         {
 128                 for (j=0; j<(h->len[i]-1); j++)
 129                 {
 130                         o = h->pos[i];     /* current position in buffer */
 131                         o += j;            /* time step through buffer */
 132                         o %= h->len[i];    /* wrap around on section length */
 133                         o += h->lim[i];    /* offset to start of section */
 134                         t[o] = h->t - j * h->dt;
 135                 }
 136                 j = h->len[i] - 1;
 137                 o = h->pos[i];     /* current position in buffer */
 138                 o += j;            /* time step through buffer */
 139                 o %= h->len[i];    /* wrap around on section length */
 140                 o += h->lim[i];    /* offset to start of section */
 141                 t[o] = h->t + h->dt; /* last point is next grid point */
 142         }
 143
 144         filler(fill_data, n, t, vi, h->buf);
 145
 146         sk_free(t);
 147         sk_free(vi);
 148 }
 149
 150 void sk_hist_get(sk_hist * restrict h, double t, double * restrict aff)
 151 {
 152         int i;
 153
 154         if (h==NULL)
 155                 return;
 156
 157         for (i = 0; i < h->nd; i++)
 158         {
 159                 int ui, i0, i1, i0_, i1_, p, l, o;
 160                 double dt, y0, y1, m, dx;
 161
 162                 ui = h->vi2i[h->vi[i]];
 163                 p = h->pos[ui];
 164                 l = h->len[ui];
 165                 o = h->lim[ui];
 166
 167                 dt = (t - h->del[i]) - (h->t - (l - 2)*h->dt);
 168
 169                 i0 = (int)ceil((l - 2) - dt / h->dt);
 170                 if (i0 < 0)
 171                         i0 += l;
 172                 i1 = i0 ? i0 - 1 : l - 1;
 173
 174                 i0_ = (p + i0) % l + o;
 175                 i1_ = (p + i1) % l + o;
 176
 177 #ifdef SKDEBUG
 178                 if ((i0_ < h->lim[ui]) || (i0_ >= h->lim[ui+1]))
 179                         fprintf(stderr, "[sk_hist_get] oob: i0_=%d not in [%d, %d) at %s:%d\n", i0_, h->lim[ui], h->lim[ui+1], __FILE__, __LINE__);
 180                 if ((i1_ < h->lim[ui]) || (i1_ >= h->lim[ui+1]))
 181                         fprintf(stderr, "[sk_hist_get] oob: i1_=%d not in [%d, %d) at %s:%d\n", i1_, h->lim[ui], h->lim[ui+1], __FILE__, __LINE__);
 182 #endif
 183
 184                 /* bottleneck */
 185                 y0 = h->buf[i0_];
 186                 y1 = h->buf[i1_];
 187
 188                 m = (y1 - y0) / h->dt;
 189
 190                 dx = (t - h->del[i]) - (h->t - i0*h->dt);
 191
 192                 aff[i] = m * dx + y0;
 193                 continue;
 194         }
 195 }
 196
 197 static void update_time(sk_hist *h, double new_t)
 198 {
 199         int i, n_steps;
 200
 201         /* the current time must always be contained between
 202          * pos and pos-1 in the buffer
 203          *
 204          * TODO handle case of irregular time step where new_t could
 205          * be several h->dt ahead, and we need to fill in the buffer
 206          * correctly.
 207          */
 208         n_steps = 0;
 209         while ((h->t + h->dt) <= new_t)
 210         {
 211                 h->t += h->dt;
 212                 n_steps++;
 213         }
 214
 215         /* if no steps to take, return early */
 216         if (n_steps==0)
 217                 return;
 218
 219         /* update positions */
 220         for (i=0; i<h->nu; i++)
 221         {
 222                 h->pos[i] -= n_steps;
 223                 if (h->pos[i] < 0)
 224                         h->pos[i] += h->len[i];
 225         }
 226 }
 227
 228 void sk_hist_set(sk_hist * restrict h, double t, double * restrict eff)
 229 {
 230         int i, i0, i1;
 231         double x0, dx, dt;
 232
 233         if (h==NULL)
 234                 return;
 235
 236         update_time(h, t);
 237
 238         for (i=0; i<h->nu; i++)
 239         {
 240                 i0 = h->pos[i];
 241                 i1 = i0 ? i0 - 1 : h->len[i] - 1;
 242                 i0 += h->lim[i];
 243                 i1 += h->lim[i];
 244 #ifdef SKDEBUG
 245         if ((i0 < h->lim[i]) || (i0 > h->lim[i+1]))
 246             fprintf(stderr, "[sk_hist_set] t=%.3f ui=%d, i0=%d not in [%d,%d) %s:%d\n",
 247                     t, i, i0, h->lim[i], h->lim[i+1], __FILE__, __LINE__);
 248         if ((i1 < h->lim[i]) || (i1 > h->lim[i+1]))
 249             fprintf(stderr, "[sk_hist_set] t=%.3f ui=%d, i1=%d not in [%d,%d) %s:%d\n",
 250                     t, i, i1, h->lim[i], h->lim[i+1], __FILE__, __LINE__);
 251 #endif
 252                 x0 = h->buf[i0];
 253                 dt = t - h->t;
 254
 255 #ifdef SKDEBUG
 256                 if (dt < 0)
 257                         fprintf(stderr, "[sk_hist] unhandled dt<0 at %s:%d\n", __FILE__, __LINE__);
 258 #endif
 259
 260                 if (dt > 0) {
 261                         /* extrapolate from (x(h->t), x(t)) to next grid point*/
 262                         dx = eff[h->vi2i[h->uvi[i]]] - h->buf[i0];
 263                         h->buf[i1] = (dx / dt) * h->dt + x0;
 264                 } else {
 265                         /* reset grid point value */
 266                         h->buf[i0] = eff[h->vi2i[h->uvi[i]]];
 267                 }
 268         }
 269 }
 270
 271 int sk_hist_nbytes(sk_hist *h)
 272 {
 273         int nb;
 274         nb = sizeof(sk_hist);
 275         nb += sizeof(int) * ((h->nu+1) + 3*h->nu + h->nd + h->maxvi);
 276         nb += sizeof(double) * (h->lim[h->nu] + h->nu + h->nd);
 277         return nb;
 278 }
 279
 280 void sk_hist_zero_filler(void * restrict data, int n, double * restrict t, int * restrict indices, double * restrict buf)
 281 {
 282         int i;
 283         /* suppress unused parameter warnings */
 284         (void) data; (void) n; (void) t; (void) indices;
 285         for (i=0; i<n; i++)
 286                 buf[i] = 0.0;
 287 }
 288
 289 int sk_hist_get_nu(sk_hist *h) {
 290         return h->nu;
 291 }
 292
 293 double sk_hist_get_buf_lin(sk_hist *h, int index) {
 294 #ifdef SKDEBUG
 295         if ((index < 0) || (index >= h->lim[h->nu]))
 296                 fprintf(stderr, "[sk_hist_get_buf_lin] oob index=%d not in [0, %d)\n", index, h->lim[h->nu]);
 297 #endif
 298         return h->buf[index];
 299 }
 300
 301 double sk_hist_get_t(sk_hist *h) {
 302         return h->t;
 303 }
 304
 305 double sk_hist_get_dt(sk_hist *h) {
 306         return h->dt;
 307 }
 308
 309 int sk_hist_get_lim(sk_hist *h, int i) {
 310         return h->lim[i];
 311 }
 312
 313 int sk_hist_get_len(sk_hist *h, int i) {
 314         return h->len[i];
 315 }
 316
 317 int sk_hist_get_nd(sk_hist *h) {
 318         return h->nd;
 319 }
 320
 321 int sk_hist_get_pos(sk_hist *h, int i) {
 322         return h->pos[i];
 323 }
 324
 325 int sk_hist_get_uvi(sk_hist *h, int i) {
 326         return h->uvi[i];
 327 }
 328
 329 int sk_hist_get_vi(sk_hist *h, int i) {
 330         return h->vi[i];
 331 }
 332
 333 double sk_hist_get_maxd(sk_hist *h, int i) {
 334         return h->maxd[i];
 335 }
 336
 337 int sk_hist_buf_is_null(sk_hist *h) {
 338         return h->buf == NULL;
 339 }
 340
 341 double sk_hist_get_vd(sk_hist *h, int i) {
 342         return h->del[i];
 343 }