class Dobjects::Function

Function is a class that embeds two Dvectors, one for X data and one for Y data. It provides

facilities for sorting the X while keeping the Y matching, with sort and ::joint_sort;
to check if X data is sorted: sorted?, is_sorted;
interpolation, with compute_spline, compute_spline_data and interpolate
some functions for data access : x, y, point;
some utiliy functions: split_monotonic, strip_nan, reverse!
data inspection: min, max;
some computational functions: integrate, primitive, derivative, and now 4th-order accurate first and second derivatives: diff_5p and diff2_5p
utility for fuzzy operations, when the X values of two functions differ, but only slightly, of when points are missing: fuzzy_sub!
linear regression reglin
a function to approximate data using a low-order spline: spline_approximation

And getting bigger (almost) everyday…

Public Class Methods

joint_sort(x,y) click to toggle source

Sorts x, while ensuring that the corresponding y values keep matching. Should be pretty fast, as it is derived from glibc's quicksort.

a = Dvector[3,2,1]
b = a * 2                 -> [6,4,2]
Function.joint_sort(a,b)  -> [[1,2,3], [2,4,6]]

static VALUE function_joint_sort(VALUE self, VALUE x, VALUE y)
{
  long x_len, y_len;
  double * x_values = Dvector_Data_for_Write(x, &x_len);
  double * y_values = Dvector_Data_for_Write(y, &y_len);
  if(x_len != y_len)
    rb_raise(rb_eArgError,"both vectors must have the same size");
  else 
    {
      /* we temporarily freeze both Dvectors before sorting */
      FL_SET(x, DVEC_TMPLOCK);
      FL_SET(y, DVEC_TMPLOCK);
      joint_quicksort(x_values, y_values, (size_t) x_len);
      /* and unfreeze them */
      FL_UNSET(x, DVEC_TMPLOCK);
      FL_UNSET(y, DVEC_TMPLOCK);
    }
  /* we return the array of both Dvectors */
  return rb_ary_new3(2,x,y); 
}

new(x,y) click to toggle source

Creates a Function object with given x and y values.

static VALUE function_initialize(VALUE self, VALUE x, VALUE y)
{
  if(IS_A_DVECTOR(x) && IS_A_DVECTOR(y)) 
    {
      if(DVECTOR_SIZE(x) == DVECTOR_SIZE(y)) {
        set_x_vector(self, x);
        set_y_vector(self, y);
        /* fine, this could have been written in pure Ruby...*/
        set_spline_vector(self,Qnil);
        /* We initialize the @spline_cache var */
      }
      else
        rb_raise(rb_eArgError,"both vectors must have the same size");
    }
  else 
    rb_raise(rb_eArgError,"both arguments must be Dvector");
  return self;
}

Public Instance Methods

[](p1) click to toggle source

Returns a Dvector with two elements: the X and Y values of the point at the given index.

static VALUE function_point(VALUE self, VALUE index)
{
  if(! NUMERIC(index))
    rb_raise(rb_eArgError, "index has to be numeric");
  else
    {
      long i = NUM2LONG(index);
      long size = function_sanity_check(self);
      if(size > 0 && i < size)
        {
          VALUE point = rb_funcall(cDvector, idNew, 1, INT2NUM(2));
          double * dat = Dvector_Data_for_Write(point, NULL);
          double *x = Dvector_Data_for_Read(get_x_vector(self),NULL);
          double *y = Dvector_Data_for_Read(get_y_vector(self),NULL);
          dat[0] = x[i];
          dat[1] = y[i];
          return point;
        }
      else
        return Qnil;
    }
  return Qnil;
}

bound_values(xmin, xmax, ymin, ymax) click to toggle source

This function browses the points inside the Function and stores in the resulting new function only points which are within boundaries, and the points just next to them (so the general direction on the sides looks fine).

Make sure xmin < xmax and ymin < ymax, else you simply won't get any output.

