[thirdparty/chrony.git] / samplefilt.c

/*
  chronyd/chronyc - Programs for keeping computer clocks accurate.

 **********************************************************************
 * Copyright (C) Miroslav Lichvar  2009-2011, 2014, 2016, 2018
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of version 2 of the GNU General Public License as
 * published by the Free Software Foundation.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License along
 * with this program; if not, write to the Free Software Foundation, Inc.,
 * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 * 
 **********************************************************************

  =======================================================================

  Routines implementing a median sample filter.

  */

#include "config.h"

#include "local.h"
#include "logging.h"
#include "memory.h"
#include "regress.h"
#include "samplefilt.h"
#include "util.h"

#define MIN_SAMPLES 1
#define MAX_SAMPLES 256

struct SPF_Instance_Record {
  int min_samples;
  int max_samples;
  int index;
  int used;
  int last;
  int avg_var_n;
  double avg_var;
  double max_var;
  double combine_ratio;
  NTP_Sample *samples;
  int *selected;
  double *x_data;
  double *y_data;
  double *w_data;
};

/* ================================================== */

SPF_Instance
SPF_CreateInstance(int min_samples, int max_samples, double max_dispersion, double combine_ratio)
{
  SPF_Instance filter;

  filter = MallocNew(struct SPF_Instance_Record);

  min_samples = CLAMP(MIN_SAMPLES, min_samples, MAX_SAMPLES);
  max_samples = CLAMP(MIN_SAMPLES, max_samples, MAX_SAMPLES);
  max_samples = MAX(min_samples, max_samples);
  combine_ratio = CLAMP(0.0, combine_ratio, 1.0);

  filter->min_samples = min_samples;
  filter->max_samples = max_samples;
  filter->index = -1;
  filter->used = 0;
  filter->last = -1;
  /* Set the first estimate to the system precision */
  filter->avg_var_n = 0;
  filter->avg_var = SQUARE(LCL_GetSysPrecisionAsQuantum());
  filter->max_var = SQUARE(max_dispersion);
  filter->combine_ratio = combine_ratio;
  filter->samples = MallocArray(NTP_Sample, filter->max_samples);
  filter->selected = MallocArray(int, filter->max_samples);
  filter->x_data = MallocArray(double, filter->max_samples);
  filter->y_data = MallocArray(double, filter->max_samples);
  filter->w_data = MallocArray(double, filter->max_samples);

  return filter;
}

/* ================================================== */

void
SPF_DestroyInstance(SPF_Instance filter)
{
  Free(filter->samples);
  Free(filter->selected);
  Free(filter->x_data);
  Free(filter->y_data);
  Free(filter->w_data);
  Free(filter);
}

/* ================================================== */

/* Check that samples times are strictly increasing */

static int
check_sample(SPF_Instance filter, NTP_Sample *sample)
{
  if (filter->used <= 0)
    return 1;

  if (UTI_CompareTimespecs(&filter->samples[filter->last].time, &sample->time) >= 0) {
    DEBUG_LOG("filter non-increasing sample time %s", UTI_TimespecToString(&sample->time));
    return 0;
  }

  return 1;
}

/* ================================================== */

int
SPF_AccumulateSample(SPF_Instance filter, NTP_Sample *sample)
{
  if (!check_sample(filter, sample))
      return 0;

  filter->index++;
  filter->index %= filter->max_samples;
  filter->last = filter->index;
  if (filter->used < filter->max_samples)
    filter->used++;

  filter->samples[filter->index] = *sample;

  DEBUG_LOG("filter sample %d t=%s offset=%.9f peer_disp=%.9f",
            filter->index, UTI_TimespecToString(&sample->time),
            sample->offset, sample->peer_dispersion);
  return 1;
}

/* ================================================== */

int
SPF_GetLastSample(SPF_Instance filter, NTP_Sample *sample)
{
  if (filter->last < 0)
    return 0;

  *sample = filter->samples[filter->last];
  return 1;
}

/* ================================================== */

int
SPF_GetNumberOfSamples(SPF_Instance filter)
{
  return filter->used;
}

/* ================================================== */

int
SPF_GetMaxSamples(SPF_Instance filter)
{
  return filter->max_samples;
}

/* ================================================== */

double
SPF_GetAvgSampleDispersion(SPF_Instance filter)
{
  return sqrt(filter->avg_var);
}

