field.c 15.2 KB
Newer Older
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1
2
3
4
/*
  This file is part of CDO. CDO is a collection of Operators to
  manipulate and analyse Climate model Data.

Uwe Schulzweida's avatar
Uwe Schulzweida committed
5
  Copyright (C) 2003-2015 Uwe Schulzweida, <uwe.schulzweida AT mpimet.mpg.de>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
6
7
8
9
10
11
12
13
14
15
16
17
18
19
  See COPYING file for copying and redistribution conditions.

  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; version 2 of the License.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
*/

#include "cdo.h"
#include "cdo_int.h"
Ralf Mueller's avatar
Ralf Mueller committed
20
#include <cdi.h>
21
#include "percentiles.h"
22
#include "merge_sort2.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
23

Uwe Schulzweida's avatar
Uwe Schulzweida committed
24
double crps_det_integrate(double *a, const double d, const size_t n);
25

Uwe Schulzweida's avatar
Uwe Schulzweida committed
26
double _FADD_(const double x, const double y, const double missval1, const double missval2)
27
28
29
30
{
  return FADD(x,y);
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
31
double _FSUB_(const double x, const double y, const double missval1, const double missval2)
32
33
34
35
{
  return FSUB(x, y);
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
36
double _FMUL_(const double x, const double y, const double missval1, const double missval2)
37
38
39
40
{
  return FMUL(x, y);
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
41
double _FDIV_(const double x, const double y, const double missval1, const double missval2)
42
43
44
45
{
  return FDIV(x, y);
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
46
double _FPOW_(const double x, const double y, const double missval1, const double missval2)
47
{
48
49
50
  return FPOW(x, y);
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
51
double _FSQRT_(const double x, const double missval1)
52
53
{
  return FSQRT(x);
54
55
56
}


Uwe Schulzweida's avatar
Uwe Schulzweida committed
57
double fldfun(field_t field, const int function)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
58
59
60
{
  double rval = 0;

61
  // Why is this not a switch-case statement?
Uwe Schulzweida's avatar
Uwe Schulzweida committed
62
  if      ( function == func_min  )  rval = fldmin(field);
63
64
65
66
67
  else if ( function == func_max  )  rval = fldmax(field);
  else if ( function == func_sum  )  rval = fldsum(field);
  else if ( function == func_mean )  rval = fldmean(field);
  else if ( function == func_avg  )  rval = fldavg(field);
  else if ( function == func_std  )  rval = fldstd(field);
68
  else if ( function == func_std1 )  rval = fldstd1(field);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
69
  else if ( function == func_var  )  rval = fldvar(field);
70
  else if ( function == func_var1 )  rval = fldvar1(field);
71
  
72
  else if ( function == func_crps )  rval = fldcrps(field);
73
  else if ( function == func_brs )   rval = fldbrs(field);
74
75
76

  else if ( function == func_rank )  rval = fldrank(field);
  else if ( function == func_roc )   rval = fldroc(field);
77
  else cdoAbort("%s: function %d not implemented!", __func__, function);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
78
79
80
81

  return rval;
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
82

83
84
double fldrank(field_t field) 
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
85
  double res = 0;
86
87
  // Using first value as reference (observation)
  double *array  =  &(field.ptr[1]);
88
  double val     = array[-1];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
89
  const double missval = field.missval;
90
  int nmiss      = field.nmiss;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
91
92
  const size_t len       = field.size-1;
  size_t j;
93
  
Uwe Schulzweida's avatar
Uwe Schulzweida committed
94
  if ( nmiss )  return(missval);
95

Uwe Schulzweida's avatar
Uwe Schulzweida committed
96
  sort_iter_single(len,array, 1);
97

98
99
  if ( val > array[len-1] ) 
    res=(double)len;
100
101
  else 
    for ( j=0; j<len; j++ )
102
103
104
105
      if ( array[j] >= val ) {
	res=(double)j; 
	break;
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
106

107
  return res;
108
109
110
111
112
113
114
115
}


double fldroc(field_t field) 
{
  return field.missval;
}

116
117
double fldcrps(field_t field)
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
118
119
  const size_t len     = field.size;
  const int    nmiss   = field.nmiss;
120
121
122
123
  double *array  = field.ptr;

  if ( nmiss > 0 ) 
    cdoAbort("Missing values not implemented in crps calculation");
124
125
126
127
128
  // possible handling of missing values:
  // (1) strip them off, and sort array without missing values
  //     using only (len - 1 - nmiss) values
  // (2) modify merge_sort in a way, that missing values will
  //     always go to the end of the list
129
130
131

  // Use first value as reference
  sort_iter_single(len-1,&array[1],ompNumThreads);
132

Uwe Schulzweida's avatar
Uwe Schulzweida committed
133
  return crps_det_integrate(&array[1],array[0],len-1);
134
135
}

136

137
138
double fldbrs(field_t field) 
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
139
140
  const size_t    len   = field.size;
  const int     nmiss   = field.nmiss;
141
  double *array   = field.ptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
142
  const double missval  = field.missval;
143

144
  double brs = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
145
  size_t i, count=0;
146

147
148
149
150
151
152
153
154
155
  // Using first value as reference
  if ( nmiss == 0 ) 
    {
      for ( i=1; i<len; i++ )
	brs += (array[i] - array[0]) * (array[i] - array[0]);
      count = i-1;
    }
  else 
    {
156
157
      if ( DBL_IS_EQUAL(array[0], missval) ) return missval;

158
      for ( i=1; i<len; i++ )
159
	if ( !DBL_IS_EQUAL(array[i], missval) ) 
160
161
162
163
164
	  {
	    brs += (array[i] - array[0]) * (array[i] - array[0]);
	    count ++;
	  }
    }
165

166
  return brs/count;
167
168
}

169

Uwe Schulzweida's avatar
Uwe Schulzweida committed
170
double fldmin(field_t field)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
171
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
172
173
174
175
  size_t   i;
  const size_t len     = field.size;
  const int    nmiss   = field.nmiss;
  const double missval = field.missval;
176
  const double *restrict array  = field.ptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
177
178
179
180
181
182
183
184
185
  double rmin = 0;

  if ( nmiss > 0 )
    {
      rmin = DBL_MAX;
      for ( i = 0; i < len; i++ ) 
	if ( !DBL_IS_EQUAL(array[i], missval) )
	  if ( array[i] < rmin ) rmin = array[i];

186
      if ( IS_EQUAL(rmin, DBL_MAX) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
187
188
189
190
191
	rmin = missval;
    }
  else
    {
      rmin = array[0];
192
      //#pragma simd reduction(min:rmin) 
Uwe Schulzweida's avatar
Uwe Schulzweida committed
193
194
195
196
197
198
199
200
      for ( i = 1; i < len; i++ ) 
	if ( array[i] < rmin )  rmin = array[i];
    }

  return (rmin);
}


Uwe Schulzweida's avatar
Uwe Schulzweida committed
201
double fldmax(field_t field)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
202
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
203
204
205
206
  size_t   i;
  const size_t len     = field.size;
  const int    nmiss   = field.nmiss;
  const double missval = field.missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
207
208
209
210
211
  double *array  = field.ptr;
  double rmax = 0;

  if ( nmiss > 0 )
    {
212
213
      rmax = -DBL_MAX;
      for ( i = 0; i < len; i++ )
214
215
216
        if ( !DBL_IS_EQUAL(array[i], missval) )
          if ( array[i] > rmax ) rmax = array[i];
      
217
      if ( IS_EQUAL(rmax, -DBL_MAX) )
218
        rmax = missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
219
220
221
222
223
    }
  else
    {
      rmax = array[0];
      for ( i = 1; i < len; i++ ) 
224
        if ( array[i] > rmax )  rmax = array[i];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
225
226
227
228
229
230
    }

  return (rmax);
}


Uwe Schulzweida's avatar
Uwe Schulzweida committed
231
double fldsum(field_t field)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
232
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
233
234
235
236
237
  size_t   i;
  size_t   nvals   = 0;
  const size_t len     = field.size;
  const int    nmiss   = field.nmiss;
  const double missval = field.missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
238
239
240
241
242
243
244
  double *array  = field.ptr;
  double rsum = 0;

  if ( nmiss > 0 )
    {
      for ( i = 0; i < len; i++ ) 
	if ( !DBL_IS_EQUAL(array[i], missval) )
245
246
247
248
249
250
	  {
	    rsum += array[i];
	    nvals++;
	  }

      if ( !nvals ) rsum = missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
251
252
253
254
255
256
257
258
259
260
261
    }
  else
    {
      for ( i = 0; i < len; i++ ) 
	rsum += array[i];
    }

  return (rsum);
}


Uwe Schulzweida's avatar
Uwe Schulzweida committed
262
double fldmean(field_t field)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
263
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
264
265
266
267
268
  size_t   i;
  const size_t len      = field.size;
  const int    nmiss    = field.nmiss;
  const double missval1 = field.missval;
  const double missval2 = field.missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
  double *array   = field.ptr;
  double *w       = field.weight;
  double rsum = 0, rsumw = 0, ravg = 0;

  if ( nmiss > 0 )
    {
      for ( i = 0; i < len; i++ ) 
	if ( !DBL_IS_EQUAL(array[i], missval1) && !DBL_IS_EQUAL(w[i], missval1) )
	  {
	    rsum  += w[i] * array[i];
	    rsumw += w[i];
	  }
    }
  else
    {
      for ( i = 0; i < len; i++ ) 
	{
	  rsum  += w[i] * array[i];
	  rsumw += w[i];
	}
    }

  ravg = DIV(rsum, rsumw);

  return (ravg);
}


Uwe Schulzweida's avatar
Uwe Schulzweida committed
297
double fldavg(field_t field)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
298
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
299
300
301
302
303
  size_t   i;
  const size_t len      = field.size;
  const int    nmiss    = field.nmiss;
  const double missval1 = field.missval;
  const double missval2 = field.missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
  double *array   = field.ptr;
  double *w       = field.weight;
  double rsum = 0, rsumw = 0, ravg = 0;

  if ( nmiss > 0 )
    {
      for ( i = 0; i < len; i++ ) 
	if ( !DBL_IS_EQUAL(w[i], missval1) )
	  {
	    rsum  = ADD(rsum, MUL(w[i], array[i]));
	    rsumw = ADD(rsumw, w[i]);
	  }
    }
  else
    {
      for ( i = 0; i < len; i++ ) 
	{
	  rsum  += w[i] * array[i];
	  rsumw += w[i];
	}
    }

  ravg = DIV(rsum, rsumw);

  return (ravg);
}

