Statistics
| Branch: | Revision:

root / src / sim / cvarhist.cc @ 08285dff

History | View | Annotate | Download (13.4 KB)

1 01873262 Georg Kunz
//=========================================================================
2
//  CVARHIST.CC - part of
3
//
4
//                  OMNeT++/OMNEST
5
//           Discrete System Simulation in C++
6
//
7
//   Member functions of
8
//    cVarHistogram     : Variable bin size histogram
9
//
10
//   Authors: Gabor Lencse
11
//   Adapted by: Andras Varga
12
//
13
//=========================================================================
14
15
/*--------------------------------------------------------------*
16
  Copyright (C) 1992-2008 Andras Varga
17
  Copyright (C) 2006-2008 OpenSim Ltd.
18

19
  This file is distributed WITHOUT ANY WARRANTY. See the file
20
  `license' for details on this and other legal matters.
21
*--------------------------------------------------------------*/
22
23
#include <stdio.h>
24
#include <stdlib.h>
25
#include <string.h>
26
#include <math.h>
27
#include "globals.h"
28
#include "random.h"
29
#include "distrib.h"
30
#include "cvarhist.h"
31
#include "cexception.h"
32
33
#ifdef WITH_PARSIM
34
#include "ccommbuffer.h"
35
#endif
36
37
USING_NAMESPACE
38
39
#define MIN(a,b) ((a)<(b) ? (a) : (b))
40
#define MAX(a,b) ((a)>(b) ? (a) : (b))
41
42
Register_Class(cVarHistogram);
43
44
//----
45
// cVarHistogram - member functions
46
47
cVarHistogram::cVarHistogram(const char *name, int maxnumcells, int transformtype) :
48
cHistogramBase(name, -1) //--LG
49
{
50
    // num_cells==-1 means that no bin boundaries are defined (num_cells+1 is 0)
51
    range_mode = RANGE_NOTSET;
52
    transform_type = transformtype;
53
    max_num_cells = maxnumcells;
54
    bin_bounds = NULL;
55
    ASSERT(firstvals); // base class must have allocated it for RANGE_AUTO
56
57
    if ((transform_type == HIST_TR_AUTO_EPC_DBL ||
58
         transform_type == HIST_TR_AUTO_EPC_INT) && max_num_cells<2)
59
        throw cRuntimeError(this, "constructor: the maximal number of cells should be >=2");
60
}
61
62
cVarHistogram::~cVarHistogram()
63
{
64
    delete [] bin_bounds;
65
}
66
67
void cVarHistogram::parsimPack(cCommBuffer *buffer)
68
{
69
#ifndef WITH_PARSIM
70
    throw cRuntimeError(this, eNOPARSIM);
71
#else
72
    cHistogramBase::parsimPack(buffer);
73
74
    buffer->pack(max_num_cells);
75
    if (buffer->packFlag(bin_bounds!=NULL))
76
        buffer->pack(bin_bounds, max_num_cells + 1);
77
#endif
78
}
79
80
void cVarHistogram::parsimUnpack(cCommBuffer *buffer)
81
{
82
#ifndef WITH_PARSIM
83
    throw cRuntimeError(this, eNOPARSIM);
84
#else
85
    cHistogramBase::parsimUnpack(buffer);
86
87
    buffer->unpack(max_num_cells);
88
    if (buffer->checkFlag())
89
    {
90
        bin_bounds = new double[max_num_cells + 1];
91
        buffer->unpack(bin_bounds, max_num_cells + 1);
92
    }
93
#endif
94
}
95
96
void cVarHistogram::addBinBound(double x) //--LG
97
{
98
    if (isTransformed())
99
        throw cRuntimeError(this, "cannot add bin bound after transform()");
100
101
    // create bin_bounds if not exists
102
    if (bin_bounds == NULL)
103
        bin_bounds = new double [max_num_cells+1];
104
105
    // expand if full
106
    if (num_cells == max_num_cells)
107
    {
108
        double * temp = new double [max_num_cells*2+1];
109
        memcpy(temp, bin_bounds, (max_num_cells+1)*sizeof(double));
110
        delete [] bin_bounds;
111
        bin_bounds = temp;
112
        max_num_cells*=2;
113
    }
114
115
    // insert bound
116
    int i;
117
    for (i = num_cells+1; bin_bounds[i]>x; i--)
118
        bin_bounds[i] = bin_bounds[i-1]; // shift up bin boundaries
119
    bin_bounds[i]=x;
120
121
    num_cells++;
122
}
123
124
cVarHistogram& cVarHistogram::operator=(const cVarHistogram& res) //--LG
125
{
126
    if (this==&res) return *this;
127
128
    cHistogramBase::operator=(res);
129
    // hack: as this ^ uses num_cells instead of max_num_cells, we must correct it:
130
    if (res.cellv)
131
    {
132
        delete [] cellv;
133
        cellv = new unsigned [max_num_cells];
134
        memcpy(cellv, res.cellv, num_cells*sizeof(unsigned));
135
    }
136
137
    max_num_cells = res.max_num_cells;
138
    transform_type = res.transform_type;
139
140
    delete [] bin_bounds;
141
    bin_bounds = NULL;
142
    if (res.bin_bounds)
143
    {
144
        bin_bounds = new double [max_num_cells+1];
145
        memcpy(bin_bounds, res.bin_bounds, (num_cells+1)*sizeof(double));
146
    }
147
148
    return *this;
149
}
150
151
static int double_compare_function(const void *p1, const void *p2) //--LG
152
{
153
    double x1 = * (double *) p1;
154
    double x2 = * (double *) p2;
155
156
    if (x1 == x2)
157
        return 0;
158
    else if (x1 < x2)
159
        return -1;
160
    else
161
        return 1;
162
}
163
164
void cVarHistogram::createEquiprobableCells()
165
{
166
    // this method is called from transform() if equi-probable cells (automatic setup) was requested
167
    if (num_cells>0)
168
        throw cRuntimeError(this, "some bin bounds already present when making equi-probable cells");
169
170
    // setRange() methods must not be used with cVarHistogram's equi-probable cell auto-setup mode,
171
    // so range_mode should still be the RANGE_NOTSET that we set in the ctor
172
    if (range_mode != RANGE_NOTSET)
173
        throw cRuntimeError(this, "setRange..() only supported with HIST_TR_NO_TRANSFORM mode");
174
175
    // this version automatically sets the cell boundaries...
176
    ASSERT(max_num_cells>=2); // maybe 1 is enough...
177
178
    // allocate cellv and bin_bounds
179
    cellv = new unsigned [max_num_cells];
180
    bin_bounds = new double [max_num_cells+1];
181
182
    qsort(firstvals, num_vals, sizeof(double), double_compare_function);
183
184
    // expected sample number per cell
185
    double esnpc = num_vals/(double)max_num_cells;
186
187
    int cell;       // index of cell being constructed
188
    int prev_index; // index of first observation in firstvals[] that will go into cellv[cell]
189
    int index;      // index of first observation in firstvals[] that will be left for the next cell (cellv[cell+1])
190
191
    double prev_boundary; // previous value of boundary
192
    double boundary;      // firstvals[index]
193
194
    // construct cells; last cell will be handled separately
195
    for (cell = 0, prev_index = 0, prev_boundary = firstvals[prev_index],
196
         rangemin = bin_bounds[0]=firstvals[0], index = prev_index+(int)esnpc;
197
198
         cell<max_num_cells-1 && index<num_vals;
199
200
         cell++, prev_index = index, prev_boundary = boundary,
201
         index = (int)MAX(prev_index+esnpc, (cell+1)*esnpc))
202
    {
203
        boundary = firstvals[index];
204
        if (boundary == prev_boundary)
205
        {
206
            // try to find a greater one
207
            int j;
208
            for (j=index; j<num_vals && firstvals[j] == prev_boundary; j++)
209
                 ;
210
            // remark: either j == num_vals or
211
            //  prev_boundary == firstvals[j-1] < firstvals[j] holds
212
            if (j == num_vals)
213
                 break; // the cell-th cell/bin will be the last cell/bin
214
            else
215
            {
216
                 index = j; // a greater one was found
217
                 boundary = firstvals[index];
218
            }
219
        }
220
        else
221
        {
222
            // go backwards in firstvals[] to find first observation that
223
            // equals `boundary' (that is, check if there's a j:
224
            // j<index, firstvals[j] == firstvals[index])
225
            int j;
226
            // sure: prev_boundary==firstvals[prev_index] < boundary
227
            //       AND index > prev_index (otherwise   ^^^ here '=' would be)
228
            //        ==> j >= prev_index when firstvals[j] is evaluated
229
            // for this reason we do not need to check j>=0
230
            for (j=index-1; firstvals[j] == boundary; j--)
231
                ; // may run 0 or more times
232
            index = j+1; // unnecessary if cycle ran 0 times
233
        }
234
        bin_bounds[cell+1]=boundary;
235
        cellv[cell]=index-prev_index;
236
    }
237
238
    // the last cell/bin:
239
    cellv[cell] = num_vals-prev_index;
240
241
    // the last boundary:
242
    rangemax = bin_bounds[cell+1]=firstvals[num_vals-1];
243
244
    // correction of the last boundary (depends on DBL/INT)
245
    if (transform_type == HIST_TR_AUTO_EPC_DBL)
246
    {
247
        double range = firstvals[num_vals-1]-firstvals[0];
248
        double epsilon = range*1e-6;   // hack: value < boundary; not '<='
249
        rangemax = bin_bounds[cell+1] += epsilon;
250
    }
251
    else if (transform_type == HIST_TR_AUTO_EPC_INT)
252
    {
253
        rangemax = bin_bounds[cell+1] += 1;  // hack: take the next integer
254
    }
255
256
    // remark: cellv[0]...cellv[cell] are the valid cells
257
    num_cells = cell+1; // maybe num_cells < max_num_cells
258
}
259
260
void cVarHistogram::transform() //--LG
261
{
262
    if (isTransformed())
263
        throw cRuntimeError(this, "transform(): histogram already transformed");
264
265
    setupRange();
266
267
    if (transform_type == HIST_TR_AUTO_EPC_DBL || transform_type == HIST_TR_AUTO_EPC_INT)
268
    {
269
        // create bin bounds based on firstvals[]
270
        createEquiprobableCells();
271
    }
272
    else
273
    {
274
        ASSERT(transform_type == HIST_TR_NO_TRANSFORM);
275
276
        // all manually added bin bounds must be in the range
277
        if (range_mode != RANGE_NOTSET)
278
        {
279
            if (rangemin>bin_bounds[0] || rangemax<bin_bounds[num_cells])
280
                throw cRuntimeError(this, "some bin bounds out of preset range");
281
282
            if (rangemin<bin_bounds[0])
283
                addBinBound(rangemin);
284
            if (rangemax>bin_bounds[num_cells])
285
                addBinBound(rangemax);
286
        }
287
288
        // create cell vector and insert observations
289
        cellv = new unsigned [num_cells];
290
        for (int i=0; i<num_cells; i++)
291
            cellv[i] = 0;
292
293
        for (int i=0; i<num_vals; i++)
294
            collectTransformed(firstvals[i]);
295
    }
296
297
    delete [] firstvals;
298
    firstvals = NULL;
299
300
    transfd = true;
301
}
302
303
void cVarHistogram::collectTransformed(double val)
304
{
305
    if (val < rangemin) // rangemin == bin_bounds[0]
306
    {
307
       cell_under++;
308
    }
309
    else if (val >= rangemax) // rangemax == bin_bounds[num_cells]
310
    {
311
       cell_over++;
312
    }
313
    else // sample falls in the range of ordinary cells/bins
314
    {
315
       // rangemin <= val < rangemax
316
       // binary search
317
       int lower_index, upper_index, index;
318
       for (lower_index = 0, upper_index = num_cells,
319
             index = (lower_index+upper_index)/2;
320
321
             lower_index<index;
322
323
             index = (lower_index+upper_index)/2)
324
       {
325
          // cycle invariant: bin_bound[lower_index]<=val<bin_bounds[upper_index]
326
           if (val < bin_bounds[index])
327
               upper_index = index;
328
           else
329
               lower_index = index;
330
       }
331
       // here, bin_bound[lower_index]<=val<bin_bounds[lower_index+1]
332
333
       // increment the appropriate counter
334
       cellv[lower_index]++;
335
    }
336
}
337
338
// clear results
339
void cVarHistogram::clearResult() //--LG
340
{
341
    cHistogramBase::clearResult();
342
343
    delete [] bin_bounds;
344
    bin_bounds = NULL;
345
}
346
347
// return kth basepoint
348
double cVarHistogram::getBasepoint(int k) const
349
{
350
    if (k<num_cells+1)
351
        return bin_bounds[k];
352
    else
353
        throw cRuntimeError(this, "invalid basepoint index %u", k);
354
}
355
356
double cVarHistogram::getCellValue(int k) const
357
{
358
    if (k<num_cells)
359
        return cellv[k];
360
    else
361
        throw cRuntimeError(this, "invalid cell index %u", k);
362
}
363
364
double cVarHistogram::random() const //--LG
365
{
366
    if (num_vals == 0)
367
        return 0;
368
369
    if (num_vals < num_firstvals)
370
    {
371
        // randomly select a sample from the stored ones
372
        return firstvals[genk_intrand(genk, num_vals)];
373
    }
374
    else
375
    {
376
        double lower, upper;
377
378
        // generate in [lower, upper)
379
        double m = genk_intrand(genk, num_vals-cell_under-cell_over);
380
381
        // select a random interval (k-1) and return a random number from
382
        // that interval generated according to uniform distribution.
383
        m -= cell_under;
384
        int k;
385
        for (k=0; m>=0; k++)
386
            m -= cellv[k];
387
        lower = bin_bounds[k-1];
388
        upper = bin_bounds[k];
389
390
        return lower + genk_dblrand(genk)*(upper-lower);
391
    }
392
}
393
394
double cVarHistogram::getPDF(double x) const // --LG
395
{
396
    if (!num_vals)
397
        return 0.0;
398
399
    if (!isTransformed())
400
        throw cRuntimeError(this, "getPDF(x) cannot be called before histogram is transformed");
401
402
    if (x < rangemin || x >= rangemax)
403
        return 0.0;
404
405
    // rangemin <= x < rangemax
406
    // binary search
407
    int lower_index, upper_index, index;
408
    for (lower_index = 0, upper_index = num_cells,
409
         index = (lower_index+upper_index)/2;
410
411
         lower_index<index;
412
413
         index = (lower_index+upper_index)/2)
414
    {
415
        // cycle invariant: bin_bound[lower_index]<=x<bin_bounds[upper_index]
416
        if (x < bin_bounds[index])
417
            upper_index = index;
418
        else
419
            lower_index = index;
420
    }
421
422
    // here, bin_bound[lower_index]<=x<bin_bounds[lower_index+1]
423
    return cellv[lower_index]/(bin_bounds[lower_index+1]-bin_bounds[lower_index])/num_vals;
424
}
425
426
double cVarHistogram::getCDF(double) const
427
{
428
    throw cRuntimeError(this, "getCDF(x) not implemented");
429
}
430
431
void cVarHistogram::saveToFile(FILE *f) const //--LG
432
{
433
    cHistogramBase::saveToFile(f);
434
435
    fprintf(f, "%d\t #= transform_type\n", transform_type);
436
    fprintf(f, "%u\t #= max_num_cells\n", max_num_cells);
437
438
    fprintf(f, "%d\t #= bin_bounds[] exists\n", bin_bounds!=NULL);
439
    if (bin_bounds) for (int i=0; i<max_num_cells+1; i++) fprintf(f, " %g\n", bin_bounds[i]);
440
}
441
442
void cVarHistogram::loadFromFile(FILE *f)
443
{
444
    cHistogramBase::loadFromFile(f);
445
446
    freadvarsf(f, "%d\t #= transform_type", &transform_type);
447
    freadvarsf(f, "%u\t #= max_num_cells", &max_num_cells);
448
449
    // increase allocated size of cellv[] to max_num_cells
450
    if (cellv && max_num_cells>num_cells)
451
    {
452
        unsigned int *new_cellv = new unsigned [max_num_cells];
453
        memcpy(new_cellv, cellv, num_cells*sizeof(unsigned));
454
        delete [] cellv; cellv = new_cellv;
455
    }
456
457
    int binbounds_exists;
458
    freadvarsf(f, "%d\t #= bin_bounds[] exists", &binbounds_exists);
459
    delete [] bin_bounds; bin_bounds = NULL;
460
    if (binbounds_exists)
461
    {
462
        bin_bounds = new double[max_num_cells+1];
463
        for (int i=0; i<max_num_cells+1; i++) freadvarsf(f, " %g", bin_bounds+i);
464
    }
465
}
466
467