## root / src / sim / cksplit.cc @ fbe00e73

History | View | Annotate | Download (20.5 KB)

1 | 01873262 | Georg Kunz | ```
//=========================================================================
``` |
---|---|---|---|

2 | ```
// CKSPLIT.CC - part of
``` |
||

3 | ```
//
``` |
||

4 | ```
// OMNeT++/OMNEST
``` |
||

5 | ```
// Discrete System Simulation in C++
``` |
||

6 | ```
//
``` |
||

7 | ```
// Member functions of
``` |
||

8 | ```
// cKSplit : implements the k-split algorithm in 1 dimension
``` |
||

9 | ```
//
``` |
||

10 | ```
// Original version by Babak Fakhamzadeh, TU Delft, Mar-Jun 1996
``` |
||

11 | ```
// This version written by Andras Varga, 1997-98
``` |
||

12 | ```
//
``` |
||

13 | ```
//=========================================================================
``` |
||

14 | ```
/*--------------------------------------------------------------*
``` |
||

15 | ```
Copyright (C) 1992-2008 Andras Varga
``` |
||

16 | ```
Copyright (C) 2006-2008 OpenSim Ltd.
``` |
||

17 | |||

18 | ```
This file is distributed WITHOUT ANY WARRANTY. See the file
``` |
||