/* ================================================== */

static void
drop_samples(SPF_Instance filter, int keep_last)
{
  filter->index = -1;
  filter->used = 0;
  if (!keep_last)
    filter->last = -1;
}

/* ================================================== */

void
SPF_DropSamples(SPF_Instance filter)
{
  drop_samples(filter, 0);
}

/* ================================================== */

static const NTP_Sample *tmp_sort_samples;

static int
compare_samples(const void *a, const void *b)
{
  const NTP_Sample *s1, *s2;

  s1 = &tmp_sort_samples[*(int *)a];
  s2 = &tmp_sort_samples[*(int *)b];

  if (s1->offset < s2->offset)
    return -1;
  else if (s1->offset > s2->offset)
    return 1;
  return 0;
}

/* ================================================== */

static int
select_samples(SPF_Instance filter)
{
  int i, j, k, o, from, to, *selected;
  double min_dispersion;

  if (filter->used < filter->min_samples)
    return 0;

  selected = filter->selected;

  /* With 4 or more samples, select those that have peer dispersion smaller
     than 1.5x of the minimum dispersion */
  if (filter->used > 4) {
    for (i = 1, min_dispersion = filter->samples[0].peer_dispersion; i < filter->used; i++) {
      if (min_dispersion > filter->samples[i].peer_dispersion)
        min_dispersion = filter->samples[i].peer_dispersion;
    }

    for (i = j = 0; i < filter->used; i++) {
      if (filter->samples[i].peer_dispersion <= 1.5 * min_dispersion)
        selected[j++] = i;
    }
  } else {
    j = 0;
  }

  if (j < 4) {
    /* Select all samples */

    for (j = 0; j < filter->used; j++)
      selected[j] = j;
  }

  /* And sort their indices by offset */
  tmp_sort_samples = filter->samples;
  qsort(selected, j, sizeof (int), compare_samples);

  /* Select samples closest to the median */
  if (j > 2) {
    from = j * (1.0 - filter->combine_ratio) / 2.0;
    from = CLAMP(1, from, (j - 1) / 2);
  } else {
    from = 0;
  }

  to = j - from;

  /* Mark unused samples and sort the rest by their time */

  o = filter->used - filter->index - 1;

  for (i = 0; i < from; i++)
    selected[i] = -1;
  for (; i < to; i++)
    selected[i] = (selected[i] + o) % filter->used;
  for (; i < filter->used; i++)
    selected[i] = -1;

  for (i = from; i < to; i++) {
    j = selected[i];
    selected[i] = -1;
    while (j != -1 && selected[j] != j) {
      k = selected[j];
      selected[j] = j;
      j = k;
    }
  }

  for (i = j = 0; i < filter->used; i++) {
    if (selected[i] != -1)
      selected[j++] = (selected[i] + filter->used - o) % filter->used;
  }

  assert(j > 0 && j <= filter->max_samples);

  return j;
}

/* ================================================== */

