/* This is -*- C -*- */
/* $Id: guppi-seq-scalar.c,v 1.42 2002/01/14 05:01:17 trow Exp $ */

/*
 * guppi-seq-scalar.c
 *
 * Copyright (C) 2000 EMC Capital Management, Inc.
 * Copyright (C) 2001 The Free Software Foundation
 *
 * Developed by Jon Trowbridge <trow@gnu.org> and
 * Havoc Pennington <hp@pobox.com>.
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License as
 * published by the Free Software Foundation; either version 2 of the
 * License, or (at your option) any later version.
 *
 * 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 * USA
 */

#include <config.h>
#include "guppi-seq-scalar.h"

#include <math.h>
#include <stdlib.h>
#include <string.h>

#include <gtk/gtklabel.h>
#include <gtk/gtksignal.h>
#include <gtk/gtktable.h>

#include <libgnome/gnome-defs.h>
#include <libgnome/gnome-config.h>
#include <libgnome/gnome-i18n.h>

#include <guppi-convenient.h>
#include <guppi-string.h>

#include "guppi-seq-boolean.h"
/* #include <gnome.h> */

typedef struct _GuppiSeqScalarPrivate GuppiSeqScalarPrivate;
struct _GuppiSeqScalarPrivate {
  double min, max, sum, sum_abs, var, q1, median, q3;
  gint ordering;
  gpointer sorted_copy;

  guint have_ordering : 1;
  guint have_minmax : 1;
  guint have_sum : 1;
  guint have_sum_abs : 1;
  guint have_var : 1;
  guint have_quartiles : 1;
  guint save_ordering : 1;
  guint save_minmax : 1;
  guint save_sum : 1;
  guint save_sum_abs : 1;
};

static GtkObjectClass *parent_class = NULL;

static void
guppi_seq_scalar_finalize (GtkObject *obj)
{
  GuppiSeqScalar *seq = GUPPI_SEQ_SCALAR (obj);

  guppi_free0 (seq->priv->sorted_copy);
  guppi_free0 (seq->priv);

  if (parent_class->finalize)
    parent_class->finalize (obj);
}

static void
changed (GuppiData *data)
{
  GuppiSeqScalar *ss = GUPPI_SEQ_SCALAR (data);

  if (GUPPI_DATA_CLASS (parent_class)->changed)
    GUPPI_DATA_CLASS (parent_class)->changed (data);

  ss->priv->have_ordering = ss->priv->save_ordering;
  ss->priv->have_minmax = ss->priv->save_minmax;
  ss->priv->have_sum = ss->priv->save_sum;
  ss->priv->have_sum_abs = ss->priv->save_sum_abs;
  ss->priv->have_var = FALSE;
  ss->priv->have_quartiles = FALSE;

  ss->priv->save_ordering = FALSE;
  ss->priv->save_minmax = FALSE;
  ss->priv->save_sum = FALSE;

  guppi_free (ss->priv->sorted_copy);
  ss->priv->sorted_copy = NULL;
}

/* ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** */

/* Default vfunctions for handling missing values. */

static void
set (GuppiSeqScalar *ss, gint i, double x)
{
  GUPPI_SEQ_CLASS (GTK_OBJECT (ss)->klass)->set_missing (GUPPI_SEQ (ss), i, FALSE);
}

static void
set_many (GuppiSeqScalar *ss, gint i, const double *ptr, gint stride, gsize N)
{
  GUPPI_SEQ_CLASS (GTK_OBJECT (ss)->klass)->set_range_missing (GUPPI_SEQ (ss), i, i-1+N, FALSE);
}

static void
insert (GuppiSeqScalar *ss, gint i, double x)
{
  GUPPI_SEQ_CLASS (GTK_OBJECT (ss)->klass)->insert_missing (GUPPI_SEQ (ss), i, FALSE, 1);
}

static void
insert_many (GuppiSeqScalar *ss, gint i, const double *ptr, gint stride, gsize N)
{
  GUPPI_SEQ_CLASS (GTK_OBJECT (ss)->klass)->insert_missing (GUPPI_SEQ (ss), i, FALSE, N);
}

/* ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** */

/* Sequence Operation Stuff */

typedef struct _GuppiDataOp_Scalar GuppiDataOp_Scalar;
struct _GuppiDataOp_Scalar {
  GuppiDataOp op;

  gint i;
  gsize N;

  double x;

  const double *in_ptr;
  gint in_stride;

  gpointer out_ptr;
  gint out_stride;
};

static void
op_set (GuppiData *d, GuppiDataOp *op)
{
  GuppiSeqScalar *ss = GUPPI_SEQ_SCALAR (d);
  GuppiSeqScalarClass *klass = GUPPI_SEQ_SCALAR_CLASS (GTK_OBJECT (d)->klass);
  GuppiDataOp_Scalar *scalar_op = (GuppiDataOp_Scalar *) op;

  gboolean was_missing;
  double old_value;
  double x = scalar_op->x;

  was_missing = guppi_seq_missing (GUPPI_SEQ (d), scalar_op->i);
  old_value = was_missing ? 0 : guppi_seq_scalar_get (ss, scalar_op->i);

  g_assert (klass->set);
  klass->set (ss, scalar_op->i, x);

  /*
   * Try to incrementally update some of our stats, so that we can
   * possibly avoid future recalcs.
   */

  if (ss->priv->have_sum) {
    ss->priv->sum += x - old_value;
    ss->priv->save_sum = TRUE;
  }
  if (ss->priv->have_sum_abs) {
    ss->priv->sum_abs += fabs(x) - fabs(old_value);
    ss->priv->save_sum_abs = TRUE;
  }

  if (!was_missing && ss->priv->min == ss->priv->max) {

    /* We could probably do something for this case */

  } else if (!was_missing && old_value == ss->priv->min) {
    if (x < old_value) {
      ss->priv->min = x;
      ss->priv->save_minmax = TRUE;
    }
  } else if (!was_missing && old_value == ss->priv->max) {
    if (x > old_value) {
      ss->priv->max = x;
      ss->priv->save_minmax = TRUE;
    }
  } else if (x < ss->priv->min) {
    ss->priv->min = x;
    ss->priv->save_minmax = TRUE;
  } else if (x > ss->priv->max) {
    ss->priv->max = x;
    ss->priv->save_minmax = TRUE;
  } else if (ss->priv->min < x && x < ss->priv->max) {
    ss->priv->save_minmax = TRUE;
  }
}

