## root / src / sim / cpsquare.cc @ 3e29b8a0

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 | } |