Ensstat3.cc 15.8 KB
Newer Older
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-2017 Uwe Schulzweida, <uwe.schulzweida AT mpimet.mpg.de>
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
  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.
*/

/*
   This module contains the following operators:
   Ensstat3       ensrkhist_space  Ensemble ranked histogram averaged over time
   Ensstat3       ensrkhist_time   Ensemble ranked histogram averaged over space
   Ensstat3       ensroccurve      Ensamble Receiver Operating Characteristics
*/


#include <cdi.h>
#include "cdo.h"
#include "cdo_int.h"
#include "pstream.h"
#include "util.h"

Cedrick Ansorge's avatar
Cedrick Ansorge committed
32
// Defines for rank histogram
33
34
35
36
enum TDATA_TYPE {TIME, SPACE};
#define time_data TIME
#define space_data SPACE

Cedrick Ansorge's avatar
Cedrick Ansorge committed
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51

// Defines for Receiver Operating Characteristics (ROC)
#define DEBUG_ROC 0
enum CONTINGENCY_TYPE {TP, FP, FN, TN};
/* TP - True positive  ( event     forecast and     occured)  HIT               */
/* FP - False positive ( event     forecast and not occured)  FALSE ALARM       */
/* TN - True negative  ( event not forecast and not ocurred)  CORRECT REJECTION */
/* FN - False negative ( event not forecast and     ocurred)  MISSED            */

enum ROC_ENUM_TYPE {TPR, FPR};
/* TPR = True Positive Rate = TP / ( TP + FN ) */
/* FNR = False Negtive Rate = FN / ( FP + TN ) */

double roc_curve_integrate(const double **roc, const int n);