static void
op_set_many (GuppiData *d, GuppiDataOp *op)
{
  GuppiSeqScalar *ss = GUPPI_SEQ_SCALAR (d);
  GuppiSeqScalarClass *klass = GUPPI_SEQ_SCALAR_CLASS (GTK_OBJECT (d)->klass);
  GuppiDataOp_Scalar *scalar_op = (GuppiDataOp_Scalar *) op;

  if (klass->set_many) {
    klass->set_many (ss,
		     scalar_op->i,
		     scalar_op->in_ptr,
		     scalar_op->in_stride, scalar_op->N);
  } else {
    /* Do our insert one piece at a time --- ugh */
    const guchar *byte_ptr = (const guchar *) scalar_op->in_ptr;
    gint count = scalar_op->N;
    gint j = scalar_op->i;

    g_assert (klass->set);
    for (; count-- > 0; ++j) {
      klass->set (ss, j, *(const double *) byte_ptr);
      byte_ptr += scalar_op->in_stride;
    }
  }
}
static void
op_insert (GuppiData *d, GuppiDataOp *op)
{
  GuppiSeqScalar *ss = GUPPI_SEQ_SCALAR (d);
  GuppiSeqScalarClass *klass = GUPPI_SEQ_SCALAR_CLASS (GTK_OBJECT (d)->klass);
  GuppiDataOp_Scalar *scalar_op = (GuppiDataOp_Scalar *) op;
  double x = scalar_op->x;

  g_assert (klass->insert);
  klass->insert (ss, scalar_op->i, x);

  if (ss->priv->have_sum) {
    ss->priv->sum += x;
    ss->priv->save_sum = TRUE;
  }
  if (ss->priv->have_sum_abs) {
    ss->priv->sum += fabs(x);
    ss->priv->save_sum_abs = TRUE;
  }

  if (ss->priv->have_minmax) {
    ss->priv->min = MIN (ss->priv->min, x);
    ss->priv->max = MAX (ss->priv->max, x);
    ss->priv->save_minmax = TRUE;
  }
}

static void
op_insert_many (GuppiData *d, GuppiDataOp *op)
{
  GuppiSeqScalar *ss = GUPPI_SEQ_SCALAR (d);
  GuppiSeqScalarClass *klass = GUPPI_SEQ_SCALAR_CLASS (GTK_OBJECT (d)->klass);
  GuppiDataOp_Scalar *scalar_op = (GuppiDataOp_Scalar *) op;

  if (klass->insert_many) {
    klass->insert_many (ss,
			scalar_op->i,
			scalar_op->in_ptr,
			scalar_op->in_stride, scalar_op->N);
  } else {
    /* Do our insert one piece at a time --- ugh */
    const guchar *byte_ptr = (const guchar *) scalar_op->in_ptr;
    gint j;

    g_assert (klass->insert);
    for (j = 0; j < scalar_op->N; ++j) {
      klass->insert (ss, scalar_op->i + j, *(const double *) byte_ptr);
      byte_ptr += scalar_op->in_stride;
    }
  }
}

/***************************************************************************/

/* Stuff for the sorted copy */

typedef struct _SortPair SortPair;
struct _SortPair {
  double x;
  gint i;
};

static int
sorted_pair_compare (const void *a, const void *b)
{
  double x = ((const SortPair *) a)->x;
  double y = ((const SortPair *) b)->x;

  return (x > y) - (x < y);
}

static void
make_sorted_copy (GuppiSeqScalar *seq)
{
  gint i, i0, i1, j, stride;
  gsize N;
  SortPair *sc;
  gboolean has_missing;
  const double *ptr;

  if (seq->priv->sorted_copy)
    return;

  has_missing = guppi_seq_has_missing (GUPPI_SEQ (seq));
  N = guppi_seq_count (GUPPI_SEQ (seq));
  ptr = guppi_seq_scalar_raw (seq, &stride);

  sc = guppi_new (SortPair, N);

  j = 0;
  guppi_seq_bounds (GUPPI_SEQ (seq), &i0, &i1);
  for (i = i0; i <= i1; ++i) {
    if (!has_missing || guppi_seq_available (GUPPI_SEQ (seq), i)) {
      double x;
      if (ptr)
	x = guppi_seq_scalar_raw_get (ptr, stride, i);
      else
	x = guppi_seq_scalar_get (seq, i);

      sc[j].x = x;
      sc[j].i = i;
      ++j;
    }
  }

  qsort (sc, N, sizeof (SortPair), sorted_pair_compare);

  seq->priv->sorted_copy = sc;
}

static SortPair *
get_sorted_copy (GuppiSeqScalar *seq)
{
  g_return_val_if_fail (seq != NULL && GUPPI_IS_SEQ_SCALAR (seq), NULL);

  if (seq->priv->sorted_copy == NULL)
    make_sorted_copy (seq);

  return (SortPair *) seq->priv->sorted_copy;
}


/***************************************************************************/

/**
 * guppi_seq_scalar_get
 * @seq: The scalar sequence
 * @i: The position index
 *
 * Returns the value of the @i-th element of @seq.
 * @i must lie within @seq's bounds.
 */
double
guppi_seq_scalar_get (const GuppiSeqScalar *seq, gint i)
{
  GuppiSeqScalarClass *klass;

  g_return_val_if_fail (GUPPI_IS_SEQ_SCALAR (seq), 0);
  g_return_val_if_fail (guppi_seq_in_bounds (GUPPI_SEQ (seq), i), 0);

  klass = GUPPI_SEQ_SCALAR_CLASS (GTK_OBJECT (seq)->klass);

  g_assert (klass->get);
  return klass->get ((GuppiSeqScalar *) seq, i);
}

/**
 * guppi_seq_scalar_set
 * @seq: The scalar sequence
 * @i: The position index
 * @val: A scalar value
 *
 * Sets the @i-th element of @seq equal to @val.
 * @seq must be writeable, and @i must lie within the
 * sequence's bounds.
 */
void
guppi_seq_scalar_set (GuppiSeqScalar *seq, gint i, double val)
{
  GuppiDataOp_Scalar op;

  g_return_if_fail (GUPPI_IS_SEQ_SCALAR (seq));
  g_return_if_fail (guppi_data_can_change (GUPPI_DATA (seq)));
  g_return_if_fail (guppi_seq_in_bounds (GUPPI_SEQ (seq), i));

  if (guppi_seq_missing (GUPPI_SEQ (seq), i) ||
      guppi_seq_scalar_get (seq, i) != val) {

    op.op.op = op_set;
    op.i = i;
    op.x = val;

    guppi_seq_changed_set (GUPPI_SEQ (seq), i, i, (GuppiDataOp *) & op);
  }
}

/**
 * guppi_seq_scalar_set_many
 * @seq: The scalar sequence
 * @start: The first position to be set in @seq
 * @vals: A pointer to a scalar value
 * @stride: The stride to use while iterating across the memory
 * @N: The number of scalar values to copy from @vals to @seq
 *
 * Copies scalar values from memory into @seq, overwriting the
 * existing values.  If all @N values don't fit within the sequence's
 * existing bounds, it will grow to accomodate the data.
 * @sec must be writeable.
 */
void
guppi_seq_scalar_set_many (GuppiSeqScalar *seq, gint start,
			   const double *vals, gint stride, gsize N)
{
  GuppiDataOp_Scalar op;

  g_return_if_fail (seq != NULL && GUPPI_IS_SEQ_SCALAR (seq));
  g_return_if_fail (guppi_data_can_change (GUPPI_DATA (seq)));

  if (N == 0)
    return;

  g_return_if_fail (vals != NULL);

  op.op.op = op_set_many;
  op.i = start;
  op.in_ptr = vals;
  op.in_stride = stride;
  op.N = N;

  guppi_seq_size_hint (GUPPI_SEQ (seq),
		       MAX (guppi_seq_size (GUPPI_SEQ (seq)), start + N));
  guppi_seq_changed_set (GUPPI_SEQ (seq), start, start + N,
			 (GuppiDataOp *) & op);
}

/**
 * guppi_seq_scalar_prepend
 * @seq: The scalar sequence
 * @val: A scalar value
 *
 * Prepend the value @val to @seq.
 * The sequence's lower bound will grow downward.
 * @seq must be writeable.
 */