331
static
332
void prevarsum(const double *restrict array, const double *restrict w, size_t len, int nmiss, 
333
	       double missval, double *restrict rsum, double *restrict rsumw, double *restrict rsumq, double *restrict rsumwq)
334
{ 
Uwe Schulzweida's avatar
Uwe Schulzweida committed
335
  size_t i;
336
337
  *rsum = *rsumw = 0;
  *rsumq = *rsumwq = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
338
339
340
341

  if ( nmiss > 0 )
    {
      for ( i = 0; i < len; i++ ) 
342
343
        if ( !DBL_IS_EQUAL(array[i], missval) && !DBL_IS_EQUAL(w[i], missval) )
          {
344
345
346
347
            *rsum   += w[i] * array[i];
            *rsumq  += w[i] * array[i] * array[i];
            *rsumw  += w[i];
            *rsumwq += w[i] * w[i];
348
          }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
349
350
351
352
    }
  else
    {
      for ( i = 0; i < len; i++ ) 
353
        {
354
355
356
357
          *rsum   += w[i] * array[i];
          *rsumq  += w[i] * array[i] * array[i];
          *rsumw  += w[i];
          *rsumwq += w[i] * w[i];
358
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
359
    }
360
361
362
363
364
}


double fldvar(field_t field)
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
365
366
367
  const size_t len     = field.size;
  const int    nmiss   = field.nmiss;
  const double missval = field.missval;
368
369
  double *array  = field.ptr;
  double *w      = field.weight;
370
  double rsum, rsumw;
371
372
373
  double rsumq, rsumwq;

  prevarsum(array, w, len, nmiss, missval, &rsum, &rsumw, &rsumq, &rsumwq);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
374

375
  double rvar = IS_NOT_EQUAL(rsumw, 0) ? (rsumq*rsumw - rsum*rsum) / (rsumw*rsumw) : missval;
376
  if ( rvar < 0 && rvar > -1.e-5 ) rvar = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
377
378
379
380
381

  return (rvar);
}


