Statistics
| Branch: | Revision:

root / src / sim / cdetect.cc @ e26d3d25

History | View | Annotate | Download (6.91 KB)

1
//==NL=============================================================NL======
2
//  CDETECT.CC - part of
3
//
4
//                  OMNeT++/OMNEST
5
//           Discrete System Simulation in C++
6
//
7
//         File designed and written by the Hollandiaba Team
8
//
9
//   Member functions of
10
//     cTransientDetection  :  virtual base class for transient detection
11
//     cAccuracyDetection  :  virtual base class for result accuracy detection
12
//
13
//     cTDExpandingWindows    :  an algorithm for transient detection
14
//     cADByStddev    :  an algorithm for result accuracy detection
15
//
16
//   Authors: Jar Heijmans et al.
17
//=========================================================================
18

    
19
/*--------------------------------------------------------------*
20
  Copyright (C) 1992-2008 Andras Varga
21
  Copyright (C) 2006-2008 OpenSim Ltd.
22

23
  This file is distributed WITHOUT ANY WARRANTY. See the file
24
  `license' for details on this and other legal matters.
25
*--------------------------------------------------------------*/
26

    
27
#include <math.h>
28
#include <stdio.h>
29
#include "csimulation.h"  // simTime()
30
#include "cdetect.h"
31

    
32
USING_NAMESPACE
33

    
34

    
35
cTDExpandingWindows::cTDExpandingWindows(const cTDExpandingWindows& r) : cTransientDetection()
36
{
37
    func=0;
38
    setName(r.getName());
39
    operator=(r);
40
}
41

    
42
cTDExpandingWindows::cTDExpandingWindows(const char *name,
43
    int reps, int minw, double wind, double acc,
44
    PostTDFunc f, void *p) :
45
  cTransientDetection(name)
46
{
47
    pdf = f; pdfdata = p;
48
    setParameters(reps, minw, wind, acc);
49
    func = NULL; size = 0;
50
    go = true;
51
}
52

    
53
cTDExpandingWindows::~cTDExpandingWindows()
54
{
55
    reset();
56
}
57

    
58
cTDExpandingWindows& cTDExpandingWindows::operator=(const cTDExpandingWindows& res)
59
{
60
    if (this==&res) return *this;
61

    
62
    cOwnedObject::operator=(res);
63

    
64
    // setHostObject(res.getHostObject());
65
    go=res.go;
66
    transval=res.transval;
67
    accuracy=res.accuracy;
68
    minwinds=res.minwinds;
69
    windexp=res.windexp;
70
    repeats=res.repeats;
71
    detreps=res.detreps;
72
    pdf=res.pdf;
73
    pdfdata=res.pdfdata;
74

    
75
    // copy actual contents  --VA
76
    size=res.size;
77
    if (res.func==NULL)
78
    {
79
        func = NULL;
80
    }
81
    else
82
    {
83
        func = new xy;
84
        func->x = res.func->x;
85
        func->y = res.func->y;
86
        xy *res_p = res.func, *p = func;
87
        while (res_p->next)
88
        {
89
            p->next = new xy;
90
            p = p->next;
91
            res_p = res_p->next;
92
            p->x = res_p->x;
93
            p->y = res_p->y;
94
        }
95
        p->next = NULL;
96
    }
97
    return *this;
98
}
99

    
100
void cTDExpandingWindows::reset()
101
{
102
    detreps = repeats;
103
    transval = false;
104

    
105
    while (func)
106
    {
107
       xy *p = func;
108
       func = func->next;
109
       delete p;
110
    }
111
    size = 0;
112
}
113

    
114
// set transient detection parameters
115
void cTDExpandingWindows::setParameters(int reps, int minw, double wind, double acc)
116
{
117
    repeats = detreps = reps;
118
    minwinds = minw; windexp = wind;
119
    accuracy = acc;
120
}
121

    
122
// collect value
123
void cTDExpandingWindows::collect(double val)
124
{
125
  if (go)
126
  {
127
     int maxsize = minwinds + (int)(minwinds*windexp + 0.5);
128
     xy* pair = new xy;         //create new element
129
     pair->x = SIMTIME_DBL(simulation.getSimTime());
130
     pair->y = val;
131
     pair->next = NULL;
132
     if (size == 0)
133
     {                          //and add it at the end
134
        func = pair;
135
        size++;
136
     }
137
     else
138
     {
139
        xy *dum = func;
140
        while (dum->next)
141
           dum = dum->next;  // not very efficient! --VA
142
        dum->next = pair;
143
        if (size < maxsize)
144
           size++;
145
        else
146
        {                   //if size == maxsize
147
           dum = func;      //delete first element
148
           func = func->next;
149
           delete dum;
150
        }
151
     }
152
     detectTransient();  //compute the new value of transval
153
  }
154
}
155

    
156

    
157
// the actual algorithm
158
void cTDExpandingWindows::detectTransient()
159
{
160
   int ws1 = minwinds;
161
   int ws2 = (int)(0.5 + ws1 * windexp);      //round or something
162
   int maxsize = ws1 + ws2;
163
   if (size == maxsize)
164
   {
165
      xy *dum = func->next;
166
      double minval = func->y;
167
      //Get the minimum in order to move the axes
168
      while ((dum=dum->next)!=NULL)  // chged from while(dum=dum->next) --VA
169
         if (dum->y < minval) minval = dum->y;
170

    
171
      //Computation of the average of the values in the first window
172
      double avg1 = 0.0;
173
      double prevx = func->x;
174
      int i;
175
      dum = func->next;
176
      for (i=2; i<=ws1; i++)
177
      {
178
         avg1 += (dum->y - minval) * (dum->x - prevx);
179
         prevx = dum->x;
180
         dum = dum->next;
181
      }
182
      avg1 = avg1/(ws1-1);
183

    
184
      //Computation of the average of the values in the next window
185
      double avg2 = 0.0;
186
      prevx = dum->x;
187
      dum = dum->next;
188
      for (i=2; i<=ws2; i++)
189
      {
190
         avg2 += (dum->y - minval) * (dum->x - prevx);
191
         prevx = dum->x;
192
         dum = dum->next;
193
      }
194
      avg2 = avg2/(ws2-1);
195

    
196
      //Comparison of the two averages.
197
      double value = fabs(avg2!=0 ? (avg1-avg2)/avg2 : avg1);
198
      if (value < accuracy)
199
          detreps--;
200
      else
201
          detreps = repeats;
202
      transval = (detreps <= 0);
203

    
204
      // If the transient period over, call registered function
205
      if (transval && (pdf!=NULL))
206
      {
207
         (*pdf)(this,pdfdata);
208
      }
209
  }
210
}
211

    
212
//----
213

    
214
cADByStddev::cADByStddev(const cADByStddev& r) : cAccuracyDetection()
215
{
216
    go=resaccval=false;
217
    detreps=sctr=0;
218
    ssum=sqrsum=0.0;
219

    
220
    setName(r.getName());
221
    operator=(r);
222
}
223

    
224
cADByStddev::cADByStddev(const char *name,
225
                         double acc, int reps, PostADFunc f, void *p) :