void
guppi_seq_scalar_prepend (GuppiSeqScalar *seq, double val)
{
  gint first;

  g_return_if_fail (seq != NULL);
  g_return_if_fail (GUPPI_IS_SEQ (seq));

  first = guppi_seq_min_index (GUPPI_SEQ (seq));
  guppi_seq_scalar_insert (seq, first, val);
}

/**
 * guppi_seq_scalar_prepend_many
 * @seq: The scalar sequence
 * @vals: A pointer to the first value to prepend
 * @stride: The stride to use while interating across the memory
 * @N: The number of scalar values to read from memory and prepend to @seq.
 *
 * Prepend a sequences of @N scalar values starting at @ptr to @seq.
 * The sequence's lower bound will grow downward.
 * @seq must be writeable.
 */
void
guppi_seq_scalar_prepend_many (GuppiSeqScalar *seq,
			       const double *vals, gint stride, gsize N)
{
  gint first;

  g_return_if_fail (seq != NULL);
  g_return_if_fail (GUPPI_IS_SEQ (seq));

  first = guppi_seq_min_index (GUPPI_SEQ (seq));
  guppi_seq_scalar_insert_many (seq, first, vals, stride, N);
}

/**
 * guppi_seq_scalar_prepend_repeating
 * @seq: The scalar sequence
 * @val: A scalar value
 * @N: A count
 *
 * Prepend the sequences of @N occurances of @val to @seq.
 * The sequence's lower bound will grow downward.
 * @seq must be writeable.
 */
void
guppi_seq_scalar_prepend_repeating (GuppiSeqScalar *seq, double val, gsize N)
{
  guppi_seq_scalar_prepend_many (seq, &val, 0, N);
}

/**
 * guppi_seq_scalar_append
 * @seq: The scalar sequence
 * @val: A scalar value
 *
 * Append the value @val to @seq.
 * The sequence's upper bound will grow upward.
 * @seq must be writeable.
 */
void
guppi_seq_scalar_append (GuppiSeqScalar *seq, double val)
{
  gint last;
  last = guppi_seq_max_index (GUPPI_SEQ (seq));
  guppi_seq_scalar_insert (seq, last + 1, val);
}

/**
 * guppi_seq_scalar_append_many
 * @seq: The scalar sequence
 * @vals: A pointer to the first value to append
 * @stride: The stride to use while interating across the memory
 * @N: The number of scalar values to read from memory and prepend to @seq.
 *
 * Append the sequences of @N scalar values starting at @ptr to @seq.
 * The sequence's upper bound will grow upward.
 * @seq must be writeable.
 */
void
guppi_seq_scalar_append_many (GuppiSeqScalar *seq,
			      const double *vals, gint stride, gsize N)
{
  gint last;

  g_return_if_fail (seq != NULL);
  g_return_if_fail (GUPPI_IS_SEQ (seq));

  last = guppi_seq_max_index (GUPPI_SEQ (seq));
  guppi_seq_scalar_insert_many (seq, last + 1, vals, stride, N);
}

/**
 * guppi_seq_scalar_append_repeating
 * @seq: The scalar sequence
 * @val: A scalar value
 * @N: A count
 *
 * Append a sequences of @N occurances of @val to @seq.
 * The sequence's upper bound will grow upward.
 * @seq must be writeable.
 */
void
guppi_seq_scalar_append_repeating (GuppiSeqScalar *seq, double val, gsize N)
{
  guppi_seq_scalar_append_many (seq, &val, 0, N);
}

/**
 * guppi_seq_scalar_insert
 * @seq: The scalar sequence
 * @i: The position index to insert at
 * @val: The value to insert
 *
 * Inserts the element @val into @seq at position @i.
 */
void
guppi_seq_scalar_insert (GuppiSeqScalar *seq, gint i, double val)
{
  GuppiDataOp_Scalar op;

  g_return_if_fail (seq != NULL && GUPPI_IS_SEQ_SCALAR (seq));
  g_return_if_fail (guppi_data_can_change (GUPPI_DATA (seq)));

  op.op.op = op_insert;
  op.i = i;
  op.x = val;

  guppi_seq_changed_insert (GUPPI_SEQ (seq), i, 1, (GuppiDataOp *) & op);
}

/**
 * guppi_seq_scalar_insert_many
 * @seq: The scalar sequence
 * @i: The position the first element is to be inserted at
 * @vals: A pointer to the first value to insert
 * @stride: The stride to use while interating across the memory
 * @N: The number of scalar values to read from memory and inserted into @seq.
 *
 * Insert the sequence of @N scalar values starting at @ptr into @seq.
 * The first element of the inserted sequence becomes the @i-th element of @seq.
 * @seq must be writeable.
 */
void
guppi_seq_scalar_insert_many (GuppiSeqScalar *seq, gint i,
			      const double *vals, gint stride, gsize N)
{
  GuppiDataOp_Scalar op;

  g_return_if_fail (seq != NULL && GUPPI_IS_SEQ_SCALAR (seq));
  g_return_if_fail (guppi_data_can_change (GUPPI_DATA (seq)));

  if (N == 0)
    return;

  g_return_if_fail (vals != NULL);

  op.op.op = op_insert_many;
  op.i = i;
  op.in_ptr = vals;
  op.in_stride = stride;
  op.N = N;

  guppi_seq_size_hint (GUPPI_SEQ (seq), guppi_seq_size (GUPPI_SEQ (seq)) + N);
  guppi_seq_changed_insert (GUPPI_SEQ (seq), i, N, (GuppiDataOp *) & op);
}

/**
 * guppi_seq_scalar_insert_repeating
 * @seq: The scalar sequence
 * @i: The position to begin inserting @val at
 * @val: A scalar value
 * @N: A count
 *
 * Insert @N occurances of @val into @seq,
 * starting at position @i.
 * @seq must be writeable.
 */
void
guppi_seq_scalar_insert_repeating (GuppiSeqScalar *seq, gint i,
				   double x, gsize N)
{
  guppi_seq_scalar_insert_many (seq, i, &x, 0, N);
}

