« back to all changes in this revision
Viewing changes to src/sdengine.cpp
- browse files at revision 12
- compare with another revision
- download diff
- download tarball
- view history from revision 12
- Committer: Package Import Robot
- Author(s): Daniel Leidert
- Date: 2014-01-18 20:07:16 UTC
- mfrom: (1.1.8)
- Revision ID: package-import@ubuntu.com-20140118200716-whsmcg7fa1eyqecq
Tags: 140117-1
New upstream release.
- files added:
- src/plproj.cpp
- files modified:
- Makefile
added
removed
src/sdengine.cpp
1
/*****************************************************************************
2
TRAVIS - Trajectory Analyzer and Visualizer
3
http://www.travis-analyzer.de/
4
5
Copyright (c) 2009-2013 Martin Brehm
6
2012-2013 Martin Thomas
7
8
This file written by Martin Brehm.
9
10
This program is free software: you can redistribute it and/or modify
11
it under the terms of the GNU General Public License as published by
12
the Free Software Foundation, either version 3 of the License, or
13
(at your option) any later version.
14
15
This program is distributed in the hope that it will be useful,
16
but WITHOUT ANY WARRANTY; without even the implied warranty of
17
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18
GNU General Public License for more details.
19
20
You should have received a copy of the GNU General Public License
21
along with this program. If not, see <http://www.gnu.org/licenses/>.
22
*****************************************************************************/
23
24
25
#include "sdengine.h"
26
#include "lmwrapper.h"
27
#include "tools.h"
28
29
30
CSDEngine::CSDEngine()
31
{
32
m_oaGroups.SetName("CSDEngine::m_oaGroups");
33
m_iaBlockSizes.SetName("CSDEngine::m_iaBlockSizes");
34
m_faTempBuffer.SetName("CSDEngine::m_faTempBuffer");
35
m_faSum.SetName("CSDEngine::m_faSum");
36
m_faSqSum.SetName("CSDEngine::m_faSqSum");
37
m_faMean.SetName("CSDEngine::m_faMean");
38
m_faVari.SetName("CSDEngine::m_faVari");
39
m_faCorLen.SetName("CSDEngine::m_faCorLen");
40
}
41
42
43
CSDEngine::~CSDEngine()
44
{
45
}
46
47
48
CSDGroup::CSDGroup()
49
{
50
m_oaBlockSets.SetName("CSDGroup::m_oaBlockSets");
51
}
52
53
54
CSDGroup::~CSDGroup()
55
{
56
}
57
58
59
CSDBlockSet::CSDBlockSet()
60
{
61
m_iCounter = 0;
62
m_faBlocks.SetName("CSDBlockSet::m_faBlocks");
63
}
64
65
66
CSDBlockSet::~CSDBlockSet()
67
{
68
}
69
70
71
void CSDEngine::Init(int bins)
72
{
73
int z, z2;
74
CSDGroup *gr;
75
CSDBlockSet *bs;
76
77
m_iSteps = 0;
78
79
m_faTempBuffer.SetSize(bins);
80
m_faSum.SetSize(bins);
81
m_faSqSum.SetSize(bins);
82
m_faMean.SetSize(bins);
83
m_faVari.SetSize(bins);
84
m_faCorLen.SetSize(bins);
85
86
for (z=0;z<bins;z++)
87
{
88
m_faTempBuffer[z] = 0;
89
m_faSum[z] = 0;
90
m_faSqSum[z] = 0;
91
m_faMean[z] = 0;
92
m_faVari[z] = 0;
93
m_faCorLen[z] = 0;
94
95
gr = new CSDGroup();
96
m_oaGroups.Add(gr);
97
gr->m_oaBlockSets.SetSize(m_iaBlockSizes.GetSize());
98
gr->m_faBlockVari.SetSize(m_iaBlockSizes.GetSize());
99
for (z2=0;z2<m_iaBlockSizes.GetSize();z2++)
100
{
101
bs = new CSDBlockSet();
102
gr->m_oaBlockSets[z2] = bs;
103
gr->m_faBlockVari[z2] = 0;
104
105
bs->m_faBlocks.Add(0);
106
}
107
}
108
}
109
110
111
void CSDEngine::FinishStep()
112
{
113
int z, z2;
114
CSDGroup *gr;
115
CSDBlockSet *bs;
116
117
m_iSteps++;
118
119
for (z=0;z<m_oaGroups.GetSize();z++)
120
{
121
gr = (CSDGroup*)m_oaGroups[z];
122
123
m_faSum[z] += m_faTempBuffer[z];
124
m_faSqSum[z] += m_faTempBuffer[z]*m_faTempBuffer[z];
125
126
for (z2=0;z2<m_iaBlockSizes.GetSize();z2++)
127
{
128
bs = (CSDBlockSet*)gr->m_oaBlockSets[z2];
129
bs->m_faBlocks[bs->m_faBlocks.GetSize()-1] += m_faTempBuffer[z];
130
131
bs->m_iCounter++;
132
if (bs->m_iCounter >= m_iaBlockSizes[z2])
133
{
134
bs->m_iCounter = 0;
135
bs->m_faBlocks[bs->m_faBlocks.GetSize()-1] /= m_iaBlockSizes[z2];
136
bs->m_faBlocks.Add(0);
137
}
138
}
139
140
m_faTempBuffer[z] = 0;
141
}
142
}
143
144
145
void CSDEngine::DumpData(CDF *df, float timessigma, bool verbose)
146
{
147
int z, z2, ti;
148
CSDGroup *gr;
149
// FILE *a;
150
FILE *b;
151
char buf[256];
152
double *x;
153
double *y, *y2, ya;
154
double da, ds, dr, tf, dsa, dsmi, dsma;
155
CLMWrapper *lm;
156
double par[3], *curve;
157
158
// a = fopen("block.csv","wt");
159
b = NULL;
160
161
mprintf(" Calculating correlation time of histogram...\n");
162
163
try { df->m_pAdditionalSets = new double*[4]; } catch(...) { df->m_pAdditionalSets = NULL; }
164
if (df->m_pAdditionalSets == NULL) NewException((double)4*sizeof(double*),__FILE__,__LINE__,__PRETTY_FUNCTION__);
165
166
try { df->m_pAdditionalSetLabels = new char*[4]; } catch(...) { df->m_pAdditionalSetLabels = NULL; }
167
if (df->m_pAdditionalSetLabels == NULL) NewException((double)4*sizeof(char*),__FILE__,__LINE__,__PRETTY_FUNCTION__);
168
169
df->m_iAdditionalSets = 4;
170
171
for (z=0;z<4;z++)
172
{
173
try { df->m_pAdditionalSets[z] = new double[df->m_iResolution]; } catch(...) { df->m_pAdditionalSets[z] = NULL; }
174
if (df->m_pAdditionalSets[z] == NULL) NewException((double)df->m_iResolution*sizeof(double),__FILE__,__LINE__,__PRETTY_FUNCTION__);
175
}
176
177
sprintf(buf,"%.2f%c confidence range",NormalDistIntegral(timessigma)*2.0-1.0,'%');
178
try { df->m_pAdditionalSetLabels[0] = new char[strlen(buf)+1]; } catch(...) { df->m_pAdditionalSetLabels[0] = NULL; }
179
if (df->m_pAdditionalSetLabels[0] == NULL) NewException((double)(strlen(buf)+1)*sizeof(char),__FILE__,__LINE__,__PRETTY_FUNCTION__);
180
strcpy(df->m_pAdditionalSetLabels[0],buf);
181
182
sprintf(buf,"%.2f%c confidence range",NormalDistIntegral(timessigma)*2.0-1.0,'%');
183
try { df->m_pAdditionalSetLabels[1] = new char[strlen(buf)+1]; } catch(...) { df->m_pAdditionalSetLabels[1] = NULL; }
184
if (df->m_pAdditionalSetLabels[1] == NULL) NewException((double)(strlen(buf)+1)*sizeof(char),__FILE__,__LINE__,__PRETTY_FUNCTION__);
185
strcpy(df->m_pAdditionalSetLabels[1],buf);
186
187
sprintf(buf,"Standard deviation");
188
try { df->m_pAdditionalSetLabels[2] = new char[strlen(buf)+1]; } catch(...) { df->m_pAdditionalSetLabels[2] = NULL; }
189
if (df->m_pAdditionalSetLabels[2] == NULL) NewException((double)(strlen(buf)+1)*sizeof(char),__FILE__,__LINE__,__PRETTY_FUNCTION__);
190
strcpy(df->m_pAdditionalSetLabels[2],buf);
191
192
sprintf(buf,"Correlation Length");
193
try { df->m_pAdditionalSetLabels[3] = new char[strlen(buf)+1]; } catch(...) { df->m_pAdditionalSetLabels[4] = NULL; }
194
if (df->m_pAdditionalSetLabels[3] == NULL) NewException((double)(strlen(buf)+1)*sizeof(char),__FILE__,__LINE__,__PRETTY_FUNCTION__);
195
strcpy(df->m_pAdditionalSetLabels[3],buf);
196
197
try { lm = new CLMWrapper(); } catch(...) { lm = NULL; }
198
if (lm == NULL) NewException((double)sizeof(CLMWrapper),__FILE__,__LINE__,__PRETTY_FUNCTION__);
199
200
try { x = new double[m_iaBlockSizes.GetSize()]; } catch(...) { x = NULL; }
201
if (x == NULL) NewException((double)m_iaBlockSizes.GetSize()*sizeof(double),__FILE__,__LINE__,__PRETTY_FUNCTION__);
202
203
try { y = new double[m_iaBlockSizes.GetSize()]; } catch(...) { y = NULL; }
204
if (y== NULL) NewException((double)m_iaBlockSizes.GetSize()*sizeof(double),__FILE__,__LINE__,__PRETTY_FUNCTION__);
205
206
try { y2 = new double[m_iaBlockSizes.GetSize()]; } catch(...) { y2 = NULL; }
207
if (y2 == NULL) NewException((double)m_iaBlockSizes.GetSize()*sizeof(double),__FILE__,__LINE__,__PRETTY_FUNCTION__);
208
209
try { curve = new double[m_iaBlockSizes.GetSize()]; } catch(...) { curve = NULL; }
210
if (curve == NULL) NewException((double)m_iaBlockSizes.GetSize()*sizeof(double),__FILE__,__LINE__,__PRETTY_FUNCTION__);
211
212
mprintf(" Fitting trajectory block decomposition...\n");
213
mprintf(WHITE," [");
214
215
tf = m_oaGroups.GetSize() / 60.0;
216
217
dsa = 0;
218
dsmi = 1e20;
219
dsma = 0;
220
ti = 0;
221
for (z=0;z<m_oaGroups.GetSize();z++)
222
{
223
if (fmod(z,tf) < 1.0)
224
mprintf(WHITE,"#");
225
226
if (df->m_pBin[z] == 0)
227
continue;
228
229
gr = (CSDGroup*)m_oaGroups[z];
230
231
// mprintf(" * Gruppe %d:\n",z);
232
233
if (verbose)
234
{
235
sprintf(buf,"extra%03d.csv",z);
236
b = fopen(buf,"wt");
237
fprintf(b,"# Block; Block size; sqrt(size); log(size); 1/size; Vari; Size*Vari; S; log(S); log(S2); Fit\n");
238
}
239
240
ya = 0;
241
for (z2=0;z2<m_iaBlockSizes.GetSize();z2++)
242
{
243
x[z2] = log10(m_iaBlockSizes[z2]);
244
y[z2] = log10(m_iaBlockSizes[z2]*gr->m_faBlockVari[z2]/m_faVari[z]);
245
if (z2 == 0)
246
{
247
y2[z2] = y[z2];
248
ya = y[z2];
249
} else
250
{
251
if (y[z2] > ya)
252
{
253
y2[z2] = y[z2];
254
ya = y[z2];
255
} else y2[z2] = ya;
256
}
257
}
258
259
par[0] = 50.0;
260
par[1] = -0.01;
261
par[2] = 1.0;
262
263
lm->Fit_SD(m_iaBlockSizes.GetSize(),x,y2,par,curve,&dr,200);
264
265
if (verbose)
266
fprintf(b,"# A=%G; B=%G; C=%G; R=%G\n",par[0],par[1],par[2],dr);
267
268
ds = pow(10.0,par[0]*(1.0-exp(par[1]*pow(log10(m_iSteps),par[2]))));
269
270
if (ds < dsmi)
271
dsmi = ds;
272
if (ds > dsma)
273
dsma = ds;
274
dsa += ds;
275
ti++;
276
277
m_faCorLen[z] = ds;
278
279
// mprintf(" Fit done. y(x) = %10G * (1 - exp( %10G * x^%10G ) ). R = %10G. S = %10G\n",par[0],par[1],par[2],dr,ds);
280
281
if (verbose)
282
{
283
for (z2=0;z2<m_iaBlockSizes.GetSize();z2++)
284
fprintf(b,"%d; %d; %G; %G; %G; %G; %G; %G; %G; %G; %G\n",z2,(int)m_iaBlockSizes[z2],sqrt(m_iaBlockSizes[z2]),log10(m_iaBlockSizes[z2]),1.0/m_iaBlockSizes[z2],gr->m_faBlockVari[z2],m_iaBlockSizes[z2]*gr->m_faBlockVari[z2],m_iaBlockSizes[z2]*gr->m_faBlockVari[z2]/m_faVari[z],y[z2],y2[z2],curve[z2]);
285
}
286
287
/* x = new double[m_iaBlockSizes.GetSize()];
288
y = new double[m_iaBlockSizes.GetSize()];
289
290
xa = 0;
291
ya = 0;
292
for (z2=0;z2<m_iaBlockSizes.GetSize();z2++)
293
{
294
bs = (CSDBlockSet*)gr->m_oaBlockSets[z2];
295
// mprintf(" - Blockgroesse %d (%d Schritte)\n",z2,m_iaBlockSizes[z2]);
296
// mprintf(" %d Bloecke, Counter %d.\n",bs->m_faBlocks.GetSize(),bs->m_iCounter);
297
298
if (z == 100)
299
{
300
for (z3=0;z3<bs->m_faBlocks.GetSize()-1;z3++)
301
{
302
fprintf(a,"%d; %d; %G\n",z2,z3,bs->m_faBlocks[z3]);
303
}
304
}
305
306
307
x[z2] = log10(m_iaBlockSizes[z2]);
308
y[z2] = log10(m_iaBlockSizes[z2]*gr->m_faBlockVari[z2]/m_faVari[z]);
309
xa += x[z2];
310
ya += y[z2];
311
}
312
313
xa /= m_iaBlockSizes.GetSize();
314
ya /= m_iaBlockSizes.GetSize();
315
316
dxy = 0;
317
dxx = 0;
318
dyy = 0;
319
320
for (z2=0;z2<m_iaBlockSizes.GetSize();z2++)
321
{
322
dxy += (x[z2] - xa) * (y[z2] - ya);
323
dxx += (x[z2] - xa) * (x[z2] - xa);
324
dyy += (y[z2] - ya) * (y[z2] - ya);
325
}
326
327
db = dxy / dxx;
328
da = ya - db * xa;
329
330
ds = pow(10.0, da + db * log10(m_iSteps));
331
332
dr = sqrt( dxy*dxy / dxx / dyy );
333
334
mprintf(" y(x) = %10G + %10G * x; S = %10G; R = %10G\n",da,db,ds,dr);
335
336
for (z2=0;z2<m_iaBlockSizes.GetSize();z2++)
337
fprintf(b,"%d; %d; %G; %G; %G; %G; %G; %G; %G; %G\n",z2,(int)m_iaBlockSizes[z2],sqrt(m_iaBlockSizes[z2]),log10(m_iaBlockSizes[z2]),1.0/m_iaBlockSizes[z2],gr->m_faBlockVari[z2],m_iaBlockSizes[z2]*gr->m_faBlockVari[z2],m_iaBlockSizes[z2]*gr->m_faBlockVari[z2]/m_faVari[z],log10(m_iaBlockSizes[z2]*gr->m_faBlockVari[z2]/m_faVari[z]),da+db*x[z2]);
338
339
340
delete[] x;
341
delete[] y;*/
342
343
if (verbose)
344
fclose(b);
345
346
}
347
delete[] x;
348
delete[] y;
349
delete[] y2;
350
delete[] curve;
351
delete lm;
352
353
mprintf(WHITE,"]\n");
354
355
dsa /= ti;
356
357
mprintf(" Correlation length: Min %.3f, Max %.3f, Avg %.3f time steps.\n",dsmi,dsma,dsa);
358
359
for (z=0;z<m_oaGroups.GetSize();z++)
360
{
361
ds = 0;
362
da = 0;
363
for (z2=z-10;z2<=z+10;z2++)
364
{
365
if (z2 < 0)
366
continue;
367
if (z2 >= m_oaGroups.GetSize())
368
continue;
369
if (df->m_pBin[z2] == 0)
370
continue;
371
ds += m_faCorLen[z2];
372
da += 1.0;
373
}
374
ds /= da;
375
376
if (df->m_pBin[z] != 0)
377
{
378
// df->m_pAdditionalSets[2][z] = sqrt(ds / m_iSteps) * df->m_pBin[z];
379
df->m_pAdditionalSets[2][z] = sqrt(ds / m_iSteps * m_faVari[z]);
380
df->m_pAdditionalSets[3][z] = ds;
381
382
df->m_pAdditionalSets[0][z] = df->m_pBin[z] - timessigma*df->m_pAdditionalSets[2][z];
383
df->m_pAdditionalSets[1][z] = df->m_pBin[z] + timessigma*df->m_pAdditionalSets[2][z];
384
} else
385
{
386
df->m_pAdditionalSets[0][z] = 0;
387
df->m_pAdditionalSets[1][z] = 0;
388
df->m_pAdditionalSets[2][z] = 0;
389
df->m_pAdditionalSets[3][z] = 0;
390
}
391
}
392
393
// fclose(a);
394
}
395
396
397
void CSDEngine::FinishAnalysis()
398
{
399
int z, z2, z3;
400
CSDGroup *gr;
401
CSDBlockSet *bl;
402
403
mprintf(" Calculating block averages...\n");
404
405
for (z=0;z<m_oaGroups.GetSize();z++)
406
{
407
m_faMean[z] = m_faSum[z] / m_iSteps;
408
m_faVari[z] = (m_faSqSum[z] / m_iSteps) - (m_faMean[z]*m_faMean[z]);
409
410
// mprintf(" # Bin %3d: Sum %10G; SqSum %10G; Mean %10G; Vari %10G\n",z,m_faSum[z],m_faSqSum[z],m_faMean[z],m_faVari[z]);
411
412
gr = (CSDGroup*)m_oaGroups[z];
413
414
for (z2=0;z2<gr->m_oaBlockSets.GetSize();z2++)
415
{
416
bl = (CSDBlockSet*)gr->m_oaBlockSets[z2];
417
418
for (z3=0;z3<bl->m_faBlocks.GetSize()-1;z3++)
419
gr->m_faBlockVari[z2] += pow(bl->m_faBlocks[z3] - m_faMean[z], 2);
420
421
gr->m_faBlockVari[z2] /= ((double)bl->m_faBlocks.GetSize()-1);
422
}
423
}
424
}
425
426
1
/*****************************************************************************
2
TRAVIS - Trajectory Analyzer and Visualizer
3
http://www.travis-analyzer.de/
4
5
Copyright (c) 2009-2014 Martin Brehm
6
2012-2014 Martin Thomas
7
8
This file written by Martin Brehm.
9
10
This program is free software: you can redistribute it and/or modify
11
it under the terms of the GNU General Public License as published by
12
the Free Software Foundation, either version 3 of the License, or
13
(at your option) any later version.
14
15
This program is distributed in the hope that it will be useful,
16
but WITHOUT ANY WARRANTY; without even the implied warranty of
17
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18
GNU General Public License for more details.
19
20
You should have received a copy of the GNU General Public License
21
along with this program. If not, see <http://www.gnu.org/licenses/>.
22
*****************************************************************************/
23
24
25
#include "sdengine.h"
26
#include "lmwrapper.h"
27
#include "tools.h"
28
29
30
CSDEngine::CSDEngine()
31
{
32
m_oaGroups.SetName("CSDEngine::m_oaGroups");
33
m_iaBlockSizes.SetName("CSDEngine::m_iaBlockSizes");
34
m_faTempBuffer.SetName("CSDEngine::m_faTempBuffer");
35
m_faSum.SetName("CSDEngine::m_faSum");
36
m_faSqSum.SetName("CSDEngine::m_faSqSum");
37
m_faMean.SetName("CSDEngine::m_faMean");
38
m_faVari.SetName("CSDEngine::m_faVari");
39
m_faCorLen.SetName("CSDEngine::m_faCorLen");
40
}
41
42
43
CSDEngine::~CSDEngine()
44
{
45
}
46
47
48
CSDGroup::CSDGroup()
49
{
50
m_oaBlockSets.SetName("CSDGroup::m_oaBlockSets");
51
}
52
53
54
CSDGroup::~CSDGroup()
55
{
56
}
57
58
59
CSDBlockSet::CSDBlockSet()
60
{
61
m_iCounter = 0;
62
m_faBlocks.SetName("CSDBlockSet::m_faBlocks");
63
}
64
65
66
CSDBlockSet::~CSDBlockSet()
67
{
68
}
69
70
71
void CSDEngine::Init(int bins)
72
{
73
int z, z2;
74
CSDGroup *gr;
75
CSDBlockSet *bs;
76
77
m_iSteps = 0;
78
79
m_faTempBuffer.SetSize(bins);
80
m_faSum.SetSize(bins);
81
m_faSqSum.SetSize(bins);
82
m_faMean.SetSize(bins);
83
m_faVari.SetSize(bins);
84
m_faCorLen.SetSize(bins);
85
86
for (z=0;z<bins;z++)
87
{
88
m_faTempBuffer[z] = 0;
89
m_faSum[z] = 0;
90
m_faSqSum[z] = 0;
91
m_faMean[z] = 0;
92
m_faVari[z] = 0;
93
m_faCorLen[z] = 0;
94
95
gr = new CSDGroup();
96
m_oaGroups.Add(gr);
97
gr->m_oaBlockSets.SetSize(m_iaBlockSizes.GetSize());
98
gr->m_faBlockVari.SetSize(m_iaBlockSizes.GetSize());
99
for (z2=0;z2<m_iaBlockSizes.GetSize();z2++)
100
{
101
bs = new CSDBlockSet();
102
gr->m_oaBlockSets[z2] = bs;
103
gr->m_faBlockVari[z2] = 0;
104
105
bs->m_faBlocks.Add(0);
106
}
107
}
108
}
109
110
111
void CSDEngine::FinishStep()
112
{
113
int z, z2;
114
CSDGroup *gr;
115
CSDBlockSet *bs;
116
117
m_iSteps++;
118
119
for (z=0;z<m_oaGroups.GetSize();z++)
120
{
121
gr = (CSDGroup*)m_oaGroups[z];
122
123
m_faSum[z] += m_faTempBuffer[z];
124
m_faSqSum[z] += m_faTempBuffer[z]*m_faTempBuffer[z];
125
126
for (z2=0;z2<m_iaBlockSizes.GetSize();z2++)
127
{
128
bs = (CSDBlockSet*)gr->m_oaBlockSets[z2];
129
bs->m_faBlocks[bs->m_faBlocks.GetSize()-1] += m_faTempBuffer[z];
130
131
bs->m_iCounter++;
132
if (bs->m_iCounter >= m_iaBlockSizes[z2])
133
{
134
bs->m_iCounter = 0;
135
bs->m_faBlocks[bs->m_faBlocks.GetSize()-1] /= m_iaBlockSizes[z2];
136
bs->m_faBlocks.Add(0);
137
}
138
}
139
140
m_faTempBuffer[z] = 0;
141
}
142
}
143
144
145
void CSDEngine::DumpData(CDF *df, float timessigma, bool verbose)
146
{
147
int z, z2, ti;
148
CSDGroup *gr;
149
// FILE *a;
150
FILE *b;
151
char buf[256];
152
double *x;
153
double *y, *y2, ya;
154
double da, ds, dr, tf, dsa, dsmi, dsma;
155
CLMWrapper *lm;
156
double par[3], *curve;
157
158
// a = fopen("block.csv","wt");
159
b = NULL;
160
161
mprintf(" Calculating correlation time of histogram...\n");
162
163
try { df->m_pAdditionalSets = new double*[4]; } catch(...) { df->m_pAdditionalSets = NULL; }
164
if (df->m_pAdditionalSets == NULL) NewException((double)4*sizeof(double*),__FILE__,__LINE__,__PRETTY_FUNCTION__);
165
166
try { df->m_pAdditionalSetLabels = new char*[4]; } catch(...) { df->m_pAdditionalSetLabels = NULL; }
167
if (df->m_pAdditionalSetLabels == NULL) NewException((double)4*sizeof(char*),__FILE__,__LINE__,__PRETTY_FUNCTION__);
168
169
df->m_iAdditionalSets = 4;
170
171
for (z=0;z<4;z++)
172
{
173
try { df->m_pAdditionalSets[z] = new double[df->m_iResolution]; } catch(...) { df->m_pAdditionalSets[z] = NULL; }
174
if (df->m_pAdditionalSets[z] == NULL) NewException((double)df->m_iResolution*sizeof(double),__FILE__,__LINE__,__PRETTY_FUNCTION__);
175
}
176
177
sprintf(buf,"%.2f%c confidence range",NormalDistIntegral(timessigma)*2.0-1.0,'%');
178
try { df->m_pAdditionalSetLabels[0] = new char[strlen(buf)+1]; } catch(...) { df->m_pAdditionalSetLabels[0] = NULL; }
179
if (df->m_pAdditionalSetLabels[0] == NULL) NewException((double)(strlen(buf)+1)*sizeof(char),__FILE__,__LINE__,__PRETTY_FUNCTION__);
180
strcpy(df->m_pAdditionalSetLabels[0],buf);
181
182
sprintf(buf,"%.2f%c confidence range",NormalDistIntegral(timessigma)*2.0-1.0,'%');
183
try { df->m_pAdditionalSetLabels[1] = new char[strlen(buf)+1]; } catch(...) { df->m_pAdditionalSetLabels[1] = NULL; }
184
if (df->m_pAdditionalSetLabels[1] == NULL) NewException((double)(strlen(buf)+1)*sizeof(char),__FILE__,__LINE__,__PRETTY_FUNCTION__);
185
strcpy(df->m_pAdditionalSetLabels[1],buf);
186
187
sprintf(buf,"Standard deviation");
188
try { df->m_pAdditionalSetLabels[2] = new char[strlen(buf)+1]; } catch(...) { df->m_pAdditionalSetLabels[2] = NULL; }
189
if (df->m_pAdditionalSetLabels[2] == NULL) NewException((double)(strlen(buf)+1)*sizeof(char),__FILE__,__LINE__,__PRETTY_FUNCTION__);
190
strcpy(df->m_pAdditionalSetLabels[2],buf);
191
192
sprintf(buf,"Correlation Length");
193
try { df->m_pAdditionalSetLabels[3] = new char[strlen(buf)+1]; } catch(...) { df->m_pAdditionalSetLabels[4] = NULL; }
194
if (df->m_pAdditionalSetLabels[3] == NULL) NewException((double)(strlen(buf)+1)*sizeof(char),__FILE__,__LINE__,__PRETTY_FUNCTION__);
195
strcpy(df->m_pAdditionalSetLabels[3],buf);
196
197
try { lm = new CLMWrapper(); } catch(...) { lm = NULL; }
198
if (lm == NULL) NewException((double)sizeof(CLMWrapper),__FILE__,__LINE__,__PRETTY_FUNCTION__);
199
200
try { x = new double[m_iaBlockSizes.GetSize()]; } catch(...) { x = NULL; }
201
if (x == NULL) NewException((double)m_iaBlockSizes.GetSize()*sizeof(double),__FILE__,__LINE__,__PRETTY_FUNCTION__);
202
203
try { y = new double[m_iaBlockSizes.GetSize()]; } catch(...) { y = NULL; }
204
if (y== NULL) NewException((double)m_iaBlockSizes.GetSize()*sizeof(double),__FILE__,__LINE__,__PRETTY_FUNCTION__);
205
206
try { y2 = new double[m_iaBlockSizes.GetSize()]; } catch(...) { y2 = NULL; }
207
if (y2 == NULL) NewException((double)m_iaBlockSizes.GetSize()*sizeof(double),__FILE__,__LINE__,__PRETTY_FUNCTION__);
208
209
try { curve = new double[m_iaBlockSizes.GetSize()]; } catch(...) { curve = NULL; }
210
if (curve == NULL) NewException((double)m_iaBlockSizes.GetSize()*sizeof(double),__FILE__,__LINE__,__PRETTY_FUNCTION__);
211
212
mprintf(" Fitting trajectory block decomposition...\n");
213
mprintf(WHITE," [");
214
215
tf = m_oaGroups.GetSize() / 60.0;
216
217
dsa = 0;
218
dsmi = 1e20;
219
dsma = 0;
220
ti = 0;
221
for (z=0;z<m_oaGroups.GetSize();z++)
222
{
223
if (fmod(z,tf) < 1.0)
224
mprintf(WHITE,"#");
225
226
if (df->m_pBin[z] == 0)
227
continue;
228
229
gr = (CSDGroup*)m_oaGroups[z];
230
231
// mprintf(" * Gruppe %d:\n",z);
232
233
if (verbose)
234
{
235
sprintf(buf,"extra%03d.csv",z);
236
b = fopen(buf,"wt");
237
fprintf(b,"# Block; Block size; sqrt(size); log(size); 1/size; Vari; Size*Vari; S; log(S); log(S2); Fit\n");
238
}
239
240
ya = 0;
241
for (z2=0;z2<m_iaBlockSizes.GetSize();z2++)
242
{
243
x[z2] = log10(m_iaBlockSizes[z2]);
244
y[z2] = log10(m_iaBlockSizes[z2]*gr->m_faBlockVari[z2]/m_faVari[z]);
245
if (z2 == 0)
246
{
247
y2[z2] = y[z2];
248
ya = y[z2];
249
} else
250
{
251
if (y[z2] > ya)
252
{
253
y2[z2] = y[z2];
254
ya = y[z2];
255
} else y2[z2] = ya;
256
}
257
}
258
259
par[0] = 50.0;
260
par[1] = -0.01;
261
par[2] = 1.0;
262
263
lm->Fit_SD(m_iaBlockSizes.GetSize(),x,y2,par,curve,&dr,200);
264
265
if (verbose)
266
fprintf(b,"# A=%G; B=%G; C=%G; R=%G\n",par[0],par[1],par[2],dr);
267
268
ds = pow(10.0,par[0]*(1.0-exp(par[1]*pow(log10(m_iSteps),par[2]))));
269
270
if (ds < dsmi)
271
dsmi = ds;
272
if (ds > dsma)
273
dsma = ds;
274
dsa += ds;
275
ti++;
276
277
m_faCorLen[z] = ds;
278
279
// mprintf(" Fit done. y(x) = %10G * (1 - exp( %10G * x^%10G ) ). R = %10G. S = %10G\n",par[0],par[1],par[2],dr,ds);
280
281
if (verbose)
282
{
283
for (z2=0;z2<m_iaBlockSizes.GetSize();z2++)
284
fprintf(b,"%d; %d; %G; %G; %G; %G; %G; %G; %G; %G; %G\n",z2,(int)m_iaBlockSizes[z2],sqrt(m_iaBlockSizes[z2]),log10(m_iaBlockSizes[z2]),1.0/m_iaBlockSizes[z2],gr->m_faBlockVari[z2],m_iaBlockSizes[z2]*gr->m_faBlockVari[z2],m_iaBlockSizes[z2]*gr->m_faBlockVari[z2]/m_faVari[z],y[z2],y2[z2],curve[z2]);
285
}
286
287
/* x = new double[m_iaBlockSizes.GetSize()];
288
y = new double[m_iaBlockSizes.GetSize()];
289
290
xa = 0;
291
ya = 0;
292
for (z2=0;z2<m_iaBlockSizes.GetSize();z2++)
293
{
294
bs = (CSDBlockSet*)gr->m_oaBlockSets[z2];
295
// mprintf(" - Blockgroesse %d (%d Schritte)\n",z2,m_iaBlockSizes[z2]);
296
// mprintf(" %d Bloecke, Counter %d.\n",bs->m_faBlocks.GetSize(),bs->m_iCounter);
297
298
if (z == 100)
299
{
300
for (z3=0;z3<bs->m_faBlocks.GetSize()-1;z3++)
301
{
302
fprintf(a,"%d; %d; %G\n",z2,z3,bs->m_faBlocks[z3]);
303
}
304
}
305
306
307
x[z2] = log10(m_iaBlockSizes[z2]);
308
y[z2] = log10(m_iaBlockSizes[z2]*gr->m_faBlockVari[z2]/m_faVari[z]);
309
xa += x[z2];
310
ya += y[z2];
311
}
312
313
xa /= m_iaBlockSizes.GetSize();
314
ya /= m_iaBlockSizes.GetSize();
315
316
dxy = 0;
317
dxx = 0;
318
dyy = 0;
319
320
for (z2=0;z2<m_iaBlockSizes.GetSize();z2++)
321
{
322
dxy += (x[z2] - xa) * (y[z2] - ya);
323
dxx += (x[z2] - xa) * (x[z2] - xa);
324
dyy += (y[z2] - ya) * (y[z2] - ya);
325
}
326
327
db = dxy / dxx;
328
da = ya - db * xa;
329
330
ds = pow(10.0, da + db * log10(m_iSteps));
331
332
dr = sqrt( dxy*dxy / dxx / dyy );
333
334
mprintf(" y(x) = %10G + %10G * x; S = %10G; R = %10G\n",da,db,ds,dr);
335
336
for (z2=0;z2<m_iaBlockSizes.GetSize();z2++)
337
fprintf(b,"%d; %d; %G; %G; %G; %G; %G; %G; %G; %G\n",z2,(int)m_iaBlockSizes[z2],sqrt(m_iaBlockSizes[z2]),log10(m_iaBlockSizes[z2]),1.0/m_iaBlockSizes[z2],gr->m_faBlockVari[z2],m_iaBlockSizes[z2]*gr->m_faBlockVari[z2],m_iaBlockSizes[z2]*gr->m_faBlockVari[z2]/m_faVari[z],log10(m_iaBlockSizes[z2]*gr->m_faBlockVari[z2]/m_faVari[z]),da+db*x[z2]);
338
339
340
delete[] x;
341
delete[] y;*/
342
343
if (verbose)
344
fclose(b);
345
346
}
347
delete[] x;
348
delete[] y;
349
delete[] y2;
350
delete[] curve;
351
delete lm;
352
353
mprintf(WHITE,"]\n");
354
355
dsa /= ti;
356
357
mprintf(" Correlation length: Min %.3f, Max %.3f, Avg %.3f time steps.\n",dsmi,dsma,dsa);
358
359
for (z=0;z<m_oaGroups.GetSize();z++)
360
{
361
ds = 0;
362
da = 0;
363
for (z2=z-10;z2<=z+10;z2++)
364
{
365
if (z2 < 0)
366
continue;
367
if (z2 >= m_oaGroups.GetSize())
368
continue;
369
if (df->m_pBin[z2] == 0)
370
continue;
371
ds += m_faCorLen[z2];
372
da += 1.0;
373
}
374
ds /= da;
375
376
if (df->m_pBin[z] != 0)
377
{
378
// df->m_pAdditionalSets[2][z] = sqrt(ds / m_iSteps) * df->m_pBin[z];
379
df->m_pAdditionalSets[2][z] = sqrt(ds / m_iSteps * m_faVari[z]);
380
df->m_pAdditionalSets[3][z] = ds;
381
382
df->m_pAdditionalSets[0][z] = df->m_pBin[z] - timessigma*df->m_pAdditionalSets[2][z];
383
df->m_pAdditionalSets[1][z] = df->m_pBin[z] + timessigma*df->m_pAdditionalSets[2][z];
384
} else
385
{
386
df->m_pAdditionalSets[0][z] = 0;
387
df->m_pAdditionalSets[1][z] = 0;
388
df->m_pAdditionalSets[2][z] = 0;
389
df->m_pAdditionalSets[3][z] = 0;
390
}
391
}
392
393
// fclose(a);
394
}
395
396
397
void CSDEngine::FinishAnalysis()
398
{
399
int z, z2, z3;
400
CSDGroup *gr;
401
CSDBlockSet *bl;
402
403
mprintf(" Calculating block averages...\n");
404
405
for (z=0;z<m_oaGroups.GetSize();z++)
406
{
407
m_faMean[z] = m_faSum[z] / m_iSteps;
408
m_faVari[z] = (m_faSqSum[z] / m_iSteps) - (m_faMean[z]*m_faMean[z]);
409
410
// mprintf(" # Bin %3d: Sum %10G; SqSum %10G; Mean %10G; Vari %10G\n",z,m_faSum[z],m_faSqSum[z],m_faMean[z],m_faVari[z]);
411
412
gr = (CSDGroup*)m_oaGroups[z];
413
414
for (z2=0;z2<gr->m_oaBlockSets.GetSize();z2++)
415
{
416
bl = (CSDBlockSet*)gr->m_oaBlockSets[z2];
417
418
for (z3=0;z3<bl->m_faBlocks.GetSize()-1;z3++)
419
gr->m_faBlockVari[z2] += pow(bl->m_faBlocks[z3] - m_faMean[z], 2);
420
421
gr->m_faBlockVari[z2] /= ((double)bl->m_faBlocks.GetSize()-1);
422
}
423
}
424
}
425
426
Loggerhead is a web-based interface for Breezy
Version: 2.0.1