static int
combine_selected_samples(SPF_Instance filter, int n, NTP_Sample *result)
{
  double mean_peer_dispersion, mean_root_dispersion, mean_peer_delay, mean_root_delay;
  double mean_x, mean_y, disp, var, prev_avg_var;
  NTP_Sample *sample, *last_sample;
  int i, dof;

  last_sample = &filter->samples[filter->selected[n - 1]];

  /* Prepare data */
  for (i = 0; i < n; i++) {
    sample = &filter->samples[filter->selected[i]];

    filter->x_data[i] = UTI_DiffTimespecsToDouble(&sample->time, &last_sample->time);
    filter->y_data[i] = sample->offset;
    filter->w_data[i] = sample->peer_dispersion;
  }

  /* Calculate mean offset and interval since the last sample */
  for (i = 0, mean_x = mean_y = 0.0; i < n; i++) {
    mean_x += filter->x_data[i];
    mean_y += filter->y_data[i];
  }
  mean_x /= n;
  mean_y /= n;

  if (n >= 4) {
    double b0, b1, s2, sb0, sb1;

    /* Set y axis to the mean sample time */
    for (i = 0; i < n; i++)
      filter->x_data[i] -= mean_x;

    /* Make a linear fit and use the estimated standard deviation of the
       intercept as dispersion */
    RGR_WeightedRegression(filter->x_data, filter->y_data, filter->w_data, n,
                           &b0, &b1, &s2, &sb0, &sb1);
    var = s2;
    disp = sb0;
    dof = n - 2;
  } else if (n >= 2) {
    for (i = 0, disp = 0.0; i < n; i++)
      disp += (filter->y_data[i] - mean_y) * (filter->y_data[i] - mean_y);
    var = disp / (n - 1);
    disp = sqrt(var);
    dof = n - 1;
  } else {
    var = filter->avg_var;
    disp = sqrt(var);
    dof = 1;
  }

  /* Avoid working with zero dispersion */
  if (var < 1e-20) {
    var = 1e-20;
    disp = sqrt(var);
  }

  /* Drop the sample if the variance is larger than the maximum */
  if (filter->max_var > 0.0 && var > filter->max_var) {
    DEBUG_LOG("filter dispersion too large disp=%.9f max=%.9f",
              sqrt(var), sqrt(filter->max_var));
    return 0;
  }

  prev_avg_var = filter->avg_var;

  /* Update the exponential moving average of the variance */
  if (filter->avg_var_n > 50) {
    filter->avg_var += dof / (dof + 50.0) * (var - filter->avg_var);
  } else {
    filter->avg_var = (filter->avg_var * filter->avg_var_n + var * dof) /
      (dof + filter->avg_var_n);
    if (filter->avg_var_n == 0)
      prev_avg_var = filter->avg_var;
    filter->avg_var_n += dof;
  }

  /* Use the long-term average of variance instead of the estimated value
     unless it is significantly smaller in order to reduce the noise in
     sourcestats weights */
  if (var * dof / RGR_GetChi2Coef(dof) < prev_avg_var)
    disp = sqrt(filter->avg_var) * disp / sqrt(var);

  mean_peer_dispersion = mean_root_dispersion = mean_peer_delay = mean_root_delay = 0.0;

  for (i = 0; i < n; i++) {
    sample = &filter->samples[filter->selected[i]];

    mean_peer_dispersion += sample->peer_dispersion;
    mean_root_dispersion += sample->root_dispersion;
    mean_peer_delay += sample->peer_delay;
    mean_root_delay += sample->root_delay;
  }

  mean_peer_dispersion /= n;
  mean_root_dispersion /= n;
  mean_peer_delay /= n;
  mean_root_delay /= n;

  UTI_AddDoubleToTimespec(&last_sample->time, mean_x, &result->time);
  result->offset = mean_y;
  result->peer_dispersion = MAX(disp, mean_peer_dispersion);
  result->root_dispersion = MAX(disp, mean_root_dispersion);
  result->peer_delay = mean_peer_delay;
  result->root_delay = mean_root_delay;

  return 1;
}

/* ================================================== */

int
SPF_GetFilteredSample(SPF_Instance filter, NTP_Sample *sample)
{
  int n;

  n = select_samples(filter);

  DEBUG_LOG("selected %d from %d samples", n, filter->used);

  if (n < 1)
    return 0;

  if (!combine_selected_samples(filter, n, sample))
    return 0;

  drop_samples(filter, 1);

  return 1;
}

/* ================================================== */

static int
get_first_last(SPF_Instance filter, int *first, int *last)
{
  if (filter->last < 0)
    return 0;

  /* Always slew the last sample as it may be returned even if no new
     samples were accumulated */
  if (filter->used > 0) {
    *first = 0;
    *last = filter->used - 1;
  } else {
    *first = *last = filter->last;
  }

  return 1;
}


/* ================================================== */

void
SPF_SlewSamples(SPF_Instance filter, struct timespec *when, double dfreq, double doffset)
{
  int i, first, last;
  double delta_time;

  if (!get_first_last(filter, &first, &last))
    return;

  for (i = first; i <= last; i++) {
    UTI_AdjustTimespec(&filter->samples[i].time, when, &filter->samples[i].time,
                       &delta_time, dfreq, doffset);
    filter->samples[i].offset -= delta_time;
  }
}

/* ================================================== */

void
SPF_CorrectOffset(SPF_Instance filter, double doffset)
{
  int i, first, last;

  if (!get_first_last(filter, &first, &last))
    return;

  for (i = first; i <= last; i++)
    filter->samples[i].offset -= doffset;
}

/* ================================================== */