static gint
do_range_query (GuppiSeqScalar *seq,
		GuppiSeqBoolean *bseq,
		double min, double max, gboolean do_and)
{
  GuppiSeqScalarClass *klass;

  gint hits;
  SortPair *sorted_copy;

  double seq_min, seq_max;
  gint count = 0;

  g_return_val_if_fail (GUPPI_IS_SEQ_SCALAR (seq), 0);
  g_return_val_if_fail (GUPPI_IS_SEQ_BOOLEAN (bseq), 0);

  klass = GUPPI_SEQ_SCALAR_CLASS (GTK_OBJECT (seq)->klass);

  guppi_2sort (&min, &max);
  seq_min = guppi_seq_scalar_min (seq);
  seq_max = guppi_seq_scalar_max (seq);

  if (min <= seq_min && seq_max <= max) {
    if (do_and) {
      /* do nothing */
    } else {
      guppi_seq_boolean_set_all (bseq, TRUE);
    }
    return guppi_seq_size (GUPPI_SEQ (seq));
  }

  if (klass->range_query) {
    hits = (klass->range_query) (seq, bseq, min, max, do_and);
    if (hits >= 0)
      return hits;
  }

  sorted_copy = get_sorted_copy (seq);

  if (sorted_copy) {
    gint k0, k1, k, a, b;
    gint set_buffer[500];
    gint sb_N = 500, sb_i;
    gint seq_count, N;
    GuppiSeqBoolean *op_bool;

    N = guppi_seq_size (GUPPI_SEQ (seq));
    seq_count = guppi_seq_count (GUPPI_SEQ (seq));

    if (do_and) {

      op_bool = GUPPI_SEQ_BOOLEAN (guppi_data_new ("GuppiSeqBooleanCore"));
      guppi_seq_size_hint (GUPPI_SEQ (op_bool), N);
      guppi_seq_boolean_append_many (op_bool, FALSE, N);
      guppi_seq_set_min_index (GUPPI_SEQ (op_bool),
			       guppi_seq_min_index (GUPPI_SEQ (seq)));

    } else {
      guppi_seq_boolean_clear (bseq);
      op_bool = bseq;
    }

    if (min <= sorted_copy[0].x)
      k0 = 0;
    else if (min > sorted_copy[seq_count - 1].x)
      k0 = seq_count;
    else {

      a = 0;
      b = seq_count - 1;
      k0 = (a + b) / 2;

      while (!(sorted_copy[k0 - 1].x < min && min <= sorted_copy[k0].x)) {
	gint k0_old = k0;
	if (min <= sorted_copy[k0 - 1].x)
	  b = k0 + 1;
	else			/* if (min > sorted_copy[m].x) */
	  a = k0 - 1;
	k0 = (a + b) / 2;
	if (k0 == k0_old)
	  ++k0;
      }

    }

    if (sorted_copy[seq_count - 1].x <= max)
      k1 = seq_count - 1;
    else if (max < sorted_copy[0].x)
      k1 = -1;
    else {
      a = k0;
      b = seq_count - 1;
      k1 = (a + b) / 2;

      while (!(sorted_copy[k1].x <= max && max < sorted_copy[k1 + 1].x)) {
	gint k1_old = k1;
	if (sorted_copy[k1].x > max)
	  b = k1 - 1;
	else			/* if (max >= sorted_copy[k1+1].x) */
	  a = k1 + 1;
	k1 = (a + b) / 2;
	if (k1 == k1_old)
	  --k1;
      }
    }

    if (k0 == 0 && k1 == seq_count - 1) {

      guppi_seq_boolean_set_all (op_bool, TRUE);

    } else {

      sb_i = 0;
      for (k = k0; k <= k1; ++k) {
	set_buffer[sb_i] = sorted_copy[k].i;
	++sb_i;
	if (sb_i == sb_N) {
	  guppi_seq_boolean_set_many (op_bool, set_buffer, sb_N, TRUE);
	  sb_i = 0;
	}
      }

      /* flush the set buffer */
      if (sb_i > 0)
	guppi_seq_boolean_set_many (op_bool, set_buffer, sb_i, TRUE);

    }

    if (do_and) {
      guppi_seq_boolean_bitwise_and (bseq, op_bool);
      guppi_unref0 (op_bool);
    }

    return guppi_seq_boolean_true_count (bseq);

  } else {
    /* no sorted copy available */
    g_assert_not_reached ();
  }

  return count;
}

/**
 * guppi_seq_scalar_range_query
 * @seq: A scalar sequence
 * @min: The lower bound of the query
 * @max: The upper bound of the query
 *
 * The range query interates across @seq and creates a
 * #GuppiSeqBoolean whose values act like a bit-mask.  The return
 * value's elements are set to %TRUE or %FALSE depending on if the
 * corresponding values of @seq lie between @min and @max. 
 *
 * This function is used to optimize some of Guppi's graphical
 * computations.
 *
 * Returns: a boolean sequence indicating which elements of @seq lie
 * between @min and @max.
 */
GuppiSeqBoolean *
guppi_seq_scalar_range_query (GuppiSeqScalar *seq, double min, double max)
{
  GuppiSeqBoolean *bseq;

  guppi_2sort (&min, &max);

  g_return_val_if_fail (seq != NULL, NULL);

  bseq = GUPPI_SEQ_BOOLEAN (guppi_data_new ("GuppiSeqBooleanCore"));
  guppi_seq_boolean_append_many (bseq, FALSE, guppi_seq_size (GUPPI_SEQ (seq)));
  guppi_seq_set_min_index (GUPPI_SEQ (bseq), guppi_seq_min_index (GUPPI_SEQ (seq)));

  do_range_query (seq, bseq, min, max, FALSE);

  return bseq;
}

/**
 * guppi_seq_scalar_in_place_range_query
 * @seq: A scalar sequence
 * @seq_bool: An existing boolean sequence
 * @min: The lower bound of the query
 * @max: The upper bound of the query
 *
 * Performs a range query similar to guppi_seq_scalar_range_query(),
 * but the results are stored in a pre-existing #GuppiSeqBoolean.
 * In some cases this can be more efficient that repeatedly allocating
 * and freeing boolean sequences to contain the results of multiple
 * successive queries.
 *
 * This function is used to optimize some of Guppi's graphical
 * computations.
 *
 * Returns: the number of elements of @seq that lie between @min and @max.
 */
gint
guppi_seq_scalar_in_place_range_query (GuppiSeqScalar *seq,
				       GuppiSeqBoolean *seq_bool,
				       double min, double max)
{
  return do_range_query (seq, seq_bool, min, max, FALSE);
}

/**
 * guppi_seq_scalar_bitwise_and_range_query
 * @seq: A scalar sequences
 * @seq_bool: An existing boolean sequences
 * @min: The lower bound of the query
 * @max: The upper bound of the query
 *
 * This function performs a range query that stores the results in
 * an existing boolean sequence, much like guppi_seq_scalar_in_place_range_query().
 * However, instead of storing the results directly in @seq_bool, the queries
 * results bitmask is applied to @seq_bool as a "bitwise-AND".
 *
 * In other words, an element of @seq_bool is %TRUE after calling this
 * function if and only if the corresponding value of @seq lies between
 * @min and @max AND if that element of @seq_bool was %TRUE before the
 * function call.
 *
 * This function is used to optimize some of Guppi's graphical
 * computations.
 *
 * Returns: the number of %TRUE values set in @seq_bool upon return.
 */
gint
guppi_seq_scalar_bitwise_and_range_query (GuppiSeqScalar *seq,
					  GuppiSeqBoolean *seq_bool,
					  double min, double max)
{
  return do_range_query (seq, seq_bool, min, max, TRUE);
}

static void
find_range (SortPair *sp, gint N, double min, double max, gint *i0, gint *i1)
{
  gint k0, k1, a, b;

  if (min <= sp[0].x)
    k0 = 0;
  else if (min > sp[N-1].x)
    k0 = N;
  else {
    a = 0;
    b = N-1;
    k0 = (a + b) / 2;

    while (!(sp[k0 - 1].x < min && min <= sp[k0].x)) {
      gint k0_old = k0;
      if (min <= sp[k0 - 1].x)
	b = k0 + 1;
      else			/* if (min > sorted_copy[m].x) */
	a = k0 - 1;
      k0 = (a + b) / 2;
      if (k0 == k0_old)
	++k0;
    }
  }
  
  if (sp[N - 1].x <= max)
    k1 = N - 1;
  else if (max < sp[0].x)
    k1 = -1;
  else {
    a = k0;
    b = N - 1;
    k1 = (a + b) / 2;

    while (!(sp[k1].x <= max && max < sp[k1 + 1].x)) {
      gint k1_old = k1;
      if (sp[k1].x > max)
	b = k1 - 1;
      else			/* if (max >= sorted_copy[k1+1].x) */
	a = k1 + 1;
      k1 = (a + b) / 2;
      if (k1 == k1_old)
	--k1;
    }
  }

  if (i0) *i0 = k0;
  if (i1) *i1 = k1;
}

