Change client_check() to calculate the best offset and the best frequency
[dragonfly.git] / usr.sbin / dntpd / client.c
1 /*
2  * Copyright (c) 2005 The DragonFly Project.  All rights reserved.
3  * 
4  * This code is derived from software contributed to The DragonFly Project
5  * by Matthew Dillon <dillon@backplane.com>
6  * 
7  * Redistribution and use in source and binary forms, with or without
8  * modification, are permitted provided that the following conditions
9  * are met:
10  * 
11  * 1. Redistributions of source code must retain the above copyright
12  *    notice, this list of conditions and the following disclaimer.
13  * 2. Redistributions in binary form must reproduce the above copyright
14  *    notice, this list of conditions and the following disclaimer in
15  *    the documentation and/or other materials provided with the
16  *    distribution.
17  * 3. Neither the name of The DragonFly Project nor the names of its
18  *    contributors may be used to endorse or promote products derived
19  *    from this software without specific, prior written permission.
20  * 
21  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
22  * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
23  * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
24  * FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE
25  * COPYRIGHT HOLDERS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
26  * INCIDENTAL, SPECIAL, EXEMPLARY OR CONSEQUENTIAL DAMAGES (INCLUDING,
27  * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
28  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
29  * AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
30  * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
31  * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
32  * SUCH DAMAGE.
33  * 
34  * $DragonFly: src/usr.sbin/dntpd/client.c,v 1.4 2005/04/24 09:39:27 dillon Exp $
35  */
36
37 #include "defs.h"
38
39 void
40 client_init(void)
41 {
42 }
43
44 int
45 client_main(struct server_info **info_ary, int count)
46 {
47     struct server_info *best_off;
48     struct server_info *best_freq;
49     double freq;
50     double offset;
51     int i;
52
53     for (;;) {
54         /*
55          * Poll clients and accumulate data
56          */
57         for (i = 0; i < count; ++i)
58             client_poll(info_ary[i]);
59
60         /*
61          * Find the best client (or synthesize one).  A different client
62          * can be chosen for frequency and offset.  Note in particular 
63          * that offset counters and averaging code gets reset when an
64          * offset correction is made (otherwise the averaging history will
65          * cause later corrections to overshoot).  
66          * 
67          * The regression used to calculate the frequency is a much 
68          * longer-term entity and is NOT reset, so it is still possible
69          * for the offset correction code to make minor adjustments to
70          * the frequency if it so desires.
71          *
72          * client_check may replace the server_info pointer with a new
73          * one.
74          */
75         best_off = NULL;
76         best_freq = NULL;
77         for (i = 0; i < count; ++i)
78             client_check(&info_ary[i], &best_off, &best_freq);
79
80         /*
81          * Offset correction.
82          *
83          * XXX it might not be a good idea to issue an offset correction if
84          * a prior offset correction is still in progress as this will skew
85          * the offset calculation.  XXX either figure out how to correct the
86          * skew or do not issue a correction.
87          */
88         if (best_off) {
89             offset = best_off->lin_sumoffset / best_off->lin_countoffset;
90             lin_resetalloffsets(info_ary, count);
91             freq = sysntp_correct_offset(offset);
92         } else {
93             freq = 0.0;
94         }
95
96         /*
97          * Frequency correction (throw away minor freq adjusts from the
98          * offset code if we can't do a frequency correction here).
99          */
100         if (best_freq) {
101             sysntp_correct_freq(best_freq->lin_cache_freq + freq);
102         }
103
104         if (debug_sleep)
105             usleep(debug_sleep * 1000000 + random() % 100000);
106         else
107             usleep(60 * 1000000 + random() % 100000);
108     }
109 }
110
111 void
112 client_poll(server_info_t info)
113 {
114     struct timeval rtv;
115     struct timeval ltv;
116     struct timeval lbtv;
117     double offset;
118
119     if (debug_opt) {
120         fprintf(stderr, "%s: poll, ", info->target);
121         fflush(stderr);
122     }
123     if (udp_ntptimereq(info->fd, &rtv, &ltv, &lbtv) < 0) {
124         if (debug_opt)
125             fprintf(stderr, "no response\n");
126         return;
127     }
128
129     /*
130      * Figure out the offset (the difference between the reported
131      * time and our current time) for linear regression purposes.
132      */
133     offset = tv_delta_double(&rtv, &ltv);
134
135     while (info) {
136         /*
137          * Linear regression
138          */
139         if (debug_opt) {
140             struct tm *tp;
141             char buf[64];
142             time_t t;
143
144             t = rtv.tv_sec;
145             tp = localtime(&t);
146             strftime(buf, sizeof(buf), "%d-%b-%Y %H:%M:%S", tp);
147             fprintf(stderr, "%s.%03ld ", 
148                     buf, rtv.tv_usec / 1000);
149         }
150         lin_regress(info, &ltv, &lbtv, offset);
151         info = info->altinfo;
152         if (info && debug_opt) {
153             fprintf(stderr, "%*.*s: poll, ", 
154                 (int)strlen(info->target), 
155                 (int)strlen(info->target),
156                 "(alt)");
157             fflush(stderr);
158         }
159     }
160 }
161
162 /*
163  * Find the best client (or synthesize a fake info structure to return).
164  * We can find separate best clients for offset and frequency.
165  */
166 void
167 client_check(struct server_info **checkp, 
168              struct server_info **best_off,
169              struct server_info **best_freq)
170 {
171     struct server_info *check = *checkp;
172     struct server_info *info;
173
174     /*
175      * Start an alternate linear regression once our current one
176      * has passed a certain point.
177      */
178     if (check->lin_count >= LIN_RESTART / 2 && check->altinfo == NULL) {
179         info = malloc(sizeof(*info));
180         assert(info != NULL);
181         bcopy(check, info, sizeof(*info));
182         check->altinfo = info;
183         lin_reset(info);
184     }
185
186     /*
187      * Replace our current linear regression with the alternate once
188      * the current one has hit its limit (beyond a certain point the
189      * linear regression starts to work against us, preventing us from
190      * reacting to changing conditions).
191      *
192      * Report any significant change in the offset or ppm.
193      */
194     if (check->lin_count >= LIN_RESTART) {
195         if ((info = check->altinfo) && info->lin_count >= LIN_RESTART / 2) {
196             double freq_diff;
197
198             if (debug_opt) {
199                 freq_diff = info->lin_cache_freq - check->lin_cache_freq;
200                 printf("%s: Switching to alternate, Frequence difference is %6.3f ppm\n",
201                         info->target, freq_diff * 1.0E+6);
202             }
203             *checkp = info;
204             free(check);
205             check = info;
206         }
207     }
208
209     /*
210      * BEST CLIENT FOR FREQUENCY CORRECTION:
211      *
212      *  8 samples and a correllation > 0.99, or
213      * 16 samples and a correllation > 0.96
214      */
215     info = *best_freq;
216     if ((check->lin_count >= 8 && fabs(check->lin_cache_corr) >= 0.99) ||
217         (check->lin_count >= 16 && fabs(check->lin_cache_corr) >= 0.96)
218     ) {
219         if (info == NULL || 
220             fabs(check->lin_cache_corr) > fabs(info->lin_cache_corr)
221         ) {
222             info = check;
223             *best_freq = info;
224         }
225
226     }
227
228     /*
229      * BEST CLIENT FOR OFFSET CORRECTION:
230      *
231      * Use the standard-deviation and require at least 4 samples.  An
232      * offset correction is valid if the standard deviation is less then
233      * the average offset divided by 4.
234      */
235     info = *best_off;
236     if (check->lin_countoffset >= 4 && 
237         check->lin_cache_stddev < fabs(check->lin_sumoffset / check->lin_countoffset / 4)) {
238         if (info == NULL || 
239             fabs(check->lin_cache_stddev) < fabs(info->lin_cache_stddev)
240         ) {
241             info = check;
242             *best_off = info;
243         }
244     }
245 }
246
247 /*
248  * Linear regression.
249  *
250  *      ltv     local time as of when the offset error was calculated between
251  *              local time and remote time.
252  *
253  *      lbtv    base time as of when local time was obtained.  Used to
254  *              calculate the cumulative corrections made to the system's
255  *              real time clock so we can de-correct the offset for the
256  *              linear regression.
257  *
258  * X is the time axis, in seconds.
259  * Y is the uncorrected offset, in seconds.
260  */
261 void
262 lin_regress(server_info_t info, struct timeval *ltv, struct timeval *lbtv,
263             double offset)
264 {
265     double time_axis;
266     double uncorrected_offset;
267
268     /*
269      * De-correcting the offset:
270      *
271      *  The passed offset is (our_real_time - remote_real_time).  To remove
272      *  corrections from our_real_time we take the difference in the basetime
273      *  (new_base_time - old_base_time) and subtract that from the offset.
274      *  That is, if the basetime goesup, the uncorrected offset goes down.
275      */
276     if (info->lin_count == 0) {
277         info->lin_tv = *ltv;
278         info->lin_btv = *lbtv;
279         time_axis = 0;
280         uncorrected_offset = offset;
281     } else {
282         time_axis = tv_delta_double(&info->lin_tv, ltv);
283         uncorrected_offset = offset - tv_delta_double(&info->lin_btv, lbtv);
284     }
285
286     /*
287      * We have to use the uncorrected offset for frequency calculations.
288      */
289     ++info->lin_count;
290     info->lin_sumx += time_axis;
291     info->lin_sumx2 += time_axis * time_axis;
292     info->lin_sumy += uncorrected_offset;
293     info->lin_sumy2 += uncorrected_offset * uncorrected_offset;
294     info->lin_sumxy += time_axis * uncorrected_offset;
295
296     /*
297      * We have to use the corrected offset for offset calculations.
298      */
299     ++info->lin_countoffset;
300     info->lin_sumoffset += offset;
301     info->lin_sumoffset2 += offset * offset;
302
303     /*
304      * Calculate various derived values.   This gets us slope, y-intercept,
305      * and correllation from the linear regression.
306      */
307     if (info->lin_count > 1) {
308         info->lin_cache_slope = 
309          (info->lin_count * info->lin_sumxy - info->lin_sumx * info->lin_sumy) /
310          (info->lin_count * info->lin_sumx2 - info->lin_sumx * info->lin_sumx);
311
312         info->lin_cache_yint = 
313          (info->lin_sumy - info->lin_cache_slope * info->lin_sumx) /
314          (info->lin_count);
315
316         info->lin_cache_corr =
317          (info->lin_count * info->lin_sumxy - info->lin_sumx * info->lin_sumy) /
318          sqrt((info->lin_count * info->lin_sumx2 - 
319                       info->lin_sumx * info->lin_sumx) *
320              (info->lin_count * info->lin_sumy2 - 
321                       info->lin_sumy * info->lin_sumy)
322          );
323     }
324
325     /*
326      * Calculate more derived values.  This gets us the standard-deviation
327      * of offsets.  The standard deviation approximately means that 68%
328      * of the samples fall within the calculated stddev of the mean.
329      */
330     if (info->lin_countoffset > 1) {
331          info->lin_cache_stddev = 
332              sqrt((info->lin_sumoffset2 - 
333                  ((info->lin_sumoffset * info->lin_sumoffset / 
334                    info->lin_countoffset))) /
335                  (info->lin_countoffset - 1.0));
336     }
337
338     /*
339      * Save the most recent offset, we might use it in the future.
340      * Save the frequency correction (we might scale the slope later so
341      * we have a separate field for the actual frequency correction in
342      * seconds per second).
343      */
344     info->lin_cache_offset = offset;
345     info->lin_cache_freq = info->lin_cache_slope;
346
347     if (debug_opt) {
348         fprintf(stderr, "iter=%2d time=%7.3f off=%.6f uoff=%.6f",
349             (int)info->lin_count,
350             time_axis, offset, uncorrected_offset);
351         if (info->lin_count > 1) {
352             fprintf(stderr, " slope %7.6f"
353                             " yint %3.2f corr %7.6f freq_ppm %4.2f", 
354                 info->lin_cache_slope,
355                 info->lin_cache_yint,
356                 info->lin_cache_corr,
357                 info->lin_cache_freq * 1000000.0);
358         }
359         if (info->lin_countoffset > 1) {
360             fprintf(stderr, " stddev %7.6f", info->lin_cache_stddev);
361         }
362         fprintf(stderr, "\n");
363     }
364 }
365
366 /*
367  * Reset the linear regression data.  The info structure will not again be
368  * a candidate for frequency or offset correction until sufficient data
369  * has been accumulated to make a decision.
370  */
371 void
372 lin_reset(server_info_t info)
373 {
374     info->lin_count = 0;
375     info->lin_sumx = 0;
376     info->lin_sumy = 0;
377     info->lin_sumxy = 0;
378     info->lin_sumx2 = 0;
379     info->lin_sumy2 = 0;
380
381     info->lin_countoffset = 0;
382     info->lin_sumoffset = 0;
383     info->lin_sumoffset2 = 0;
384
385     info->lin_cache_slope = 0;
386     info->lin_cache_yint = 0;
387     info->lin_cache_corr = 0;
388     info->lin_cache_offset = 0;
389     info->lin_cache_freq = 0;
390 }
391
392 /*
393  * Sometimes we want to clean out the offset calculations without
394  * destroying the linear regression used to figure out the frequency
395  * correction.  This usually occurs whenever we issue an offset
396  * adjustment to the system, which invalidates any offset data accumulated
397  * up to that point.
398  */
399 void
400 lin_resetalloffsets(struct server_info **info_ary, int count)
401 {
402     server_info_t info;
403     int i;
404
405     for (i = 0; i < count; ++i) {
406         for (info = info_ary[i]; info; info = info->altinfo)
407             lin_resetoffsets(info);
408     }
409 }
410
411 void
412 lin_resetoffsets(server_info_t info)
413 {
414     info->lin_countoffset = 0;
415     info->lin_sumoffset = 0;
416     info->lin_sumoffset2 = 0;
417 }
418