From e591e3622ba5bbe88e797e9c70572ec8390a0829 Mon Sep 17 00:00:00 2001 From: Miroslav Lichvar Date: Tue, 17 Aug 2010 16:40:25 +0200 Subject: [PATCH] Extend runs test Double the number of samples that are used in the runs test. E.g. with 64 samples in regression the runs test will be tried over the 64 samples and up to 64 previous samples. The minimum number of samples in now 4. This improves the response with low-mid jitters by about 50%. --- regress.c | 27 +++--- regress.h | 7 ++ sourcestats.c | 233 ++++++++++++++++++++++++++------------------------ 3 files changed, 146 insertions(+), 121 deletions(-) diff --git a/regress.c b/regress.c index 9f85090..6913fed 100644 --- a/regress.c +++ b/regress.c @@ -226,6 +226,9 @@ RGR_FindBestRegression less reliable) */ int n, /* number of data points */ + int m, /* number of extra samples in x and y arrays + (negative index) which can be used to + extend runs test */ /* And now the results */ @@ -249,15 +252,15 @@ RGR_FindBestRegression ) { double P, Q, U, V, W; /* total */ - double resid[MAX_POINTS]; + double resid[MAX_POINTS * REGRESS_RUNS_RATIO]; double ss; double a, b, u, ui, aa; - int start, nruns, npoints, npoints_left; + int start, resid_start, nruns, npoints; int i; assert(n <= MAX_POINTS); - assert(MAX_POINTS < sizeof (critical_runs10) / sizeof (critical_runs10[0])); + assert(n * REGRESS_RUNS_RATIO < sizeof (critical_runs10) / sizeof (critical_runs10[0])); if (n < MIN_SAMPLES_FOR_REGRESS) { return 0; @@ -282,20 +285,22 @@ RGR_FindBestRegression V += ui * ui / w[i]; } - npoints = n - start; b = Q / V; a = (P / W) - (b * u); - for (i=start; i critical_runs10[npoints]) || (npoints_left < MIN_SAMPLES_FOR_REGRESS)) { + if (nruns > critical_runs10[n - resid_start] || n - start <= MIN_SAMPLES_FOR_REGRESS) { break; } else { /* Try dropping one sample at a time until the runs test passes. */ @@ -310,7 +315,7 @@ RGR_FindBestRegression ss = 0.0; for (i=start; i local fast) at a particular time */ @@ -101,7 +104,7 @@ struct SST_Stats_Record { /* This array contains the sample epochs, in terms of the local clock. */ - struct timeval sample_times[MAX_SAMPLES]; + struct timeval sample_times[MAX_SAMPLES * REGRESS_RUNS_RATIO]; /* This is an array of offsets, in seconds, corresponding to the sample times. In this module, we use the convention that @@ -109,7 +112,7 @@ struct SST_Stats_Record { means it is SLOW. This is contrary to the convention in the NTP stuff; that part of the code is written to correspond with RFC1305 conventions. */ - double offsets[MAX_SAMPLES]; + double offsets[MAX_SAMPLES * REGRESS_RUNS_RATIO]; /* This is an array of the offsets as originally measured. Local clock fast of real time is indicated by positive values. This @@ -167,6 +170,7 @@ SST_CreateInstance(unsigned long refid, IPAddr *addr) inst->refid = refid; inst->ip_addr = addr; inst->n_samples = 0; + inst->runs_samples = 0; inst->last_sample = 0; inst->estimated_frequency = 0; inst->skew = 2000.0e-6; @@ -199,6 +203,11 @@ prune_register(SST_Stats inst, int new_oldest) { assert(inst->n_samples >= new_oldest); inst->n_samples -= new_oldest; + inst->runs_samples += new_oldest; + if (inst->runs_samples > inst->n_samples * (REGRESS_RUNS_RATIO - 1)) + inst->runs_samples = inst->n_samples * (REGRESS_RUNS_RATIO - 1); + + assert(inst->n_samples + inst->runs_samples <= MAX_SAMPLES * REGRESS_RUNS_RATIO); } /* ================================================== */ @@ -210,33 +219,47 @@ SST_AccumulateSample(SST_Stats inst, struct timeval *sample_time, double root_delay, double root_dispersion, int stratum) { - int n; + int n, m; if (inst->n_samples == MAX_SAMPLES) { prune_register(inst, 1); } - n = inst->last_sample = (inst->last_sample + 1) % MAX_SAMPLES; + n = inst->last_sample = (inst->last_sample + 1) % + (MAX_SAMPLES * REGRESS_RUNS_RATIO); + m = n % MAX_SAMPLES; inst->sample_times[n] = *sample_time; inst->offsets[n] = offset; - inst->orig_offsets[n] = offset; - inst->peer_delays[n] = peer_delay; - inst->peer_dispersions[n] = peer_dispersion; - inst->root_delays[n] = root_delay; - inst->root_dispersions[n] = root_dispersion; - inst->strata[n] = stratum; + inst->orig_offsets[m] = offset; + inst->peer_delays[m] = peer_delay; + inst->peer_dispersions[m] = peer_dispersion; + inst->root_delays[m] = root_delay; + inst->root_dispersions[m] = root_dispersion; + inst->strata[m] = stratum; ++inst->n_samples; } /* ================================================== */ -/* Return index of the i-th sample in the circular buffer */ +/* Return index of the i-th sample in the sample_times and offset buffers, + i can be negative down to -runs_samples */ + +static int +get_runsbuf_index(SST_Stats inst, int i) +{ + return (unsigned int)(inst->last_sample + 2 * MAX_SAMPLES * REGRESS_RUNS_RATIO - + inst->n_samples + i + 1) % (MAX_SAMPLES * REGRESS_RUNS_RATIO); +} + +/* ================================================== */ +/* Return index of the i-th sample in the other buffers */ static int get_buf_index(SST_Stats inst, int i) { - return (unsigned int)(inst->last_sample + MAX_SAMPLES - inst->n_samples + i + 1) % MAX_SAMPLES; + return (unsigned int)(inst->last_sample + MAX_SAMPLES * REGRESS_RUNS_RATIO - + inst->n_samples + i + 1) % MAX_SAMPLES; } /* ================================================== */ @@ -251,10 +274,10 @@ convert_to_intervals(SST_Stats inst, double *times_back) int i; newest_tv = &(inst->sample_times[inst->last_sample]); - for (i=0; in_samples; i++) { + for (i = -inst->runs_samples; i < inst->n_samples; i++) { /* The entries in times_back[] should end up negative */ UTI_DiffTimevalsToDouble(×_back[i], - &inst->sample_times[get_buf_index(inst, i)], newest_tv); + &inst->sample_times[get_runsbuf_index(inst, i)], newest_tv); } } @@ -268,42 +291,25 @@ find_best_sample_index(SST_Stats inst, double *times_back) double root_distance, best_root_distance; double elapsed; - int i, j, l, n, best_index; + int i, j, best_index; - n = inst->n_samples - 1; - best_index = l = inst->last_sample; - best_root_distance = inst->root_dispersions[l] + 0.5 * fabs(inst->root_delays[l]); -#if 0 - LOG(LOGS_INFO, LOGF_SourceStats, "n=%d brd=%f", n, best_root_distance); -#endif - for (i=0; in_samples; i++) { j = get_buf_index(inst, i); + elapsed = -times_back[i]; -#if 0 - LOG(LOGS_INFO, LOGF_SourceStats, "n=%d i=%d latest=[%s] doing=[%s] elapsed=%f", n, i, - UTI_TimevalToString(&inst->sample_times[l]), - UTI_TimevalToString(&inst->sample_times[j]), - elapsed); -#endif - - /* Because the loop does not consider the most recent sample, this assertion must hold */ - if (elapsed <= 0.0) { - LOG(LOGS_ERR, LOGF_SourceStats, "Elapsed<0! n=%d i=%d latest=[%s] doing=[%s] elapsed=%f", - n, i, - UTI_TimevalToString(&inst->sample_times[l]), - UTI_TimevalToString(&inst->sample_times[j]), - elapsed); - - elapsed = fabs(elapsed); - } + assert(elapsed >= 0.0); root_distance = inst->root_dispersions[j] + elapsed * inst->skew + 0.5 * fabs(inst->root_delays[j]); if (root_distance < best_root_distance) { best_root_distance = root_distance; - best_index = j; + best_index = i; } } - + + assert(best_index >= 0); inst->best_single_sample = best_index; #if 0 @@ -331,13 +337,13 @@ find_best_sample_index(SST_Stats inst, double *times_back) void SST_DoNewRegression(SST_Stats inst) { - double times_back[MAX_SAMPLES]; - double offsets[MAX_SAMPLES]; + double times_back[MAX_SAMPLES * REGRESS_RUNS_RATIO]; + double offsets[MAX_SAMPLES * REGRESS_RUNS_RATIO]; double peer_distances[MAX_SAMPLES]; double weights[MAX_SAMPLES]; int degrees_of_freedom; - int best_start; + int best_start, times_back_start; double est_intercept, est_slope, est_var, est_intercept_sd, est_slope_sd; int i, j, nruns; double min_distance; @@ -346,17 +352,16 @@ SST_DoNewRegression(SST_Stats inst) int regression_ok; - convert_to_intervals(inst, times_back); + convert_to_intervals(inst, times_back + inst->runs_samples); if (inst->n_samples > 0) { - for (i=0; in_samples; i++) { - j = get_buf_index(inst, i); - offsets[i] = inst->offsets[j]; - peer_distances[i] = 0.5 * fabs(inst->peer_delays[j]) + inst->peer_dispersions[j]; + for (i = -inst->runs_samples; i < inst->n_samples; i++) { + offsets[i + inst->runs_samples] = inst->offsets[get_runsbuf_index(inst, i)]; } - min_distance = peer_distances[0]; - for (i=1; in_samples; i++) { + for (i = 0, min_distance = DBL_MAX; i < inst->n_samples; i++) { + j = get_buf_index(inst, i); + peer_distances[i] = 0.5 * fabs(inst->peer_delays[j]) + inst->peer_dispersions[j]; if (peer_distances[i] < min_distance) { min_distance = peer_distances[i]; } @@ -370,8 +375,9 @@ SST_DoNewRegression(SST_Stats inst) } } - regression_ok = RGR_FindBestRegression(times_back, offsets, weights, - inst->n_samples, + regression_ok = RGR_FindBestRegression(times_back + inst->runs_samples, + offsets + inst->runs_samples, weights, + inst->n_samples, inst->runs_samples, &est_intercept, &est_slope, &est_var, &est_intercept_sd, &est_slope_sd, &best_start, &nruns, °rees_of_freedom); @@ -417,18 +423,18 @@ SST_DoNewRegression(SST_Stats inst) best_start, nruns); } + times_back_start = inst->runs_samples + best_start; prune_register(inst, best_start); - } else { #if 0 LOG(LOGS_INFO, LOGF_SourceStats, "too few points (%d) for regression", inst->n_samples); #endif inst->estimated_frequency = 0.0; inst->skew = WORST_CASE_FREQ_BOUND; - best_start = 0; + times_back_start = 0; } - find_best_sample_index(inst, times_back + best_start); + find_best_sample_index(inst, times_back + times_back_start); } @@ -442,22 +448,23 @@ SST_GetReferenceData(SST_Stats inst, struct timeval *now, { double elapsed; - int n; + int i, j; *frequency = inst->estimated_frequency; *skew = inst->skew; - n = inst->best_single_sample; + i = get_runsbuf_index(inst, inst->best_single_sample); + j = get_buf_index(inst, inst->best_single_sample); - UTI_DiffTimevalsToDouble(&elapsed, now, &(inst->sample_times[n])); - *root_delay = inst->root_delays[n]; - *root_dispersion = inst->root_dispersions[n] + elapsed * inst->skew; - *offset = inst->offsets[n] + elapsed * inst->estimated_frequency; - *stratum = inst->strata[n]; + UTI_DiffTimevalsToDouble(&elapsed, now, &inst->sample_times[i]); + *root_delay = inst->root_delays[j]; + *root_dispersion = inst->root_dispersions[j] + elapsed * inst->skew; + *offset = inst->offsets[i] + elapsed * inst->estimated_frequency; + *stratum = inst->strata[j]; #ifdef TRACEON LOG(LOGS_INFO, LOGF_SourceStats, "n=%d freq=%f skew=%f del=%f disp=%f ofs=%f str=%d", - n, *frequency, *skew, *root_delay, *root_dispersion, *offset, *stratum); + inst->n_samples, *frequency, *skew, *root_delay, *root_dispersion, *offset, *stratum); #endif return; @@ -493,20 +500,22 @@ SST_GetSelectionData(SST_Stats inst, struct timeval *now, double average_offset; double sample_elapsed; double elapsed; - int n; + int i, j; double peer_distance; - n = inst->best_single_sample; - *stratum = inst->strata[n]; + i = get_runsbuf_index(inst, inst->best_single_sample); + j = get_buf_index(inst, inst->best_single_sample); + + *stratum = inst->strata[j]; *variance = inst->variance; - peer_distance = inst->peer_dispersions[n] + 0.5 * fabs(inst->peer_delays[n]); + peer_distance = inst->peer_dispersions[j] + 0.5 * fabs(inst->peer_delays[j]); UTI_DiffTimevalsToDouble(&elapsed, now, &(inst->offset_time)); - UTI_DiffTimevalsToDouble(&sample_elapsed, now, &(inst->sample_times[n])); - *best_offset = inst->offsets[n] + sample_elapsed * inst->estimated_frequency; - *best_root_delay = inst->root_delays[n]; - *best_root_dispersion = inst->root_dispersions[n] + sample_elapsed * inst->skew; + UTI_DiffTimevalsToDouble(&sample_elapsed, now, &inst->sample_times[i]); + *best_offset = inst->offsets[i] + sample_elapsed * inst->estimated_frequency; + *best_root_delay = inst->root_delays[j]; + *best_root_dispersion = inst->root_dispersions[j] + sample_elapsed * inst->skew; average_offset = inst->estimated_offset + inst->estimated_frequency * elapsed; if (fabs(average_offset - *best_offset) <= peer_distance) { @@ -517,7 +526,7 @@ SST_GetSelectionData(SST_Stats inst, struct timeval *now, #ifdef TRACEON LOG(LOGS_INFO, LOGF_SourceStats, "n=%d off=%f del=%f dis=%f var=%f pdist=%f avoff=%f avok=%d", - n, *best_offset, *best_root_delay, *best_root_dispersion, *variance, + inst->n_samples, *best_offset, *best_root_delay, *best_root_dispersion, *variance, peer_distance, average_offset, *average_ok); #endif @@ -532,26 +541,27 @@ SST_GetTrackingData(SST_Stats inst, struct timeval *now, double *accrued_dispersion, double *frequency, double *skew) { - int n; + int i, j; double peer_distance; double elapsed_offset, elapsed_sample; - n = inst->best_single_sample; + i = get_runsbuf_index(inst, inst->best_single_sample); + j = get_buf_index(inst, inst->best_single_sample); *frequency = inst->estimated_frequency; *skew = inst->skew; - peer_distance = inst->peer_dispersions[n] + 0.5 * fabs(inst->peer_delays[n]); + peer_distance = inst->peer_dispersions[j] + 0.5 * fabs(inst->peer_delays[j]); UTI_DiffTimevalsToDouble(&elapsed_offset, now, &(inst->offset_time)); *average_offset = inst->estimated_offset + inst->estimated_frequency * elapsed_offset; *offset_sd = inst->estimated_offset_sd + elapsed_offset * inst->skew; - UTI_DiffTimevalsToDouble(&elapsed_sample, now, &(inst->sample_times[n])); + UTI_DiffTimevalsToDouble(&elapsed_sample, now, &inst->sample_times[i]); *accrued_dispersion = inst->skew * elapsed_sample; #ifdef TRACEON LOG(LOGS_INFO, LOGF_SourceStats, "n=%d freq=%f (%.3fppm) skew=%f (%.3fppm) pdist=%f avoff=%f offsd=%f accrdis=%f", - n, *frequency, 1.0e6* *frequency, *skew, 1.0e6* *skew, peer_distance, *average_offset, *offset_sd, *accrued_dispersion); + inst->n_samples, *frequency, 1.0e6* *frequency, *skew, 1.0e6* *skew, peer_distance, *average_offset, *offset_sd, *accrued_dispersion); #endif } @@ -561,15 +571,13 @@ SST_GetTrackingData(SST_Stats inst, struct timeval *now, void SST_SlewSamples(SST_Stats inst, struct timeval *when, double dfreq, double doffset) { - int m, n, i; + int m, i; double delta_time; struct timeval *sample, prev; double prev_offset, prev_freq; - n = inst->n_samples; - - for (m = 0; m < n; m++) { - i = get_buf_index(inst, m); + for (m = -inst->runs_samples; m < inst->n_samples; m++) { + i = get_runsbuf_index(inst, m); sample = &(inst->sample_times[i]); prev = *sample; UTI_AdjustTimeval(sample, when, sample, &delta_time, dfreq, doffset); @@ -661,24 +669,25 @@ SST_MinRoundTripDelay(SST_Stats inst) void SST_SaveToFile(SST_Stats inst, FILE *out) { - int m, i; + int m, i, j; fprintf(out, "%d\n", inst->n_samples); for(m = 0; m < inst->n_samples; m++) { - i = get_buf_index(inst, m); + i = get_runsbuf_index(inst, m); + j = get_buf_index(inst, m); fprintf(out, "%08lx %08lx %.6e %.6e %.6e %.6e %.6e %.6e %.6e %d\n", (unsigned long) inst->sample_times[i].tv_sec, (unsigned long) inst->sample_times[i].tv_usec, inst->offsets[i], - inst->orig_offsets[i], - inst->peer_delays[i], - inst->peer_dispersions[i], - inst->root_delays[i], - inst->root_dispersions[i], + inst->orig_offsets[j], + inst->peer_delays[j], + inst->peer_dispersions[j], + inst->root_delays[j], + inst->root_dispersions[j], 1.0, /* used to be inst->weights[i] */ - inst->strata[i]); + inst->strata[j]); } } @@ -733,6 +742,7 @@ SST_LoadFromFile(SST_Stats inst, FILE *in) } inst->last_sample = inst->n_samples - 1; + inst->runs_samples = 0; return 1; @@ -743,17 +753,18 @@ SST_LoadFromFile(SST_Stats inst, FILE *in) void SST_DoSourceReport(SST_Stats inst, RPT_SourceReport *report, struct timeval *now) { - int n; + int i, j; struct timeval ago; if (inst->n_samples > 0) { - n = inst->last_sample; - report->orig_latest_meas = inst->orig_offsets[n]; - report->latest_meas = inst->offsets[n]; - report->latest_meas_err = 0.5*inst->root_delays[n] + inst->root_dispersions[n]; - report->stratum = inst->strata[n]; + i = get_runsbuf_index(inst, inst->n_samples - 1); + j = get_buf_index(inst, inst->n_samples - 1); + report->orig_latest_meas = inst->orig_offsets[j]; + report->latest_meas = inst->offsets[i]; + report->latest_meas_err = 0.5*inst->root_delays[j] + inst->root_dispersions[j]; + report->stratum = inst->strata[j]; - UTI_DiffTimevals(&ago, now, &inst->sample_times[n]); + UTI_DiffTimevals(&ago, now, &inst->sample_times[i]); report->latest_meas_ago = ago.tv_sec; } else { report->latest_meas_ago = 86400 * 365 * 10; @@ -779,28 +790,30 @@ SST_DoSourcestatsReport(SST_Stats inst, RPT_SourcestatsReport *report, struct ti { double dspan; double elapsed, sample_elapsed; - int n, nb; + int li, lj, bi, bj; report->n_samples = inst->n_samples; report->n_runs = inst->nruns; if (inst->n_samples > 1) { - n = inst->last_sample; - UTI_DiffTimevalsToDouble(&dspan, &inst->sample_times[n], - &inst->sample_times[get_buf_index(inst, 0)]); + li = get_runsbuf_index(inst, inst->n_samples - 1); + lj = get_buf_index(inst, inst->n_samples - 1); + UTI_DiffTimevalsToDouble(&dspan, &inst->sample_times[li], + &inst->sample_times[get_runsbuf_index(inst, 0)]); report->span_seconds = (unsigned long) (dspan + 0.5); if (inst->n_samples > 3) { UTI_DiffTimevalsToDouble(&elapsed, now, &inst->offset_time); - nb = inst->best_single_sample; - UTI_DiffTimevalsToDouble(&sample_elapsed, now, &(inst->sample_times[nb])); + bi = get_runsbuf_index(inst, inst->best_single_sample); + bj = get_buf_index(inst, inst->best_single_sample); + UTI_DiffTimevalsToDouble(&sample_elapsed, now, &inst->sample_times[bi]); report->est_offset = inst->estimated_offset + elapsed * inst->estimated_frequency; report->est_offset_err = (inst->estimated_offset_sd + sample_elapsed * inst->skew + - (0.5*inst->root_delays[nb] + inst->root_dispersions[nb])); + (0.5*inst->root_delays[bj] + inst->root_dispersions[bj])); } else { - report->est_offset = inst->offsets[n]; - report->est_offset_err = 0.5*inst->root_delays[n] + inst->root_dispersions[n]; + report->est_offset = inst->offsets[li]; + report->est_offset_err = 0.5*inst->root_delays[lj] + inst->root_dispersions[lj]; } } else { report->span_seconds = 0;