void
SPF_AddDispersion(SPF_Instance filter, double dispersion)
{
  int i;

  for (i = 0; i < filter->used; i++) {
    filter->samples[i].peer_dispersion += dispersion;
    filter->samples[i].root_dispersion += dispersion;
  }
}
Commit	Line	Data
c498c21f ML	1	/*
	2	chronyd/chronyc - Programs for keeping computer clocks accurate.
	3
	4	**********************************************************************
	5	* Copyright (C) Miroslav Lichvar 2009-2011, 2014, 2016, 2018
	6	*
	7	* This program is free software; you can redistribute it and/or modify
	8	* it under the terms of version 2 of the GNU General Public License as
	9	* published by the Free Software Foundation.
	10	*
	11	* This program is distributed in the hope that it will be useful, but
	12	* WITHOUT ANY WARRANTY; without even the implied warranty of
	13	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
	14	* General Public License for more details.
	15	*
	16	* You should have received a copy of the GNU General Public License along
	17	* with this program; if not, write to the Free Software Foundation, Inc.,
	18	* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
	19	*
	20	**********************************************************************
	21
	22	=======================================================================
	23
	24	Routines implementing a median sample filter.
	25
	26	*/
	27
	28	#include "config.h"
	29
	30	#include "local.h"
	31	#include "logging.h"
	32	#include "memory.h"
	33	#include "regress.h"
	34	#include "samplefilt.h"
	35	#include "util.h"
	36
	37	#define MIN_SAMPLES 1
	38	#define MAX_SAMPLES 256
	39
	40	struct SPF_Instance_Record {
	41	int min_samples;
	42	int max_samples;
	43	int index;
	44	int used;
	45	int last;
	46	int avg_var_n;
	47	double avg_var;
	48	double max_var;
	49	double combine_ratio;
	50	NTP_Sample *samples;
	51	int *selected;
	52	double *x_data;
	53	double *y_data;
	54	double *w_data;
	55	};
	56
	57	/* ================================================== */
	58
	59	SPF_Instance
	60	SPF_CreateInstance(int min_samples, int max_samples, double max_dispersion, double combine_ratio)
	61	{
	62	SPF_Instance filter;
	63
	64	filter = MallocNew(struct SPF_Instance_Record);
65
66	min_samples = CLAMP(MIN_SAMPLES, min_samples, MAX_SAMPLES);
67	max_samples = CLAMP(MIN_SAMPLES, max_samples, MAX_SAMPLES);
68	max_samples = MAX(min_samples, max_samples);
69	combine_ratio = CLAMP(0.0, combine_ratio, 1.0);
70
71	filter->min_samples = min_samples;
72	filter->max_samples = max_samples;
73	filter->index = -1;
74	filter->used = 0;
75	filter->last = -1;
76	/* Set the first estimate to the system precision */
77	filter->avg_var_n = 0;
be3c1b52	78	filter->avg_var = SQUARE(LCL_GetSysPrecisionAsQuantum());
0b709ab1	79	filter->max_var = SQUARE(max_dispersion);
c498c21f ML	80	filter->combine_ratio = combine_ratio;
	81	filter->samples = MallocArray(NTP_Sample, filter->max_samples);
	82	filter->selected = MallocArray(int, filter->max_samples);
	83	filter->x_data = MallocArray(double, filter->max_samples);
	84	filter->y_data = MallocArray(double, filter->max_samples);
	85	filter->w_data = MallocArray(double, filter->max_samples);
	86
	87	return filter;
	88	}
	89
	90	/* ================================================== */
	91
	92	void
	93	SPF_DestroyInstance(SPF_Instance filter)
	94	{
	95	Free(filter->samples);
	96	Free(filter->selected);
	97	Free(filter->x_data);
	98	Free(filter->y_data);
	99	Free(filter->w_data);
	100	Free(filter);
	101	}
	102
	103	/* ================================================== */
	104
bba29a0e ML	105	/* Check that samples times are strictly increasing */
	106
	107	static int
	108	check_sample(SPF_Instance filter, NTP_Sample *sample)
	109	{
	110	if (filter->used <= 0)
	111	return 1;
	112
	113	if (UTI_CompareTimespecs(&filter->samples[filter->last].time, &sample->time) >= 0) {
	114	DEBUG_LOG("filter non-increasing sample time %s", UTI_TimespecToString(&sample->time));
	115	return 0;
	116	}
	117
	118	return 1;
	119	}
	120
	121	/* ================================================== */
	122
	123	int