gint
guppi_seq_scalar_gather_pairs (GuppiSeqScalar *a, GuppiSeqScalar *b, 
			       double x0, double y0,
			       double x1, double y1, 
			       GuppiIndexedPairFn fn, gpointer closure)
{
  gint count;
  GuppiSeqScalarClass *klass;

  const gint N = 1000;
  GuppiIndexedPair buffer[1000];
  gint buf_i = 0;

  SortPair *sa, *sb;
  gint ak0, ak1, aN;
  gint bk0, bk1, bN;
  gint k, k0, k1;
  double z0, z1;

  SortPair *sorted;
  const double *raw;
  gint stride;
  gboolean transpose;

  g_return_val_if_fail (GUPPI_IS_SEQ_SCALAR (a), -1);
  g_return_val_if_fail (GUPPI_IS_SEQ_SCALAR (b), -1);
  
  guppi_2sort (&x0, &x1);
  guppi_2sort (&y0, &y1);

  klass = GUPPI_SEQ_SCALAR_CLASS (GTK_OBJECT (a)->klass);
  if (GTK_OBJECT_TYPE (a) == GTK_OBJECT_TYPE (b) && klass->gather_pairs) {
    count = klass->gather_pairs (a, b, x0, y0, x1, y1, fn, closure);
    if (count >= 0)
      return count;
  }

  sa = get_sorted_copy (a);
  aN = guppi_seq_count (GUPPI_SEQ (a));
  find_range (sa, aN, x0, x1, &ak0, &ak1);

  /* If either one of our ranges is empty,
     there is nothing more to do. */
  if (ak1 < ak0)
    return 0;

  sb = get_sorted_copy (b);
  bN = guppi_seq_count (GUPPI_SEQ (b));
  find_range (sb, bN, y0, y1, &bk0, &bk1);

  if (bk1 < bk0)
    return 0;

  if (ak1 - ak0 < bk1 - bk0) {
    sorted = sa;
    raw = guppi_seq_scalar_raw (b, &stride);
    transpose = FALSE;
    k0 = ak0, k1 = ak1;
    z0 = y0, z1 = y1;
  } else {
    sorted = sb;
    raw = guppi_seq_scalar_raw (a, &stride);
    transpose = TRUE;
    k0 = bk0, k1 = bk1;
    z0 = x0, z1 = x1;
  }

  count = 0;
  for (k = k0; k <= k1; ++k) {
    gint i = sorted[k].i;
    double x = sorted[k].x;
    double y = guppi_seq_scalar_raw_get (raw, stride, i);
    if (z0 <= y && y <= z1) {
      buffer[buf_i].i = i;
      buffer[buf_i].x = transpose ? y : x;
      buffer[buf_i].y = transpose ? x : y;
      ++buf_i;
      ++count;
    }

    if (buf_i == N || (k == k1 && buf_i > 0)) {
      if (fn)
	fn (buffer, buf_i, closure);
      buf_i = 0;
    }
  }
      
  return count;
}


/**
 * guppi_seq_scalar_ordering
 * @seq: A scalar sequence
 *
 * Analyzes the sequence @seq to determine if its values
 * are arranged in an increasing or decreasing pattern.
 *
 * Returns 1 if @seq is constant or increasing (non-decreasing),
 * -1 if @seq is decreasing (non-increasing), and 0 otherwise.
 */

/*
    2 = constant
    1 = increasing
   -1 = decreasing
*/

#define ord_val(x,y) ((x)==(y) ? 2 : ((x) < (y) ? 1 : -1))
gint
guppi_seq_scalar_ordering (GuppiSeqScalar *seq)
{
  gint i, i0, i1;

  g_return_val_if_fail (GUPPI_IS_SEQ_SCALAR (seq), 0);

  if (! seq->priv->have_ordering) {
    gint stride;
    const double *raw; 
    double x0, x1;
    gint ord, val;
    
    if (guppi_seq_empty (GUPPI_SEQ (seq)))
      return 0;

    /* Missing values mean unknown ordering. */
    if (guppi_seq_size (GUPPI_SEQ (seq)) != guppi_seq_count (GUPPI_SEQ (seq)))
      return 0;

    /* Singleton sequences are constant, of course. */
    if (guppi_seq_size (GUPPI_SEQ (seq)) == 1) {
      return seq->priv->ordering = 1;
    }

    guppi_seq_bounds (GUPPI_SEQ (seq), &i0, &i1);
    raw = guppi_seq_scalar_raw (seq, &stride);

    x0 = guppi_seq_scalar_get (seq, i0);
    x1 = guppi_seq_scalar_get (seq, i0+1);

    ord = ord_val (x0, x1);

    for (i = i0+2; i <= i1 && ord; ++i) {
      x0 = x1;
      x1 = raw ? guppi_seq_scalar_raw_get (seq, stride, i) : guppi_seq_scalar_get (seq, i);
      val = ord_val (x0, x1);
      if (ord == 2)
	ord = val;
      else if (ord != val)
	ord = 0;
    }
      
    seq->priv->ordering = ord;
    seq->priv->have_ordering = TRUE;
  }
  
  /* We keep the fact that the sequence is constant internal,
     and identify the ordering as increasing. */
  return seq->priv->ordering == 2 ? 1 : seq->priv->ordering;
}

static void
calc_minmax (GuppiSeqScalar *seq)
{
  GuppiSeqScalarClass *klass;

  klass = GUPPI_SEQ_SCALAR_CLASS (GTK_OBJECT (seq)->klass);

  if (klass->range) {

    klass->range (seq, &seq->priv->min, &seq->priv->max);

  } else {
    gint i0 = guppi_seq_min_index (GUPPI_SEQ (seq));
    gint i1 = guppi_seq_max_index (GUPPI_SEQ (seq));
    gint i;
    double x0 = 0, x1 = 0, y;
    gboolean has_missing = guppi_seq_has_missing (GUPPI_SEQ (seq));
    const double *ptr;
    gint stride;

    ptr = guppi_seq_scalar_raw (seq, &stride);

    while (has_missing && guppi_seq_missing (GUPPI_SEQ (seq), i0) && i0 <= i1)
      ++i0;

    if (i0 <= i1) {
      if (ptr)
	x0 = x1 = guppi_seq_scalar_raw_get (ptr, stride, i0);
      else
	x0 = x1 = guppi_seq_scalar_get (seq, i0);
    }

    for (i = i0 + 1; i <= i1; ++i) {
      if (!(has_missing && guppi_seq_missing (GUPPI_SEQ (seq), i))) {
	if (ptr)
	  y = guppi_seq_scalar_raw_get (ptr, stride, i);
	else
	  y = guppi_seq_scalar_get (seq, i);
	if (y < x0)
	  x0 = y;
	if (y > x1)
	  x1 = y;
      }
    }
    seq->priv->min = x0;
    seq->priv->max = x1;
  }

  seq->priv->have_minmax = TRUE;
}

/**
 * guppi_seq_scalar_min
 * @seq: A scalar sequence
 *
 * Finds the minimum value of the sequence.
 * Missing values are ignored.
 * The value is cached, so it is not inefficient to call 
 * guppi_seq_scalar_min() repeatedly.
 * @seq must be non-empty.
 *
 * Returns the smallest scalar element of @seq.
 */
