1 /*************************************************************************\
2 *
3 * Package: TestU01
4 * File: sspacings.c
5 * Environment: ANSI C
6 *
7 * Copyright (c) 2002 Pierre L'Ecuyer, DIRO, Université de Montréal.
8 * e-mail: lecuyer@iro.umontreal.ca
9 * All rights reserved.
10 *
11 * Redistribution and use in source and binary forms, with or without
12 * modification, are permitted without a fee for private, research,
13 * academic, or other non-commercial purposes.
14 * Any use of this software in a commercial environment requires a
15 * written licence from the copyright owner.
16 *
17 * Any changes made to this package must be clearly identified as such.
18 *
19 * In scientific publications which used this software, a reference to it
20 * would be appreciated.
21 *
22 * Redistributions of source code must retain this copyright notice
23 * and the following disclaimer.
24 *
25 * THIS PACKAGE IS PROVIDED "AS IS" AND WITHOUT ANY EXPRESS OR
26 * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
27 * WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
28 *
29 \*************************************************************************/
30
31 #include "util.h"
32 #include "chrono.h"
33 #include "num.h"
34 #include "tables.h"
35
36 #include "sspacings.h"
37 #include "unif01.h"
38 #include "wdist.h"
39 #include "swrite.h"
40
41 #include "statcoll.h"
42 #include "gofw.h"
43 #include "fbar.h"
44
45 #include <math.h>
46 #include <string.h>
47 #include <stdio.h>
48
49
50
51
52
53
54 /*------------------------------- Constants -------------------------------*/
55
56 /* String length */
57 #define LEN 50
58
59 /* Euler's constant */
60 static const double EULER = 0.577215664901533;
61
62 /* Necessary so that the multiplication of a large number of U01 will not
63 underflow; the precise value is not so important: it could be bigger. */
64 static const double Epsilon = 1.0E-200;
65
66
67 #define MAXM 201
68
69
70 /*--------------------------------- Types -------------------------------*/
71
72 typedef enum
73 {
74 LOG_EXACT_CIRC, /* Log of spacings, exact, circular */
75 LOG_EXACT_LIN, /* Log of spacings, exact, linear */
76 LOG_ASYMP_CIRC, /* Log of spacings, asymptotic, circular */
77 LOG_ASYMP_LIN, /* Log of spacings, asymptotic, linear */
78 SQUA_EXACT_CIRC, /* Square of spacings, exact, circular */
79 SQUA_EXACT_LIN, /* Square of spacings, exact, linear */
80 SQUA_ASYMP_CIRC, /* Square of spacings, asymptotic, circular */
81 SQUA_ASYMP_LIN, /* Square of spacings, asymptotic, linear */
82 Stats_N /* Number of statistics in this enum */
83 } Stats;
84
85
86 typedef struct
87 {
88 int NbColl; /* Number of collectors used for each m */
89 int Nbm; /* Number of values of m to consider */
90 int Loc[MAXM]; /* The statistic for m is in Loc [m] */
91 double Mu[MAXM][Stats_N]; /* Mean */
92 double Sig[MAXM][Stats_N]; /* Standard deviation */
93 double HMu[MAXM][Stats_N]; /* Empirical mean */
94 double HSig[MAXM][Stats_N]; /* Empirical standard deviation */
95 } Param;
96
97
98
99 /*------------------------------- Functions -------------------------------*/
100
101
sspacings_CreateRes(void)102 sspacings_Res * sspacings_CreateRes (void)
103 {
104 sspacings_Res *res;
105 res = util_Malloc (sizeof (sspacings_Res));
106 memset (res, 0, sizeof (sspacings_Res));
107 res->name = util_Calloc (1, sizeof (char));
108 res->smax = -1;
109 return res;
110 }
111
112
113 /*-------------------------------------------------------------------------*/
114
sspacings_DeleteRes(sspacings_Res * res)115 void sspacings_DeleteRes (sspacings_Res *res)
116 {
117 int j;
118 if (res == NULL)
119 return;
120 for (j = 0; j <= res->smax; j += 2)
121 res->Collectors[j] = statcoll_Delete (res->Collectors[j]);
122 util_Free (res->Collectors);
123
124 for (j = 0; j <= res->imax; j++) {
125 sres_DeleteBasic (res->LogCAMu[j]);
126 sres_DeleteBasic (res->LogCEMu[j]);
127 sres_DeleteBasic (res->SquareCAMu[j]);
128 sres_DeleteBasic (res->SquareCEMu[j]);
129 }
130
131 util_Free (res->LogCEMu);
132 util_Free (res->LogCAMu);
133 util_Free (res->SquareCEMu);
134 util_Free (res->SquareCAMu);
135
136 util_Free (res->LogCESig_sVal);
137 util_Free (res->LogCESig_pVal);
138 util_Free (res->LogCASig_sVal);
139 util_Free (res->LogCASig_pVal);
140 util_Free (res->SquareCESig_sVal);
141 util_Free (res->SquareCESig_pVal);
142 util_Free (res->SquareCASig_sVal);
143 util_Free (res->SquareCASig_pVal);
144
145 util_Free (res->name);
146 util_Free (res);
147 }
148
149
150 /*-------------------------------------------------------------------------*/
151
InitRes(sspacings_Res * res,long N,int Nbm,char * nam)152 static void InitRes (
153 sspacings_Res *res, /* Results holder */
154 long N, /* Number of replications */
155 int Nbm, /* Number of values of m to consider */
156 char *nam /* Test name */
157 )
158 /*
159 * Initializes res
160 */
161 {
162 char nom[LEN + 1];
163 char spindex[LEN + 1];
164 int j;
165
166 if (res->smax < 0) {
167 res->Collectors = util_Calloc ((size_t) 8 * Nbm,
168 sizeof (statcoll_Collector *));
169 for (j = 0; j < 8 * Nbm; j += 2)
170 res->Collectors[j] = statcoll_Create (N, "");
171
172 res->LogCAMu = util_Calloc ((size_t) Nbm, sizeof (sres_Basic *));
173 res->LogCEMu = util_Calloc ((size_t) Nbm, sizeof (sres_Basic *));
174 res->SquareCAMu = util_Calloc ((size_t) Nbm, sizeof (sres_Basic *));
175 res->SquareCEMu = util_Calloc ((size_t) Nbm, sizeof (sres_Basic *));
176 for (j = 0; j < Nbm; j++) {
177 res->LogCAMu[j] = sres_CreateBasic ();
178 res->LogCEMu[j] = sres_CreateBasic ();
179 res->SquareCAMu[j] = sres_CreateBasic ();
180 res->SquareCEMu[j] = sres_CreateBasic ();
181 }
182
183 res->LogCESig_sVal = util_Calloc ((size_t) Nbm, sizeof (double));
184 res->LogCESig_pVal = util_Calloc ((size_t) Nbm, sizeof (double));
185 res->LogCASig_sVal = util_Calloc ((size_t) Nbm, sizeof (double));
186 res->LogCASig_pVal = util_Calloc ((size_t) Nbm, sizeof (double));
187 res->SquareCESig_sVal = util_Calloc ((size_t) Nbm, sizeof (double));
188 res->SquareCESig_pVal = util_Calloc ((size_t) Nbm, sizeof (double));
189 res->SquareCASig_sVal = util_Calloc ((size_t) Nbm, sizeof (double));
190 res->SquareCASig_pVal = util_Calloc ((size_t) Nbm, sizeof (double));
191
192 } else {
193 for (j = 8 * Nbm; j <= res->smax; j += 2)
194 res->Collectors[j] = statcoll_Delete (res->Collectors[j]);
195 res->Collectors = util_Realloc (res->Collectors,
196 8 * Nbm * sizeof (statcoll_Collector *));
197 for (j = res->smax + 2; j < 8 * Nbm; j += 2)
198 res->Collectors[j] = statcoll_Create (N, "");
199
200 for (j = Nbm; j <= res->imax; j++) {
201 sres_DeleteBasic (res->LogCAMu[j]);
202 sres_DeleteBasic (res->LogCEMu[j]);
203 sres_DeleteBasic (res->SquareCAMu[j]);
204 sres_DeleteBasic (res->SquareCEMu[j]);
205 }
206 res->LogCAMu = util_Realloc (res->LogCAMu,
207 Nbm * sizeof (sres_Basic *));
208 res->LogCEMu = util_Realloc (res->LogCEMu,
209 Nbm * sizeof (sres_Basic *));
210 res->SquareCAMu = util_Realloc (res->SquareCAMu,
211 Nbm * sizeof (sres_Basic *));
212 res->SquareCEMu = util_Realloc (res->SquareCEMu,
213 Nbm * sizeof (sres_Basic *));
214
215 for (j = res->imax + 1; j < Nbm; j++) {
216 res->LogCAMu[j] = sres_CreateBasic ();
217 res->LogCEMu[j] = sres_CreateBasic ();
218 res->SquareCAMu[j] = sres_CreateBasic ();
219 res->SquareCEMu[j] = sres_CreateBasic ();
220 }
221
222 res->LogCESig_sVal =
223 util_Realloc (res->LogCESig_sVal, Nbm * sizeof (double));
224 res->LogCESig_pVal =
225 util_Realloc (res->LogCESig_pVal, Nbm * sizeof (double));
226 res->LogCASig_sVal =
227 util_Realloc (res->LogCASig_sVal, Nbm * sizeof (double));
228 res->LogCASig_pVal =
229 util_Realloc (res->LogCASig_pVal, Nbm * sizeof (double));
230 res->SquareCESig_sVal =
231 util_Realloc (res->SquareCESig_sVal, Nbm * sizeof (double));
232 res->SquareCESig_pVal =
233 util_Realloc (res->SquareCESig_pVal, Nbm * sizeof (double));
234 res->SquareCASig_sVal =
235 util_Realloc (res->SquareCASig_sVal, Nbm * sizeof (double));
236 res->SquareCASig_pVal =
237 util_Realloc (res->SquareCASig_pVal, Nbm * sizeof (double));
238 }
239
240 for (j = 0; j < 8 * Nbm; j += 2)
241 statcoll_Init (res->Collectors[j], N);
242 res->smax = 8 * Nbm - 2;
243
244 for (j = 0; j < Nbm; j++) {
245 strncpy (nom, "LogCEMu[", (size_t) LEN);
246 sprintf (spindex, "%1d", j);
247 strncat (nom, spindex, (size_t) 5);
248 strncat (nom, "]", (size_t) 2);
249 sres_InitBasic (res->LogCEMu[j], N, nom);
250
251 strncpy (nom, "LogCAMu[", (size_t) LEN);
252 sprintf (spindex, "%1d", j);
253 strncat (nom, spindex, (size_t) 5);
254 strncat (nom, "]", (size_t) 2);
255 sres_InitBasic (res->LogCAMu[j], N, nom);
256
257 strncpy (nom, "SquareCEMu[", (size_t) LEN);
258 sprintf (spindex, "%1d", j);
259 strncat (nom, spindex, (size_t) 5);
260 strncat (nom, "]", (size_t) 2);
261 sres_InitBasic (res->SquareCEMu[j], N, nom);
262
263 strncpy (nom, "SquareCAMu[", (size_t) LEN);
264 sprintf (spindex, "%1d", j);
265 strncat (nom, spindex, (size_t) 5);
266 strncat (nom, "]", (size_t) 2);
267 sres_InitBasic (res->SquareCAMu[j], N, nom);
268
269 res->imax = j;
270 }
271 res->name = util_Realloc (res->name, 1 + strlen (nam) * sizeof (char));
272 strcpy (res->name, nam);
273
274 memset (res->LogCESig_sVal, 0, sizeof (res->LogCESig_sVal));
275 memset (res->LogCESig_pVal, 0, sizeof (res->LogCESig_pVal));
276 memset (res->LogCASig_sVal, 0, sizeof (res->LogCASig_sVal));
277 memset (res->LogCASig_pVal, 0, sizeof (res->LogCASig_pVal));
278 memset (res->SquareCESig_sVal, 0, sizeof (res->SquareCESig_sVal));
279 memset (res->SquareCESig_pVal, 0, sizeof (res->SquareCESig_pVal));
280 memset (res->SquareCASig_sVal, 0, sizeof (res->SquareCASig_sVal));
281 memset (res->SquareCASig_pVal, 0, sizeof (res->SquareCASig_pVal));
282 }
283
284
285 /*=========================================================================*/
286
287
sspacings_SumLogsSpacings(unif01_Gen * gen,sspacings_Res * res,long N,long n,int r,int m)288 void sspacings_SumLogsSpacings (unif01_Gen * gen, sspacings_Res * res,
289 long N, long n, int r, int m)
290 {
291 util_Error ("sspacings_SumLogsSpacings is implemented inside "
292 "sspacings_AllSpacings2");
293 }
294
295
296 /*=========================================================================*/
297
sspacings_SumSquaresSpacings(unif01_Gen * gen,sspacings_Res * res,long N,long n,int r,int m)298 void sspacings_SumSquaresSpacings (unif01_Gen * gen, sspacings_Res * res,
299 long N, long n, int r, int m)
300 {
301 util_Error ("sspacings_SumSquaresSpacings is implemented inside "
302 "sspacings_AllSpacings2");
303 }
304
305
306 /*=========================================================================*/
307
sspacings_ScanSpacings(unif01_Gen * gen,sspacings_Res * res,long N,long n,int r,double d)308 void sspacings_ScanSpacings (unif01_Gen * gen, sspacings_Res * res,
309 long N, long n, int r, double d)
310 {
311 util_Error ("sspacings_ScanSpacings not implemented");
312 }
313
314
315 /*=========================================================================*/
316
WrRes(char S1[],long N,Param * par,int m,int i,statcoll_Collector * Collectors[],gofw_TestArray sv,gofw_TestArray pv)317 static void WrRes (char S1[], long N, Param * par, int m, int i,
318 statcoll_Collector * Collectors[], gofw_TestArray sv, gofw_TestArray pv)
319 {
320 printf ("%s", S1);
321 printf ("\n Mu = ");
322 num_WriteD (par->Mu[m][i], 12, 8, 7);
323 printf ("\n Sigma = ");
324 num_WriteD (par->Sig[m][i], 12, 8, 7);
325 printf ("\n\nEmpirical mean of standardized values :");
326 num_WriteD (par->HMu[m][i] / N, 12, 8, 7);
327 printf ("\n");
328 gofw_Writep1 (fbar_Normal1 (par->HMu[m][i] / N));
329 printf ("Second empirical moment of standardized values:");
330 num_WriteD (par->HSig[m][i] / N, 12, 8, 7);
331 printf ("\n");
332 gofw_Writep1 (fbar_ChiSquare2 (N, 12, par->HSig[m][i]));
333 i += 8 * par->Loc[m];
334 if (N > 1)
335 gofw_WriteActiveTests0 (N, sv, pv);
336
337 if (swrite_Collectors) {
338 statcoll_Write (Collectors[i], 5, 14, 4, 3);
339 printf ("\n");
340 }
341 printf ("\n");
342 }
343
344 /*-------------------------------------------------------------------------*/
345
UpdateStat(Param * par,int m,int i,double v,statcoll_Collector * Collectors[])346 static void UpdateStat (Param * par, int m, int i, double v,
347 statcoll_Collector * Collectors[])
348 {
349 double x;
350
351 x = (v - par->Mu[m][i]) / par->Sig[m][i];
352 par->HMu[m][i] += x;
353 par->HSig[m][i] += x * x;
354 i += 8 * par->Loc[m];
355 statcoll_AddObs (Collectors[i], x);
356 }
357
358 /*-------------------------------------------------------------------------*/
359
InitAllSpacings(unif01_Gen * gen,char * TestName,Param * par,long N,long n0,int r,int M0,int M1,int D,int LgEps)360 static void InitAllSpacings (unif01_Gen * gen, char *TestName, Param * par,
361 long N, long n0, int r, int M0, int M1, int D, int LgEps)
362 {
363 double Rmn[MAXM]; /* R(m,n) and Q(m,n) from Cressie (1976) */
364 double Qmn[MAXM];
365 double Rmm1[MAXM];
366 double Qmm1[MAXM];
367 double Rm1;
368 double Qm1; /* R(1,m-1) and Q(1,m-1) from Cressie */
369 double Fact; /* n^2 / (n+2) */
370 double Mu0;
371 double x, y;
372 double I;
373 double R;
374 double Q;
375 double nLR;
376 double mLR;
377 int m;
378 int i;
379
380 if (swrite_Basic) {
381 swrite_Head (gen, TestName, N, n0, r);
382 printf (", M0 = %1d, M1 = %1d, D = %1d\n", M0, M1, D);
383 printf (" LgEps = %1d\n\n\n", LgEps);
384 }
385 util_Assert (M1 < MAXM, "InitAllSpacings: M1 is too large");
386
387 par->Nbm = 1 + (M1 - M0) / D;
388 for (i = 0; i < par->Nbm; i++)
389 par->Loc[M0 + i * D] = i;
390 if (M0 == 0)
391 par->Loc[1] = 0;
392
393 /* Initialize constants and arrays */
394 nLR = n0;
395 Fact = nLR * nLR / (nLR + 2.0);
396 R = 0.0;
397 Q = 0.0;
398
399 for (i = n0; i >= M1; i--) {
400 I = 1.0 / i;
401 R += I;
402 Q += I * I;
403 }
404
405 Rmn[M1] = R;
406 Qmn[M1] = Q;
407 for (m = M1 - 1; m >= 1; m--) {
408 I = 1.0 / m;
409 Rmn[m] = Rmn[m + 1] + I;
410 Qmn[m] = Qmn[m + 1] + I * I;
411 }
412
413 m = 1;
414 Rmm1[1] = 0.0;
415 Qmm1[1] = 0.0;
416 for (m = 2; m <= M1; m++) {
417 I = 1.0 / (m - 1);
418 Rmm1[m] = Rmm1[m - 1] + I;
419 Qmm1[m] = Qmm1[m - 1] + I * I;
420 }
421
422 /* Compute theoretical means and variances */
423 if (M0 == 0)
424 m = 1;
425 else
426 m = M0;
427
428 while (m <= M1) {
429 mLR = m;
430 Rm1 = Rmm1[m];
431 Qm1 = Qmm1[m];
432 par->Mu[m][LOG_EXACT_CIRC] = -(nLR + 1.0) * Rmn[m];
433 par->Mu[m][LOG_EXACT_LIN] =
434 par->Mu[m][LOG_EXACT_CIRC] * (nLR + 2.0 - mLR) / (nLR + 1.0);
435 par->Mu[m][LOG_ASYMP_CIRC] =
436 -(nLR + 1.0) * (log (nLR + 1.0) + EULER - Rm1);
437 par->Mu[m][LOG_ASYMP_LIN] =
438 par->Mu[m][LOG_ASYMP_CIRC] * (nLR + 2.0 - mLR) / (nLR + 1.0);
439 Mu0 = mLR * (mLR + 1.0);
440 par->Mu[m][SQUA_EXACT_CIRC] = Mu0 * Fact;
441 par->Mu[m][SQUA_EXACT_LIN] =
442 par->Mu[m][SQUA_EXACT_CIRC] * (nLR - mLR + 2.0) / (nLR + 1.0);
443 par->Mu[m][SQUA_ASYMP_CIRC] = Mu0 * (nLR + 1.0);
444 par->Mu[m][SQUA_ASYMP_LIN] = Mu0 * (nLR - mLR + 2.0);
445
446 x = (2 * m * (m - 1) + 1) * ((num_Pi * num_Pi) / 6.0 - Qm1) +
447 (-2 * m + 1);
448 util_Assert (x > 0.0, "Negative Sig [m, 2]");
449 par->Sig[m][LOG_ASYMP_CIRC] = sqrt (nLR * x);
450 par->Sig[m][LOG_ASYMP_LIN] = par->Sig[m][LOG_ASYMP_CIRC]; /* See Holst
451 1979 */
452
453 x = Qmn[m] + nLR * Qmn[1] - 2.0 * (mLR - 1.0) * (mLR * Qm1 + 1.0) +
454 (2.0 * mLR * (mLR - 1.0) - nLR) * num_Pi * num_Pi / 6.0;
455 util_Assert (x > 0.0, "Negative Sig [m, 0] ...");
456 par->Sig[m][LOG_EXACT_CIRC] = sqrt ((nLR + 1.0) * x);
457
458 y = 2 * m * (m + 1) * (2 * m + 1) / 3.0;
459 par->Sig[m][SQUA_ASYMP_CIRC] = sqrt (nLR * y);
460 par->Sig[m][SQUA_ASYMP_LIN] = par->Sig[m][SQUA_ASYMP_CIRC];
461
462 x = 2.0 * mLR * (1.0 + mLR) * (2.0 + mLR * (1.0 - 3.0 * mLR) +
463 nLR * (1.0 + 2.0 * mLR)) / 3.0;
464 x = x / ((nLR + 3.0) * (nLR + 4.0));
465 util_Assert (x > 0.0, "Negative Sig [m, 4]");
466 par->Sig[m][SQUA_EXACT_CIRC] = sqrt (x) * Fact;
467
468 x = 20.0 + mLR * (-54.0 + mLR * (6.0 + mLR * (58.0 - 30.0 * mLR)));
469 x += nLR * (34.0 + mLR * (-37.0 + mLR * (-27.0 + mLR * (48.0 -
470 12.0 * mLR))));
471 x += nLR * nLR * (16.0 + mLR * (3.0 + mLR * (-15.0 + 8.0 * mLR)));
472 x += nLR * nLR * nLR * 2.0 * (1.0 + 2.0 * mLR);
473 x = x * mLR * (1.0 + mLR) / 3.0;
474 x = x / ((nLR + 3.0) * (nLR + 4.0));
475 util_Assert (x > 0.0, "Negative Sig [m, 5]");
476 par->Sig[m][SQUA_EXACT_LIN] = sqrt (x) * Fact / (nLR + 1.0);
477
478 /* Initalization of tables for the tests */
479 for (i = 0; i < Stats_N; i++) {
480 par->HMu[m][i] = 0.0;
481 par->HSig[m][i] = 0.0;
482 }
483 if (M0 == 0 && m == 1)
484 m = D;
485 else
486 m += D;
487 }
488 /*
489 if (swrite_Basic)
490 printf (" Number of collectors initialized to N: %1d\n\n",
491 par->NbColl * par->Nbm);
492 */
493 }
494
495
496 /*=========================================================================*/
497
CopyResults(sspacings_Res * res,Param * par,long N,int M0,int M1,int D,int flag)498 static void CopyResults (sspacings_Res * res, Param * par, long N,
499 int M0, int M1, int D, int flag)
500 {
501 int m, s, i, k;
502
503 if (M0 == 0)
504 m = 1;
505 else
506 m = M0;
507 s = 0;
508
509 while (m <= M1) {
510 i = LOG_EXACT_CIRC;
511 k = i + 8 * par->Loc[m];
512 tables_CopyTabD (res->Collectors[k]->V, res->LogCEMu[s]->sVal1->V, 1,
513 N);
514 res->LogCEMu[s]->sVal1->NObs = N;
515 gofw_ActiveTests2 (res->Collectors[k]->V, res->LogCEMu[s]->pVal1->V,
516 N, wdist_Normal, (double *) NULL, res->LogCEMu[s]->sVal2,
517 res->LogCEMu[s]->pVal2);
518
519 res->LogCEMu[s]->sVal2[gofw_Mean] = par->HMu[m][i] / N;
520 res->LogCEMu[s]->pVal2[gofw_Mean] = fbar_Normal1 (par->HMu[m][i] / N);
521 res->LogCESig_sVal[s] = par->HSig[m][i] / N;
522 res->LogCESig_pVal[s] = fbar_ChiSquare2 (N, 12, par->HSig[m][i]);
523
524 i = LOG_ASYMP_CIRC;
525 k = i + 8 * par->Loc[m];
526 tables_CopyTabD (res->Collectors[k]->V, res->LogCAMu[s]->sVal1->V, 1,
527 N);
528 res->LogCAMu[s]->sVal1->NObs = N;
529 if (flag) {
530 gofw_ActiveTests2 (res->Collectors[k]->V, res->LogCAMu[s]->pVal1->V,
531 N, wdist_Normal, (double *) NULL, res->LogCAMu[s]->sVal2,
532 res->LogCAMu[s]->pVal2);
533
534 res->LogCAMu[s]->sVal2[gofw_Mean] = par->HMu[m][i] / N;
535 res->LogCAMu[s]->pVal2[gofw_Mean] =
536 fbar_Normal1 (par->HMu[m][i] / N);
537 res->LogCASig_sVal[s] = par->HSig[m][i] / N;
538 res->LogCASig_pVal[s] = fbar_ChiSquare2 (N, 12, par->HSig[m][i]);
539 } else {
540 res->LogCAMu[s]->sVal2[gofw_Mean] = 0.0;
541 res->LogCAMu[s]->pVal2[gofw_Mean] = 0.0;
542 res->LogCASig_sVal[s] = 0.0;
543 res->LogCASig_pVal[s] = 0.0;
544 memset (res->LogCAMu[s]->sVal2, 0, sizeof (res->LogCAMu[s]->sVal2));
545 memset (res->LogCAMu[s]->pVal2, 0, sizeof (res->LogCAMu[s]->pVal2));
546 }
547
548 i = SQUA_EXACT_CIRC;
549 k = i + 8 * par->Loc[m];
550 tables_CopyTabD (res->Collectors[k]->V, res->SquareCEMu[s]->sVal1->V,
551 1, N);
552 res->SquareCEMu[s]->sVal1->NObs = N;
553 gofw_ActiveTests2 (res->Collectors[k]->V, res->SquareCEMu[s]->pVal1->V,
554 N, wdist_Normal, (double *) NULL, res->SquareCEMu[s]->sVal2,
555 res->SquareCEMu[s]->pVal2);
556
557 res->SquareCEMu[s]->sVal2[gofw_Mean] = par->HMu[m][i] / N;
558 res->SquareCEMu[s]->pVal2[gofw_Mean] =
559 fbar_Normal1 (par->HMu[m][i] / N);
560 res->SquareCESig_sVal[s] = par->HSig[m][i] / N;
561 res->SquareCESig_pVal[s] = fbar_ChiSquare2 (N, 12, par->HSig[m][i]);
562
563 i = SQUA_ASYMP_CIRC;
564 k = i + 8 * par->Loc[m];
565 tables_CopyTabD (res->Collectors[k]->V, res->SquareCAMu[s]->sVal1->V,
566 1, N);
567 res->SquareCAMu[s]->sVal1->NObs = N;
568 if (flag) {
569 gofw_ActiveTests2 (res->Collectors[k]->V,
570 res->SquareCAMu[s]->pVal1->V, N, wdist_Normal, (double *) NULL,
571 res->SquareCAMu[s]->sVal2, res->SquareCAMu[s]->pVal2);
572 res->SquareCAMu[s]->sVal2[gofw_Mean] = par->HMu[m][i] / N;
573 res->SquareCAMu[s]->pVal2[gofw_Mean] =
574 fbar_Normal1 (par->HMu[m][i] / N);
575 res->SquareCASig_sVal[s] = par->HSig[m][i] / N;
576 res->SquareCASig_pVal[s] = fbar_ChiSquare2 (N, 12, par->HSig[m][i]);
577
578 } else {
579 res->SquareCAMu[s]->sVal2[gofw_Mean] = 0.0;
580 res->SquareCAMu[s]->pVal2[gofw_Mean] = 0.0;
581 res->SquareCASig_sVal[s] = 0.0;
582 res->SquareCASig_pVal[s] = 0.0;
583 memset (res->SquareCAMu[s]->sVal2, 0,
584 sizeof (res->SquareCAMu[s]->sVal2));
585 memset (res->SquareCAMu[s]->pVal2, 0,
586 sizeof (res->SquareCAMu[s]->pVal2));
587 }
588
589 if (M0 == 0 && m == 1)
590 m = D;
591 else
592 m += D;
593 s++;
594 }
595 }
596
597
598 /*=========================================================================*/
599
sspacings_AllSpacings(unif01_Gen * gen,sspacings_Res * res,long N,long n0,int r,int M0,int M1,int D,int LgEps)600 void sspacings_AllSpacings (unif01_Gen * gen, sspacings_Res * res,
601 long N, long n0, int r, int M0, int M1, int D, int LgEps)
602 {
603 long i;
604 int m, s;
605 long Seq;
606 double Eps; /* Minimal spacing for SumLogsSpacings */
607 double *U;
608 double Prod, x;
609 double LnProd;
610 double SumSq;
611 int NbMinus[MAXM]; /* Number of spacings < Eps */
612 Param par;
613 lebool localRes = FALSE;
614 chrono_Chrono *Timer;
615 char *TestName = "sspacings_AllSpacings test";
616
617 Timer = chrono_Create ();
618 memset (&par, 0, sizeof (Param));
619 par.NbColl = 4;
620 InitAllSpacings (gen, TestName, &par, N, n0, r, M0, M1, D, LgEps);
621 Eps = 1.0 / num_TwoExp[LgEps];
622
623 if (res == NULL) {
624 localRes = TRUE;
625 res = sspacings_CreateRes ();
626 }
627 InitRes (res, N, par.Nbm, "sspacings_AllSpacings");
628 res->step = 2;
629
630 U = util_Calloc ((size_t) n0 + M1 + 2, sizeof (double));
631 U[0] = 0.0;
632
633 /* Beginning of tests */
634 for (Seq = 1; Seq <= N; Seq++) {
635 for (i = 1; i <= n0; i++)
636 U[i] = unif01_StripD (gen, r);
637 tables_QuickSortD (U, 1, n0);
638 util_Assert (U[1] >= 0.0, "sspacings_AllSpacings: U[1] < 0.0");
639 util_Assert (U[n0] <= 1.0, "sspacings_AllSpacings: U[n] > 1.0");
640
641 for (i = 1; i <= M1; i++)
642 U[n0 + i] = 1.0 + U[i - 1];
643 if (M0 == 0)
644 m = 1;
645 else
646 m = M0;
647
648 while (m <= M1) {
649 NbMinus[m] = 0;
650 Prod = 1.0;
651 LnProd = 0.0;
652 SumSq = 0.0;
653
654 for (i = 0; i <= n0; i++) {
655 x = U[i + m] - U[i];
656 SumSq += x * x;
657 /* In case a spacing is zero */
658 if (x < Eps) {
659 x = Eps;
660 ++NbMinus[m];
661 }
662 /* Compute log of product instead of sum of logs; it is faster */
663 Prod *= x;
664 if (Prod < Epsilon) {
665 /* Take log once in a while to avoid underflow */
666 LnProd += log (Prod);
667 Prod = 1.0;
668 }
669 }
670 LnProd += log (Prod);
671
672 UpdateStat (&par, m, LOG_EXACT_CIRC, LnProd, res->Collectors);
673 UpdateStat (&par, m, LOG_ASYMP_CIRC, LnProd, res->Collectors);
674 UpdateStat (&par, m, SQUA_EXACT_CIRC, SumSq * n0 * n0,
675 res->Collectors);
676 UpdateStat (&par, m, SQUA_ASYMP_CIRC, SumSq * n0 * n0,
677 res->Collectors);
678 if (M0 == 0 && m == 1)
679 m = D;
680 else
681 m += D;
682 }
683 }
684
685 CopyResults (res, &par, N, M0, M1, D, 1);
686 if (swrite_Basic) {
687 printf ("\nResults:");
688 if (M0 == 0)
689 m = 1;
690 else
691 m = M0;
692 s = 0;
693 while (m <= M1) {
694 printf ("\n----------------------------------------------------\n");
695 printf ("m = %1d\n\n", m);
696 if (NbMinus[m] > 0)
697 printf ("%1d spacings < 1 / 2^%1d\n\n", NbMinus[m], LgEps);
698
699 printf ("Logs of spacings:\n-----------------\n\n");
700 WrRes ("Exact mean and standard deviation, circular:",
701 N, &par, m, LOG_EXACT_CIRC, res->Collectors,
702 res->LogCEMu[s]->sVal2, res->LogCEMu[s]->pVal2);
703 WrRes ("Asymptotic mean and standard deviation, circular:",
704 N, &par, m, LOG_ASYMP_CIRC, res->Collectors,
705 res->LogCAMu[s]->sVal2, res->LogCAMu[s]->pVal2);
706
707 printf ("\nSquares of spacings:\n--------------------\n\n");
708 WrRes ("Exact mean and standard deviation, circular:",
709 N, &par, m, SQUA_EXACT_CIRC, res->Collectors,
710 res->SquareCEMu[s]->sVal2, res->SquareCEMu[s]->pVal2);
711 WrRes ("Asymptotic mean and standard deviation, circular:",
712 N, &par, m, SQUA_ASYMP_CIRC, res->Collectors,
713 res->SquareCAMu[s]->sVal2, res->SquareCAMu[s]->pVal2);
714
715 if (M0 == 0 && m == 1)
716 m = D;
717 else
718 m += D;
719 s++;
720 }
721 printf ("\n");
722 swrite_Final (gen, Timer);
723 }
724
725 U = util_Free (U);
726 if (localRes)
727 sspacings_DeleteRes (res);
728 chrono_Delete (Timer);
729 }
730
731
732 /*=========================================================================*/
733
sspacings_AllSpacings2(unif01_Gen * gen,sspacings_Res * res,long N,long n0,int r,int M0,int M1,int D,int LgEps)734 void sspacings_AllSpacings2 (unif01_Gen * gen, sspacings_Res * res,
735 long N, long n0, int r, int M0, int M1, int D, int LgEps)
736 {
737 long i;
738 int m, j, s;
739 long Seq;
740 double Eps; /* Minimal spacing for SumLogsSpacings */
741 double *U;
742 double Prod, x;
743 double LnProd;
744 double SumSq;
745 int NbMinus[MAXM]; /* Number of spacings < Eps */
746 Param par;
747 lebool localRes = FALSE;
748 chrono_Chrono *Timer;
749 char *TestName = "sspacings_AllSpacings2 test";
750
751 Timer = chrono_Create ();
752 memset (&par, 0, sizeof (Param));
753 par.NbColl = 2;
754 InitAllSpacings (gen, TestName, &par, N, n0, r, M0, M1, D, LgEps);
755 Eps = 1.0 / num_TwoExp[LgEps];
756
757 if (res == NULL) {
758 localRes = TRUE;
759 res = sspacings_CreateRes ();
760 }
761 InitRes (res, N, par.Nbm, "sspacings_AllSpacings2");
762 res->step = 4;
763
764 U = util_Calloc ((size_t) n0 + M1 + 2, sizeof (double));
765 U[0] = 0.0;
766
767 /* Beginning of tests */
768 for (Seq = 1; Seq <= N; Seq++) {
769 for (i = 1; i <= n0; i++)
770 U[i] = unif01_StripD (gen, r);
771 tables_QuickSortD (U, 1, n0);
772 for (j = 1; j <= M1; j++)
773 U[n0 + j] = 1.0 + U[j - 1];
774 if (M0 == 0)
775 m = 1;
776 else
777 m = M0;
778
779 while (m <= M1) {
780 NbMinus[m] = 0;
781 Prod = 1.0;
782 LnProd = 0.0;
783 SumSq = 0.0;
784
785 for (i = 0; i <= n0; i++) {
786 x = U[i + m] - U[i];
787 SumSq += x * x;
788 /* In case a spacing is zero */
789 if (x < Eps) {
790 x = Eps;
791 ++NbMinus[m];
792 }
793 /* Compute log of product instead of sum of logs; it is faster */
794 Prod *= x;
795 if (Prod < Epsilon) {
796 /* Take log once in a while to avoid underflow */
797 LnProd += log (Prod);
798 Prod = 1.0;
799 }
800 }
801 LnProd += log (Prod);
802
803 UpdateStat (&par, m, LOG_EXACT_CIRC, LnProd, res->Collectors);
804 UpdateStat (&par, m, SQUA_EXACT_CIRC, SumSq * n0 * n0,
805 res->Collectors);
806 if (M0 == 0 && m == 1)
807 m = D;
808 else
809 m += D;
810 }
811 }
812
813 CopyResults (res, &par, N, M0, M1, D, 0);
814 if (swrite_Basic) {
815 printf ("\nResults:");
816 if (M0 == 0)
817 m = 1;
818 else
819 m = M0;
820 s = 0;
821
822 while (m <= M1) {
823 printf ("\n----------------------------------------------------\n");
824 printf ("m = %1d\n\n", m);
825 if (NbMinus[m] > 0)
826 printf ("%1d spacings < 1 / 2^%1d\n\n", NbMinus[m], LgEps);
827 printf ("Logs of spacings:\n-----------------\n\n");
828 WrRes ("Exact mean and standard deviation, circular:",
829 N, &par, m, LOG_EXACT_CIRC, res->Collectors,
830 res->LogCEMu[s]->sVal2, res->LogCEMu[s]->pVal2);
831 printf ("\nSquares of spacings:\n--------------------\n\n");
832 WrRes ("Exact mean and standard deviation, circular:",
833 N, &par, m, SQUA_EXACT_CIRC, res->Collectors,
834 res->SquareCEMu[s]->sVal2, res->SquareCEMu[s]->pVal2);
835 if (M0 == 0 && m == 1)
836 m = D;
837 else
838 m += D;
839 s++;
840 }
841 printf ("\n");
842 swrite_Final (gen, Timer);
843 }
844
845 U = util_Free (U);
846 if (localRes)
847 sspacings_DeleteRes (res);
848 chrono_Delete (Timer);
849 }
850