c498c21f ML	124	SPF_AccumulateSample(SPF_Instance filter, NTP_Sample *sample)
c498c21f ML	125	{
bba29a0e ML	126	if (!check_sample(filter, sample))
	127	return 0;
	128
c498c21f ML	129	filter->index++;
	130	filter->index %= filter->max_samples;
	131	filter->last = filter->index;
	132	if (filter->used < filter->max_samples)
	133	filter->used++;
	134
	135	filter->samples[filter->index] = *sample;
	136
	137	DEBUG_LOG("filter sample %d t=%s offset=%.9f peer_disp=%.9f",
	138	filter->index, UTI_TimespecToString(&sample->time),
	139	sample->offset, sample->peer_dispersion);
bba29a0e	140	return 1;
c498c21f ML	141	}
	142
	143	/* ================================================== */
	144
	145	int
	146	SPF_GetLastSample(SPF_Instance filter, NTP_Sample *sample)
	147	{
	148	if (filter->last < 0)
	149	return 0;
	150
	151	*sample = filter->samples[filter->last];
	152	return 1;
	153	}
	154
	155	/* ================================================== */
	156
	157	int
	158	SPF_GetNumberOfSamples(SPF_Instance filter)
	159	{
	160	return filter->used;
	161	}
	162
	163	/* ================================================== */
	164
a4349b13 ML	165	int
	166	SPF_GetMaxSamples(SPF_Instance filter)
	167	{
	168	return filter->max_samples;
	169	}
	170
	171	/* ================================================== */
	172
c498c21f ML	173	double
	174	SPF_GetAvgSampleDispersion(SPF_Instance filter)
	175	{
	176	return sqrt(filter->avg_var);
	177	}
	178
	179	/* ================================================== */
	180