double
guppi_seq_scalar_min (GuppiSeqScalar *seq)
{
  g_return_val_if_fail (GUPPI_IS_SEQ_SCALAR (seq), 0);
  g_return_val_if_fail (guppi_seq_nonempty (GUPPI_SEQ (seq)), 0);

  if (!seq->priv->have_minmax) {
    calc_minmax (seq);
  } 

  return seq->priv->min;
}

/**
 * guppi_seq_scalar_max
 * @seq: A scalar sequence
 *
 * Finds the maximum value of the sequence.
 * Missing values are ignored.
 * The value is cached, so it is not inefficient to call 
 * guppi_seq_scalar_max() repeatedly.
 * @seq must be non-empty.
 *
 * Returns the largest scalar element of @seq.
 */
double
guppi_seq_scalar_max (GuppiSeqScalar *seq)
{
  g_return_val_if_fail (GUPPI_IS_SEQ_SCALAR (seq), 0);
  g_return_val_if_fail (guppi_seq_nonempty (GUPPI_SEQ (seq)), 0);

  if (!seq->priv->have_minmax) {
    calc_minmax (seq);
  }

  return seq->priv->max;
}

/**
 * guppi_seq_scalar_sum
 * @seq: A scalar sequence
 *
 * Finds the sum of all of the values in the sequence.
 * Missing values are ignored.
 * The value is cached, so it is not inefficient to call 
 * guppi_seq_scalar_sum() repeatedly.
 * @seq must be non-empty.
 *
 * Returns the sum of the elements of @seq.
 */
double
guppi_seq_scalar_sum (GuppiSeqScalar *seq)
{
  g_return_val_if_fail (GUPPI_IS_SEQ_SCALAR (seq), 0);

  if (guppi_seq_empty (GUPPI_SEQ (seq)))
    return 0;

  if (!seq->priv->have_sum) {

    GuppiSeqScalarClass *klass;

    klass = GUPPI_SEQ_SCALAR_CLASS (GTK_OBJECT (seq)->klass);

    if (klass->sum) {

      seq->priv->sum = klass->sum (seq);

    } else {

      gint i0 = guppi_seq_min_index (GUPPI_SEQ (seq));
      gint i1 = guppi_seq_max_index (GUPPI_SEQ (seq));
      gint i;
      double y = 0;
      gboolean has_missing = guppi_seq_has_missing (GUPPI_SEQ (seq));

      for (i = i0; i <= i1; ++i)
	if (!(has_missing && guppi_seq_missing (GUPPI_SEQ (seq), i)))
	  y += guppi_seq_scalar_get (seq, i);

      seq->priv->sum = y;
    }

    seq->priv->have_sum = TRUE;
  }

  return seq->priv->sum;
}

/**
 * guppi_seq_scalar_sum_abs
 * @seq: A scalar sequence
 *
 * Finds the sum of the absolute values of all of the elements of the sequence.
 * Missing values are ignored.
 * The value is cached, so it is not inefficient to call 
 * guppi_seq_scalar_sum_abs() repeatedly.
 * @seq must be non-empty.
 *
 * Returns the sum of the absolute values of the elements of @seq.
 */
double
guppi_seq_scalar_sum_abs (GuppiSeqScalar *seq)
{
  g_return_val_if_fail (GUPPI_IS_SEQ_SCALAR (seq), 0);

  if (guppi_seq_empty (GUPPI_SEQ (seq)))
    return 0;

  if (!seq->priv->have_sum_abs) {
    gint i0 = guppi_seq_min_index (GUPPI_SEQ (seq));
    gint i1 = guppi_seq_max_index (GUPPI_SEQ (seq));
    gint i;
    double y = 0;
    gboolean has_missing = guppi_seq_has_missing (GUPPI_SEQ (seq));

    for (i = i0; i <= i1; ++i)
      if (!(has_missing && guppi_seq_missing (GUPPI_SEQ (seq), i)))
	y += fabs (guppi_seq_scalar_get (seq, i));

    seq->priv->sum_abs = y;
    seq->priv->have_sum_abs = TRUE;
  }

  return seq->priv->sum_abs;
}

/**
 * guppi_seq_scalar_mean
 * @seq: A scalar sequence
 *
 * Finds the mean (or average) of all of the values in the sequence.
 * Missing values are ignored.
 * The value is cached, so it is not inefficient to call 
 * guppi_seq_scalar_mean() repeatedly.
 * @seq must be non-empty.
 *
 * Returns the mean of the elements of @seq.
 */
double
guppi_seq_scalar_mean (GuppiSeqScalar *seq)
{
  double sum;
  gsize n;

  g_return_val_if_fail (GUPPI_IS_SEQ_SCALAR (seq), 0);

  n = guppi_seq_count (GUPPI_SEQ (seq));
  g_return_val_if_fail (n > 0, 0);
  sum = guppi_seq_scalar_sum (seq);

  return sum / n;
}


/**
 * guppi_seq_scalar_var
 * @seq: A scalar sequence
 *
 * Finds the population variance of the sequence.
 * Missing values are ignored.
 * The value is cached, so it is not inefficient to call 
 * guppi_seq_scalar_var() repeatedly.
 * @seq must be non-empty.
 *
 * Returns the population variance of the elements of @seq.
 */
double
guppi_seq_scalar_var (GuppiSeqScalar *seq)
{
  g_return_val_if_fail (GUPPI_IS_SEQ_SCALAR (seq), 0);

  if (!seq->priv->have_var) {

    GuppiSeqScalarClass *klass;

    klass = GUPPI_SEQ_SCALAR_CLASS (GTK_OBJECT (seq)->klass);

    if (klass->var) {

      seq->priv->var = klass->var (seq);

    } else {

      double x, om, mean = 0, sumsq = 0;
      gint i, i0, i1, count = 0, stride;
      const double *ptr;
      gboolean has_missing;

      guppi_seq_indices (GUPPI_SEQ (seq), &i0, &i1);
      has_missing = guppi_seq_has_missing (GUPPI_SEQ (seq));

      ptr = guppi_seq_scalar_raw (seq, &stride);

      for (i = i0; i <= i1; ++i) {

	if (!(has_missing && guppi_seq_missing (GUPPI_SEQ (seq), i))) {

	  if (ptr)
	    x = guppi_seq_scalar_raw_get (ptr, stride, i);
	  else
	    x = guppi_seq_scalar_get (seq, i);

	  ++count;
	  om = mean;
	  mean += (x - mean) / count;
	  if (count > 1)
	    sumsq += (x - mean) * (x - om);
	}
      }

      seq->priv->var = sumsq / count;
    }

    seq->priv->have_var = TRUE;
  }

  return seq->priv->var;
}

/**
 * guppi_seq_scalar_vars
 * @seq: A scalar sequence
 *
 * Finds the sample variance of the sequence.
 * Missing values are ignored.
 * The value is cached, so it is not inefficient to call 
 * guppi_seq_scalar_vars() repeatedly.
 * @seq must contain at least two elements.
 *
 * Returns the population variance of the elements of @seq.
 */
double
guppi_seq_scalar_vars (GuppiSeqScalar *seq)
{
  double v;
  gsize n;

  g_return_val_if_fail (seq != NULL && GUPPI_IS_SEQ_SCALAR (seq), 0);

  n = guppi_seq_count (GUPPI_SEQ (seq));
  g_return_val_if_fail (n > 1, 0);

  v = guppi_seq_scalar_var (seq);
  return (n * v) / (n - 1);
}