52
53
void *Ensstat3(void *argument)
{
Cedrick Ansorge's avatar
Cedrick Ansorge committed
54
  int i,j;
55
  int nrecs = 0, nrecs0, nmiss;
Cedrick Ansorge's avatar
Cedrick Ansorge committed
56
57
  int cum;
  int chksum;                  // for check of histogram population 
58
59
  int levelID, varID, binID = 0;
  int vlistID;
60
  int gridID, gridID2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
61
62
  int have_miss = 0;
  int streamID = 0, streamID2 = 0;
Cedrick Ansorge's avatar
Cedrick Ansorge committed
63
  int **array2 = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
64
  int **ctg_tab = NULL, *hist = NULL;         // contingency table and histogram
65
  double missval;
Cedrick Ansorge's avatar
Cedrick Ansorge committed
66
  double *dat;                  // pointer to ensemble data for ROC
Uwe Schulzweida's avatar
Uwe Schulzweida committed
67
68
  double *uThresh = NULL, *lThresh = NULL;    // thresholds for histograms
  double **roc = NULL;                 // receiver operating characteristics table
69
  double val;
70
  int ival;
71
72
73
74
75
76
77
78
79
80

  typedef struct
  {
    int streamID;
    int vlistID;
    double *array;
  } ens_file_t;

  cdoInitialize(argument);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
81
  // clang-format off
82
83
84
  cdoOperatorAdd("ensroc",          func_roc,  0,          NULL);
  cdoOperatorAdd("ensrkhist_space", func_rank, space_data, NULL);
  cdoOperatorAdd("ensrkhist_time",  func_rank, time_data,  NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
85
  // clang-format on
86
  
87
88
89
  int operatorID = cdoOperatorID();
  int operfunc = cdoOperatorF1(operatorID);
  int datafunc = cdoOperatorF2(operatorID);
Cedrick Ansorge's avatar
Cedrick Ansorge committed
90

91
  int nbins = 0;
Cedrick Ansorge's avatar
Cedrick Ansorge committed
92
93
  if ( operfunc == func_roc ) {
    operatorInputArg("Number of eigen functions to write out");
94
    nbins = parameter2int(operatorArgv()[0]);
Cedrick Ansorge's avatar
Cedrick Ansorge committed
95
  }
96
  
97
98
  int nfiles = cdoStreamCnt() - 1;
  int nens = nfiles-1;
99
100
101
102

  if ( cdoVerbose )
    cdoPrint("Ensemble over %d files.", nfiles);

103
  const char *ofilename = cdoStreamName(nfiles)->args;
104

105
106
  if ( !cdoOverwriteMode && fileExists(ofilename) && !userFileOverwrite(ofilename) )
    cdoAbort("Outputfile %s already exists!", ofilename);
107

108
  ens_file_t *ef = (ens_file_t*) Malloc(nfiles*sizeof(ens_file_t));
109
110
111
112
113
114

  /* *************************************************** */
  /* should each thread be allocating memory locally???? */
  /* ("first touch strategy")                            */
  /* --> #pragma omp parallel for ...                    */
  /* *************************************************** */
115
  field_type *field = (field_type*) Malloc(ompNumThreads*sizeof(field_type));
116
  for ( i = 0; i < ompNumThreads; i++ )
117
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
118
      field_init(&field[i]);
119
      field[i].size   = nfiles;
120
121
      field[i].ptr    = (double*) Malloc(nfiles*sizeof(double));
      field[i].weight = (double*) Malloc(nfiles*sizeof(double));
122
      for ( int fileID = 0; fileID < nfiles; fileID++ )
123
124
125
	field[i].weight[fileID] = 1;
    }

126
  for ( int fileID = 0; fileID < nfiles; fileID++ )
127
    {
128
      streamID = pstreamOpenRead(cdoStreamName(fileID));
129

130
      vlistID = pstreamInqVlist(streamID);
131
132
133
134
135

      ef[fileID].streamID = streamID;
      ef[fileID].vlistID = vlistID;
    }

Cedrick Ansorge's avatar
Cedrick Ansorge committed
136
  /* check for identical contents of all ensemble members */
137
  for ( int fileID = 1; fileID < nfiles; fileID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
138
    vlistCompare(ef[0].vlistID, ef[fileID].vlistID, CMP_ALL);
139

140
141
142
143
  int vlistID1 = ef[0].vlistID;
  int vlistID2 = vlistCreate();
  int nvars = vlistNvars(vlistID1);
  int *varID2 = (int*) Malloc(nvars*sizeof(int));
144

145
146
  double *levs = (double*) Calloc(nfiles, sizeof(double));
  int zaxisID2 = zaxisCreate(ZAXIS_GENERIC, nfiles);
147
148
149
150
  for ( i=0; i<nfiles; i++ )
    levs[i] = i;
  zaxisDefLevels(zaxisID2,levs);
  zaxisDefName(zaxisID2, "histogram_binID");
Cedrick Ansorge's avatar
Cedrick Ansorge committed
151

152
  int time_mode = datafunc == TIME? TSTEP_INSTANT : TSTEP_CONSTANT;
Cedrick Ansorge's avatar
Cedrick Ansorge committed
153

154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
  for ( varID=0; varID<nvars; varID++) {

    /* **************************************************************** */
    /* nfiles includes the observation, so there are nfiles-1 ensembles */
    /* and exactly nfiles bins, in which the observation could fall     */
    /* **************************************************************** */

    if ( datafunc == TIME ) 
      {
	val = 0;
	gridID2 = gridCreate(GRID_LONLAT, 1);
	gridDefXsize(gridID2, 1);
	gridDefYsize(gridID2, 1);
	gridDefXvals(gridID2, &val);
	gridDefYvals(gridID2, &val);
      }
    else // datafunc == SPACE
      gridID2 = vlistInqVarGrid(vlistID1,varID);

    varID2[varID] = vlistDefVar(vlistID2, gridID2, zaxisID2, time_mode);
  }

176
177
  int taxisID1 = vlistInqTaxis(vlistID1);
  int taxisID2 = taxisDuplicate(taxisID1);
178
179
180
181
182
183
184
185
186
187
188
  vlistDefTaxis(vlistID2, taxisID2);
  
  for ( varID=0; varID< nvars; varID++ ){
    if ( zaxisInqSize(vlistInqVarZaxis(vlistID1, varID)) > 1 ) {
      cdoWarning("More than one level not supported when processing ranked histograms.");
      cdoWarning("Try to use `cdo splitlevel` to split the dataset into levels and apply");
      cdoWarning("the operator seperately to each level.");
      cdoAbort("Exit due to unsupported file structure");
    }
  } 

189
190
  if ( operfunc != func_roc )
    {
191
192
      streamID2 = pstreamOpenWrite(cdoStreamName(nfiles), cdoFiletype());
      pstreamDefVlist(streamID2, vlistID2);
193
    }
194

195
  int gridsize = vlistGridsizeMax(vlistID1);
196

197
  for ( int fileID = 0; fileID < nfiles; fileID++ )
198
    ef[fileID].array = (double*) Malloc(gridsize*sizeof(double));
199

Cedrick Ansorge's avatar
Cedrick Ansorge committed
200
  if ( operfunc == func_rank && datafunc == SPACE ) 
201
    { /*need to memorize data for entire grid before writing          */
202
      array2 = (int**) Malloc((nfiles+1)*sizeof(int*));
203
      for ( binID=0; binID<nfiles; binID++ ) 
204
	array2[binID] = (int*) Calloc( gridsize, sizeof(int));
205
    }
Cedrick Ansorge's avatar
Cedrick Ansorge committed
206
  else if ( operfunc == func_rank )
207
208
    {  /* can process data separately for each timestep and only need */
       /* to cumulate values over the grid                            */
209
      array2    = (int **) Malloc( (nfiles+1)*sizeof(int *));
210
      for ( binID=0; binID<nfiles; binID++ )
211
	array2[binID] = (int*) Calloc( 1, sizeof(int));
212
    }
Cedrick Ansorge's avatar
Cedrick Ansorge committed
213
214

  if ( operfunc == func_roc ) {
215
    ctg_tab = (int**) Malloc((nbins+1)*sizeof(int*));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
216
217
218
    hist =    (int*) Malloc(nbins*sizeof(int));
    uThresh = (double*) Malloc(nbins*sizeof(double));
    lThresh = (double*) Malloc(nbins*sizeof(double));
219
    roc     = (double**) Malloc((nbins+1)*sizeof(double*));
Cedrick Ansorge's avatar
Cedrick Ansorge committed
220
221
    
    for  ( i=0; i<nbins; i++ ) {
222
223
      ctg_tab[i] = (int*) Calloc( 4,sizeof(int));
      roc[i]     = (double*) Calloc( 2,sizeof(double));
Cedrick Ansorge's avatar
Cedrick Ansorge committed
224
225
226
      uThresh[i] = ((double)i+1)/nbins;
      lThresh[i] = (double)i/nbins;
    }
227
228
    ctg_tab[nbins] = (int*) Calloc(4,sizeof(int));
    roc[nbins]     = (double*) Calloc(2,sizeof(double));
Cedrick Ansorge's avatar
Cedrick Ansorge committed
229
  }
230
231
  
  
232
  int tsID = 0;
233
234
  do
    {
235
      nrecs0 = pstreamInqTimestep(ef[0].streamID, tsID);
236
      for ( int fileID = 1; fileID < nfiles; fileID++ )
237
238
	{
	  streamID = ef[fileID].streamID;
239
	  nrecs = pstreamInqTimestep(streamID, tsID);
240
	  if ( nrecs != nrecs0 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
241
242
	    {
	      if ( nrecs == 0 )
243
		cdoAbort("Inconsistent ensemble file, too few time steps in %s!", cdoStreamName(fileID)->args);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
244
245
	      else
		cdoAbort("Inconsistent ensemble file, number of records at time step %d of %s and %s differ!",
246
			   tsID+1, cdoStreamName(0)->args, cdoStreamName(fileID)->args);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
247
	    }
248
249
	}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
250
251
252
      if ( operfunc == func_rank && ( datafunc == TIME || tsID == 0 ) )
	{
	  taxisCopyTimestep(taxisID2, taxisID1);
253
	  if ( nrecs0 > 0 ) pstreamDefTimestep(streamID2, tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
254
	}
255

Cedrick Ansorge's avatar
Cedrick Ansorge committed
256
      //      fprintf(stderr,"TIMESTEP %i varID %i rec %i\n",tsID,varID,recID);
257
      
258
      for ( int recID = 0; recID < nrecs0; recID++ )
259
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
260
#if defined(_OPENMP)
261
262
263
#pragma omp parallel for default(none) shared(ef, nfiles)      \
                                      private(streamID, nmiss) \
                                  lastprivate(varID, levelID)
264
#endif
265
	  for ( int fileID = 0; fileID < nfiles; fileID++ )
266
267
	    {
	      streamID = ef[fileID].streamID;
268
269
	      pstreamInqRecord(streamID, &varID, &levelID);
	      pstreamReadRecord(streamID, ef[fileID].array, &nmiss);
270
271
272
273
274
275
276
	    }

	  gridID   = vlistInqVarGrid(vlistID1, varID);
	  gridsize = gridInqSize(gridID);
	  missval  = vlistInqVarMissval(vlistID1, varID);

	  nmiss = 0;
Cedrick Ansorge's avatar
Cedrick Ansorge committed
277
	  if ( datafunc == TIME && operfunc == func_rank) 
278
279
	    for ( binID=0;binID<nfiles;binID++ )
	      array2[binID][0] = 0;
280

Uwe Schulzweida's avatar
Uwe Schulzweida committed
281
#if defined(_OPENMP)
282
#pragma omp parallel for default(shared) private(binID)
283
284
285
#endif
	  for ( i = 0; i < gridsize; i++ )
	    {
286
287
	      int ompthID = cdo_omp_get_thread_num();

288
289
290
	      field[ompthID].missval = missval;
	      field[ompthID].nmiss = 0;
	      have_miss = 0;
291
	      for ( int fileID = 0; fileID < nfiles; fileID++ )
292
293
294
295
296
297
298
299
300
301
302
		{
		  field[ompthID].ptr[fileID] = ef[fileID].array[i];
		  if ( DBL_IS_EQUAL(field[ompthID].ptr[fileID], missval) ) 
		    {
		      have_miss = 1;
		      break;
		    }
		}
	      
	      // need to ignore all data for a gridpoint if a single ensemble
	      // has a missing value at that gridpoint.
Cedrick Ansorge's avatar
Cedrick Ansorge committed
303
304
305
306
307
308
309
310
	      if ( ! have_miss )  // only process if no missing value in ensemble
		{
		  switch( operfunc ) 
		    {
		    case ( func_rank ): 
		      /* ****************/
		      /* RANK HISTOGRAM */
		      /* ************** */
311
312
		      // for ( j=0; j<nfiles; j++ )
		      //   fprintf(stderr,"%5.2g ",field[ompthID].ptr[j]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
313
#if defined(_OPENMP)
314
315
#pragma omp critical
#endif
316
317
		      binID = (int) fldfun(field[ompthID], operfunc);
		      // fprintf(stderr,"-->%i\n",binID);
Cedrick Ansorge's avatar
Cedrick Ansorge committed
318
319
		      
		      if ( datafunc == SPACE && ! have_miss) 
320
			array2[binID][i]++;
Cedrick Ansorge's avatar
Cedrick Ansorge committed
321
		      else if ( ! have_miss ) 
322
			array2[binID][0]++;
Cedrick Ansorge's avatar
Cedrick Ansorge committed
323
		      break;
324

Cedrick Ansorge's avatar
Cedrick Ansorge committed
325
326
327
328
329
		    case ( func_roc ):
		      /* ********************************** */
		      /* RECEIVER OPERATING CHARACTERISTICS */
		      /* ********************************** */
		      dat = &field[ompthID].ptr[1];
330
		      ival = field[ompthID].ptr[0] > 0.5? 1 : 0;
Cedrick Ansorge's avatar
Cedrick Ansorge committed
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346

		      for ( binID=0; binID<nbins; binID++ )
			hist[binID] = 0;

		      for ( j=0; j<nens; j++ ) 
			for ( binID=0; binID<nbins; binID++ ) 
			  if ( dat[j] >= lThresh[binID] && dat[j] < uThresh[binID] )
			    hist[binID]++;

		      chksum = 0;
		      for ( binID=0; binID<nbins; binID++ )
			chksum += hist[binID];

		      if ( chksum != nens )  exit(1);

		      cum = 0;
347
		      if ( ival == 1 ) 
Cedrick Ansorge's avatar
Cedrick Ansorge committed
348
349
350
351
352
353
354
355
356
357
358
359
360
361
			{
			  // all true positives in first bin
			  ctg_tab[0][TP] += nens;
			  
			  cum += hist[0];
			  for ( binID=1; binID<nbins; binID++ ) 
			    {
			      ctg_tab[binID][TP] += nens-cum;
			      ctg_tab[binID][FN] += cum;
			      cum += hist[binID];
			    }
			  ctg_tab[binID][TP] += nens-cum;
			  ctg_tab[binID][FN] += cum;
			}
362
		      else if ( ival == 0 ) 
Cedrick Ansorge's avatar
Cedrick Ansorge committed
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
			{
			  // all false positives in first bin
			  ctg_tab[0][FP] += nens;
			  cum += hist[0];
			  for ( binID=1; binID<nbins; binID++ ) 
			    {
			      ctg_tab[binID][FP] += nens-cum;
			      ctg_tab[binID][TN] += cum;
			      cum += hist[binID];
			    }
			  ctg_tab[binID][FP] += nens-cum;
			  ctg_tab[binID][TN] += cum;
			}
		      break;
		      
		    }// switch ( operfunc )
		}    // if ( ! have_miss ) 
	    }        // for ( i=0; i<gridsize; i++ )
381

Cedrick Ansorge's avatar
Cedrick Ansorge committed
382
383
384
385
386
	  if ( datafunc == TIME && operfunc == func_rank ) 
	    {
	      for ( binID=0; binID<nfiles; binID++ ) {
		val = (double)array2[binID][0];
		//		fprintf(stderr,"%i ",(int)val);
387
388
		pstreamDefRecord(streamID2, varID2[varID], binID);
		pstreamWriteRecord(streamID2,&val, nmiss);
Cedrick Ansorge's avatar
Cedrick Ansorge committed
389
	      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
390
	      //fprintf(stderr,"\n");
391
	    }
Cedrick Ansorge's avatar
Cedrick Ansorge committed
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
	  else if ( operfunc == func_roc ) 
	    {
	      if ( DEBUG_ROC ) {
		fprintf(stderr, "#             :     TP     FP     FN     TN         TPR        FPR\n");
		
		for ( binID=0; binID<= nbins; binID++ )  {
		  int p = ctg_tab[binID][TP] + ctg_tab[binID][FN];
		  int n = ctg_tab[binID][FP] + ctg_tab[binID][TN];
		  double tpr = ctg_tab[binID][TP]/(double) p;
		  double fpr = ctg_tab[binID][FP]/(double) n;
		  chksum += ctg_tab[binID][0] + ctg_tab[binID][1] + ctg_tab[binID][2] + ctg_tab[binID][3];
		  
		  roc[binID][TPR] = tpr;
		  roc[binID][FPR] = fpr;
		  
		  fprintf(stderr, "%3i %10.4g: %6i %6i %6i %6i: %10.4g %10.4g\n",
			  binID,binID<nbins?lThresh[binID]:1,
			  ctg_tab[binID][0],ctg_tab[binID][1],ctg_tab[binID][2],ctg_tab[binID][3],
			  tpr,fpr);
		}
		fprintf(stderr,"nbins %10i\n",nbins);
		fprintf(stderr,"#ROC CurveArea: %10.6f\n",
			roc_curve_integrate((const double **)roc,nbins));
	      } // if ( DEBUG_ROC )
	    }   // else if (operfunc == func_roc )
	}       // for ( recID=0; recID<nrecs; recID++ )
      tsID++;
    }           // do [...]
  while ( nrecs0 > 0 );
421
422


Cedrick Ansorge's avatar
Cedrick Ansorge committed
423
424

  if ( operfunc == func_rank )
425
426
427
428
429
    {
      double *tmpdoub;
      int osize = gridsize;

      if ( datafunc == TIME ) osize = 1;
430
      tmpdoub = (double*) Malloc(osize*sizeof(double));
431
432
433
434
435
436

      for ( binID = 0; binID < nfiles; binID++ )
	{
	  for ( i = 0; i < osize; i++ )
	    tmpdoub[i] = (double) array2[binID][i];

437
438
	  pstreamDefRecord(streamID2, varID2[varID], binID);
	  pstreamWriteRecord(streamID2, tmpdoub, nmiss);
439
440
	}

441
      Free(tmpdoub);
442
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
443
  else if ( operfunc == func_roc ) {
Cedrick Ansorge's avatar
Cedrick Ansorge committed
444
445
446
447
448
449
450
451
452
453
454
455
    fprintf(stdout, "#             :     TP     FP     FN     TN         TPR        FPR\n");
    
    for ( i=0; i<= nbins; i++ )  {
      int p = ctg_tab[i][TP] + ctg_tab[i][FN];
      int n = ctg_tab[i][FP] + ctg_tab[i][TN];
      double tpr = ctg_tab[i][TP]/(double) p;
      double fpr = ctg_tab[i][FP]/(double) n;
      chksum += ctg_tab[i][0] + ctg_tab[i][1] + ctg_tab[i][2] + ctg_tab[i][3];
      
      roc[i][TPR] = tpr;
      roc[i][FPR] = fpr;
      
456
      int sum = ctg_tab[i][TP] + ctg_tab[i][TN] + ctg_tab[i][FP] + ctg_tab[i][FN];
Cedrick Ansorge's avatar
Cedrick Ansorge committed
457
458
459
460
461
462
463
464
465

      fprintf(stdout, "%3i %10.4g: %6i %6i %6i %6i (%6i): %10.4g %10.4g\n",
	      i,i<nbins?lThresh[i]:1,
	      ctg_tab[i][0],ctg_tab[i][1],ctg_tab[i][2],ctg_tab[i][3],sum,
	      tpr,fpr);
    }
 
    fprintf(stdout,"#ROC CurveArea: %10.6f\n",
	    roc_curve_integrate((const double **)roc,nbins));
466
467
  }

468
  for ( int fileID = 0; fileID < nfiles; fileID++ )
469
470
    {
      streamID = ef[fileID].streamID;
471
      pstreamClose(streamID);
472
473
    }

Cedrick Ansorge's avatar
Cedrick Ansorge committed
474
  if ( operfunc != func_roc ) 
475
    pstreamClose(streamID2);
476

477
  for ( int fileID = 0; fileID < nfiles; fileID++ )
478
    if ( ef[fileID].array ) Free(ef[fileID].array);
479

480
  if ( ef ) Free(ef);
481
482
  if ( array2 ) {
    for (binID=0; binID<nfiles; binID++ )
483
484
      Free(array2[binID]);
    Free(array2);
485
486
487
488
  }

  for ( i = 0; i < ompNumThreads; i++ )
    {
489
490
      if ( field[i].ptr    ) Free(field[i].ptr);
      if ( field[i].weight ) Free(field[i].weight);
491
492
    }
  
493
  if ( field ) Free(field);
494
495
496
  
  cdoFinish();
  
Uwe Schulzweida's avatar
Uwe Schulzweida committed
497
  return 0;
498
}
Cedrick Ansorge's avatar
Cedrick Ansorge committed
499
500
501


double roc_curve_integrate(const double **roc, const int n) {
502
  double y1, y0, x1,x0, dx, dy;
Cedrick Ansorge's avatar
Cedrick Ansorge committed
503
  double step_area;
504
  double area = 0;
Cedrick Ansorge's avatar
Cedrick Ansorge committed
505

506
  for ( int i = 1; i <= n; ++i ) {
Cedrick Ansorge's avatar
Cedrick Ansorge committed
507
508
509
510
511
512
513
514
515
516
517
518
    x1 = roc[i][FPR]; x0 = roc[i-1][FPR];
    y1 = roc[i][TPR]; y0 = roc[i-1][TPR];
    dx = x1-x0;
    dy = y1-y0;

    step_area = -0.5*dx*dy - dx*y0;
    area += step_area;
  }

  return area-0.5;

}