1 /*
2 chronyd/chronyc - Programs for keeping computer clocks accurate.
3
4 **********************************************************************
5 * Copyright (C) Richard P. Curnow 1997-2003
6 * Copyright (C) Miroslav Lichvar 2011-2014, 2016-2018, 2021
7 *
8 * This program is free software; you can redistribute it and/or modify
9 * it under the terms of version 2 of the GNU General Public License as
10 * published by the Free Software Foundation.
11 *
12 * This program is distributed in the hope that it will be useful, but
13 * WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License along
18 * with this program; if not, write to the Free Software Foundation, Inc.,
19 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
20 *
21 **********************************************************************
22
23 =======================================================================
24
25 This file contains the routines that do the statistical
26 analysis on the samples obtained from the sources,
27 to determined frequencies and error bounds. */
28
29 #include "config.h"
30
31 #include "sysincl.h"
32
33 #include "sourcestats.h"
34 #include "memory.h"
35 #include "regress.h"
36 #include "util.h"
37 #include "conf.h"
38 #include "logging.h"
39 #include "local.h"
40
41 /* ================================================== */
42 /* Define the maxumum number of samples that we want
43 to store per source */
44 #define MAX_SAMPLES 64
45
46 /* This is the assumed worst case bound on an unknown frequency,
47 2000ppm, which would be pretty bad */
48 #define WORST_CASE_FREQ_BOUND (2000.0/1.0e6)
49
50 /* The minimum and maximum assumed skew */
51 #define MIN_SKEW 1.0e-12
52 #define MAX_SKEW 1.0e+02
53
54 /* The minimum standard deviation */
55 #define MIN_STDDEV 1.0e-9
56
57 /* The worst case bound on an unknown standard deviation of the offset */
58 #define WORST_CASE_STDDEV_BOUND 4.0
59
60 /* The asymmetry of network jitter when all jitter is in one direction */
61 #define MAX_ASYMMETRY 0.5
62
63 /* The minimum estimated asymmetry that can activate the offset correction */
64 #define MIN_ASYMMETRY 0.45
65
66 /* The minimum number of consecutive asymmetries with the same sign needed
67 to activate the offset correction */
68 #define MIN_ASYMMETRY_RUN 10
69
70 /* The maximum value of the counter */
71 #define MAX_ASYMMETRY_RUN 1000
72
73 /* ================================================== */
74
75 static LOG_FileID logfileid;
76
77 /* ================================================== */
78 /* This data structure is used to hold the history of data from the
79 source */
80
81 struct SST_Stats_Record {
82
83 /* Reference ID and IP address of source, used for logging to statistics log */
84 uint32_t refid;
85 IPAddr *ip_addr;
86
87 /* User defined minimum and maximum number of samples */
88 int min_samples;
89 int max_samples;
90
91 /* User defined minimum delay */
92 double fixed_min_delay;
93
94 /* User defined asymmetry of network jitter */
95 double fixed_asymmetry;
96
97 /* Number of samples currently stored. The samples are stored in circular
98 buffer. */
99 int n_samples;
100
101 /* Number of extra samples stored in sample_times, offsets and peer_delays
102 arrays that are used to extend the runs test */
103 int runs_samples;
104
105 /* The index of the newest sample */
106 int last_sample;
107
108 /* Flag indicating whether last regression was successful */
109 int regression_ok;
110
111 /* The best individual sample that we are holding, in terms of the minimum
112 root distance at the present time */
113 int best_single_sample;
114
115 /* The index of the sample with minimum delay in peer_delays */
116 int min_delay_sample;
117
118 /* This is the estimated offset (+ve => local fast) at a particular time */
119 double estimated_offset;
120 double estimated_offset_sd;
121 struct timespec offset_time;
122
123 /* Number of runs of the same sign amongst the residuals */
124 int nruns;
125
126 /* Number of consecutive estimated asymmetries with the same sign.
127 The sign of the number encodes the sign of the asymmetry. */
128 int asymmetry_run;
129
130 /* This is the latest estimated asymmetry of network jitter */
131 double asymmetry;
132
133 /* This value contains the estimated frequency. This is the number
134 of seconds that the local clock gains relative to the reference
135 source per unit local time. (Positive => local clock fast,
136 negative => local clock slow) */
137 double estimated_frequency;
138 double estimated_frequency_sd;
139
140 /* This is the assumed worst case bounds on the estimated frequency.
141 We assume that the true frequency lies within +/- half this much
142 about estimated_frequency */
143 double skew;
144
145 /* This is the estimated standard deviation of the data points */
146 double std_dev;
147
148 /* This array contains the sample epochs, in terms of the local
149 clock. */
150 struct timespec sample_times[MAX_SAMPLES * REGRESS_RUNS_RATIO];
151
152 /* This is an array of offsets, in seconds, corresponding to the
153 sample times. In this module, we use the convention that
154 positive means the local clock is FAST of the source and negative
155 means it is SLOW. This is contrary to the convention in the NTP
156 stuff. */
157 double offsets[MAX_SAMPLES * REGRESS_RUNS_RATIO];
158
159 /* This is an array of the offsets as originally measured. Local
160 clock fast of real time is indicated by positive values. This
161 array is not slewed to adjust the readings when we apply
162 adjustments to the local clock, as is done for the array
163 'offset'. */
164 double orig_offsets[MAX_SAMPLES];
165
166 /* This is an array of peer delays, in seconds, being the roundtrip
167 measurement delay to the peer */
168 double peer_delays[MAX_SAMPLES * REGRESS_RUNS_RATIO];
169
170 /* This is an array of peer dispersions, being the skew and local
171 precision dispersion terms from sampling the peer */
172 double peer_dispersions[MAX_SAMPLES];
173
174 /* This array contains the root delays of each sample, in seconds */
175 double root_delays[MAX_SAMPLES];
176
177 /* This array contains the root dispersions of each sample at the
178 time of the measurements */
179 double root_dispersions[MAX_SAMPLES];
180 };
181
182 /* ================================================== */
183
184 static void find_min_delay_sample(SST_Stats inst);
185 static int get_buf_index(SST_Stats inst, int i);
186
187 /* ================================================== */
188
189 void
SST_Initialise(void)190 SST_Initialise(void)
191 {
192 logfileid = CNF_GetLogStatistics() ? LOG_FileOpen("statistics",
193 " Date (UTC) Time IP Address Std dev'n Est offset Offset sd Diff freq Est skew Stress Ns Bs Nr Asym")
194 : -1;
195 }
196
197 /* ================================================== */
198
199 void
SST_Finalise(void)200 SST_Finalise(void)
201 {
202 }
203
204 /* ================================================== */
205 /* This function creates a new instance of the statistics handler */
206
207 SST_Stats
SST_CreateInstance(uint32_t refid,IPAddr * addr,int min_samples,int max_samples,double min_delay,double asymmetry)208 SST_CreateInstance(uint32_t refid, IPAddr *addr, int min_samples, int max_samples,
209 double min_delay, double asymmetry)
210 {
211 SST_Stats inst;
212 inst = MallocNew(struct SST_Stats_Record);
213
214 inst->min_samples = min_samples;
215 inst->max_samples = max_samples;
216 inst->fixed_min_delay = min_delay;
217 inst->fixed_asymmetry = asymmetry;
218
219 SST_SetRefid(inst, refid, addr);
220 SST_ResetInstance(inst);
221
222 return inst;
223 }
224
225 /* ================================================== */
226 /* This function deletes an instance of the statistics handler. */
227
228 void
SST_DeleteInstance(SST_Stats inst)229 SST_DeleteInstance(SST_Stats inst)
230 {
231 Free(inst);
232 }
233
234 /* ================================================== */
235
236 void
SST_ResetInstance(SST_Stats inst)237 SST_ResetInstance(SST_Stats inst)
238 {
239 inst->n_samples = 0;
240 inst->runs_samples = 0;
241 inst->last_sample = 0;
242 inst->regression_ok = 0;
243 inst->best_single_sample = 0;
244 inst->min_delay_sample = 0;
245 inst->estimated_frequency = 0;
246 inst->estimated_frequency_sd = WORST_CASE_FREQ_BOUND;
247 inst->skew = WORST_CASE_FREQ_BOUND;
248 inst->estimated_offset = 0.0;
249 inst->estimated_offset_sd = WORST_CASE_STDDEV_BOUND;
250 UTI_ZeroTimespec(&inst->offset_time);
251 inst->std_dev = WORST_CASE_STDDEV_BOUND;
252 inst->nruns = 0;
253 inst->asymmetry_run = 0;
254 inst->asymmetry = 0.0;
255 }
256
257 /* ================================================== */
258
259 void
SST_SetRefid(SST_Stats inst,uint32_t refid,IPAddr * addr)260 SST_SetRefid(SST_Stats inst, uint32_t refid, IPAddr *addr)
261 {
262 inst->refid = refid;
263 inst->ip_addr = addr;
264 }
265
266 /* ================================================== */
267 /* This function is called to prune the register down when it is full.
268 For now, just discard the oldest sample. */
269
270 static void
prune_register(SST_Stats inst,int new_oldest)271 prune_register(SST_Stats inst, int new_oldest)
272 {
273 if (!new_oldest)
274 return;
275
276 assert(inst->n_samples >= new_oldest);
277 inst->n_samples -= new_oldest;
278 inst->runs_samples += new_oldest;
279 if (inst->runs_samples > inst->n_samples * (REGRESS_RUNS_RATIO - 1))
280 inst->runs_samples = inst->n_samples * (REGRESS_RUNS_RATIO - 1);
281
282 assert(inst->n_samples + inst->runs_samples <= MAX_SAMPLES * REGRESS_RUNS_RATIO);
283
284 find_min_delay_sample(inst);
285 }
286
287 /* ================================================== */
288
289 void
SST_AccumulateSample(SST_Stats inst,NTP_Sample * sample)290 SST_AccumulateSample(SST_Stats inst, NTP_Sample *sample)
291 {
292 int n, m;
293
294 /* Make room for the new sample */
295 if (inst->n_samples > 0 &&
296 (inst->n_samples == MAX_SAMPLES || inst->n_samples == inst->max_samples)) {
297 prune_register(inst, 1);
298 }
299
300 /* Make sure it's newer than the last sample */
301 if (inst->n_samples &&
302 UTI_CompareTimespecs(&inst->sample_times[inst->last_sample], &sample->time) >= 0) {
303 LOG(LOGS_WARN, "Out of order sample detected, discarding history for %s",
304 inst->ip_addr ? UTI_IPToString(inst->ip_addr) : UTI_RefidToString(inst->refid));
305 SST_ResetInstance(inst);
306 }
307
308 n = inst->last_sample = (inst->last_sample + 1) %
309 (MAX_SAMPLES * REGRESS_RUNS_RATIO);
310 m = n % MAX_SAMPLES;
311
312 /* WE HAVE TO NEGATE OFFSET IN THIS CALL, IT IS HERE THAT THE SENSE OF OFFSET
313 IS FLIPPED */
314 inst->sample_times[n] = sample->time;
315 inst->offsets[n] = -sample->offset;
316 inst->orig_offsets[m] = -sample->offset;
317 inst->peer_delays[n] = sample->peer_delay;
318 inst->peer_dispersions[m] = sample->peer_dispersion;
319 inst->root_delays[m] = sample->root_delay;
320 inst->root_dispersions[m] = sample->root_dispersion;
321
322 if (inst->peer_delays[n] < inst->fixed_min_delay)
323 inst->peer_delays[n] = 2.0 * inst->fixed_min_delay - inst->peer_delays[n];
324
325 if (!inst->n_samples || inst->peer_delays[n] < inst->peer_delays[inst->min_delay_sample])
326 inst->min_delay_sample = n;
327
328 ++inst->n_samples;
329 }
330
331 /* ================================================== */
332 /* Return index of the i-th sample in the sample_times and offset buffers,
333 i can be negative down to -runs_samples */
334
335 static int
get_runsbuf_index(SST_Stats inst,int i)336 get_runsbuf_index(SST_Stats inst, int i)
337 {
338 return (unsigned int)(inst->last_sample + 2 * MAX_SAMPLES * REGRESS_RUNS_RATIO -
339 inst->n_samples + i + 1) % (MAX_SAMPLES * REGRESS_RUNS_RATIO);
340 }
341
342 /* ================================================== */
343 /* Return index of the i-th sample in the other buffers */
344
345 static int
get_buf_index(SST_Stats inst,int i)346 get_buf_index(SST_Stats inst, int i)
347 {
348 return (unsigned int)(inst->last_sample + MAX_SAMPLES * REGRESS_RUNS_RATIO -
349 inst->n_samples + i + 1) % MAX_SAMPLES;
350 }
351
352 /* ================================================== */
353 /* This function is used by both the regression routines to find the
354 time interval between each historical sample and the most recent
355 one */
356
357 static void
convert_to_intervals(SST_Stats inst,double * times_back)358 convert_to_intervals(SST_Stats inst, double *times_back)
359 {
360 struct timespec *ts;
361 int i;
362
363 ts = &inst->sample_times[inst->last_sample];
364 for (i = -inst->runs_samples; i < inst->n_samples; i++) {
365 /* The entries in times_back[] should end up negative */
366 times_back[i] = UTI_DiffTimespecsToDouble(&inst->sample_times[get_runsbuf_index(inst, i)], ts);
367 }
368 }
369
370 /* ================================================== */
371
372 static void
find_best_sample_index(SST_Stats inst,double * times_back)373 find_best_sample_index(SST_Stats inst, double *times_back)
374 {
375 /* With the value of skew that has been computed, see which of the
376 samples offers the tightest bound on root distance */
377
378 double root_distance, best_root_distance;
379 double elapsed;
380 int i, j, best_index;
381
382 if (!inst->n_samples)
383 return;
384
385 best_index = -1;
386 best_root_distance = DBL_MAX;
387
388 for (i = 0; i < inst->n_samples; i++) {
389 j = get_buf_index(inst, i);
390
391 elapsed = -times_back[i];
392 assert(elapsed >= 0.0);
393
394 root_distance = inst->root_dispersions[j] + elapsed * inst->skew + 0.5 * inst->root_delays[j];
395 if (root_distance < best_root_distance) {
396 best_root_distance = root_distance;
397 best_index = i;
398 }
399 }
400
401 assert(best_index >= 0);
402 inst->best_single_sample = best_index;
403 }
404
405 /* ================================================== */
406
407 static void
find_min_delay_sample(SST_Stats inst)408 find_min_delay_sample(SST_Stats inst)
409 {
410 int i, index;
411
412 inst->min_delay_sample = get_runsbuf_index(inst, -inst->runs_samples);
413
414 for (i = -inst->runs_samples + 1; i < inst->n_samples; i++) {
415 index = get_runsbuf_index(inst, i);
416 if (inst->peer_delays[index] < inst->peer_delays[inst->min_delay_sample])
417 inst->min_delay_sample = index;
418 }
419 }
420
421 /* ================================================== */
422 /* This function estimates asymmetry of network jitter on the path to the
423 source as a slope of offset against network delay in multiple linear
424 regression. If the asymmetry is significant and its sign doesn't change
425 frequently, the measured offsets (which are used later to estimate the
426 offset and frequency of the clock) are corrected to correspond to the
427 minimum network delay. This can significantly improve the accuracy and
428 stability of the estimated offset and frequency. */
429
430 static int
estimate_asymmetry(double * times_back,double * offsets,double * delays,int n,double * asymmetry,int * asymmetry_run)431 estimate_asymmetry(double *times_back, double *offsets, double *delays, int n,
432 double *asymmetry, int *asymmetry_run)
433 {
434 double a;
435
436 /* Reset the counter when the regression fails or the sign changes */
437 if (!RGR_MultipleRegress(times_back, delays, offsets, n, &a) ||
438 a * *asymmetry_run < 0.0) {
439 *asymmetry = 0;
440 *asymmetry_run = 0.0;
441 return 0;
442 }
443
444 if (a <= -MIN_ASYMMETRY && *asymmetry_run > -MAX_ASYMMETRY_RUN)
445 (*asymmetry_run)--;
446 else if (a >= MIN_ASYMMETRY && *asymmetry_run < MAX_ASYMMETRY_RUN)
447 (*asymmetry_run)++;
448
449 if (abs(*asymmetry_run) < MIN_ASYMMETRY_RUN)
450 return 0;
451
452 *asymmetry = CLAMP(-MAX_ASYMMETRY, a, MAX_ASYMMETRY);
453
454 return 1;
455 }
456
457 /* ================================================== */
458
459 static void
correct_asymmetry(SST_Stats inst,double * times_back,double * offsets)460 correct_asymmetry(SST_Stats inst, double *times_back, double *offsets)
461 {
462 double min_delay, delays[MAX_SAMPLES * REGRESS_RUNS_RATIO];
463 int i, n;
464
465 /* Check if the asymmetry was not specified to be zero */
466 if (inst->fixed_asymmetry == 0.0)
467 return;
468
469 min_delay = SST_MinRoundTripDelay(inst);
470 n = inst->runs_samples + inst->n_samples;
471
472 for (i = 0; i < n; i++)
473 delays[i] = inst->peer_delays[get_runsbuf_index(inst, i - inst->runs_samples)] -
474 min_delay;
475
476 if (fabs(inst->fixed_asymmetry) <= MAX_ASYMMETRY) {
477 inst->asymmetry = inst->fixed_asymmetry;
478 } else {
479 if (!estimate_asymmetry(times_back, offsets, delays, n,
480 &inst->asymmetry, &inst->asymmetry_run))
481 return;
482 }
483
484 /* Correct the offsets */
485 for (i = 0; i < n; i++)
486 offsets[i] -= inst->asymmetry * delays[i];
487 }
488
489 /* ================================================== */
490
491 /* This defines the assumed ratio between the standard deviation of
492 the samples and the peer distance as measured from the round trip
493 time. E.g. a value of 4 means that we think the standard deviation
494 is four times the fluctuation of the peer distance */
495
496 #define SD_TO_DIST_RATIO 0.7
497
498 /* ================================================== */
499 /* This function runs the linear regression operation on the data. It
500 finds the set of most recent samples that give the tightest
501 confidence interval for the frequency, and truncates the register
502 down to that number of samples */
503
504 void
SST_DoNewRegression(SST_Stats inst)505 SST_DoNewRegression(SST_Stats inst)
506 {
507 double times_back[MAX_SAMPLES * REGRESS_RUNS_RATIO];
508 double offsets[MAX_SAMPLES * REGRESS_RUNS_RATIO];
509 double peer_distances[MAX_SAMPLES];
510 double weights[MAX_SAMPLES];
511
512 int degrees_of_freedom;
513 int best_start, times_back_start;
514 double est_intercept, est_slope, est_var, est_intercept_sd, est_slope_sd;
515 int i, j, nruns;
516 double min_distance, median_distance;
517 double sd_weight, sd;
518 double old_skew, old_freq, stress;
519 double precision;
520
521 convert_to_intervals(inst, times_back + inst->runs_samples);
522
523 if (inst->n_samples > 0) {
524 for (i = -inst->runs_samples; i < inst->n_samples; i++) {
525 offsets[i + inst->runs_samples] = inst->offsets[get_runsbuf_index(inst, i)];
526 }
527
528 for (i = 0, min_distance = DBL_MAX; i < inst->n_samples; i++) {
529 j = get_buf_index(inst, i);
530 peer_distances[i] = 0.5 * inst->peer_delays[get_runsbuf_index(inst, i)] +
531 inst->peer_dispersions[j];
532 if (peer_distances[i] < min_distance) {
533 min_distance = peer_distances[i];
534 }
535 }
536
537 /* And now, work out the weight vector */
538
539 precision = LCL_GetSysPrecisionAsQuantum();
540 median_distance = RGR_FindMedian(peer_distances, inst->n_samples);
541
542 sd = (median_distance - min_distance) / SD_TO_DIST_RATIO;
543 sd = CLAMP(precision, sd, min_distance);
544 min_distance += precision;
545
546 for (i=0; i<inst->n_samples; i++) {
547 sd_weight = 1.0;
548 if (peer_distances[i] > min_distance)
549 sd_weight += (peer_distances[i] - min_distance) / sd;
550 weights[i] = SQUARE(sd_weight);
551 }
552 }
553
554 correct_asymmetry(inst, times_back, offsets);
555
556 inst->regression_ok = RGR_FindBestRegression(times_back + inst->runs_samples,
557 offsets + inst->runs_samples, weights,
558 inst->n_samples, inst->runs_samples,
559 inst->min_samples,
560 &est_intercept, &est_slope, &est_var,
561 &est_intercept_sd, &est_slope_sd,
562 &best_start, &nruns, °rees_of_freedom);
563
564 if (inst->regression_ok) {
565
566 old_skew = inst->skew;
567 old_freq = inst->estimated_frequency;
568
569 inst->estimated_frequency = est_slope;
570 inst->estimated_frequency_sd = CLAMP(MIN_SKEW, est_slope_sd, MAX_SKEW);
571 inst->skew = est_slope_sd * RGR_GetTCoef(degrees_of_freedom);
572 inst->estimated_offset = est_intercept;
573 inst->offset_time = inst->sample_times[inst->last_sample];
574 inst->estimated_offset_sd = est_intercept_sd;
575 inst->std_dev = MAX(MIN_STDDEV, sqrt(est_var));
576 inst->nruns = nruns;
577
578 inst->skew = CLAMP(MIN_SKEW, inst->skew, MAX_SKEW);
579 stress = fabs(old_freq - inst->estimated_frequency) / old_skew;
580
581 DEBUG_LOG("off=%e freq=%e skew=%e n=%d bs=%d runs=%d asym=%f arun=%d",
582 inst->estimated_offset, inst->estimated_frequency, inst->skew,
583 inst->n_samples, best_start, inst->nruns,
584 inst->asymmetry, inst->asymmetry_run);
585
586 if (logfileid != -1) {
587 LOG_FileWrite(logfileid, "%s %-15s %10.3e %10.3e %10.3e %10.3e %10.3e %7.1e %3d %3d %3d %5.2f",
588 UTI_TimeToLogForm(inst->offset_time.tv_sec),
589 inst->ip_addr ? UTI_IPToString(inst->ip_addr) : UTI_RefidToString(inst->refid),
590 inst->std_dev,
591 inst->estimated_offset, inst->estimated_offset_sd,
592 inst->estimated_frequency, inst->skew, stress,
593 inst->n_samples, best_start, inst->nruns,
594 inst->asymmetry);
595 }
596
597 times_back_start = inst->runs_samples + best_start;
598 prune_register(inst, best_start);
599 } else {
600 inst->estimated_frequency_sd = WORST_CASE_FREQ_BOUND;
601 inst->skew = WORST_CASE_FREQ_BOUND;
602 inst->estimated_offset_sd = WORST_CASE_STDDEV_BOUND;
603 inst->std_dev = WORST_CASE_STDDEV_BOUND;
604 inst->nruns = 0;
605
606 if (inst->n_samples > 0) {
607 inst->estimated_offset = inst->offsets[inst->last_sample];
608 inst->offset_time = inst->sample_times[inst->last_sample];
609 } else {
610 inst->estimated_offset = 0.0;
611 UTI_ZeroTimespec(&inst->offset_time);
612 }
613
614 times_back_start = 0;
615 }
616
617 find_best_sample_index(inst, times_back + times_back_start);
618
619 }
620
621 /* ================================================== */
622 /* Return the assumed worst case range of values that this source's
623 frequency lies within. Frequency is defined as the amount of time
624 the local clock gains relative to the source per unit local clock
625 time. */
626 void
SST_GetFrequencyRange(SST_Stats inst,double * lo,double * hi)627 SST_GetFrequencyRange(SST_Stats inst,
628 double *lo, double *hi)
629 {
630 double freq, skew;
631 freq = inst->estimated_frequency;
632 skew = inst->skew;
633 *lo = freq - skew;
634 *hi = freq + skew;
635
636 /* This function is currently used only to determine the values of delta
637 and epsilon in the ntp_core module. Limit the skew to a reasonable maximum
638 to avoid failing the dispersion test too easily. */
639 if (skew > WORST_CASE_FREQ_BOUND) {
640 *lo = -WORST_CASE_FREQ_BOUND;
641 *hi = WORST_CASE_FREQ_BOUND;
642 }
643 }
644
645 /* ================================================== */
646
647 void
SST_GetSelectionData(SST_Stats inst,struct timespec * now,double * offset_lo_limit,double * offset_hi_limit,double * root_distance,double * std_dev,double * first_sample_ago,double * last_sample_ago,int * select_ok)648 SST_GetSelectionData(SST_Stats inst, struct timespec *now,
649 double *offset_lo_limit,
650 double *offset_hi_limit,
651 double *root_distance,
652 double *std_dev,
653 double *first_sample_ago,
654 double *last_sample_ago,
655 int *select_ok)
656 {
657 double offset, sample_elapsed;
658 int i, j;
659
660 if (!inst->n_samples) {
661 *select_ok = 0;
662 return;
663 }
664
665 i = get_runsbuf_index(inst, inst->best_single_sample);
666 j = get_buf_index(inst, inst->best_single_sample);
667
668 *std_dev = inst->std_dev;
669
670 sample_elapsed = fabs(UTI_DiffTimespecsToDouble(now, &inst->sample_times[i]));
671 offset = inst->offsets[i] + sample_elapsed * inst->estimated_frequency;
672 *root_distance = 0.5 * inst->root_delays[j] +
673 inst->root_dispersions[j] + sample_elapsed * inst->skew;
674
675 *offset_lo_limit = offset - *root_distance;
676 *offset_hi_limit = offset + *root_distance;
677
678 #if 0
679 double average_offset, elapsed;
680 int average_ok;
681 /* average_ok ignored for now */
682 elapsed = UTI_DiffTimespecsToDouble(now, &inst->offset_time);
683 average_offset = inst->estimated_offset + inst->estimated_frequency * elapsed;
684 if (fabs(average_offset - offset) <=
685 inst->peer_dispersions[j] + 0.5 * inst->peer_delays[i]) {
686 average_ok = 1;
687 } else {
688 average_ok = 0;
689 }
690 #endif
691
692 i = get_runsbuf_index(inst, 0);
693 *first_sample_ago = UTI_DiffTimespecsToDouble(now, &inst->sample_times[i]);
694 i = get_runsbuf_index(inst, inst->n_samples - 1);
695 *last_sample_ago = UTI_DiffTimespecsToDouble(now, &inst->sample_times[i]);
696
697 *select_ok = inst->regression_ok;
698
699 /* If maxsamples is too small to have a successful regression, enable the
700 selection as a special case for a fast update/print-once reference mode */
701 if (!*select_ok && inst->n_samples < 3 && inst->n_samples == inst->max_samples) {
702 *std_dev = CNF_GetMaxJitter();
703 *select_ok = 1;
704 }
705
706 DEBUG_LOG("n=%d off=%f dist=%f sd=%f first_ago=%f last_ago=%f selok=%d",
707 inst->n_samples, offset, *root_distance, *std_dev,
708 *first_sample_ago, *last_sample_ago, *select_ok);
709 }
710
711 /* ================================================== */
712
713 void
SST_GetTrackingData(SST_Stats inst,struct timespec * ref_time,double * average_offset,double * offset_sd,double * frequency,double * frequency_sd,double * skew,double * root_delay,double * root_dispersion)714 SST_GetTrackingData(SST_Stats inst, struct timespec *ref_time,
715 double *average_offset, double *offset_sd,
716 double *frequency, double *frequency_sd, double *skew,
717 double *root_delay, double *root_dispersion)
718 {
719 int i, j;
720 double elapsed_sample;
721
722 assert(inst->n_samples > 0);
723
724 i = get_runsbuf_index(inst, inst->best_single_sample);
725 j = get_buf_index(inst, inst->best_single_sample);
726
727 *ref_time = inst->offset_time;
728 *average_offset = inst->estimated_offset;
729 *offset_sd = inst->estimated_offset_sd;
730 *frequency = inst->estimated_frequency;
731 *frequency_sd = inst->estimated_frequency_sd;
732 *skew = inst->skew;
733 *root_delay = inst->root_delays[j];
734
735 elapsed_sample = UTI_DiffTimespecsToDouble(&inst->offset_time, &inst->sample_times[i]);
736 *root_dispersion = inst->root_dispersions[j] + inst->skew * elapsed_sample + *offset_sd;
737
738 DEBUG_LOG("n=%d off=%f offsd=%f freq=%e freqsd=%e skew=%e delay=%f disp=%f",
739 inst->n_samples, *average_offset, *offset_sd,
740 *frequency, *frequency_sd, *skew, *root_delay, *root_dispersion);
741 }
742
743 /* ================================================== */
744
745 void
SST_SlewSamples(SST_Stats inst,struct timespec * when,double dfreq,double doffset)746 SST_SlewSamples(SST_Stats inst, struct timespec *when, double dfreq, double doffset)
747 {
748 int m, i;
749 double delta_time;
750 struct timespec *sample, prev;
751 double prev_offset, prev_freq;
752
753 if (!inst->n_samples)
754 return;
755
756 for (m = -inst->runs_samples; m < inst->n_samples; m++) {
757 i = get_runsbuf_index(inst, m);
758 sample = &inst->sample_times[i];
759 prev = *sample;
760 UTI_AdjustTimespec(sample, when, sample, &delta_time, dfreq, doffset);
761 inst->offsets[i] += delta_time;
762 }
763
764 /* Update the regression estimates */
765 prev = inst->offset_time;
766 prev_offset = inst->estimated_offset;
767 prev_freq = inst->estimated_frequency;
768 UTI_AdjustTimespec(&inst->offset_time, when, &inst->offset_time,
769 &delta_time, dfreq, doffset);
770 inst->estimated_offset += delta_time;
771 inst->estimated_frequency = (inst->estimated_frequency - dfreq) / (1.0 - dfreq);
772
773 DEBUG_LOG("n=%d m=%d old_off_time=%s new=%s old_off=%f new_off=%f old_freq=%.3f new_freq=%.3f",
774 inst->n_samples, inst->runs_samples,
775 UTI_TimespecToString(&prev), UTI_TimespecToString(&inst->offset_time),
776 prev_offset, inst->estimated_offset,
777 1.0e6 * prev_freq, 1.0e6 * inst->estimated_frequency);
778 }
779
780 /* ================================================== */
781
782 void
SST_CorrectOffset(SST_Stats inst,double doffset)783 SST_CorrectOffset(SST_Stats inst, double doffset)
784 {
785 int i;
786
787 if (!inst->n_samples)
788 return;
789
790 for (i = -inst->runs_samples; i < inst->n_samples; i++)
791 inst->offsets[get_runsbuf_index(inst, i)] += doffset;
792
793 inst->estimated_offset += doffset;
794 }
795
796 /* ================================================== */
797
798 void
SST_AddDispersion(SST_Stats inst,double dispersion)799 SST_AddDispersion(SST_Stats inst, double dispersion)
800 {
801 int m, i;
802
803 for (m = 0; m < inst->n_samples; m++) {
804 i = get_buf_index(inst, m);
805 inst->root_dispersions[i] += dispersion;
806 inst->peer_dispersions[i] += dispersion;
807 }
808 }
809
810 /* ================================================== */
811
812 double
SST_PredictOffset(SST_Stats inst,struct timespec * when)813 SST_PredictOffset(SST_Stats inst, struct timespec *when)
814 {
815 double elapsed;
816
817 if (inst->n_samples < 3) {
818 /* We don't have any useful statistics, and presumably the poll
819 interval is minimal. We can't do any useful prediction other
820 than use the latest sample or zero if we don't have any samples */
821 if (inst->n_samples > 0) {
822 return inst->offsets[inst->last_sample];
823 } else {
824 return 0.0;
825 }
826 } else {
827 elapsed = UTI_DiffTimespecsToDouble(when, &inst->offset_time);
828 return inst->estimated_offset + elapsed * inst->estimated_frequency;
829 }
830
831 }
832
833 /* ================================================== */
834
835 double
SST_MinRoundTripDelay(SST_Stats inst)836 SST_MinRoundTripDelay(SST_Stats inst)
837 {
838 if (inst->fixed_min_delay > 0.0)
839 return inst->fixed_min_delay;
840
841 if (!inst->n_samples)
842 return DBL_MAX;
843
844 return inst->peer_delays[inst->min_delay_sample];
845 }
846
847 /* ================================================== */
848
849 int
SST_GetDelayTestData(SST_Stats inst,struct timespec * sample_time,double * last_sample_ago,double * predicted_offset,double * min_delay,double * skew,double * std_dev)850 SST_GetDelayTestData(SST_Stats inst, struct timespec *sample_time,
851 double *last_sample_ago, double *predicted_offset,
852 double *min_delay, double *skew, double *std_dev)
853 {
854 if (inst->n_samples < 6)
855 return 0;
856
857 *last_sample_ago = UTI_DiffTimespecsToDouble(sample_time, &inst->offset_time);
858 *predicted_offset = inst->estimated_offset +
859 *last_sample_ago * inst->estimated_frequency;
860 *min_delay = SST_MinRoundTripDelay(inst);
861 *skew = inst->skew;
862 *std_dev = inst->std_dev;
863
864 return 1;
865 }
866
867 /* ================================================== */
868 /* This is used to save the register to a file, so that we can reload
869 it after restarting the daemon */
870
871 int
SST_SaveToFile(SST_Stats inst,FILE * out)872 SST_SaveToFile(SST_Stats inst, FILE *out)
873 {
874 int m, i, j;
875
876 if (inst->n_samples < 1)
877 return 0;
878
879 if (fprintf(out, "%d %d\n", inst->n_samples, inst->asymmetry_run) < 0)
880 return 0;
881
882 for(m = 0; m < inst->n_samples; m++) {
883 i = get_runsbuf_index(inst, m);
884 j = get_buf_index(inst, m);
885
886 if (fprintf(out, "%s %.6e %.6e %.6e %.6e %.6e %.6e\n",
887 UTI_TimespecToString(&inst->sample_times[i]),
888 inst->offsets[i], inst->orig_offsets[j],
889 inst->peer_delays[i], inst->peer_dispersions[j],
890 inst->root_delays[j], inst->root_dispersions[j]) < 0)
891 return 0;
892 }
893
894 return 1;
895 }
896
897 /* ================================================== */
898 /* This is used to reload samples from a file */
899
900 int
SST_LoadFromFile(SST_Stats inst,FILE * in)901 SST_LoadFromFile(SST_Stats inst, FILE *in)
902 {
903 int i, n_samples, arun;
904 struct timespec now;
905 double sample_time;
906 char line[256];
907
908 if (!fgets(line, sizeof (line), in) ||
909 sscanf(line, "%d %d", &n_samples, &arun) != 2 ||
910 n_samples < 1 || n_samples > MAX_SAMPLES)
911 return 0;
912
913 SST_ResetInstance(inst);
914
915 LCL_ReadCookedTime(&now, NULL);
916
917 for (i = 0; i < n_samples; i++) {
918 if (!fgets(line, sizeof (line), in) ||
919 sscanf(line, "%lf %lf %lf %lf %lf %lf %lf",
920 &sample_time, &inst->offsets[i], &inst->orig_offsets[i],
921 &inst->peer_delays[i], &inst->peer_dispersions[i],
922 &inst->root_delays[i], &inst->root_dispersions[i]) != 7)
923 return 0;
924
925 if (!UTI_IsTimeOffsetSane(&now, sample_time - UTI_TimespecToDouble(&now)))
926 return 0;
927
928 /* Some resolution is lost in the double format, but that's ok */
929 UTI_DoubleToTimespec(sample_time, &inst->sample_times[i]);
930
931 /* Make sure the samples are sane and they are in order */
932 if (!UTI_IsTimeOffsetSane(&inst->sample_times[i], -inst->offsets[i]) ||
933 !(fabs(inst->peer_delays[i]) < 1.0e6 && fabs(inst->peer_dispersions[i]) < 1.0e6 &&
934 fabs(inst->root_delays[i]) < 1.0e6 && fabs(inst->root_dispersions[i]) < 1.0e6) ||
935 (i > 0 && UTI_CompareTimespecs(&inst->sample_times[i],
936 &inst->sample_times[i - 1]) <= 0))
937 return 0;
938 }
939
940 inst->n_samples = n_samples;
941 inst->last_sample = inst->n_samples - 1;
942 inst->asymmetry_run = CLAMP(-MAX_ASYMMETRY_RUN, arun, MAX_ASYMMETRY_RUN);
943
944 find_min_delay_sample(inst);
945 SST_DoNewRegression(inst);
946
947 return 1;
948 }
949
950 /* ================================================== */
951
952 void
SST_DoSourceReport(SST_Stats inst,RPT_SourceReport * report,struct timespec * now)953 SST_DoSourceReport(SST_Stats inst, RPT_SourceReport *report, struct timespec *now)
954 {
955 int i, j;
956 struct timespec last_sample_time;
957
958 if (inst->n_samples > 0) {
959 i = get_runsbuf_index(inst, inst->n_samples - 1);
960 j = get_buf_index(inst, inst->n_samples - 1);
961 report->orig_latest_meas = inst->orig_offsets[j];
962 report->latest_meas = inst->offsets[i];
963 report->latest_meas_err = 0.5*inst->root_delays[j] + inst->root_dispersions[j];
964
965 /* Align the sample time to reduce the leak of the receive timestamp */
966 last_sample_time = inst->sample_times[i];
967 last_sample_time.tv_nsec = 0;
968 report->latest_meas_ago = UTI_DiffTimespecsToDouble(now, &last_sample_time);
969 } else {
970 report->latest_meas_ago = (uint32_t)-1;
971 report->orig_latest_meas = 0;
972 report->latest_meas = 0;
973 report->latest_meas_err = 0;
974 report->stratum = 0;
975 }
976 }
977
978 /* ================================================== */
979
980 int
SST_Samples(SST_Stats inst)981 SST_Samples(SST_Stats inst)
982 {
983 return inst->n_samples;
984 }
985
986 /* ================================================== */
987
988 void
SST_DoSourcestatsReport(SST_Stats inst,RPT_SourcestatsReport * report,struct timespec * now)989 SST_DoSourcestatsReport(SST_Stats inst, RPT_SourcestatsReport *report, struct timespec *now)
990 {
991 double dspan;
992 double elapsed, sample_elapsed;
993 int bi, bj;
994
995 report->n_samples = inst->n_samples;
996 report->n_runs = inst->nruns;
997
998 if (inst->n_samples > 0) {
999 bi = get_runsbuf_index(inst, inst->best_single_sample);
1000 bj = get_buf_index(inst, inst->best_single_sample);
1001
1002 dspan = UTI_DiffTimespecsToDouble(&inst->sample_times[inst->last_sample],
1003 &inst->sample_times[get_runsbuf_index(inst, 0)]);
1004 elapsed = UTI_DiffTimespecsToDouble(now, &inst->offset_time);
1005 sample_elapsed = UTI_DiffTimespecsToDouble(now, &inst->sample_times[bi]);
1006
1007 report->span_seconds = round(dspan);
1008 report->est_offset = inst->estimated_offset + elapsed * inst->estimated_frequency;
1009 report->est_offset_err = inst->estimated_offset_sd + sample_elapsed * inst->skew +
1010 (0.5 * inst->root_delays[bj] + inst->root_dispersions[bj]);
1011 } else {
1012 report->span_seconds = 0;
1013 report->est_offset = 0;
1014 report->est_offset_err = 0;
1015 }
1016
1017 report->resid_freq_ppm = 1.0e6 * inst->estimated_frequency;
1018 report->skew_ppm = 1.0e6 * inst->skew;
1019 report->sd = inst->std_dev;
1020 }
1021
1022 /* ================================================== */
1023
1024 double
SST_GetJitterAsymmetry(SST_Stats inst)1025 SST_GetJitterAsymmetry(SST_Stats inst)
1026 {
1027 return inst->asymmetry;
1028 }
1029
1030 /* ================================================== */
1031