/**
 * guppi_seq_scalar_sdev
 * @seq: A scalar sequence
 *
 * Finds the population standard deviation of the sequence.
 * Missing values are ignored.
 * The value is cached, so it is not inefficient to call 
 * guppi_seq_scalar_sdev() repeatedly.
 * @seq must be non-empty.
 *
 * Returns the population standard deviation of the elements of @seq.
 */
double
guppi_seq_scalar_sdev (GuppiSeqScalar *seq)
{
  g_return_val_if_fail (seq != NULL && GUPPI_IS_SEQ_SCALAR (seq), 0);
  return sqrt (guppi_seq_scalar_var (seq));
}

/**
 * guppi_seq_scalar_sdevs
 * @seq: A scalar sequence
 *
 * Finds the sample standard deviation of the sequence.
 * Missing values are ignored.
 * The value is cached, so it is not inefficient to call 
 * guppi_seq_scalar_sdevs() repeatedly.
 * @seq must contain at least two elements.
 *
 * Returns the population standard deviation of the elements of @seq.
 */
double
guppi_seq_scalar_sdevs (GuppiSeqScalar *seq)
{
  g_return_val_if_fail (seq != NULL && GUPPI_IS_SEQ_SCALAR (seq), 0);
  return sqrt (guppi_seq_scalar_vars (seq));
}

static void
calc_quartiles (GuppiSeqScalar *seq)
{
  GuppiSeqScalarClass *klass;

  if (seq->priv->have_quartiles)
    return;

  klass = GUPPI_SEQ_SCALAR_CLASS (GTK_OBJECT (seq)->klass);

  /* Try to use an implementation-specific calculation. */
  seq->priv->have_quartiles = klass->quartiles
    && klass->quartiles (seq, &seq->priv->q1, &seq->priv->median, &seq->priv->q3);

  /* If that doesn't work, use our fallback calculation. */
  if (!seq->priv->have_quartiles) {
    gsize N = guppi_seq_count (GUPPI_SEQ (seq));

    if (N > 0) {
      SortPair *sc = get_sorted_copy (seq);
      gint i;
      double t;

      g_assert (sc != NULL);

      t = 0.25 * (N - 1);
      i = (gint) t;
      seq->priv->q1 = (i + 1 - t) * sc[i].x + (t - i) * sc[i + 1].x;

      t = 0.5 * (N - 1);
      i = (gint) t;
      seq->priv->median = (i + 1 - t) * sc[i].x + (t - i) * sc[i + 1].x;

      t = 0.75 * (N - 1);
      i = (gint) t;
      seq->priv->q3 = (i + 1 - t) * sc[i].x + (t - i) * sc[i + 1].x;

    } else {

      seq->priv->q1 = seq->priv->median = seq->priv->q3 = 0;

    }
  }

  seq->priv->have_quartiles = TRUE;
}

/**
 * guppi_seq_scalar_q1
 * @seq: A scalar sequence
 *
 * Computes the first quartile (Q1) of the elements of @seq.
 * Missing values are ignored.
 * The value is cached, so it is not inefficient to call 
 * guppi_seq_scalar_q1() repeatedly.
 * @seq must be non-empty.
 *
 * Returns: Q1 of @seq.
 */
double
guppi_seq_scalar_q1 (GuppiSeqScalar *seq)
{
  g_return_val_if_fail (GUPPI_IS_SEQ_SCALAR (seq), 0);
  g_return_val_if_fail (guppi_seq_nonempty (GUPPI_SEQ (seq)), 0);

  if (!seq->priv->have_quartiles)
    calc_quartiles (seq);

  return seq->priv->q1;
}

/**
 * guppi_seq_scalar_median
 * @seq: A scalar sequence
 *
 * Computes the median of the elements of @seq.
 * Missing values are ignored.
 * The value is cached, so it is not inefficient to call 
 * guppi_seq_scalar_median() repeatedly.
 * @seq must be non-empty.
 *
 * Returns: the median of @seq.
 */
double
guppi_seq_scalar_median (GuppiSeqScalar *seq)
{
  g_return_val_if_fail (GUPPI_IS_SEQ_SCALAR (seq), 0);
  g_return_val_if_fail (guppi_seq_nonempty (GUPPI_SEQ (seq)), 0);

  if (!seq->priv->have_quartiles)
    calc_quartiles (seq);

  return seq->priv->median;
}

/**
 * guppi_seq_scalar_q3
 * @seq: A scalar sequence
 *
 * Computes the third quartile (Q3) of the elements of @seq.
 * Missing values are ignored.
 * The value is cached, so it is not inefficient to call 
 * guppi_seq_scalar_q3() repeatedly.
 * @seq must be non-empty.
 *
 * Returns: Q3 of @seq.
 */
double
guppi_seq_scalar_q3 (GuppiSeqScalar *seq)
{
  g_return_val_if_fail (GUPPI_IS_SEQ_SCALAR (seq), 0);
  g_return_val_if_fail (guppi_seq_nonempty (GUPPI_SEQ (seq)), 0);

  if (!seq->priv->have_quartiles)
    calc_quartiles (seq);

  return seq->priv->q3;
}

/**
 * guppi_seq_scalar_percentile
 * @seq: A scalar sequence
 * @p: A percentile value between 0 and 1.
 *
 * Computes the value at the @p-th percentile of the elements of @seq,
 * interpolating between elements as necessary.
 * Missing values are ignored.
 * @seq must be non-empty.
 *
 * Returns: the @p-th percentile of @seq.
 */
double
guppi_seq_scalar_percentile (GuppiSeqScalar *seq, double p)
{
  GuppiSeqScalarClass *klass;
  double x;
  gint N;

  g_return_val_if_fail (GUPPI_IS_SEQ_SCALAR (seq), 0);
  g_return_val_if_fail (guppi_seq_nonempty (GUPPI_SEQ (seq)), 0);
  g_return_val_if_fail (0 <= p && p <= 1, 0);

  klass = GUPPI_SEQ_SCALAR_CLASS (GTK_OBJECT (seq)->klass);

  if (klass->percentile && klass->percentile (seq, p, &x))
    return x;

  N = guppi_seq_count (GUPPI_SEQ (seq));
  if (N > 0) {
    SortPair *sc = get_sorted_copy (seq);
    gint i;
    double t;

    t = p * (N - 1);
    i = (gint) t;
    return (i + 1 - t) * sc[i].x + (t - i) * sc[i + 1].x;

  }

  g_assert_not_reached ();
  return 0;
}

/**
 * guppi_seq_scalar_raw
 * @seq: A scalar sequence
 * @stride: A variable to store the required memory stride
 *
 * For certain #GuppiSeqScalar implementations, the data will be
 * stored in an array in the machine's memory.  If this is the case
 * for @seq, guppi_seq_scalar_raw() provides a pointer and a stride
 * that can be used to read the data directly, which is much more
 * efficiently that using repeatedly calling guppi_seq_scalar_get().
 *
 * The returned pointer and the value in @stride can be combined and
 * used to find the value of a particular element of @seq by
 * using guppi_seq_scalar_raw_get().
 *
 * If this feature is not supported by @seq's particular implementation,
 * %NULL is returned.
 *
 * Returns: a pointer into memory that can be used to retrieve
 * sequence values directly.
 */
