root / src / sim / cvarhist.cc @ e26d3d25
History  View  Annotate  Download (13.4 KB)
1 
//=========================================================================


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) 19922008 Andras Varga

17 
Copyright (C) 20062008 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[i1]; // 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 equiprobable cells (automatic setup) was requested

167 
if (num_cells>0) 
168 
throw cRuntimeError(this, "some bin bounds already present when making equiprobable cells"); 
169  
170 
// setRange() methods must not be used with cVarHistogram's equiprobable cell autosetup 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_cells1 && 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[j1] < firstvals[j] holds

212 
if (j == num_vals)

213 
break; // the cellth 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=index1; 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]=indexprev_index; 
236 
} 
237  
238 
// the last cell/bin:

239 
cellv[cell] = num_valsprev_index; 
240  
241 
// the last boundary:

242 
rangemax = bin_bounds[cell+1]=firstvals[num_vals1]; 
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_vals1]firstvals[0]; 
248 
double epsilon = range*1e6; // 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_valscell_undercell_over);

380  
381 
// select a random interval (k1) 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[k1];

388 
upper = bin_bounds[k]; 
389  
390 
return lower + genk_dblrand(genk)*(upperlower);

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  
468 