static VALUE function_bound_values(VALUE self, 
                                   VALUE vxmin, VALUE vxmax,
                                   VALUE vymin, VALUE vymax)
{
  long ss = function_sanity_check(self);
  const double *xs = Dvector_Data_for_Read(get_x_vector(self),NULL);
  const double *ys = Dvector_Data_for_Read(get_y_vector(self),NULL);
  double xmin = NUM2DBL(vxmin);
  double xmax = NUM2DBL(vxmax);
  double ymin = NUM2DBL(vymin);
  double ymax = NUM2DBL(vymax);

  /* Now, two dvectors for writing: */
  VALUE x_out = rb_funcall(cDvector, idNew, 0);
  VALUE y_out = rb_funcall(cDvector, idNew, 0);

  /* No forward computation of the size of the targets, meaning
     memory allocation penalty.
  */
  
  int last_point_in = 0;        /* Whether the last point was in */
  long i;
  for(i = 0; i < ss; i++) {
    double x = xs[i];
    double y = ys[i];
    if( (xmin <= x) && (xmax >= x) && (ymin <= y) && (ymax >= y)) {
      if(! last_point_in) {
        last_point_in = 1;
        if(i) {                        /* Not for the first element */
          Dvector_Push_Double(x_out, xs[i-1]);
          Dvector_Push_Double(y_out, ys[i-1]);
        }
      }
      Dvector_Push_Double(x_out, x);
      Dvector_Push_Double(y_out, y);
    }
    else {                      /* Outside boundaries */
      if(last_point_in) {
        last_point_in = 0;
        Dvector_Push_Double(x_out, x);
        Dvector_Push_Double(y_out, y);
      }
    }
  }
  return Function_Create(x_out, y_out);
}

bounds() click to toggle source

Returns [xmin, ymin, xmax, ymax]

# File lib/Dobjects/Function_extras.rb, line 27
def bounds
  xmin,xmax = x.bounds
  ymin,ymax = y.bounds
  return [xmin, ymin, xmax, ymax]
end

compute_spline(p1) click to toggle source

Interpolates the value of the function at the points given. Returns a brand new Dvector. The X values must be sorted !

static VALUE function_compute_spline(VALUE self, VALUE x_values)
{
  VALUE x_vec = get_x_vector(self);
  VALUE y_vec = get_y_vector(self);
  VALUE cache;
  VALUE ret_val;
  long dat_size = function_sanity_check(self);
  long size = DVECTOR_SIZE(x_values);
  
  function_ensure_spline_data_present(self);

  cache = get_spline_vector(self);

  ret_val = rb_funcall(cDvector, rb_intern("new"),
                       1, LONG2NUM(size));
  double * x_dat = Dvector_Data_for_Read(x_vec,NULL);
  double * y_dat = Dvector_Data_for_Read(y_vec,NULL);
  double * spline = Dvector_Data_for_Read(cache,NULL);
  double * x = Dvector_Data_for_Read(x_values,NULL);
  double * y = Dvector_Data_for_Write(ret_val,NULL);
  
  function_compute_spline_interpolation(dat_size, x_dat,
                                        y_dat, spline,
                                        size, x, y);
  return ret_val;
}

compute_spline_data() click to toggle source