const double *
guppi_seq_scalar_raw (const GuppiSeqScalar *seq, gint *stride)
{
  GuppiSeqScalarClass *klass;

  g_return_val_if_fail (GUPPI_IS_SEQ_SCALAR (seq), NULL);
  g_return_val_if_fail (stride != NULL, NULL);

  klass = GUPPI_SEQ_SCALAR_CLASS (GTK_OBJECT (seq)->klass);

  return klass->raw_access ? klass->raw_access ((GuppiSeqScalar *) seq, stride) : NULL;
}

/* This is kind of funny, embedding the docs for the macro here... */

/**
 * guppi_seq_scalar_raw_get
 * @ptr: The pointer returns by guppi_seq_scalar_raw()
 * @str: The stride returned (in an argument) by guppi_seq_scalar_raw()
 * @i: The position index
 *
 * This macro can be used to extract the elements of a sequence from
 * memory very efficiently, using the data return from a call to
 * guppi_seq_scalar_raw().
 */

/* ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** */

enum {
  MIN_LABEL,
  MAX_LABEL,
  MEAN_LABEL,
  SDEV_LABEL,
  MEDIAN_LABEL,
  Q1_LABEL,
  Q3_LABEL,
  IQR_LABEL,
  LAST_LABEL
};

static const gchar *label_names [LAST_LABEL] = {
  N_("Min"), N_("Max"), N_("Mean"), N_("Std Dev"),
  N_("Median"), N_("Q1"), N_("Q3"), N_("IRQ")
};

typedef struct _info info;
struct _info {
  GuppiData *data;
  GtkLabel *labels[LAST_LABEL];
};

static void
push_state (GuppiData *data, gpointer user_data)
{
  GuppiSeqScalar *seq = GUPPI_SEQ_SCALAR (data);
  info *x = (info *)user_data;
  gchar buf[64];
  double q1, q3;

  g_snprintf (buf, 64, "%g", guppi_seq_scalar_min (seq));
  gtk_label_set_text (x->labels[MIN_LABEL], buf);

  g_snprintf (buf, 64, "%g", guppi_seq_scalar_max (seq));
  gtk_label_set_text (x->labels[MAX_LABEL], buf);

  g_snprintf (buf, 64, "%g", guppi_seq_scalar_mean (seq));
  gtk_label_set_text (x->labels[MEAN_LABEL], buf);

  g_snprintf (buf, 64, "%g", guppi_seq_scalar_sdev (seq));
  gtk_label_set_text (x->labels[SDEV_LABEL], buf);

  g_snprintf (buf, 64, "%g", guppi_seq_scalar_median (seq));
  gtk_label_set_text (x->labels[MEDIAN_LABEL], buf);

  g_snprintf (buf, 64, "%g", q1 = guppi_seq_scalar_q1 (seq));
  gtk_label_set_text (x->labels[Q1_LABEL], buf);

  g_snprintf (buf, 64, "%g", q3 = guppi_seq_scalar_q3 (seq));
  gtk_label_set_text (x->labels[Q3_LABEL], buf);

  g_snprintf (buf, 64, "%g", q3 - q1);
  gtk_label_set_text (x->labels[IQR_LABEL], buf);
}

static void
destroy_cb (GtkWidget *w, gpointer user_data)
{
  info *x = (info *)user_data;

  gtk_signal_disconnect_by_data (GTK_OBJECT (x->data), x);
  guppi_unref0 (x->data);
}

static GtkWidget *
info_display (GuppiData *data)
{
  info *x = guppi_new0 (info, 1);
  GtkTable *table = GTK_TABLE (gtk_table_new (LAST_LABEL, 2, FALSE));
  GtkWidget *w;
  gint i;

  x->data = data;
  guppi_ref (data);

  for (i=0; i<LAST_LABEL; ++i) {
    w = gtk_label_new (_(label_names[i]));
    gtk_misc_set_alignment (GTK_MISC (w), 1.0, 0.5);
    gtk_table_attach (table, w, 0, 1, i, i+1, GTK_EXPAND | GTK_FILL, 0, 6, 1);
    gtk_widget_show (w);
    x->labels[i] = GTK_LABEL (gtk_label_new (""));
    gtk_table_attach_defaults (table, w = GTK_WIDGET (x->labels[i]),
			       1, 2, i, i+1);
    gtk_widget_show (w);
  }

  push_state (data, x);

  gtk_signal_connect (GTK_OBJECT (data), "changed",
		      GTK_SIGNAL_FUNC (push_state), x);
  gtk_signal_connect (GTK_OBJECT (table), "destroy",
		      GTK_SIGNAL_FUNC (destroy_cb), x);

  return GTK_WIDGET (table);
}

static xmlNodePtr
export_xml_element (GuppiSeq *seq, gint i, GuppiXMLDocument *doc)
{
  double x = guppi_seq_scalar_get (GUPPI_SEQ_SCALAR (seq), i);
  return guppi_xml_new_text_nodef (doc, "scalar", "%g", x);
}

static gboolean
import_xml_element (GuppiSeq *seq, GuppiXMLDocument *doc, xmlNodePtr node)
{
  xmlChar *s;
  double x;

  if (strcmp (node->name, "scalar"))
    return FALSE;

  s = xmlNodeGetContent (node);
  if (sscanf (s, "%lg", &x) != 1)
    return FALSE;

  guppi_seq_scalar_append (GUPPI_SEQ_SCALAR (seq), x);

  xmlFree (s);
  return TRUE;
}

/* ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** */

static void
guppi_seq_scalar_class_init (GuppiSeqScalarClass *klass)
{
  GtkObjectClass *object_class = (GtkObjectClass *) klass;
  GuppiDataClass *data_class = GUPPI_DATA_CLASS (klass);
  GuppiSeqClass *seq_class = GUPPI_SEQ_CLASS (klass);

  parent_class = gtk_type_class (GUPPI_TYPE_SEQ);

  klass->set         = set;
  klass->set_many    = set_many;
  klass->insert      = insert;
  klass->insert_many = insert_many;

  data_class->changed      = changed;
  data_class->info_display = info_display;

  seq_class->export_xml_element = export_xml_element;
  seq_class->import_xml_element = import_xml_element;

  object_class->finalize = guppi_seq_scalar_finalize;
}

static void
guppi_seq_scalar_init (GuppiSeqScalar *obj)
{
  obj->priv = guppi_new0 (GuppiSeqScalarPrivate, 1);
}

GtkType guppi_seq_scalar_get_type (void)
{
  static GtkType guppi_seq_scalar_type = 0;
  if (!guppi_seq_scalar_type) {
    static const GtkTypeInfo guppi_seq_scalar_info = {
      "GuppiSeqScalar",
      sizeof (GuppiSeqScalar),
      sizeof (GuppiSeqScalarClass),
      (GtkClassInitFunc) guppi_seq_scalar_class_init,
      (GtkObjectInitFunc) guppi_seq_scalar_init,
      NULL, NULL, (GtkClassInitFunc) NULL
    };
    guppi_seq_scalar_type = gtk_type_unique (GUPPI_TYPE_SEQ, &guppi_seq_scalar_info);
  }
  return guppi_seq_scalar_type;
}

/* $Id: guppi-seq-scalar.c,v 1.42 2002/01/14 05:01:17 trow Exp $ */


syntax highlighted by Code2HTML, v. 0.9.1