# Monthly Archives: May 2009

In an engineering project, I have been assigned to a task constructing smooth curve for given sets of sample points. This is accomplished by implementing Cubic Splines and Thomas algorithm to solve linear algebraic equation. So, I write PLSQL stored procedure as follows:

CREATE OR REPLACE FUNCTION spline(
IN x_arr numeric[], IN y_arr numeric[], IN resolution numeric,
OUT xs numeric, OUT ys numeric)
RETURNS SETOF record AS
\$BODY\$
DECLARE
x numeric[];
y numeric[];
e numeric[] DEFAULT ‘{}’::numeric[];
f numeric[] DEFAULT ‘{}’::numeric[];
g numeric[] DEFAULT ‘{}’::numeric[];
r numeric[] DEFAULT ‘{}’::numeric[];
f2 numeric[] DEFAULT ‘{}’::numeric[];
count integer;
i integer;
j integer;
idx integer;
tmp numeric;
BEGIN
x:=x_arr;
y:=y_arr;
count:=ARRAY_UPPER(x,1);
IF count<3 THEN
RAISE EXCEPTION ‘sample size must be >=3’;
END IF;
IF ARRAY_UPPER(y,1)<>count THEN
RAISE EXCEPTION ‘x and y array size must be identical’;
END IF;

——–sort x, ascendingly:
FOR i IN 1..count-1 LOOP
idx:=i+1;
FOR j IN i+2..count LOOP
IF x[j]<x[idx] THEN
idx:=j;
END IF;
END LOOP;
IF x[idx]=x[i] THEN
RAISE EXCEPTION ‘x values must be unique’;
ELSIF x[idx]<x[i] THEN –swap
tmp=x[i];
x[i]=x[idx];
x[idx]=tmp;
tmp=y[i];
y[i]=y[idx];
y[idx]=tmp;
END IF;
END LOOP;

——–prepare variables for Thomas algorithm
FOR i IN 1..count LOOP
IF i=1 OR i=count THEN
e:= e || 0.0;
f:= f || 0.0;
g:= g || 0.0;
r:= r || 0.0;
f2:=f2 || 0.0;
CONTINUE;
END IF;
e:= e || x[i] – x[i-1];
f:= f || 2*(x[i+1]-x[i-1]);
g:= g || x[i+1] – x[i];
r:= r || ( (6/(x[i+1]-x[i])) * (y[i+1]-y[i]) )
+ ( (6/(x[i]-x[i-1])) * (y[i-1]-y[i]) );
f2:=f2 || 0.0;
END LOOP;

——–decomposition
FOR i IN 2..count LOOP
IF f[i-1]<>0 THEN
e[i]:=e[i]/f[i-1];
END IF;
f[i]:=f[i]-(e[i]*g[i-1]);
END LOOP;

——–forward substitution
FOR i IN 2..count LOOP
r[i]:=r[i]-(e[i]*r[i-1]);
END LOOP;

——–backward substitution
FOR i IN REVERSE count-1..1 LOOP
IF f[i]<>0 THEN
f2[i]:=(r[i] – (g[i]*f2[i+1])) / f[i];
END IF;
END LOOP;

——–do interpolation
FOR i IN 2..count LOOP
xs:=x[i-1]-resolution;
LOOP
xs:=xs+resolution;
EXIT WHEN i=count AND xs>x[i];
EXIT WHEN i=x[i];
ys:= ( (f2[i-1] * (x[i]-xs)^3)
+ (f2[i] * (xs-x[i-1])^3)
) / (6*(x[i]-x[i-1]));
ys:= ys + (
(y[i-1]/(x[i]-x[i-1])) –
(f2[i-1]*(x[i]-x[i-1])/6)
) * (x[i]-xs);
ys:= ys + (
(y[i]/(x[i]-x[i-1])) –
(f2[i]*(x[i]-x[i-1])/6)
) * (xs-x[i-1]);
RETURN NEXT;
EXIT WHEN xs>=x[i];
END LOOP;
END LOOP;

RETURN;
END;
\$BODY\$
LANGUAGE ‘plpgsql’ IMMUTABLE STRICT;

The stored procedure takes 3 input arguments and 2 outs. x_arr and y_arr are sample points’ x and y. With resolution, you can determine the resolution of smoothed curve’s x.

Output arguments: xs and ys contain the set of x and y of smoothed curve, respectively.

Let us use the function. At first, we have sample points:

# x y
1 3.3 1.2
2 4.2 -1.1
3 6.8 1.5
4 8.2 1.2
5 10.2 2.5

Call the function, supply x and y array as displayed on above table, set resolution to 0.1:

SELECT * FROM spline(
‘{3.3,4.2,6.8,8.2,10.2}’::numeric[],
‘{1.2,-1.1,1.5,1.2,2.5}’::numeric[],
0.1);

We will have set of records containing xs and ys.:

xs  ys
3.3  1.199999999999999970
3.4  0.8860175045157231811
3.5  0.5764170295261004851
3.6  0.2755805955257859750
3.7  -0.0121097769905662565
3.8  -0.28227206752830211620
… and so on

Let us put them in graph: The blue-y1 and red-y2 are sample and interpolation points as returned by the spline function, respectively.