382
383
double fldvar1(field_t field)
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
384
385
386
  const size_t len     = field.size;
  const int    nmiss   = field.nmiss;
  const double missval = field.missval;
387
388
  double *array  = field.ptr;
  double *w      = field.weight;
389
  double rsum, rsumw;
390
391
392
393
  double rsumq, rsumwq;

  prevarsum(array, w, len, nmiss, missval, &rsum, &rsumw, &rsumq, &rsumwq);

394
  double rvar = (rsumw*rsumw > rsumwq) ? (rsumq*rsumw - rsum*rsum) / (rsumw*rsumw - rsumwq) : missval;
395
396
397
398
399
400
  if ( rvar < 0 && rvar > -1.e-5 ) rvar = 0;

  return (rvar);
}


401
double var_to_std(double rvar, double missval)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
402
{
403
  double rstd;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
404

405
406
407
408
409
410
411
412
  if ( DBL_IS_EQUAL(rvar, missval) || rvar < 0 )
    {
      rstd = missval;
    }
  else
    {
      rstd = IS_NOT_EQUAL(rvar, 0) ? sqrt(rvar) : 0;
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
413

414
  return rstd;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
415
416
}

417
double fldstd(field_t field)
418
{
419
420
  return var_to_std(fldvar(field), field.missval);
}
421
422


423
424
425
double fldstd1(field_t field)
{
  return var_to_std(fldvar1(field), field.missval);
426
427
428
}


Uwe Schulzweida's avatar
Uwe Schulzweida committed
429
void fldrms(field_t field, field_t field2, field_t *field3)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
430
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
431
432
  size_t   i;
  size_t len;
433
  int    rnmiss = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
434
  int    grid1    = field.grid;
435
  //  int    nmiss1   = field.nmiss;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
436
437
  double *array1  = field.ptr;
  int    grid2    = field2.grid;
438
  //  int    nmiss2   = field2.nmiss;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
439
  double *array2  = field2.ptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
440
441
  const double missval1 = field.missval;
  const double missval2 = field2.missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
442
443
444
445
  double *w       = field.weight;
  double rsum = 0, rsumw = 0, ravg = 0;

  len    = gridInqSize(grid1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
446
  if ( len != (size_t) gridInqSize(grid2) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
    cdoAbort("fields have different size!");

  /*
  if ( nmiss1 > 0 )
  */
    {
      for ( i = 0; i < len; i++ ) 
	if ( !DBL_IS_EQUAL(w[i], missval1) )
	  {
	    rsum  = ADD(rsum, MUL(w[i], MUL(SUB(array2[i], array1[i]),
                                            SUB(array2[i], array1[i]))));
	    rsumw = ADD(rsumw, w[i]);
	  }
    }
    /*
  else
    {
      for ( i = 0; i < len; i++ ) 
	{
	  rsum  += w[i] * array1[i];
	  rsumw += w[i];
	}
    }
    */

472
  ravg = SQRT(DIV(rsum, rsumw));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
473
474
475
476
477
478
479
480

  if ( DBL_IS_EQUAL(ravg, missval1) ) rnmiss++;

  field3->ptr[0] = ravg;
  field3->nmiss  = rnmiss;
}


Uwe Schulzweida's avatar
Uwe Schulzweida committed
481
void varrms(field_t field, field_t field2, field_t *field3)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
482
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
483
  size_t   i, k, nlev, len;
484
  int    rnmiss = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
485
486
  int    zaxis    = field.zaxis;
  int    grid1    = field.grid;
487
  //  int    nmiss1   = field.nmiss;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
488
489
  double *array1  = field.ptr;
  int    grid2    = field2.grid;
490
  //  int    nmiss2   = field2.nmiss;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
491
  double *array2  = field2.ptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
492
493
  const double missval1 = field.missval;
  const double missval2 = field2.missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
494
495
496
497
498
  double *w       = field.weight;
  double rsum = 0, rsumw = 0, ravg = 0;

  nlev   = zaxisInqSize(zaxis);
  len    = gridInqSize(grid1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
499
  if ( len != (size_t) gridInqSize(grid2) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
    cdoAbort("fields have different size!");

  /*
  if ( nmiss1 > 0 )
  */
    {
      for ( k = 0; k < nlev; k++ )
	for ( i = 0; i < len; i++ )
	  /*	  if ( !DBL_IS_EQUAL(w[i], missval1) ) */
	    {
	      rsum  = ADD(rsum, MUL(w[i], MUL(SUB(array2[k*len+i], array1[k*len+i]),
                                              SUB(array2[k*len+i], array1[k*len+i]))));
	      rsumw = ADD(rsumw, w[i]);
	    }
    }
    /*
  else
    {
      for ( i = 0; i < len; i++ )
	{
	  rsum  += w[i] * array1[i];
	  rsumw += w[i];
	}
    }
    */

526
  ravg = SQRT(DIV(rsum, rsumw));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
527
528
529
530
531
532

  if ( DBL_IS_EQUAL(ravg, missval1) ) rnmiss++;

  field3->ptr[0] = ravg;
  field3->nmiss  = rnmiss;
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
533
534

/* RQ */
535
double fldpctl(field_t field, const double pn)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
536
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
537
538
539
  const size_t len     = field.size;
  const int    nmiss   = field.nmiss;
  const double missval = field.missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
540
541
  double *array  = field.ptr;
  double pctl = missval;
542

Uwe Schulzweida's avatar
Uwe Schulzweida committed
543
544
545
546
  if ( len - nmiss > 0 )
    {
      if ( nmiss > 0 )
        {
547
          double *array2 = (double*) Malloc((len - nmiss)*sizeof(double));
548

549
550
          size_t j = 0;
          for ( size_t i = 0; i < len; i++ ) 
Uwe Schulzweida's avatar
Uwe Schulzweida committed
551
552
553
            if ( !DBL_IS_EQUAL(array[i], missval) )
              array2[j++] = array[i];

554
          pctl = percentile(array2, j, pn);
555

556
          Free(array2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
557
558
559
        }
      else
        {
560
          pctl = percentile(array, len, pn);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
561
562
563
564
565
566
        }
    }

  return pctl;
}
/* QR */
567

Ralf Mueller's avatar
Ralf Mueller committed
568
/*  field_t UTILITIES */
569
570
571
/*  update the number non missing values */
void fldunm(field_t *field)
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
572
  size_t i;
573
574
575
576
577

  field->nmiss = 0;
  for ( i = 0; i < field->size; i++ )
    if ( DBL_IS_EQUAL(field->ptr[i], field->missval) ) field->nmiss++;
}
Ralf Mueller's avatar
Ralf Mueller committed
578
579

/*  check for non missval values */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
580
int fldhvs(field_t *fieldPtr, const size_t nlevels)
Ralf Mueller's avatar
Ralf Mueller committed
581
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
582
  size_t level;
Ralf Mueller's avatar
Ralf Mueller committed
583
584
585
586
587
  field_t field;

  for ( level = 0; level < nlevels; level++)
    {
      field = fieldPtr[level];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
588
      if ( (size_t)field.nmiss != field.size )
Ralf Mueller's avatar
Ralf Mueller committed
589
590
        return TRUE;
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
591

Ralf Mueller's avatar
Ralf Mueller committed
592
593
  return FALSE;
}
594
595


Uwe Schulzweida's avatar
Uwe Schulzweida committed
596
double crps_det_integrate(double *a, const double d, const size_t n)
597
598
599
600
601
602
603
604
605
606
607
608
609
610
{
  /* *************************************************************************** */
  /* This routine finds the area between the cdf described by the ordered array  */
  /* of doubles (double *a) and the Heavyside function H(d)                      */
  /* INPUT ARGUMENTS:                                                            */
  /*     double *a  - ordered array of doubles describing a cdf                  */
  /*                  as cdf(a[i]) = ( (double)i )/ n                            */
  /*     double d   - describing a reference value                               */
  /*     int n      - the length of array a                                      */
  /* RETURN VALUE:                                                               */
  /*     double     - area under the curve in units of a                         */
  /* *************************************************************************** */

  double area = 0; 
611
  //  double tmp;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
612
  size_t i;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
613
#if defined(_OPENMP)
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
#pragma omp parallel for if ( n>10000 ) shared(a) private(i) \
  reduction(+:area) schedule(static,10000) 
#endif                                                     /* **************************** */
  for ( i=1; i<n; i++ ) {                                  /* INTEGRATE CURVE AREA         */
    if ( a[i] < d )                                        /* left of heavyside            */
      area += (a[i]-a[i-1])*(double)i*i/n/n;                   /*                              */
    else if ( a[i-1] > d )                                 /* right of heavyside           */
      area += (a[i]-a[i-1])*(1.-(double)i/n)*(1.-(double)i/n);              /*                              */
    else if ( a[i-1] < d && a[i] > d ) {                   /* hitting jump pf heavyside    */
      area += (d-a[i-1]) * (double)i*i/n/n;                    /* (occurs exactly once!)       */
      area += (a[i]-d) * (1.-(double)i/n)*(1.-(double)i/n);                 /* **************************** */
    }
  }


  return(area);
}