## root / src / sim / cvarhist.cc @ e26d3d25

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 |