root / src / sim / cksplit.cc @ 3e29b8a0
History  View  Annotate  Download (20.5 KB)
1 
//=========================================================================


2 
// CKSPLIT.CC  part of

3 
//

4 
// OMNeT++/OMNEST

5 
// Discrete System Simulation in C++

6 
//

7 
// Member functions of

8 
// cKSplit : implements the ksplit algorithm in 1 dimension

9 
//

10 
// Original version by Babak Fakhamzadeh, TU Delft, MarJun 1996

11 
// This version written by Andras Varga, 199798

12 
//

13 
//=========================================================================

14 
/**

15 
Copyright (C) 19922008 Andras Varga

16 
Copyright (C) 20062008 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.totalg.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) ((xgridmin) / cellsize);

281  
282 
// guard against rounding errors

283 
if (i<0) i=0; 
284 
if (i>=K) i=K1; 
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) ((xgridmin) / cellsize);

315  
316 
if (i<0) i=0; 
317 
if (i>=K) i=K1; 
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 += K1;

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 += K1;

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[(K1)/2] = old_rootgrid; 
404  
405 
rangemin = (K1)/2*gridsize; 
406 
rangemax += (K1)/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 * (xrangemin) / (rangemaxrangemin)); 
417  
418 
// if it falls in the old root grid, we have to correct it

419 
if (i == (K1)/2) 
420 
{ 
421 
if (x > (rangemaxrangemin)/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.totalgrid.mother);

479  
480 
return cell_mot + (1lambda)*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_cells1; 
780 
cell = _beg ? 0 : K1; 
781 
cellsize = (ks>rangemax  ks>rangemin)/K; 
782 
gridmin = ks>rangemin; 
783  
784 
dive(_beg ? 0 : K1); 
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 = i1; // actual '' operation 
852 
} 
853 
dive(K1);

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 