e66f1df8 ML	181	static void
e66f1df8 ML	182	drop_samples(SPF_Instance filter, int keep_last)
c498c21f ML	183	{
	184	filter->index = -1;
	185	filter->used = 0;
e66f1df8 ML	186	if (!keep_last)
	187	filter->last = -1;
	188	}
	189
	190	/* ================================================== */
	191
	192	void
	193	SPF_DropSamples(SPF_Instance filter)
	194	{
	195	drop_samples(filter, 0);
c498c21f ML	196	}
	197
	198	/* ================================================== */
	199
	200	static const NTP_Sample *tmp_sort_samples;
	201
	202	static int
	203	compare_samples(const void a, const void b)
	204	{
	205	const NTP_Sample s1, s2;
	206
	207	s1 = &tmp_sort_samples[(int )a];
	208	s2 = &tmp_sort_samples[(int )b];
	209
	210	if (s1->offset < s2->offset)
	211	return -1;
	212	else if (s1->offset > s2->offset)
	213	return 1;
	214	return 0;
	215	}
	216
	217	/* ================================================== */
	218
	219	static int
	220	select_samples(SPF_Instance filter)
	221	{
	222	int i, j, k, o, from, to, *selected;
	223	double min_dispersion;
	224
	225	if (filter->used < filter->min_samples)
	226	return 0;
	227
	228	selected = filter->selected;
	229
	230	/* With 4 or more samples, select those that have peer dispersion smaller
	231	than 1.5x of the minimum dispersion */
	232	if (filter->used > 4) {
	233	for (i = 1, min_dispersion = filter->samples[0].peer_dispersion; i < filter->used; i++) {
	234	if (min_dispersion > filter->samples[i].peer_dispersion)
	235	min_dispersion = filter->samples[i].peer_dispersion;
	236	}
	237
	238	for (i = j = 0; i < filter->used; i++) {
	239	if (filter->samples[i].peer_dispersion <= 1.5 * min_dispersion)
	240	selected[j++] = i;
	241	}
	242	} else {
	243	j = 0;
	244	}
	245
	246	if (j < 4) {
	247	/* Select all samples */
	248
	249	for (j = 0; j < filter->used; j++)
	250	selected[j] = j;
	251	}
	252
	253	/* And sort their indices by offset */
	254	tmp_sort_samples = filter->samples;
	255	qsort(selected, j, sizeof (int), compare_samples);
	256
	257	/* Select samples closest to the median */
	258	if (j > 2) {
	259	from = j * (1.0 - filter->combine_ratio) / 2.0;
260	from = CLAMP(1, from, (j - 1) / 2);
261	} else {
262	from = 0;
263	}
264
265	to = j - from;
266
267	/* Mark unused samples and sort the rest by their time */
268
269	o = filter->used - filter->index - 1;
270
271	for (i = 0; i < from; i++)
272	selected[i] = -1;
273	for (; i < to; i++)
274	selected[i] = (selected[i] + o) % filter->used;
275	for (; i < filter->used; i++)
276	selected[i] = -1;
277
278	for (i = from; i < to; i++) {
279	j = selected[i];
280	selected[i] = -1;
281	while (j != -1 && selected[j] != j) {
282	k = selected[j];
283	selected[j] = j;
284	j = k;
285	}
286	}
287
105b3faa	288	for (i = j = 0; i < filter->used; i++) {
c498c21f ML	289	if (selected[i] != -1)
	290	selected[j++] = (selected[i] + filter->used - o) % filter->used;
	291	}
	292
	293	assert(j > 0 && j <= filter->max_samples);
	294
	295	return j;
	296	}
	297
	298	/* ================================================== */
	299
	300	static int
	301	combine_selected_samples(SPF_Instance filter, int n, NTP_Sample *result)
	302	{
	303	double mean_peer_dispersion, mean_root_dispersion, mean_peer_delay, mean_root_delay;
	304	double mean_x, mean_y, disp, var, prev_avg_var;
	305	NTP_Sample sample, last_sample;
	306	int i, dof;
	307
	308	last_sample = &filter->samples[filter->selected[n - 1]];
	309
	310	/* Prepare data */
	311	for (i = 0; i < n; i++) {
	312	sample = &filter->samples[filter->selected[i]];
	313
	314	filter->x_data[i] = UTI_DiffTimespecsToDouble(&sample->time, &last_sample->time);
	315	filter->y_data[i] = sample->offset;
	316	filter->w_data[i] = sample->peer_dispersion;
	317	}
	318
	319	/* Calculate mean offset and interval since the last sample */
	320	for (i = 0, mean_x = mean_y = 0.0; i < n; i++) {
	321	mean_x += filter->x_data[i];
	322	mean_y += filter->y_data[i];
	323	}
	324	mean_x /= n;
	325	mean_y /= n;
	326
	327	if (n >= 4) {
	328	double b0, b1, s2, sb0, sb1;
	329
	330	/* Set y axis to the mean sample time */
	331	for (i = 0; i < n; i++)
	332	filter->x_data[i] -= mean_x;
	333
	334	/* Make a linear fit and use the estimated standard deviation of the
	335	intercept as dispersion */
	336	RGR_WeightedRegression(filter->x_data, filter->y_data, filter->w_data, n,
	337	&b0, &b1, &s2, &sb0, &sb1);
	338	var = s2;
	339	disp = sb0;
	340	dof = n - 2;
	341	} else if (n >= 2) {
	342	for (i = 0, disp = 0.0; i < n; i++)
	343	disp += (filter->y_data[i] - mean_y) * (filter->y_data[i] - mean_y);
	344	var = disp / (n - 1);
	345	disp = sqrt(var);
	346	dof = n - 1;
	347	} else {
	348	var = filter->avg_var;
	349	disp = sqrt(var);
	350	dof = 1;
	351	}
	352
353	/* Avoid working with zero dispersion */
354	if (var < 1e-20) {
355	var = 1e-20;
356	disp = sqrt(var);
357	}
358
359	/* Drop the sample if the variance is larger than the maximum */
360	if (filter->max_var > 0.0 && var > filter->max_var) {
361	DEBUG_LOG("filter dispersion too large disp=%.9f max=%.9f",
362	sqrt(var), sqrt(filter->max_var));
363	return 0;
364	}
365
366	prev_avg_var = filter->avg_var;
367
368	/* Update the exponential moving average of the variance */
369	if (filter->avg_var_n > 50) {
370	filter->avg_var += dof / (dof + 50.0) * (var - filter->avg_var);
371	} else {
372	filter->avg_var = (filter->avg_var * filter->avg_var_n + var * dof) /
373	(dof + filter->avg_var_n);
374	if (filter->avg_var_n == 0)
375	prev_avg_var = filter->avg_var;
376	filter->avg_var_n += dof;
377	}
378
379	/* Use the long-term average of variance instead of the estimated value
380	unless it is significantly smaller in order to reduce the noise in
381	sourcestats weights */
382	if (var * dof / RGR_GetChi2Coef(dof) < prev_avg_var)
383	disp = sqrt(filter->avg_var) * disp / sqrt(var);
384
385	mean_peer_dispersion = mean_root_dispersion = mean_peer_delay = mean_root_delay = 0.0;
386
387	for (i = 0; i < n; i++) {
388	sample = &filter->samples[filter->selected[i]];
389
390	mean_peer_dispersion += sample->peer_dispersion;
391	mean_root_dispersion += sample->root_dispersion;
392	mean_peer_delay += sample->peer_delay;
393	mean_root_delay += sample->root_delay;
394	}
395
396	mean_peer_dispersion /= n;
397	mean_root_dispersion /= n;
398	mean_peer_delay /= n;
399	mean_root_delay /= n;
400
401	UTI_AddDoubleToTimespec(&last_sample->time, mean_x, &result->time);
402	result->offset = mean_y;
403	result->peer_dispersion = MAX(disp, mean_peer_dispersion);
404	result->root_dispersion = MAX(disp, mean_root_dispersion);
405	result->peer_delay = mean_peer_delay;
406	result->root_delay = mean_root_delay;
c498c21f ML	407
	408	return 1;
	409	}
	410
	411	/* ================================================== */
	412
	413	int
	414	SPF_GetFilteredSample(SPF_Instance filter, NTP_Sample *sample)
	415	{
	416	int n;
	417
	418	n = select_samples(filter);
	419
a16094ad ML	420	DEBUG_LOG("selected %d from %d samples", n, filter->used);
a16094ad ML	421
c498c21f ML	422	if (n < 1)
	423	return 0;
	424
	425	if (!combine_selected_samples(filter, n, sample))
	426	return 0;
	427
e66f1df8	428	drop_samples(filter, 1);
c498c21f ML	429
	430	return 1;
	431	}
	432
	433	/* ================================================== */
	434