19 | ```
`license' for details on this and other legal matters.
``` |
||

20 | ```
*--------------------------------------------------------------*/
``` |
||

21 | |||

22 | #include <stdio.h> |
||

23 | #include <math.h> |
||

24 | #include <stdlib.h> |
||

25 | #include "globals.h" |
||

26 | #include "cenvir.h" |
||

27 | #include "cksplit.h" |
||

28 | #include "random.h" |
||

29 | #include "distrib.h" |
||

30 | #include "cexception.h" |
||

31 | |||

32 | ```
#ifdef WITH_PARSIM
``` |
||

33 | #include "ccommbuffer.h" |
||

34 | ```
#endif
``` |
||

35 | |||

36 | ```
using std::ostream;
``` |
||

37 | ```
using std::endl;
``` |
||

38 | |||

39 | NAMESPACE_BEGIN |
||

40 | |||

41 | Register_Class(cKSplit); |
||

42 | |||

43 | |||

44 | ```
//----
``` |
||

45 | ```
// Cell division criteria - they are used to decide whether a cell should be split.
``` |
||

46 | |||

47 | double critdata_default[] = {20, 4, 2}; |
||

48 | |||

49 | int critfunc_const(const cKSplit&, cKSplit::Grid& g, int i, double *c) |
||

50 | { |
||

51 | return g.cells[i] >= c[0]; |
||

52 | } |
||

53 | |||

54 | int critfunc_depth(const cKSplit& ks, cKSplit::Grid& g, int i, double *c) |
||

55 | { |
||

56 | ```
int depth = g.reldepth - ks.getRootGrid().reldepth;
``` |
||

57 | return g.cells[i] >= c[1]*pow(c[2], depth); |
||

58 | } |
||

59 | |||

60 | double divdata_default[] = {0.5}; |
||

61 | |||

62 | double divfunc_const(const cKSplit&, cKSplit::Grid&, double, double *d) |
||

63 | { |
||

64 | return d[0]; // lambda=constant |
||

65 | } |
||

66 | |||

67 | double divfunc_babak(const cKSplit&, cKSplit::Grid& g, double mother, double *d) |
||

68 | { |
||

69 | ```
int newobs = g.total-g.mother;
``` |
||

70 | double lambda = newobs / (d[1]*mother); |
||

71 | return lambda<1.0 ? lambda : 1.0; |
||

72 | } |
||

73 | |||

74 | ```
//----
``` |
||

75 | |||

76 | ```
cKSplit::cKSplit(const cKSplit& r) : cDensityEstBase()
``` |
||

77 | { |
||

78 | setName(r.getName()); |
||

79 | gridv = NULL; iter = NULL; |
||

80 | ```
operator=(r);
``` |
||

81 | } |
||

82 | |||

83 | cKSplit::cKSplit(const char *name) : cDensityEstBase(name) |
||

84 | { |
||

85 | if (K<2 || (K>2 && K!=2*(int)(K/2)+1)) |
||

86 | throw cRuntimeError("cKSplit: K must be 2 or a >=3 odd number"); |
||

87 | |||

88 | ```
num_cells = 0;
``` |
||

89 | |||

90 | critfunc = critfunc_depth; |
||

91 | critdata = critdata_default; |
||

92 | divfunc = divfunc_const; |
||

93 | divdata = divdata_default; |
||

94 | ```
rangeext_enabled = true;
``` |
||

95 | |||

96 | ```
gridv = NULL;
``` |
||

97 | ```
gridv_size = 0;
``` |
||

98 | |||

99 | ```
iter = NULL;
``` |
||

100 | } |
||

101 | |||

102 | cKSplit::~cKSplit() |
||

103 | { |
||

104 | ```
delete [] gridv;
``` |
||

105 | ```
delete iter;
``` |
||

106 | } |
||

107 | |||

108 | ```
void cKSplit::parsimPack(cCommBuffer *buffer)
``` |
||

109 | { |
||

110 | ```
// cDensityEstBase::parsimPack(buffer);
``` |
||

111 | throw cRuntimeError(this, "parsimPack() not implemented"); |
||

112 | } |
||

113 | |||

114 | ```
void cKSplit::parsimUnpack(cCommBuffer *buffer)
``` |
||

115 | { |
||

116 | ```
// cDensityEstBase::parsimUnpack(buffer);
``` |
||

117 | throw cRuntimeError(this, "parsimUnpack() not implemented"); |
||

118 | } |
||

119 | |||

120 | cKSplit& cKSplit::operator=(const cKSplit& res) |
||

121 | { |
||

122 | if (this==&res) return *this; |
||

123 | |||

124 | ```
cDensityEstBase::operator=(res);
``` |
||

125 | |||

126 | num_cells = res.num_cells; |
||

127 | |||

128 | rootgrid = res.rootgrid; |
||

129 | lastgrid = res.lastgrid; |
||

130 | gridv_size = res.gridv_size; |
||

131 | |||

132 | ```
delete [] gridv;
``` |
||

133 | ```
if (!res.gridv)
``` |
||

134 | ```
gridv = NULL;
``` |
||

135 | ```
else
``` |
||

136 | { |
||

137 | gridv = new Grid[gridv_size + 1]; |
||

138 | for (int i = 1; i <= lastgrid; i++) |
||

139 | gridv[i] = res.gridv[i]; |
||

140 | } |
||

141 | |||

142 | critfunc = res.critfunc; |
||

143 | critdata = res.critdata; |
||

144 | divfunc = res.divfunc; |
||

145 | divdata = res.divdata; |
||

146 | rangeext_enabled = res.rangeext_enabled; |
||

147 | |||

148 | ```
delete iter;
``` |
||

149 | ```
iter = NULL;
``` |
||

150 | |||

151 | return (*this); |
||

152 | } |
||

153 | |||

154 | std::string cKSplit::detailedInfo() const |
||

155 | { |
||

156 | std::stringstream os; |
||

157 | os << cDensityEstBase::detailedInfo(); |
||

158 | |||

159 | int nn = num_vals <= num_cells+1 ? num_vals : num_cells+1; //??? |
||

160 | |||

161 | ```
os << "\n cells:\n";
``` |
||

162 | for (int i=0; i<nn; i++) |
||

163 | os << " >=" << getBasepoint(i) << ": " << getCellValue(i) << endl; |
||

164 | ```
return os.str();
``` |
||

165 | } |
||

166 | |||

167 | void cKSplit::setCritFunc(CritFunc _critfunc, double *_critdata) |
||

168 | { |
||

169 | critfunc = _critfunc; |
||

170 | critdata = _critdata; |
||

171 | } |
||

172 | |||

173 | void cKSplit::setDivFunc(DivFunc _divfunc, double *_divdata) |
||

174 | { |
||

175 | divfunc = _divfunc; |
||

176 | divdata = _divdata; |
||

177 | } |
||

178 | |||

179 | void cKSplit::rangeExtension(bool enabled) |
||

180 | { |
||

181 | rangeext_enabled = enabled; |
||

182 | } |
||

183 | |||

184 | void cKSplit::resetGrids(int grid) |
||

185 | { |
||

186 | Grid *g = &(gridv[grid]); |
||

187 | ```
g->total = g->mother = 0;
``` |
||

188 | for (int i=0; i<K; i++) |
||

189 | { |
||

190 | if (g->cells[i] < 0) |
||

191 | resetGrids(- g->cells[i]); |
||

192 | ```
else
``` |
||

193 | ```
g->cells[i] = 0;
``` |
||

194 | } |
||

195 | } |
||

196 | |||

197 | void cKSplit::merge(const cStatistic *other) |
||

198 | { |
||

199 | throw cRuntimeError(this, "The cKSplit class does not support merge()"); |
||

200 | } |
||

201 | |||

202 | void cKSplit::doMergeCellValues(const cDensityEstBase *other) |
||

203 | { |
||

204 | ASSERT(false); // never comes here, as merge() already throws an error |
||

205 | } |
||

206 | |||

207 | ```
void cKSplit::transform()
``` |
||

208 | { |
||

209 | ```
if (isTransformed())
``` |
||

210 | throw cRuntimeError(this, "transform(): histogram already transformed"); |
||

211 | |||

212 | setupRange(); |
||

213 | createRootGrid(); |
||

214 | |||

215 | ```
// Trick to avoid error from cell splits: first, insert observations
``` |
||

216 | ```
// just to create grid structure; second, reset all cells to zero and
``` |
||

217 | ```
// insert observations again, with cell splits disabled now.
``` |
||

218 | |||

219 | for (int i=0; i<num_vals; i++) |
||

220 | collectTransformed(firstvals[i]); |
||

221 | |||

222 | resetGrids(rootgrid); |
||

223 | |||

224 | for (int i=0; i<num_vals; i++) |
||

225 | ```
insertIntoGrids(firstvals[i], false);
``` |
||

226 | |||

227 | ```
delete [] firstvals;
``` |
||

228 | ```
firstvals = NULL;
``` |
||

229 | |||

230 | ```
transfd = true;
``` |
||

231 | } |
||

232 | |||

233 | ```
void cKSplit::createRootGrid()
``` |
||

234 | { |
||

235 | ```
gridv_size = 8;
``` |
||

236 | gridv = new Grid[gridv_size+1]; |
||

237 | |||

238 | ```
rootgrid = 1;
``` |
||

239 | ```
lastgrid = 1;
``` |
||

240 | num_cells = K; |
||

241 | |||

242 | gridv[1].parent = 0; |
||

243 | gridv[1].reldepth = 0; |
||

244 | gridv[1].total = 0; |
||

245 | gridv[1].mother = 0; |
||

246 | for (int i=0; i<K; i++) |
||

247 | gridv[1].cells[i] = 0; |
||

248 | } |
||

249 | |||

250 | void cKSplit::collectTransformed(double x) |
||

251 | { |
||

252 | ```
// see if x fits into current range and act accordingly
``` |
||

253 | ```
if (x >= rangemin && x < rangemax)
``` |
||

254 | ```
insertIntoGrids(x, true);
``` |
||

255 | else if (rangeext_enabled) |
||

256 | newRootGrids(x); |
||

257 | else if (x < rangemin) |
||

258 | cell_under++; |
||

259 | else if (x >= rangemax) |
||

260 | cell_over++; |
||

261 | } |
||

262 | |||

263 | void cKSplit::insertIntoGrids(double x, int enable_splits) |
||

264 | { |
||

265 | ```
int i;
``` |
||

266 | |||

267 | ```
// new observation fits in the current grid
``` |
||

268 | ```
double gridmin = rangemin;
``` |
||

269 | ```
double gridmax = rangemax;
``` |
||

270 | double cellsize = (gridmax - gridmin) / (double)K; |
||

271 | |||

272 | ```
int location = rootgrid;
``` |
||

273 | |||

274 | ```
// keep searching until the right grid has been found
``` |
||

275 | ```
for(;;)
``` |
||

276 | { |
||

277 | gridv[location].total++; |
||

278 | |||

279 | ```
// calc. the cell in which the new observation falls
``` |
||

280 | ```
i = (int) ((x-gridmin) / cellsize);
``` |
||

281 | |||

282 | ```
// guard against rounding errors
``` |
||

283 | if (i<0) i=0; |
||

284 | if (i>=K) i=K-1; |
||

285 | |||

286 | ```
// exit loop if no subgrid
``` |
||

287 | if (gridv[location].cells[i] >= 0) |
||

288 | ```
break;
``` |
||

289 | |||

290 | ```
// go down to subgrid
``` |
||

291 | location = - gridv[location].cells[i]; |
||

292 | |||

293 | gridmin += i * cellsize; |
||

294 | ```
//gridmax = gridmin + cellsize;
``` |
||

295 | cellsize /= K; |
||

296 | } |
||

297 | |||

298 | ```
// arrived at gridv[location].cells[i] -- possibly split this cell
``` |
||

299 | if (enable_splits && critfunc(*this, gridv[location], i, critdata)) |
||

300 | { |
||

301 | splitCell(location, i); |
||

302 | |||

303 | ```
// go down to new subgrid and insert the observation there
``` |
||

304 | ```
// (NOTE: code is copied here from the above loop)
``` |
||

305 | |||

306 | location = - gridv[location].cells[i]; |
||

307 | |||

308 | gridmin += i * cellsize; |
||

309 | ```
//gridmax = gridmin + cellsize;
``` |
||

310 | cellsize /= K; |
||

311 | |||

312 | gridv[location].total++; |
||

313 | |||

314 | ```
i = (int) ((x-gridmin) / cellsize);
``` |
||

315 | |||

316 | if (i<0) i=0; |
||

317 | if (i>=K) i=K-1; |
||

318 | } |
||

319 | |||

320 | ```
// increment cell value
``` |
||

321 | gridv[location].cells[i]++; |
||

322 | } |
||

323 | |||

324 | void cKSplit::splitCell(int grid, int cell) |
||

325 | { |
||

326 | ```
// make room for new grid
``` |
||

327 | ```
if (lastgrid == gridv_size)
``` |
||

328 | expandGridVector(); |
||

329 | |||

330 | ```
// create shorthands for the two grids in question
``` |
||

331 | Grid& g = gridv[grid]; |
||

332 | ```
Grid& subg = gridv[lastgrid+1];
``` |
||

333 | |||

334 | ```
// if none of the cells in the current grid has divided yet, than
``` |
||

335 | ```
// the observations from the mother cell have to be divided as well
``` |
||

336 | ```
#ifdef DISTRIBUTE_ON_CHILD_SPLIT
``` |
||

337 | if (g.mother>0) distributeMotherObservations(grid); |
||

338 | ```
#endif
``` |
||

339 | |||

340 | ```
// insert subg into g.cell[cell]
``` |
||

341 | subg.parent = grid; |
||

342 | ```
subg.reldepth = g.reldepth+1;
``` |
||

343 | subg.mother = g.cells[cell]; |
||

344 | subg.total = subg.mother; |
||

345 | for (int i=0; i<K; i++) |
||

346 | ```
subg.cells[i] = 0;
``` |
||

347 | |||

348 | ```
g.cells[cell] = -(lastgrid+1);
``` |
||

349 | |||

350 | lastgrid++; |
||

351 | ```
num_cells += K-1;
``` |
||

352 | } |
||

353 | |||

354 | ```
#ifdef DISTRIBUTE_ON_CHILD_SPLIT
``` |
||

355 | void cKSplit::distributeMotherObservations(int grid) |
||

356 | { |
||

357 | Grid& g = gridv[grid]; |
||

358 | Grid orig = g; |
||

359 | |||

360 | for(int i=0; i<K; i++) |
||

361 | g.cells[i] = (unsigned)getRealCellValue(orig, i); // WRONG!!! |
||

362 | |||

363 | ```
g.mother = 0;
``` |
||

364 | } |
||

365 | ```
#endif
``` |
||

366 | |||

367 | void cKSplit::newRootGrids(double x) |
||

368 | { |
||

369 | ```
// new "supergrid" has to be inserted until it includes x
``` |
||

370 | ```
for(;;)
``` |
||

371 | { |
||

372 | ```
// a new grid still has to be inserted around the current root grid
``` |
||

373 | ```
if (lastgrid == gridv_size)
``` |
||

374 | expandGridVector(); |
||

375 | |||

376 | ```
int old_rootgrid = rootgrid;
``` |
||

377 | ```
rootgrid = lastgrid+1;
``` |
||

378 | |||

379 | lastgrid++; |
||

380 | ```
num_cells += K-1;
``` |
||

381 | |||

382 | gridv[old_rootgrid].parent = rootgrid; |
||

383 | ```
gridv[rootgrid].parent = 0;
``` |
||

384 | ```
gridv[rootgrid].reldepth = gridv[old_rootgrid].reldepth - 1;
``` |
||

385 | gridv[rootgrid].total = gridv[old_rootgrid].total; |
||

386 | ```
gridv[rootgrid].mother = 0;
``` |
||

387 | for (int i=0; i<K; i++) |
||

388 | ```
gridv[rootgrid].cells[i] = 0;
``` |
||

389 | |||

390 | ```
double gridsize = rangemax - rangemin;
``` |
||

391 | #if K==2 |
||

392 | ```
if (x<rangemin)
``` |
||

393 | { |
||

394 | ```
gridv[rootgrid].cells[1] = -old_rootgrid;
``` |
||

395 | rangemin -= gridsize; |
||

396 | } |
||

397 | else // (x>=rangemax) |
||

398 | { |
||

399 | ```
gridv[rootgrid].cells[0] = -old_rootgrid;
``` |
||

400 | rangemax += gridsize; |
||

401 | } |
||

402 | ```
#else
``` |
||

403 | gridv[rootgrid].cells[(K-1)/2] = -old_rootgrid; |
||

404 | |||

405 | rangemin -= (K-1)/2*gridsize; |
||

406 | rangemax += (K-1)/2*gridsize; |
||

407 | ```
#endif
``` |
||

408 | |||

409 | ```
if (x >= rangemin && x < rangemax)
``` |
||

410 | ```
break;
``` |
||

411 | } |
||

412 | |||

413 | ```
// now, insert x into new root grid
``` |
||

414 | |||

415 | ```
// calc. in which cell x falls
``` |
||

416 | int i = (int)(K * (x-rangemin) / (rangemax-rangemin)); |
||

417 | |||

418 | ```
// if it falls in the old root grid, we have to correct it
``` |
||

419 | if (i == (K-1)/2) |
||

420 | { |
||

421 | if (x > (rangemax-rangemin)/2) |
||

422 | i++; |
||

423 | ```
else
``` |
||

424 | i--; |
||

425 | } |
||

426 | |||

427 | ```
// increment cell value
``` |
||

428 | gridv[rootgrid].cells[i]++; |
||

429 | gridv[rootgrid].total++; |
||

430 | } |
||

431 | |||

432 | ```
void cKSplit::expandGridVector()
``` |
||

433 | { |
||

434 | ```
gridv_size += 8;
``` |
||

435 | Grid *newgridv = new Grid[gridv_size+1]; |
||

436 | |||

437 | for (int i = 1; i <= lastgrid; i++) |
||

438 | newgridv[i] = gridv[i]; |
||

439 | |||

440 | ```
delete [] gridv;
``` |
||

441 | gridv = newgridv; |
||

442 | } |
||

443 | |||

444 | double cKSplit::getRealCellValue(Grid& grid, int i) const |
||

445 | { |
||

446 | ```
// returns the actual amount of observations in cell 'i' of 'grid'
``` |
||

447 | |||

448 | ```
#ifdef DISTRIBUTE_ON_CHILD_SPLIT
``` |
||

449 | ```
double mother = grid.mother;
``` |
||

450 | ```
#else
``` |
||

451 | ```
// might go up until the root grid to collect mother observations
``` |
||

452 | ```
double mother;
``` |
||

453 | if (grid.parent == 0) |
||

454 | { |
||

455 | ```
mother = 0;
``` |
||

456 | } |
||

457 | ```
else
``` |
||

458 | { |
||

459 | ```
// find (into k) which cell of the parent our 'grid' was
``` |
||

460 | ```
int k;
``` |
||

461 | Grid& parentgrid = gridv[grid.parent]; |
||

462 | ```
int gridnum = &grid - gridv;
``` |
||

463 | for (k=0; k<K; k++) |
||

464 | ```
if (parentgrid.cells[k] == -gridnum)
``` |
||

465 | ```
break;
``` |
||

466 | |||

467 | ```
// parent grid's undistributed observations
``` |
||

468 | mother = getRealCellValue(parentgrid, k); |
||

469 | } |
||

470 | ```
#endif
``` |
||

471 | ```
// interpolate between even and proportional division of the mother obs.
``` |
||

472 | double lambda = divfunc(*this, grid, mother, divdata); |
||

473 | |||

474 | int cell_tot = grid.cells[i]; if (cell_tot<0) cell_tot = gridv[-cell_tot].total; |
||

475 | int cell_mot = grid.cells[i]; if (cell_mot<0) cell_mot = gridv[-cell_mot].mother; |
||

476 | |||

477 | ```
double even = mother / K;
``` |
||

478 | ```
double prop = mother * cell_tot/(grid.total-grid.mother);
``` |
||

479 | |||

480 | return cell_mot + (1-lambda)*even + lambda*prop; |
||

481 | } |
||

482 | |||

483 | int cKSplit::getTreeDepth() const |
||

484 | { |
||

485 | ```
return getTreeDepth(gridv[rootgrid]);
``` |
||

486 | } |
||

487 | |||

488 | int cKSplit::getTreeDepth(Grid& grid) const |
||

489 | { |
||

490 | int d, maxd = 0; |
||

491 | ```
double c;
``` |
||

492 | for (int i=0; i<K; i++) |
||

493 | { |
||

494 | c = grid.cells[i]; |
||

495 | if (c<0) |
||

496 | { |
||

497 | ```
d = getTreeDepth(gridv[-(int)c]);
``` |
||

498 | ```
if (d>maxd) maxd = d;
``` |
||

499 | } |
||

500 | } |
||

501 | return maxd+1; |
||

502 | } |
||

503 | |||

504 | void cKSplit::printGrids() const |
||

505 | { |
||

506 | ```
if (!isTransformed())
``` |
||

507 | { |
||

508 | ```
ev << "collecting initial observations; no grids yet.\n";
``` |
||

509 | ```
return;
``` |
||

510 | } |
||

511 | ev << "Range: " << rangemin << "..." << rangemax << "\n"; |
||

512 | ev << "Root grid: #" << rootgrid << "\n"; |
||

513 | for (int i = 1; i <= lastgrid; i++) |
||

514 | { |
||

515 | ev << "grid #" << i << ": parent=#" << gridv[i].parent; |
||

516 | ```
ev << " total=" << gridv[i].total;
``` |
||

517 | ev << " mother=" << gridv[i].mother << " ("; |
||

518 | |||

519 | for (int j=0; j<K; j++) |
||

520 | if (gridv[i].cells[j] < 0) |
||

521 | ```
ev << " " << gridv[i].cells[j];
``` |
||

522 | ```
else
``` |
||

523 | ```
ev << " " << gridv[i].cells[j];
``` |
||

524 | ```
ev << ")\n";
``` |
||

525 | } |
||

526 | } |
||

527 | |||

528 | void cKSplit::iteratorToCell(int cell_nr) const |
||

529 | { |
||

530 | ```
// create iterator or reinit if it is stale
``` |
||

531 | ```
iter=0;
``` |
||

532 | ```
if (!iter)
``` |
||

533 | {iter = new Iterator(*this); iter_num_vals = num_vals;} |
||

534 | else if (num_vals!=iter_num_vals) |
||

535 | {iter->init(*this, cell_nr<num_cells/2); iter_num_vals = num_vals;} |
||

536 | |||

537 | ```
// drive iterator up or down to reach cell_nr
``` |
||

538 | ```
if (cell_nr>iter->getCellNumber())
``` |
||

539 | ```
while (cell_nr!=iter->getCellNumber())
``` |
||

540 | (*iter)++; |
||

541 | else if (cell_nr<iter->getCellNumber()) |
||

542 | ```
while (cell_nr!=iter->getCellNumber())
``` |
||

543 | (*iter)--; |
||

544 | } |
||

545 | |||

546 | int cKSplit::getNumCells() const |
||

547 | { |
||

548 | ```
if (!isTransformed())
``` |
||

549 | return 0; |
||

550 | ```
return num_cells;
``` |
||

551 | } |
||

552 | |||

553 | double cKSplit::getCellValue(int nr) const |
||

554 | { |
||

555 | ```
if (nr >= num_cells)
``` |
||

556 | return 0.0; |
||

557 | |||

558 | iteratorToCell(nr); |
||

559 | ```
return iter->getCellValue();
``` |
||

560 | } |
||

561 | |||

562 | double cKSplit::getBasepoint(int nr) const |
||

563 | { |
||

564 | ```
if (nr >= num_cells)
``` |
||

565 | ```
return rangemax;
``` |
||

566 | |||

567 | iteratorToCell(nr); |
||

568 | ```
return iter->getCellMin();
``` |
||

569 | } |
||

570 | |||

571 | |||

572 | ```
/////////////////////////////////////////////////////////////////
``` |
||

573 | ```
// random()
``` |
||

574 | ```
// This random number generator works by generating a random number
``` |
||

575 | ```
// between 0 and the total number of observations. After that, it walks
``` |
||

576 | ```
// through the cells, until the total of observations encountered
``` |
||

577 | ```
// is larger than the random number drawn.
``` |
||

578 | ```
//
``` |
||

579 | ```
//
``` |
||

580 | ```
// CODE NOT CLEANED UP BY VA YET!
``` |
||

581 | ```
/////////////////////////////////////////////////////////////////
``` |
||

582 | |||

583 | double cKSplit::random() const |
||

584 | { |
||

585 | ```
int i;
``` |
||

586 | ```
//int dp = getTreeDepth();
``` |
||

587 | int cd = 1; |
||

588 | |||

589 | ```
double x = genk_intrand(genk, num_vals);
``` |
||

590 | |||

591 | ```
int location = rootgrid;
``` |
||

592 | |||

593 | ```
double xi = rangemin, xa;
``` |
||

594 | |||

595 | short finish1 = 0; |
||

596 | |||

597 | while (finish1 == 0) |
||

598 | { |
||

599 | ```
i = 0;
``` |
||

600 | short finish2 = 0; |
||

601 | |||

602 | double sum = 0; |
||

603 | ```
double hlp, hlp4;
``` |
||

604 | |||

605 | while (finish2 == 0) |
||

606 | { |
||

607 | Grid lp = gridv[location]; |
||

608 | |||

609 | if (lp.cells[i] >= 0) |
||

610 | hlp = getRealCellValue(lp, i); |
||

611 | ```
else
``` |
||

612 | hlp = lp.cells[i]; |
||

613 | |||

614 | ```
hlp4 = 0;
``` |
||

615 | |||

616 | if (hlp < 0) |
||

617 | { |
||

618 | int hlp3 = - (int)hlp; |
||

619 | hlp = gridv[hlp3].total; |
||

620 | hlp4 = gridv[hlp3].mother; |
||

621 | } |
||

622 | |||

623 | ```
if ((sum + hlp + hlp4) > x)
``` |
||

624 | ```
finish2 = 1;
``` |
||

625 | ```
else {
``` |
||

626 | sum += hlp + hlp4; |
||

627 | i ++; |
||

628 | } |
||

629 | ```
} //while (finish2..
``` |
||

630 | x -= sum; |
||

631 | |||

632 | ```
xi += i * (rangemax - rangemin) / pow ((double)K, cd);
``` |
||

633 | |||

634 | if (gridv[location].cells[i] < 0) |
||

635 | { |
||

636 | ```
location = -(int)gridv[location].cells[i];
``` |
||

637 | cd++; |
||

638 | } |
||

639 | ```
else
``` |
||

640 | ```
finish1 = 1;
``` |
||

641 | |||

642 | ```
} //while (finish1..
``` |
||

643 | |||

644 | ```
xa = xi + (rangemax - rangemin) / pow ((double)K, cd);
``` |
||

645 | |||

646 | ```
x = rand() / (double)RAND_MAX * (xa - xi) + xi;
``` |
||

647 | ```
return x;
``` |
||

648 | } |
||

649 | |||

650 | ```
/////////////////////////////////////////////////////////////////////////////////////////////
``` |
||

651 | ```
// getPDF()
``` |
||

652 | ```
// The location of the observation in the tree of cells is found. After that, the pdf
``` |
||

653 | ```
// value of that location is calculated and returned. When the observation falls outside of
``` |
||

654 | ```
// root grid, zero is returned.
``` |
||

655 | ```
// The vars dpth and cdpth are used to calculate the actual pdf value of a possibly divided
``` |
||

656 | ```
// cell.
``` |
||

657 | ```
//
``` |
||

658 | ```
// CODE NOT CLEANED UP BY VA YET!
``` |
||

659 | ```
/////////////////////////////////////////////////////////////////////////////////////////////
``` |
||

660 | |||

661 | double cKSplit::getPDF(double x) const |
||

662 | { |
||

663 | ```
if (!isTransformed())
``` |
||

664 | return 0; |
||

665 | |||

666 | ```
int i;
``` |
||

667 | |||

668 | ```
int dpth = getTreeDepth();
``` |
||

669 | int cdpth = 1; |
||

670 | |||

671 | double xi = rangemin; //used for computing the boundary values of the current |
||

672 | double xa = rangemax; //location |
||

673 | |||

674 | int location = rootgrid; //start at the beginning of the gridv |
||

675 | short finish = 0; |
||

676 | |||

677 | while (finish == 0) //keep searching until the right grid has been found |
||

678 | ```
{ //returns the cellnr, in which the point falls
``` |
||

679 | ```
i = (int)(K * (x - xi) / (xa - xi));
``` |
||

680 | |||

681 | if (gridv[location].cells[i] < 0) |
||

682 | { |
||

683 | cdpth ++; |
||

684 | |||

685 | location = - gridv[location].cells[i]; |
||

686 | ```
double hlp = xi;
``` |
||

687 | |||

688 | xi += i * (xa - hlp) / K; |
||

689 | xa = xi + (xa - hlp) / K; |
||

690 | } |
||

691 | else //the right grid has been found |
||

692 | ```
finish = 1;
``` |
||

693 | ```
} //while...
``` |
||

694 | |||

695 | Grid lp = gridv[location]; |
||

696 | |||

697 | return getRealCellValue(lp, i) / pow ((double)K, dpth - cdpth); |
||

698 | } |
||

699 | |||

700 | double cKSplit::getCDF(double) const |
||

701 | { |
||

702 | throw cRuntimeError(this, "getCDF() not implemented"); |
||

703 | } |
||

704 | |||

705 | void cKSplit::saveToFile(FILE *f) const |
||

706 | { |
||

707 | cDensityEstBase::saveToFile(f); |
||

708 | |||

709 | ```
fprintf(f, "%d\t #= num_cells\n", num_cells);
``` |
||

710 | ```
fprintf(f, "%d\t #= gridv_size\n", gridv_size);
``` |
||

711 | ```
fprintf(f, "%d\t #= rootgrid\n", rootgrid);
``` |
||

712 | ```
fprintf(f, "%d\t #= lastgrid\n", lastgrid);
``` |
||

713 | |||

714 | fprintf(f, "%d\t #= gridv_exists\n", gridv!=NULL); |
||

715 | ```
if (gridv)
``` |
||

716 | { |
||

717 | for (int i = 1; i <= lastgrid; i++) |
||

718 | { |
||

719 | ```
fprintf(f, "# grid %d\n", i);
``` |
||

720 | ```
fprintf(f, "%d\t #= parent\n", gridv[i].parent);
``` |
||

721 | ```
fprintf(f, "%d\t #= reldepth\n", gridv[i].reldepth);
``` |
||

722 | ```
fprintf(f, "%ld\t #= total\n", gridv[i].total);
``` |
||

723 | ```
fprintf(f, "%d\t #= mother\n", gridv[i].mother);
``` |
||

724 | #if K==2 |
||

725 | fprintf(f, "%d %d\t #= cells[0], cells[1]\n", gridv[i].cells[0], gridv[i].cells[1]); |
||

726 | ```
#else
``` |
||

727 | ```
fprintf(f, "#= cells[]\n");
``` |
||

728 | for (int j=0; j<K; j++) fprintf(f, " %d\n", gridv[i].cells[j]); |
||

729 | ```
#endif
``` |
||

730 | } |
||

731 | } |
||

732 | } |
||

733 | |||

734 | ```
void cKSplit::loadFromFile(FILE *f)
``` |
||

735 | { |
||

736 | cDensityEstBase::loadFromFile(f); |
||

737 | |||

738 | ```
freadvarsf(f, "%d\t #= num_cells", &num_cells);
``` |
||

739 | ```
freadvarsf(f, "%d\t #= gridv_size", &gridv_size);
``` |
||

740 | ```
freadvarsf(f, "%d\t #= rootgrid", &rootgrid);
``` |
||

741 | ```
freadvarsf(f, "%d\t #= lastgrid", &lastgrid);
``` |
||

742 | |||

743 | ```
int gridv_exists;
``` |
||

744 | ```
freadvarsf(f, "%d\t #= gridv_exists", &gridv_exists);
``` |
||

745 | |||

746 | ```
delete [] gridv;
``` |
||

747 | ```
gridv = NULL;
``` |
||

748 | ```
if (gridv_exists)
``` |
||

749 | { |
||

750 | gridv = new Grid[gridv_size+1]; |
||

751 | for (int i = 1; i <= lastgrid; i++) |
||

752 | { |
||

753 | ```
freadvarsf(f, "");
``` |
||

754 | ```
freadvarsf(f, "%d\t #= parent", &gridv[i].parent);
``` |
||

755 | ```
freadvarsf(f, "%d\t #= reldepth", &gridv[i].reldepth);
``` |
||

756 | ```
freadvarsf(f, "%ld\t #= total", &gridv[i].total);
``` |
||

757 | ```
freadvarsf(f, "%d\t #= mother", &gridv[i].mother);
``` |
||

758 | #if K==2 |
||

759 | freadvarsf(f, "%d %d\t #= cells[0], cells[1]", gridv[i].cells+0, gridv[i].cells+1); |
||

760 | ```
#else
``` |
||

761 | ```
freadvarsf(f, "#= cells[]");
``` |
||

762 | for (int j=0; j<K; j++) freadvarsf(f, " %d", gridv[i].cells+j); |
||

763 | ```
#endif
``` |
||

764 | } |
||

765 | } |
||

766 | } |
||

767 | |||

768 | ```
//----
``` |
||

769 | |||

770 | cKSplit::Iterator::Iterator(const cKSplit& _ks, bool _beg) |
||

771 | { |
||

772 | init(_ks, _beg); |
||

773 | } |
||

774 | |||

775 | void cKSplit::Iterator::init(const cKSplit& _ks, bool _beg) |
||

776 | { |
||

777 | ```
ks = const_cast<cKSplit*>(&_ks);
``` |
||

778 | grid = ks->rootgrid; |
||

779 | cellnum = _beg ? 0 : ks->num_cells-1; |
||

780 | cell = _beg ? 0 : K-1; |
||

781 | cellsize = (ks->rangemax - ks->rangemin)/K; |
||

782 | gridmin = ks->rangemin; |
||

783 | |||

784 | dive(_beg ? 0 : K-1); |
||

785 | } |
||

786 | |||

787 | void cKSplit::Iterator::dive(int where) |
||

788 | { |
||

789 | ```
// go into subgrids as deep as possible, along cells[where]
``` |
||

790 | while (ks->gridv[grid].cells[cell]<0) |
||

791 | { |
||

792 | gridmin += cell*cellsize; |
||

793 | grid = - ks->gridv[grid].cells[cell]; |
||

794 | |||

795 | cell = where; |
||

796 | cellsize /= K; |
||

797 | } |
||

798 | } |
||

799 | |||

800 | void cKSplit::Iterator::operator++(int) |
||

801 | { |
||

802 | if (grid == 0) {init(*ks, 1); return;} |
||

803 | |||

804 | cellnum++; |
||

805 | |||

806 | cell++; |
||

807 | ```
while (cell == K)
``` |
||

808 | { |
||

809 | ```
// step back to parent
``` |
||

810 | ```
int i, oldgrid = grid;
``` |
||

811 | grid = ks->gridv[grid].parent; |
||

812 | if (grid == 0) return; |
||

813 | |||

814 | ```
// find (into i) which cell of the parent 'grid' was
``` |
||

815 | for (i=0; i<K; i++) |
||

816 | ```
if (ks->gridv[grid].cells[i] == -oldgrid)
``` |
||

817 | ```
break;
``` |
||

818 | |||

819 | ```
// calc. cellsize and left edge of parent grid
``` |
||

820 | cellsize *= K; |
||

821 | gridmin -= i*cellsize; |
||

822 | |||

823 | cell = i+1; // actual '++' operation |
||

824 | } |
||

825 | ```
dive(0);
``` |
||

826 | } |
||

827 | |||

828 | void cKSplit::Iterator::operator--(int) |
||

829 | { |
||

830 | if (grid == 0) {init(*ks, 0); return;} |
||

831 | |||

832 | cellnum--; |
||

833 | |||

834 | cell--; |
||

835 | while (cell == -1) |
||

836 | { |
||

837 | ```
// step to parent
``` |
||

838 | ```
int i, oldgrid = grid;
``` |
||

839 | grid = ks->gridv[grid].parent; |
||

840 | if (grid == 0) return; |
||

841 | |||

842 | ```
// find (into i) which cell of the parent 'grid' was
``` |
||

843 | for (i=0; i<K; i++) |
||

844 | ```
if (ks->gridv[grid].cells[i] == -oldgrid)
``` |
||

845 | ```
break;
``` |
||

846 | |||

847 | ```
// calc. cellsize and left edge of parent grid
``` |
||

848 | cellsize *= K; |
||

849 | gridmin -= i*cellsize; |
||

850 | |||

851 | cell = i-1; // actual '--' operation |
||

852 | } |
||

853 | ```
dive(K-1);
``` |
||

854 | } |
||

855 | |||

856 | double cKSplit::Iterator::getCellValue() const |
||

857 | { |
||

858 | cKSplit::Grid& g = ks->gridv[grid]; |
||

859 | ```
return ks->getRealCellValue(g, cell);
``` |
||

860 | } |
||

861 | |||

862 | NAMESPACE_END |