Computes spline data and caches it inside the object. Both X and Y vectors are cleared (see Dobjects::Dvector#clear) to make sure the cache is kept up-to-date. If the function is not sorted, sorts it.

static VALUE function_compute_spline_data(VALUE self)
{
  VALUE x_vec = get_x_vector(self);
  VALUE y_vec = get_y_vector(self);
  VALUE cache = get_spline_vector(self);
  long size = DVECTOR_SIZE(x_vec);

  if(DVECTOR_SIZE(y_vec) != size)
    rb_raise(rb_eRuntimeError, 
             "x and y should have the same size !");
  if(! IS_A_DVECTOR(cache))    /* create it -- and silently ignores
                                  its previous values */
      cache = rb_funcall(cDvector, idNew,
                         1, LONG2NUM(size));
  if(DVECTOR_SIZE(cache) != size) /* switch to the required size for cache */
    Dvector_Data_Resize(cache, size);

  /* we make sure that the X values are sorted */
  if(! RTEST(function_is_sorted(self)))
     function_sort(self);
  
  double * x, *y, *spline;
  x = Dvector_Data_for_Read(x_vec, NULL);
  y = Dvector_Data_for_Read(y_vec, NULL);
  spline = Dvector_Data_for_Write(cache, NULL);

  function_fill_second_derivatives(size, x, y, spline,1.0/0.0, 1.0/0.0);
  set_spline_vector(self, cache);

  /* now, we clear both X and Y */
  DVECTOR_CLEAR(x_vec);
  DVECTOR_CLEAR(y_vec);
  return self;
}

derivative() click to toggle source

Computes the derivative of the Function and returns it as a new Function. The newly created function shares the X vector with the previous one.

WARNING: this is a very naive 3-points algorithm; you should consider using #diff_5p

static VALUE function_derivative(VALUE self)
{
  long size = function_sanity_check(self);
  const double *x = Dvector_Data_for_Read(get_x_vector(self),NULL);
  const double *y = Dvector_Data_for_Read(get_y_vector(self),NULL);
  VALUE derivative = Dvector_Create();
  long i = 0;
  /* First value */
  Dvector_Push_Double(derivative, (y[i+1] - y[i]) /(x[i+1] - x[i]));
  i++;
  while(i < (size - 1))
    {
      Dvector_Push_Double(derivative, 
                          .5 * (
                                (y[i+1] - y[i]) /(x[i+1] - x[i]) + 
                                (y[i] - y[i-1]) /(x[i] - x[i-1])
                                ));
      i++;
    }
  Dvector_Push_Double(derivative, (y[i] - y[i-1]) /(x[i] - x[i-1]));
  return Function_Create(get_x_vector(self), derivative);
}

diff2_5p() click to toggle source

Computes a 4th order accurate second derivative of the Function.

This function requires that there are at the very least 5 data points!

static VALUE function_diff2_5p(VALUE self)
{
  long size = function_sanity_check(self);
  const double *x = Dvector_Data_for_Read(get_x_vector(self),NULL);
  const double *y = Dvector_Data_for_Read(get_y_vector(self),NULL);
  VALUE derivative = Dvector_Create();
  long i = 0;
  double delta_1, delta_2, delta_3, delta_4;
  double alpha_1, alpha_2, alpha_3, alpha_4;
  double v0,v1,v2,v3,v4;

  for(i = 0; i < size; i++) {
    /* First initialize values, though this is very suboptimal */
    v0 = y[i];
    if(i == 0) {
      delta_1 = x[1] - x[0]; v1 = y[1];
      delta_2 = x[2] - x[0]; v2 = y[2];
      delta_3 = x[3] - x[0]; v3 = y[3];
      delta_4 = x[4] - x[0]; v4 = y[4];
    } else if(i == 1) {
      delta_1 = x[0] - x[1]; v1 = y[0];
      delta_2 = x[2] - x[1]; v2 = y[2];
      delta_3 = x[3] - x[1]; v3 = y[3];
      delta_4 = x[4] - x[1]; v4 = y[4];
    } else if(i == size - 2) {
      delta_1 = x[size-1] - x[size-2]; v1 = y[size-1];
      delta_2 = x[size-3] - x[size-2]; v2 = y[size-3];
      delta_3 = x[size-4] - x[size-2]; v3 = y[size-4];
      delta_4 = x[size-5] - x[size-2]; v4 = y[size-5];
    } else if(i == size - 1) {
      delta_1 = x[size-2] - x[size-1]; v1 = y[size-2];
      delta_2 = x[size-3] - x[size-1]; v2 = y[size-3];
      delta_3 = x[size-4] - x[size-1]; v3 = y[size-4];
      delta_4 = x[size-5] - x[size-1]; v4 = y[size-5];
    } else {
      delta_1 = x[i-2] - x[i]; v1 = y[i-2];
      delta_2 = x[i-1] - x[i]; v2 = y[i-1];
      delta_3 = x[i+2] - x[i]; v3 = y[i+2];
      delta_4 = x[i+1] - x[i]; v4 = y[i+1];
    }
    alpha_1 = -2 * (delta_2*delta_3 + delta_2*delta_4 + delta_3*delta_4)/
      (delta_1 * (delta_2 - delta_1) * (delta_3 - delta_1) 
       * (delta_4 - delta_1));
    alpha_2 = -2 * (delta_1*delta_3 + delta_1*delta_4 + delta_3*delta_4)/
      (delta_2 * (delta_1 - delta_2) * (delta_3 - delta_2) 
       * (delta_4 - delta_2));
    alpha_3 = -2 * (delta_2*delta_1 + delta_2*delta_4 + delta_1*delta_4)/
      (delta_3 * (delta_1 - delta_3) * (delta_2 - delta_3) 
       * (delta_4 - delta_3));
    alpha_4 = -2 * (delta_2*delta_3 + delta_2*delta_1 + delta_3*delta_1)/
      (delta_4 * (delta_1 - delta_4) * (delta_2 - delta_4) 
       * (delta_3 - delta_4));
    Dvector_Push_Double(derivative,
                        -(alpha_1 + alpha_2 + alpha_3 + alpha_4) * v0 +
                        alpha_1 * v1 + alpha_2 * v2 + 
                        alpha_3 * v3 + alpha_4 * v4);
  }
  return Function_Create(get_x_vector(self), derivative);
}

diff_5p() click to toggle source

Computes a 4th order accurate derivative of the Function.

This function requires that there are at the very least 5 data points !

static VALUE function_diff_5p(VALUE self)
{
  long size = function_sanity_check(self);
  const double *x = Dvector_Data_for_Read(get_x_vector(self),NULL);
  const double *y = Dvector_Data_for_Read(get_y_vector(self),NULL);
  VALUE derivative = Dvector_Create();
  long i = 0;
  double delta_1, delta_2, delta_3, delta_4;
  double alpha_1, alpha_2, alpha_3, alpha_4;
  double v0,v1,v2,v3,v4;
  /* TODO: what happens when there are less than 5 points ? */

  for(i = 0; i < size; i++) {
    /* First initialize values, though this is very suboptimal */
    v0 = y[i];
    if(i == 0) {
      delta_1 = x[1] - x[0]; v1 = y[1];
      delta_2 = x[2] - x[0]; v2 = y[2];
      delta_3 = x[3] - x[0]; v3 = y[3];
      delta_4 = x[4] - x[0]; v4 = y[4];
    } else if(i == 1) {
      delta_1 = x[0] - x[1]; v1 = y[0];
      delta_2 = x[2] - x[1]; v2 = y[2];
      delta_3 = x[3] - x[1]; v3 = y[3];
      delta_4 = x[4] - x[1]; v4 = y[4];
    } else if(i == size - 2) {
      delta_1 = x[size-1] - x[size-2]; v1 = y[size-1];
      delta_2 = x[size-3] - x[size-2]; v2 = y[size-3];
      delta_3 = x[size-4] - x[size-2]; v3 = y[size-4];
      delta_4 = x[size-5] - x[size-2]; v4 = y[size-5];
    } else if(i == size - 1) {
      delta_1 = x[size-2] - x[size-1]; v1 = y[size-2];
      delta_2 = x[size-3] - x[size-1]; v2 = y[size-3];
      delta_3 = x[size-4] - x[size-1]; v3 = y[size-4];
      delta_4 = x[size-5] - x[size-1]; v4 = y[size-5];
    } else {
      delta_1 = x[i-2] - x[i]; v1 = y[i-2];
      delta_2 = x[i-1] - x[i]; v2 = y[i-1];
      delta_3 = x[i+2] - x[i]; v3 = y[i+2];
      delta_4 = x[i+1] - x[i]; v4 = y[i+1];
    }
    alpha_1 = delta_2*delta_3*delta_4/
      (delta_1 * (delta_2 - delta_1) * (delta_3 - delta_1) 
       * (delta_4 - delta_1));
    alpha_2 = delta_1*delta_3*delta_4/
      (delta_2 * (delta_1 - delta_2) * (delta_3 - delta_2) 
       * (delta_4 - delta_2));
    alpha_3 = delta_1*delta_2*delta_4/
      (delta_3 * (delta_1 - delta_3) * (delta_2 - delta_3) 
       * (delta_4 - delta_3));
    alpha_4 = delta_1*delta_2*delta_3/
      (delta_4 * (delta_1 - delta_4) * (delta_2 - delta_4) 
       * (delta_3 - delta_4));
    Dvector_Push_Double(derivative,
                        -(alpha_1 + alpha_2 + alpha_3 + alpha_4) * v0 +
                        alpha_1 * v1 + alpha_2 * v2 + 
                        alpha_3 * v3 + alpha_4 * v4);
  }
  return Function_Create(get_x_vector(self), derivative);
}

distance(x,y) → a_number click to toggle source

distance(x,y, xscale, yscale) → a_number

Returns the distance of the function to the given point. Optionnal xscale and yscale says by how much we should divide the x and y coordinates before computing the distance. Use it if the distance is not homogeneous.

static VALUE function_distance(int argc, VALUE *argv, VALUE self)
{
  switch(argc)
    {
    case 2:
      return rb_float_new(private_function_distance(self, 
                                                    NUM2DBL(argv[0]),
                                                    NUM2DBL(argv[1]),
                                                    1.0,1.0,NULL));
    case 4:
      return rb_float_new(private_function_distance(self, 
                                                    NUM2DBL(argv[0]),
                                                    NUM2DBL(argv[1]),
                                                    NUM2DBL(argv[2]),
                                                    NUM2DBL(argv[3]),
                                                    NULL));
    default:
      rb_raise(rb_eArgError, "distance should have 2 or 4 parameters");
    }
  return Qnil;
}

each do |x,y| _code_ end click to toggle source

Iterates over all the points in the Function, yielding X and Y for each point.

static VALUE function_each(VALUE self) /* :yields: x,y */
{

  long x_len, y_len;
  VALUE x = get_x_vector(self);
  VALUE y = get_y_vector(self);
  double * x_values = Dvector_Data_for_Write(x, &x_len);
  double * y_values = Dvector_Data_for_Write(y, &y_len);
  if(x_len != y_len)
    rb_raise(rb_eRuntimeError,"X and Y must have the same size");
  else 
    {
      /* we temporarily freeze both Dvectors during iteration */
      FL_SET(x, DVEC_TMPLOCK);
      FL_SET(y, DVEC_TMPLOCK);
      while(x_len--)
        {
          VALUE flt_x = rb_float_new(*x_values++);
          VALUE flt_y = rb_float_new(*y_values++);
          rb_yield_values(2, flt_x, flt_y);
        }
      /* and unfreeze them */
      FL_UNSET(x, DVEC_TMPLOCK);
      FL_UNSET(y, DVEC_TMPLOCK);
    }
  return self; /* nothing interesting */
  
}

fuzzy_sub!(p1) click to toggle source

Fuzzy substraction of two curves. Substracts the Y values of op to the current Function, by making sure that the Y value substracted to a given point corresponds to the closest X_ value of the point in op. This function somehow assumes that the data is reasonably organised, and will never go backwards to find a matching X value in op.

In any case, you really should consider using #split_monotonic on it first.

static VALUE function_fuzzy_substract(VALUE self, VALUE op)
{
  long ss = function_sanity_check(self);
  const double *xs = Dvector_Data_for_Read(get_x_vector(self),NULL);
  double *ys = Dvector_Data_for_Write(get_y_vector(self),NULL);
  long so = function_sanity_check(op);
  const double *xo = Dvector_Data_for_Read(get_x_vector(op),NULL);
  const double *yo = Dvector_Data_for_Read(get_y_vector(op),NULL);
  long i,j = 0;
  double diff;
  double fuzz = 0;              /* The actual sum of the terms */
  
  for(i = 0; i < ss; i++) 
    {
      /* We first look for the closest point */
      diff = fabs(xs[i] - xo[j]);
      while((j < (so - 1)) && (fabs(xs[i] - xo[j+1]) <  diff))
        diff = fabs(xs[i] - xo[++j]);
      fuzz += diff;
      ys[i] -= yo[j];
    }
  return rb_float_new(fuzz);
}

integrate() → value click to toggle source

integrate(start_index, end_index) → value

Returns the value of the integral of the function between the two indexes given, or over the whole function if no indexes are specified.

static VALUE function_integrate(int argc, VALUE *argv, VALUE self)
{
  long start,end;
  switch(argc) 
    {
    case 0:
      start = 0;
      end = function_sanity_check(self) - 1; 
      break;
    case 2:
      start = NUM2LONG(argv[0]);
      end = NUM2LONG(argv[1]);
      break;
    default:
      rb_raise(rb_eArgError, "integrate should have 0 or 2 parameters");
    }
  return rb_float_new(private_function_integrate(self,start,end));
}

interpolate(x_values) click to toggle source

interpolate(a_number)

Computes interpolated values of the data contained in f and returns a Function object holding both x_values and the computed Y values. x_values will be sorted if necessary.

With the second form, specify only the number of points, and the function will construct the appropriate vector with equally spaced points within the function range.

static VALUE function_interpolate(VALUE self, VALUE x_values)
{
  if(NUMERIC(x_values))
    {
      /* we're in the second case, although I sincerely doubt it would
         come useful 
      */
      long size,i;
      /* we make sure the function is sorted */
      function_ensure_sorted(self);
      double * data;
      double x_min;
      double x_max;
      data = Dvector_Data_for_Read(get_x_vector(self), &size);
      x_min = *data;
      x_max = *(data + size -1);
      x_values = rb_funcall(cDvector, idNew, 1, x_values);
      data = Dvector_Data_for_Write(x_values, &size);
      for(i = 0;i < size; i++)
        data[i] = x_min + ((x_max - x_min)/((double) (size-1))) * i;
    }
  if(! IS_A_DVECTOR(x_values))
    rb_raise(rb_eArgError, "x_values should be a Dvector or a number");
  else 
    {
      /* sort x_values */
      if(! dvector_is_sorted(x_values))
        rb_funcall(x_values, idSort,0);
      VALUE y_values = function_compute_spline(self, x_values);
      return rb_funcall(cFunction, idNew, 2, x_values, y_values);
    }
  return Qnil;
}

is_sorted()

Alias for: sorted?

length()

Alias for: size

make_interpolant() click to toggle source

Returns an interpolant that can be fed to Special_Paths#append_interpolant_to_path to make nice splines.

Can be used this way:

f = Function.new(x,y)
t.append_interpolant_to_path(f.make_interpolant)
t.stroke

static VALUE function_make_interpolant(VALUE self)
{
  VALUE x_vec = get_x_vector(self);
  VALUE y_vec = get_y_vector(self);
  VALUE cache;
  VALUE a_vec,b_vec,c_vec;
  VALUE ret_val;
  double *x, *y, *a, *b, *c, *y2;
  double delta_x;
  long size = function_sanity_check(self);
  long i;
  
  function_ensure_spline_data_present(self);

  cache = get_spline_vector(self);
  x = Dvector_Data_for_Read(x_vec,NULL);
  y = Dvector_Data_for_Read(y_vec,NULL);
  y2 = Dvector_Data_for_Read(cache,NULL);

  a_vec  = rb_funcall(cDvector, idNew, 1, LONG2NUM(size));
  a = Dvector_Data_for_Write(a_vec, NULL);
  b_vec  = rb_funcall(cDvector, idNew, 1, LONG2NUM(size));
  b = Dvector_Data_for_Write(b_vec, NULL);
  c_vec  = rb_funcall(cDvector, idNew, 1, LONG2NUM(size));
  c = Dvector_Data_for_Write(c_vec, NULL);

  /* from my computations, the formula is the following:
     A = (y_2n+1 - y_2n)/(6 * delta_x)
     B = 0.5 * y_2n
     C = (y_n+1 - y_n)/delta_x - (2 * y_2n + y_2n+1) * delta_x/6
  */

  for(i = 0; i < size - 1; i++)
    {
      delta_x = x[i+1] - x[i];
      a[i] = (y2[i+1] - y2[i]) / (6.0 * delta_x);
      b[i] = 0.5 * y2[i];
      c[i] = (y[i+1] - y[i])/delta_x - 
        (2 * y2[i] + y2[i+1]) * (delta_x / 6.0);
    }
  a[i] = b[i] = c[i] = 0.0;
  ret_val = rb_ary_new();
  rb_ary_push(ret_val, x_vec);
  rb_ary_push(ret_val, y_vec);
  rb_ary_push(ret_val, a_vec);
  rb_ary_push(ret_val, b_vec);
  rb_ary_push(ret_val, c_vec);

  return ret_val;
}

max() click to toggle source

Returns the point where Y is the maximum

# File lib/Dobjects/Function_extras.rb, line 39
def max
  return point(y.where_max)
end

min() click to toggle source

Returns the point where Y is the minimum

# File lib/Dobjects/Function_extras.rb, line 34
def min
  return point(y.where_min)
end

point(p1) click to toggle source

Returns a Dvector with two elements: the X and Y values of the point at the given index.

static VALUE function_point(VALUE self, VALUE index)
{
  if(! NUMERIC(index))
    rb_raise(rb_eArgError, "index has to be numeric");
  else
    {
      long i = NUM2LONG(index);
      long size = function_sanity_check(self);
      if(size > 0 && i < size)
        {
          VALUE point = rb_funcall(cDvector, idNew, 1, INT2NUM(2));
          double * dat = Dvector_Data_for_Write(point, NULL);
          double *x = Dvector_Data_for_Read(get_x_vector(self),NULL);
          double *y = Dvector_Data_for_Read(get_y_vector(self),NULL);
          dat[0] = x[i];
          dat[1] = y[i];
          return point;
        }
      else
        return Qnil;
    }
  return Qnil;
}

primitive() click to toggle source

Computes the primitive of the Function (whose value for the first point is 0) and returns it as a new Function. The newly created function shares the X vector with the previous one.

static VALUE function_primitive(VALUE self)
{
  long size = function_sanity_check(self);
  const double *x = Dvector_Data_for_Read(get_x_vector(self),NULL);
  const double *y = Dvector_Data_for_Read(get_y_vector(self),NULL);
  VALUE primitive = Dvector_Create();
  long i = 0;
  double val = 0;
  while(i < (size - 1))
    {
      Dvector_Push_Double(primitive, val);
      val += (y[i] + y[i+1]) * (x[i+1] - x[i]) * 0.5;
      i++;
    }
  Dvector_Push_Double(primitive, val);
  return Function_Create(get_x_vector(self), primitive);
}

reglin(*args) click to toggle source

Performs a linear regression of the Function; returns the pair

[ a, b]

where f(x) = a*x + b

if the optional arguments first and last are provided, they represent the indices of the first and last elements.

static VALUE function_reglin(int argc, VALUE *argv, VALUE self)
{
  long len = function_sanity_check(self);
  const double *x = Dvector_Data_for_Read(get_x_vector(self),NULL);
  const double *y = Dvector_Data_for_Read(get_y_vector(self),NULL);
  VALUE ret = rb_funcall(cDvector, idNew, 1, INT2NUM(2));
  double * dat = Dvector_Data_for_Write(ret, NULL);
  long nb;
  if(argc == 2) {
    long f = NUM2LONG(argv[0]);
    long l = NUM2LONG(argv[1]);
    if(f < 0)
      f = len + f;
    if(l < 0)
      l = len + l;
    x += f;
    y += f;
    nb = l - f;
  }
  else if(argc == 0) {
    nb = len;
  }
  else {
    rb_raise(rb_eArgError, "reglin should have 0 or 2 parameters");
  }
  reglin(x,y,nb,dat,dat+1);
  return ret;
}

reverse!() click to toggle source

Reverses the function. Equivalent to doing 

x.reverse!
y.reverse!

excepted that it is faster (though not much faster).

static VALUE function_reverse(VALUE self)
{
  long len = function_sanity_check(self);
  double *xs = Dvector_Data_for_Write(get_x_vector(self),NULL);
  double *ys = Dvector_Data_for_Write(get_y_vector(self),NULL);
  
  double *xe = xs+len-1;
  double *ye = ys+len-1;
  double tmp;
  long i;
  for(i = 0; i < len/2; i++, xs++, ys++, xe--, ye--) {
    tmp = *xe; *xe = *xs; *xs = tmp;
    tmp = *ye; *ye = *ys; *ys = tmp;
  }
  return self;
}

size() click to toggle source

Returns the number of points inside the function.

static VALUE function_size(VALUE self)
{
  long size = function_sanity_check(self);
  return LONG2NUM(size);
}

Also aliased as: length

smooth_pick(*args) click to toggle source

Attempts to pick a smooth value for a point, according to the algorithm implented for “smooth” markers in Soas. See DOI: 10.1016/j.bioelechem.2009.02.010

Warning: be wary of this function as it will return a correct value only for rather noisy data !

static VALUE function_smooth_pick(int argc, VALUE *argv, VALUE self)
{
  long len = function_sanity_check(self);
  const double *x = Dvector_Data_for_Read(get_x_vector(self),NULL);
  const double *y = Dvector_Data_for_Read(get_y_vector(self),NULL);
  long idx;
  long range;
  switch(argc) {
  case 2:
    range = NUM2LONG(argv[1]);
    break;
  case 1:
    range = len > 500 ? 50 : len/10;
    break;
  default:
    rb_raise(rb_eArgError, "smooth_a=t should have 1 or 2 parameters");
  }
  idx = NUM2LONG(argv[0]);
  if(idx < 0)
    idx = len + idx;
  return rb_float_new(smooth_pick(x,y,len,idx,range));
}

sort() click to toggle source

Sorts the X values while keeping the matching Y values.

static VALUE function_sort(VALUE self)
{
  return function_joint_sort(self,get_x_vector(self), get_y_vector(self));
}

sorted?() click to toggle source

Checks if the X values of the Function are sorted.

static VALUE function_is_sorted(VALUE self)
{
  if(dvector_is_sorted(get_x_vector(self)))
    return Qtrue;
  else
    return Qfalse;
}

Also aliased as: is_sorted

spline_approximation(p1) click to toggle source

Filters the Function through interpolation. params holds a hash with the following values:

??

It returns a hash.

static VALUE function_spline_approximation(VALUE self, VALUE params)
{
  long len = function_sanity_check(self);
  const double *x = Dvector_Data_for_Read(get_x_vector(self),NULL);
  const double *y = Dvector_Data_for_Read(get_y_vector(self),NULL);
  VALUE xiret, yiret, y2iret, yintret,ret;
  double * xi, *yi, *y2i, *yint;
  long nbavg = 9;  
  long nbmax = 20;
  if(RTEST(rb_hash_aref(params, rb_str_new2("number"))))
    nbmax = NUM2LONG(rb_hash_aref(params, rb_str_new2("number")));
  if(RTEST(rb_hash_aref(params, rb_str_new2("average"))))
    nbavg = NUM2LONG(rb_hash_aref(params, rb_str_new2("average")));

  /* TODO: add checks that monotonic and growing. */
  
  xiret = rb_funcall(cDvector, idNew, 1, INT2NUM(nbmax)); 
  xi = Dvector_Data_for_Write(xiret, NULL);
  yiret = rb_funcall(cDvector, idNew, 1, INT2NUM(nbmax)); 
  yi = Dvector_Data_for_Write(yiret, NULL);
  y2iret = rb_funcall(cDvector, idNew, 1, INT2NUM(nbmax)); 
  y2i = Dvector_Data_for_Write(y2iret, NULL);
  yintret = rb_funcall(cDvector, idNew, 1, INT2NUM(len)); 
  yint = Dvector_Data_for_Write(yintret, NULL);

  internal_spline_approximation(x, y, len, xi, yi, y2i,
                                nbmax, nbavg, yint);
  ret = rb_hash_new();
  rb_hash_aset(ret, rb_str_new2("xi"), xiret);
  rb_hash_aset(ret, rb_str_new2("yi"), yiret);
  rb_hash_aset(ret, rb_str_new2("y2i"), y2iret);
  rb_hash_aset(ret, rb_str_new2("y"), yintret);
  return ret;
}

split_monotonic() click to toggle source

Splits the function into strictly monotonic sub-functions. Returns the array of the subfunctions. The returned values are necessarily new values.

static VALUE function_split_monotonic(VALUE self)
{
  VALUE ret = rb_ary_new();
  VALUE cur_x = Dvector_Create();
  VALUE cur_y = Dvector_Create();

  long size = function_sanity_check(self);
  long i;
  if(size < 2)
    rb_raise(rb_eRuntimeError, "Function needs to have at least 2 points");

  double *x = Dvector_Data_for_Read(get_x_vector(self),NULL);
  double *y = Dvector_Data_for_Read(get_y_vector(self),NULL);

  double last_x;
  double direction; /* -1 if down, +1 if up, so that the product of 
                       (x - last_x) with direction should always be positive
                    */
  VALUE f;
                     
                       
  /* bootstrap */
  if(x[1] > x[0])
    direction = 1;
  else
    direction = -1;
  last_x = x[1];
  for(i = 0; i < 2; i++)
    {
      Dvector_Push_Double(cur_x, x[i]);
      Dvector_Push_Double(cur_y, y[i]);
    }

  for(i = 2; i < size; i++) 
    {
      if(direction * (x[i] - last_x) <= 0) 
        {
          /* we need to add a new set of Dvectors */
          f = Function_Create(cur_x, cur_y);
          rb_ary_push(ret, f);
          cur_x = Dvector_Create();
          cur_y = Dvector_Create();
          /* We don't store the previous point if 
           the X value is the same*/
          if(x[i] != last_x) 
            {
              Dvector_Push_Double(cur_x, x[i-1]);
              Dvector_Push_Double(cur_y, y[i-1]);
            }
          direction *= -1;
        }
      /* store the current point */
      Dvector_Push_Double(cur_x, x[i]);
      Dvector_Push_Double(cur_y, y[i]);
      last_x = x[i];
    }
  f = Function_Create(cur_x, cur_y);
  rb_ary_push(ret, f);
  return ret;
}

split_on_nan(p1) click to toggle source

Splits the function on NaN values for x, y or xy, depending on whether sym is :x, :y or :xy (or, as a matter of fact, anything else than :x or :y).

This returns an array of new Function objects.

This function will return empty Function objects between consecutive NaN values.

static VALUE function_split_on_nan(VALUE self, VALUE sym)
{
  VALUE ret = rb_ary_new();
  VALUE cur_x = Dvector_Create();
  VALUE cur_y = Dvector_Create();
  int on_x = 1;
  int on_y = 1;
  long size = function_sanity_check(self);
  long cur_size = 0;
  long i;
  if(size < 2)
    rb_raise(rb_eRuntimeError, "Function needs to have at least 2 points");

  double *x = Dvector_Data_for_Read(get_x_vector(self),NULL);
  double *y = Dvector_Data_for_Read(get_y_vector(self),NULL);

  VALUE f;
  
  if(sym == ID2SYM(rb_intern("x")))
    on_y = 0;
  else if(sym == ID2SYM(rb_intern("y")))
    on_x = 0;


  for(i = 0; i < size; i++) {
    if((on_x && isnan(x[i])) ||
       (on_y && isnan(y[i]))) {
      /* We split */
      f = Function_Create(cur_x, cur_y);
      rb_ary_push(ret, f);
      cur_x = Dvector_Create();
      cur_y = Dvector_Create();
    }
    else {
      Dvector_Push_Double(cur_x, x[i]);
      Dvector_Push_Double(cur_y, y[i]);
    }
  }
  f = Function_Create(cur_x, cur_y);
  rb_ary_push(ret, f);
  return ret;
}

strip_nan() click to toggle source

Strips all the points containing NaN values from the function, and returns the number of points stripped.

static VALUE function_strip_nan(VALUE self)
{
  long size = function_sanity_check(self);
  long nb_stripped = 0;
  long i;

  double *x = Dvector_Data_for_Write(get_x_vector(self),NULL);
  double *y = Dvector_Data_for_Write(get_y_vector(self),NULL);
  for( i = 0; i < size; i++)
    {
      if(isnan(x[i]) || isnan(y[i]))
        nb_stripped ++;
      else
        {
          x[i - nb_stripped] = x[i];
          y[i - nb_stripped] = y[i];
        }
    }
  if(nb_stripped)
    {
      Dvector_Data_Resize(get_x_vector(self), size - nb_stripped);
      Dvector_Data_Resize(get_y_vector(self), size - nb_stripped);
    }
  return INT2NUM(nb_stripped);
}

x() click to toggle source

The X vector.

static VALUE get_x_vector(VALUE self) 
{
  return rb_iv_get(self, X_VAL);
}

y() click to toggle source

The Y vector.

static VALUE get_y_vector(VALUE self) 
{
  return rb_iv_get(self, Y_VAL);
}