Statistics
| Branch: | Revision:

root / src / sim / cpsquare.cc @ e1750c09

History | View | Annotate | Download (9.89 KB)

1 01873262 Georg Kunz
//=========================================================================
2
//  CPSQUARE.CC - part of
3
//
4
//                  OMNeT++/OMNEST
5
//           Discrete System Simulation in C++
6
//
7
//   Member functions of
8
//     cPSquare: calculates quantile values without storing the observations
9
//
10
//   Written by Babak Fakhamzadeh, TU Delft, Dec 1996
11
//   Modifications by Andras Varga
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 "cpsquare.h"
27
#include "cexception.h"
28
#include "random.h"
29
#include "distrib.h"
30
31
#ifdef WITH_PARSIM
32
#include "ccommbuffer.h"
33
#endif
34
35
USING_NAMESPACE
36
37
using std::ostream;
38
using std::endl;
39
40
41
Register_Class(cPSquare);
42
43
44
cPSquare::cPSquare(const cPSquare& r) : cDensityEstBase()
45
{
46
    setName( r.getName() );
47
    operator=(r);
48
}
49
50
51
cPSquare::cPSquare (const char *name, int b) : cDensityEstBase(name)
52
{
53
    transfd = true;
54
55
    numcells=b;
56
    numobs=0;
57
    n = new int[numcells+2];
58
    q = new double[numcells+2];
59
60
    for (int i=0; i<=numcells+1; i++)
61
    {
62
      n[i]=i;
63
      q[i]=-1e50;       //this should be -(max(double))
64
    }
65
}
66
67
cPSquare::~cPSquare()
68
{
69
    delete [] q;
70
    delete [] n;
71
}
72
73
74
void cPSquare::parsimPack(cCommBuffer *buffer)
75
{
76
#ifndef WITH_PARSIM
77
    throw cRuntimeError(this,eNOPARSIM);
78
#else
79
    cDensityEstBase::parsimPack(buffer);
80
81
    buffer->pack(numcells);
82
    buffer->pack(numobs);
83
84
    if (buffer->packFlag(n!=NULL))
85
        buffer->pack(n, numcells + 2);
86
    if (buffer->packFlag(q!=NULL))
87
        buffer->pack(q, numcells + 2);
88
#endif
89
}
90
91
void cPSquare::parsimUnpack(cCommBuffer *buffer)
92
{
93
#ifndef WITH_PARSIM
94
    throw cRuntimeError(this,eNOPARSIM);
95
#else
96
    cDensityEstBase::parsimUnpack(buffer);
97
98
    buffer->unpack(numcells);
99
    buffer->unpack(numobs);
100
101
    if (buffer->checkFlag())
102
    {
103
        n = new int[numcells + 2];
104
        buffer->unpack(n, numcells + 2);
105
    }
106
107
    if (buffer->checkFlag())
108
    {
109
        q = new double[numcells + 2];
110
        buffer->unpack(q, numcells + 2);
111
    }
112
#endif
113
}
114
115
cPSquare& cPSquare::operator=(const cPSquare& res)
116
{
117
    if (this==&res) return *this;
118
119
    cDensityEstBase::operator=(res);
120
121
    numobs=res.numobs;
122
    numcells=res.numcells;
123
    delete [] n;
124
    delete [] q;
125
    n=new int[numcells+2];
126
    q=new double[numcells+2];
127
    for (int i=0; i<=numcells+1; i++)
128
    {
129
      n[i]=res.n[i];
130
      q[i]=res.q[i];
131
    }
132
    return (*this);
133
}
134
135
void cPSquare::giveError()
136
{
137
    throw cRuntimeError(this, "setRange..() and setNumFirstVals() makes no sense with cPSquare");
138
}
139
140
#ifdef NEW_CODE_FROM_PUPPIS_WITH_FLOATING_POINT_EXCEPTION
141
void cPSquare::collectTransformed(double val)
142
{
143
    int i, j;
144
145
    numobs++;
146
147
    if (numobs <= numcells + 1)
148
    {
149
       q[numcells+2-numobs] = val;
150
    }
151
    else
152
    {
153
      // now numobs==b+1, number of observations is enough to produce 'b' cells,
154
      // estimation has to be done
155
      if (numobs == numcells+2) {
156
         for (i=1; i<numcells+1; i++) {
157
            for (j=i+1; j<numcells+2; j++) {
158
               if (q[j] < q[i])
159
               {
160
                  double temp = q[i];
161
                  q[i] = q[j];
162
                  q[j] = temp;
163
               }
164
            }
165
         }
166
      }
167
168
      int k = 0;        //the cell number in which 'val' falls
169
170
      for (i=1; i<=numcells+1; i++)
171
      {
172
         if (val <= q[i]) {
173
            if (i==1)
174
            {
175
               q[1] = val;
176
               k = 1;
177
            }
178
            else {
179
               k = i-1;
180
            }
181
            break;
182
         }
183
      }
184
185
      if (k==0) //the new value falls outside of the current cells
186
      {
187
         q[numcells+1] = val;
188
         k = numcells;
189
      }
190
191
      for (i=k+1; i<=numcells+1; i++)
192
         n[i] = n[i]+1;
193
194
      double d, qd;
195
      int e;
196
      for (i=2; i<=numcells; i++)
197
      {
198
         d = 1 + (i - 1) * (numobs - 1) / ((double)numcells) - n[i];
199
200
         if ((d>=1 && n[i+1]-n[i]>1) || (d<=-1 && n[i-1]-n[i]<-1))
201
         //if it is possible to adjust the marker position
202
         {
203
            d < 0 ? e = -1 : e = 1;
204
205
            //try the parabolic formula
206
            qd = q[i] + e / ((double)(n[i+1] - n[i-1])) * ((n[i] - n[i-1] + e)
207
                 * (q[i+1] - q[i]) / ((double)(n[i+1] - n[i])) + (n[i+1]
208
                 - n[i] - e) * (q[i] - q[i-1]) / ((double)(n[i] - n[i-1])));
209
210
            if ((qd>q[i-1]) && (qd<q[i+1]))       // if it is possible to fit the new height
211
               q[i] = qd;                         // then do so (in between the neigbouring heights)
212
            else                                  // else use the linear formula
213
               q[i] += e * (q[i+e] - q[i]) / (double)(n[i+e] - n[i]);
214
            n[i] += e;
215
         }
216
      }
217
    }
218
}
219
#else /* OLDER_WORKING_CODE */
220
void cPSquare::collectTransformed(double val)
221
{
222
    int i;
223
224
    numobs++;          //an extra observation is added
225
226
    if (numobs <= numcells + 1)
227
    {
228
      // old code:
229
      //q[numcells+2-numobs] = val;
230
      // places val in front, because qsort puts everything at the end,
231
      // because initialized with q[i]=-max
232
      //qsort(q, numcells+2, sizeof(*q), CompDouble);
233
234
      // insert val into q[]; q[0] is not used!
235
      for(i=numobs; i>=2 && val<=q[i-1]; i--)
236
         q[i]=q[i-1];
237
      q[i]=val;
238
    }
239
    else
240
    // now numobs==b+1, number of observations is enough to produce 'b' cells,
241
    // estimation has to be done
242
    { int k = 0;        //the cell number in which 'val' falls
243
244
      for (i=1; i<=numcells+1; i++)
245
      { if (val <= q[i])
246
        { if (i==1)
247
          {  q [1] = val;
248
             k = 1;
249
          }
250
          else
251
          { k = i - 1;}
252
          break;
253
        }
254
      }
255
      if (k==0) //the new value falls outside of the current cells
256
      { q[numcells+1]=val;
257
        k = numcells;
258
      }
259
      for (i=k+1; i<=numcells+1; i++)
260
      { n[i]=n[i]+1;}
261
262
      double d, qd;
263
      int e;
264
      for (i=2; i<=numcells; i++)
265
      { d = 1 + (i - 1) * (numobs - 1) / ((double)numcells) - n[i];
266
267
        if ((d>=1 && n[i+1]-n[i]>1) || (d<=-1 && n[i-1]-n[i]<-1))
268
        //if it is possible to adjust the marker position
269
        { if (d < 0)
270
          { e = - 1;}
271
          else
272
          { e = 1;}
273
          //try the parabolic formula
274
          qd = q[i] + e / ((double)(n [i + 1] - n [i - 1])) * ((n [i] - n [i - 1] + e)
275
               * (q [i + 1] - q [i]) / ((double)(n [i + 1] - n [i])) + (n [i + 1]
276
               - n [i] - e) * (q [i] - q [i - 1]) / ((double)(n [i] - n [i - 1])));
277
278
          if ((qd>q[i-1]) && (qd<q[i+1]))       //if it is possible to fit the new height
279
280
          { q [i] = qd;}                        //then do so (in between the neigbouring heights)
281
          else                                  //else
282
          { q [i] += e * (q [i + e] - q [i])    //use the linear formula
283
                     / ((double)(n [i + e]
284
                     - n [i]));
285
          }
286
          n [i] += e;
287
        }
288
      }
289
    }
290
}
291
#endif
292
293
void cPSquare::merge(const cStatistic *other)
294
{
295
    throw cRuntimeError(this, "The cPSquare class does not support merge()");
296
}
297
298
void cPSquare::doMergeCellValues(const cDensityEstBase *other)
299
{
300
    ASSERT(false); // never comes here, as merge() already throws an error
301
}
302
303
double cPSquare::random() const
304
{
305
    double s;
306
    int k = 0;
307
    int l = 0;
308
309
    //if (num_obs==0)   // newer, from PUPPIS
310
    if (numobs<numcells+1)
311
        throw cRuntimeError(this,"must collect at least num_cells values before random() can be used");
312
313
    s = numobs * genk_dblrand(genk);
314
315
    for (int i=1; i<=numcells+1; i++)
316
    {
317
       if (s < n[i])
318
       {
319
          k=i;
320
          l=k-1;
321
          break;
322
       }
323
    }
324
325
    if (k==1)
326
       l=k;
327
328
    if (numobs<numcells+1)
329
    {
330
       k += numcells-numobs+1;
331
       l += numcells-numobs+1;
332
    }
333
334
    return dblrand()*(q[k]-q[l])+q[l];
335
}
336
337
int cPSquare::getNumCells() const
338
{
339
    if (numobs<2)
340
       return 0;
341
    else if (numobs<numcells)
342
       return numobs-1;
343
    else
344
       return numcells;
345
}
346
347
double cPSquare::getBasepoint(int k) const
348
{
349
    return q[k+1];
350
}
351
352
double cPSquare::getCellValue(int k) const
353
{
354
    return n[k+2] - n[k+1] + (k==0);
355
}
356
357
std::string cPSquare::detailedInfo() const
358
{
359
    std::stringstream os;
360
    os << cDensityEstBase::detailedInfo();
361
362
    int nn = numobs<=numcells+1 ? numobs : numcells+1;
363
364
    os << "\n  The quantiles (#(observations: observation<=marker)):\n";
365
    os << "      #observations\t<=marker\n";
366
    for (int i=1; i<=nn; i++)
367
       os << "       " << n[i] << "\t " << q[i] << endl;
368
    return os.str();
369
}
370
371
double cPSquare::getCDF(double x) const
372
{
373
   // returns 0..1; uses linear approximation between two markers
374
   for (int i=1; i<numcells+2 ; i++)
375
      if (q[i]>x)
376
         return ((x-q[i-1]) / (q[i]-q[i-1]) * (n[i]-n[i-1] + (i==1)) + n[i-1] + (i==1)) / numobs;
377
   return 1.0;
378
}
379
380
double cPSquare::getPDF(double x) const
381
{
382
   // returns 0..1; assumes constant PDF within a cell
383
   for (int i=1 ; i<numcells+2 ; i++)
384
      if (q[i]>x)
385
        return (n[i]-n[i-1] + (i==1))/(q[i]-q[i-1])/numobs;
386
   return 0;
387
}
388
389
void cPSquare::saveToFile(FILE *f) const
390
{
391
   cDensityEstBase::saveToFile(f);
392
393
   fprintf(f,"%u\t #= numcells\n",numcells);
394
   fprintf(f,"%ld\t #= numobs\n",numobs);
395
396
   int i;
397
   fprintf(f,"#= n[]\n");
398
   for (i=0; i<numcells+2; i++)  fprintf(f," %d\n",n[i]);
399
400
   fprintf(f,"#= q[]\n");
401
   for (i=0; i<numcells+2; i++)  fprintf(f," %g\n",q[i]);
402
}
403
404
void cPSquare::loadFromFile(FILE *f)
405
{
406
   cDensityEstBase::loadFromFile(f);
407
408
   freadvarsf(f,"%u\t #= numcells",&numcells);
409
   freadvarsf(f,"%ld\t #= numobs",&numobs);
410
411
   int i;
412
   freadvarsf(f,"#= n[]");
413
   for (i=0; i<numcells+2; i++)  freadvarsf(f," %d",n+i);
414
415
   freadvarsf(f,"#= q[]");
416
   for (i=0; i<numcells+2; i++)  freadvarsf(f," %g",q+i);
417
}