d5e645eb ML	435	static int
d5e645eb ML	436	get_first_last(SPF_Instance filter, int first, int last)
c498c21f	437	{
c498c21f	438	if (filter->last < 0)
d5e645eb	439	return 0;
c498c21f ML	440
	441	/* Always slew the last sample as it may be returned even if no new
	442	samples were accumulated */
	443	if (filter->used > 0) {
d5e645eb ML	444	*first = 0;
d5e645eb ML	445	*last = filter->used - 1;
c498c21f	446	} else {
d5e645eb	447	first = last = filter->last;
c498c21f ML	448	}
c498c21f ML	449
d5e645eb ML	450	return 1;
	451	}
	452
	453
	454	/* ================================================== */
	455
	456	void
	457	SPF_SlewSamples(SPF_Instance filter, struct timespec *when, double dfreq, double doffset)
	458	{
	459	int i, first, last;
	460	double delta_time;
	461
	462	if (!get_first_last(filter, &first, &last))
	463	return;
	464
c498c21f ML	465	for (i = first; i <= last; i++) {
	466	UTI_AdjustTimespec(&filter->samples[i].time, when, &filter->samples[i].time,
	467	&delta_time, dfreq, doffset);
	468	filter->samples[i].offset -= delta_time;
	469	}
	470	}
	471
	472	/* ================================================== */
	473
d5e645eb ML	474	void
	475	SPF_CorrectOffset(SPF_Instance filter, double doffset)
	476	{
	477	int i, first, last;
	478
	479	if (!get_first_last(filter, &first, &last))
	480	return;
	481
	482	for (i = first; i <= last; i++)
	483	filter->samples[i].offset -= doffset;
	484	}
	485
	486	/* ================================================== */
	487
c498c21f ML	488	void
	489	SPF_AddDispersion(SPF_Instance filter, double dispersion)
	490	{
	491	int i;
	492
	493	for (i = 0; i < filter->used; i++) {
	494	filter->samples[i].peer_dispersion += dispersion;
	495	filter->samples[i].root_dispersion += dispersion;
	496	}
	497	}