226
  cAccuracyDetection(name)
227
{
228
   pdf = f; pdfdata = p;
229
   accuracy = acc;
230
   repeats=detreps=reps;
231
   go=resaccval=false;
232
   sctr=0;
233
   ssum=sqrsum=0.0;
234
}
235

    
236
cADByStddev& cADByStddev::operator=(const cADByStddev& res)
237
{
238
   if (this==&res) return *this;
239

    
240
   cOwnedObject::operator=(res);
241
   // setHostObject(res.getHostObject());
242
   go=res.go; resaccval=res.resaccval;
243
   accuracy=res.accuracy; sctr=res.sctr;
244
   ssum=res.ssum; sqrsum=sqrsum;
245
   repeats=res.repeats; detreps=res.detreps;
246
   pdf=res.pdf; pdfdata=res.pdfdata;
247
   return *this;
248
}
249

    
250
double cADByStddev::getStddev() const
251
{
252
   if (!sctr)
253
      return 0.0;
254
   else
255
   {
256
      double devsqr = (sqrsum - ssum*ssum/sctr)/sctr;
257
      if (devsqr<=0)
258
          return 0.0;
259
      else
260
          return sqrt(devsqr);
261
   }
262
}
263

    
264
void cADByStddev::reset()
265
{
266
   ssum = sqrsum = 0.0; sctr = 0;
267
   resaccval = false;
268
}
269

    
270
// collect a value
271
void cADByStddev::collect(double val)
272
{
273
   if (go)
274
   {
275
      ssum+=val;  sqrsum+=val*val; sctr++;
276
      if (sctr>6) detectAccuracy();
277
   }
278
}
279

    
280

    
281
void cADByStddev::detectAccuracy()
282
{
283
   // the actual algorithm: divide the standard deviation by the square of
284
   // the number of values and check if this is small enough
285

    
286
   double sdvmean = (getStddev()/(sctr*sctr));
287
   if (sdvmean <= accuracy) detreps--;
288
   else detreps = repeats;
289
   resaccval = (detreps <= 0);
290

    
291
   // If we have the accuracy, call registered function
292
   if (resaccval && pdf!=NULL)
293
   {
294
      (*pdf)(this, pdfdata);
295
   }
296
}
297