Class: Dobjects::Function
- Inherits:
-
Object
- Object
- Dobjects::Function
- Defined in:
- ext/Dobjects/Function/function.c,
lib/Dobjects/Function_extras.rb,
ext/Dobjects/Function/function.c
Overview
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 Function.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…
Class Method Summary collapse
-
.joint_sort(x, y) ⇒ Object
Sorts
x
, while ensuring that the correspondingy
values keep matching.
Instance Method Summary collapse
-
#[](index) ⇒ Object
Returns a Dvector with two elements: the X and Y values of the point at the given index.
-
#bound_values(xmin, xmax, ymin, ymax) ⇒ Object
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).
-
#bounds ⇒ Object
Returns [xmin, ymin, xmax, ymax].
-
#compute_spline(x_values) ⇒ Object
Interpolates the value of the function at the points given.
-
#compute_spline_data ⇒ Object
Computes spline data and caches it inside the object.
-
#derivative ⇒ Object
Computes the derivative of the Function and returns it as a new Function.
-
#diff2_5p ⇒ Object
Computes a 4th order accurate second derivative of the Function.
-
#diff_5p ⇒ Object
Computes a 4th order accurate derivative of the Function.
-
#distance(*args) ⇒ Object
Returns the distance of the function to the given point.
-
#each ⇒ Object
:yields: x,y.
-
#fuzzy_sub!(op) ⇒ Object
Fuzzy substraction of two curves.
-
#new(x, y) ⇒ Object
constructor
Creates a Function object with given
x
andy
values. -
#integrate(*args) ⇒ Object
:call-seq: f.integrate() -> value f.integrate(start_index, end_index) -> value.
-
#interpolate(x_values) ⇒ Object
Computes interpolated values of the data contained in
f
and returns a Function object holding bothx_values
and the computed Y values. -
#make_interpolant ⇒ Object
Returns an interpolant that can be fed to Special_Paths#append_interpolant_to_path to make nice splines.
-
#max ⇒ Object
Returns the point where Y is the maximum.
-
#min ⇒ Object
Returns the point where Y is the minimum.
-
#point(index) ⇒ Object
Returns a Dvector with two elements: the X and Y values of the point at the given index.
-
#primitive ⇒ Object
Computes the primitive of the Function (whose value for the first point is 0) and returns it as a new Function.
-
#reglin(*args) ⇒ Object
Performs a linear regression of the Function; returns the pair [ a, b] where f(x) = a*x + b.
-
#reverse! ⇒ Object
Reverses the function.
-
#size ⇒ Object
(also: #length)
Returns the number of points inside the function.
-
#smooth_pick(*args) ⇒ Object
Attempts to pick a smooth value for a point, according to the algorithm implented for “smooth” markers in Soas.
-
#sort ⇒ Object
Sorts the X values while keeping the matching Y values.
-
#sorted? ⇒ Boolean
(also: #is_sorted)
Checks if the X values of the Function are sorted.
-
#spline_approximation(params) ⇒ Object
Filters the Function through interpolation.
-
#split_monotonic ⇒ Object
Splits the function into strictly monotonic sub-functions.
-
#split_on_nan(sym) ⇒ Object
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
). -
#strip_nan ⇒ Object
Strips all the points containing NaN values from the function, and returns the number of points stripped.
- #x ⇒ Object
- #y ⇒ Object
Constructor Details
#new(x, y) ⇒ Object
Creates a Function object with given x
and y
values.
156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 |
# File 'ext/Dobjects/Function/function.c', line 156
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;
}
|
Class Method Details
.joint_sort(x, y) ⇒ Object
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]]
500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 |
# File 'ext/Dobjects/Function/function.c', line 500
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);
}
|
Instance Method Details
#[](index) ⇒ Object
Returns a Dvector with two elements: the X and Y values of the point at the given index.
781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 |
# File 'ext/Dobjects/Function/function.c', line 781
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) ⇒ Object
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.
1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 1210 1211 1212 1213 1214 1215 1216 1217 1218 1219 1220 |
# File 'ext/Dobjects/Function/function.c', line 1175
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 ⇒ Object
Returns [xmin, ymin, xmax, ymax]
27 28 29 30 31 |
# 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(x_values) ⇒ Object
Interpolates the value of the function at the points given.
Returns a brand new Dvector. The X values must be sorted !
391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 |
# File 'ext/Dobjects/Function/function.c', line 391
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 ⇒ Object
Computes spline data and caches it inside the object. Both X and Y vectors are cleared (see Dvector#clear) to make sure the cache is kept up-to-date. If the function is not sorted, sorts it.
282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 |
# File 'ext/Dobjects/Function/function.c', line 282
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 ⇒ Object
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
964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 |
# File 'ext/Dobjects/Function/function.c', line 964
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 ⇒ Object
Computes a 4th order accurate second derivative of the Function.
This function requires that there are at the very least 5 data points!
1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 |
# File 'ext/Dobjects/Function/function.c', line 1060
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 ⇒ Object
Computes a 4th order accurate derivative of the Function.
This function requires that there are at the very least 5 data points !
993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 |
# File 'ext/Dobjects/Function/function.c', line 993
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) ⇒ Numeric #distance(x, y, xscale, yscale) ⇒ Numeric
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.
861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 |
# File 'ext/Dobjects/Function/function.c', line 861
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 ⇒ Object
:yields: x,y
528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 |
# File 'ext/Dobjects/Function/function.c', line 528
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!(op) ⇒ Object
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.
1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 |
# File 'ext/Dobjects/Function/function.c', line 1139
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(*args) ⇒ Object
:call-seq:
f.integrate() -> value
f.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.
915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 |
# File 'ext/Dobjects/Function/function.c', line 915
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) ⇒ Object #interpolate(a_number) ⇒ Object
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.
582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 |
# File 'ext/Dobjects/Function/function.c', line 582
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;
}
|
#make_interpolant ⇒ Object
429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 |
# File 'ext/Dobjects/Function/function.c', line 429
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 ⇒ Object
Returns the point where Y is the maximum
39 40 41 |
# File 'lib/Dobjects/Function_extras.rb', line 39 def max return point(y.where_max) end |
#min ⇒ Object
Returns the point where Y is the minimum
34 35 36 |
# File 'lib/Dobjects/Function_extras.rb', line 34 def min return point(y.where_min) end |
#point(index) ⇒ Object
Returns a Dvector with two elements: the X and Y values of the point at the given index.
781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 |
# File 'ext/Dobjects/Function/function.c', line 781
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 ⇒ Object
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.
939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 |
# File 'ext/Dobjects/Function/function.c', line 939
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) ⇒ Object
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.
1282 1283 1284 1285 1286 1287 1288 1289 1290 1291 1292 1293 1294 1295 1296 1297 1298 1299 1300 1301 1302 1303 1304 1305 1306 1307 1308 1309 |
# File 'ext/Dobjects/Function/function.c', line 1282
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! ⇒ Object
Reverses the function. Equivalent to doing
x.reverse!
y.reverse!
excepted that it is faster (though not *much* faster).
1229 1230 1231 1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 |
# File 'ext/Dobjects/Function/function.c', line 1229
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 ⇒ Object Also known as: length
Returns the number of points inside the function.
1123 1124 1125 1126 1127 |
# File 'ext/Dobjects/Function/function.c', line 1123
static VALUE function_size(VALUE self)
{
long size = function_sanity_check(self);
return LONG2NUM(size);
}
|
#smooth_pick(*args) ⇒ Object
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 !
1381 1382 1383 1384 1385 1386 1387 1388 1389 1390 1391 1392 1393 1394 1395 1396 1397 1398 1399 1400 1401 1402 |
# File 'ext/Dobjects/Function/function.c', line 1381
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 ⇒ Object
Sorts the X values while keeping the matching Y values.
772 773 774 775 |
# File 'ext/Dobjects/Function/function.c', line 772
static VALUE function_sort(VALUE self)
{
return function_joint_sort(self,get_x_vector(self), get_y_vector(self));
}
|
#sorted? ⇒ Boolean Also known as: is_sorted
Checks if the X values of the Function are sorted.
201 202 203 204 205 206 207 |
# File 'ext/Dobjects/Function/function.c', line 201
static VALUE function_is_sorted(VALUE self)
{
if(dvector_is_sorted(get_x_vector(self)))
return Qtrue;
else
return Qfalse;
}
|
#spline_approximation(params) ⇒ Object
Filters the Function through interpolation. params holds a hash with the following values: ??
It returns a hash.
1598 1599 1600 1601 1602 1603 1604 1605 1606 1607 1608 1609 1610 1611 1612 1613 1614 1615 1616 1617 1618 1619 1620 1621 1622 1623 1624 1625 1626 1627 1628 1629 1630 1631 |
# File 'ext/Dobjects/Function/function.c', line 1598
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 ⇒ Object
Splits the function into strictly monotonic sub-functions. Returns the array of the subfunctions. The returned values are necessarily new values.
653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 |
# File 'ext/Dobjects/Function/function.c', line 653
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(sym) ⇒ Object
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.
725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 |
# File 'ext/Dobjects/Function/function.c', line 725
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 ⇒ Object
Strips all the points containing NaN values from the function, and returns the number of points stripped.
622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 |
# File 'ext/Dobjects/Function/function.c', line 622
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);
}
|