Statistics
| Branch: | Revision:

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

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