line source
+ − (*
+ −
+ − /home/ghh/math/polyreg/polyregz.m 17.01.2008 Gerhard Heinzel AEI
+ −
+ − This is version 1.1
+ −
+ − history:
+ − version 1.0 17.1.08 initial release
+ − version 1.1 18.1.08 fixed two typos in comment
+ −
+ − Mathematica code to produce C code for efficient polynomial detrending
+ −
+ − The polynomials are defined in terms of z = 2i/(n-1)-1, i=0...n
+ − such that the range 0...n-1 is mapped to -1...+1
+ −
+ − Thereby the high dynamic range of i^n is reduced to -1..1,
+ − and furthermore every second entry in the matrix of
+ − the normal equations and its inverse vanishes due to symmetry.
+ −
+ − The regression is done via explicit solution of the normal equations.
+ −
+ − These are of dimension (p+1) for degree p:
+ −
+ − ( s_0 s_1 s_2 .. s_p ) ( a_0 ) = ( b_0 )
+ − ( s_1 s_2 s_3 .. s_p+1 ) ( a_1 ) ( b_1 )
+ − ( ... ) ( ... ) ( ... )
+ − ( s_p s_p+1 .. s_2p ) ( a_p ) ( b_p )
+ −
+ − with
+ − s_k = Sum(i=0...n-1) z^k,
+ − b_k = Sum(i=0...n-1) x_i z^k,
+ − a_k = The unknown coefficient of z^k in the best-fit regression.
+ −
+ − The Matrix with s depends only on p and n, not on the data.
+ −
+ − The Mathematica code below finds exact explicit expressions for
+ − the Matrix and its Inverse and writes the corresponding C code.
+ −
+ − If any modifications of the C code are desired apart from "indent"ing,
+ − this file below should be altered but not the C code.
+ −
+ − *)
+ −
+ − <<clear.m
+ −
+ − (*
+ − <<polyregz.m
+ − *)
+ −
+ − z[i_,n_] := 2 i /(n-1) -1;
+ −
+ − s[k_,n_] := Sum[z[i,n]^k,{i,0,n-1}];
+ −
+ − m[n_,p_] := Table[s[i+j,n],{i,0,p},{j,0,p}];
+ −
+ − sol[p_] := sol[p] = Simplify[Inverse[m[n, p]]];
+ −
+ − cln[x_] := CoefficientList[Numerator[x],n];
+ − cld[x_] := CoefficientList[Denominator[x],n];
+ −
+ − max[x_] := Max[cln[x],cld[x]];
+ −
+ − pn[x_] := "0" /; x== 0;
+ − pn[x_] := ToString[CForm[x]] <> "L";
+ −
+ − ncln[x_] := Map[pn, N[cln[x] / max[x], 22]];
+ − ncld[x_] := Map[pn, N[cld[x] / max[x], 22]];
+ −
+ − polyprint[p_, x_] := p[[1]] /; Length[p]==1;
+ − polyprint[p_, x_] := "(" <> polyprint[Drop[p,1], x] <> " * " <> x <> " + " <> p[[1]] <> ")" /; Length[p]>1 ;
+ −
+ − cform[x_] := "0" /; FreeQ[x,n];
+ − cform[x_] := polyprint[ncln[x], "n"] <> " / " <> polyprint[ncld[x], "n"];
+ −
+ − wlist[sbeg_, list_, sep_, send_] :=
+ − Write[fp, sbeg <> Apply[StringJoin, Riffle[list, sep]] <> send];
+ −
+ − out[p_] := Module[ {fn, s, i, j, s1, list},
+ − fn="polyreg"<>ToString[p] <> ".c";
+ − fp=OpenWrite[fn,FormatType->OutputForm, PageWidth->1024];
+ − Write[fp,"void"];
+ − Write[fp,"polyreg",p," (double *x, int nn, double *y, double *a)"];
+ − Write[fp,"{"];
+ − Write[fp,""];
+ − Write[fp,"/*"];
+ − Write[fp,"polynominal detrending of a time series, order ",p];
+ − Write[fp,"machine-generated file, do not edit!"];
+ − Write[fp,"made by polyregz.m (Mathematica )"];
+ − Write[fp,"Gerhard Heinzel AEI 17.01.2008"];
+ − Write[fp,""];
+ − Write[fp,"x[]: input, read-only: time series to be detrended"];
+ − Write[fp,"nn: input, read-only: length of x[]"];
+ − Write[fp,"y[]: output: time series with trend subtracted"];
+ − Write[fp,"a[]: fitting coefficients in for z^0, z^1, z^2,... with"];
+ − Write[fp," z = 2*i/(nn-1)-1 ; i=0,...,nn-1"];
+ − Write[fp,"*/"];
+ − Write[fp,""];
+ − Write[fp," long double n = nn, n1 = 2.L / (n - 1), z;"];
+ − For[i=0, i<=p, i++,
+ − Write[fp," long double a",i,", temp",i,", sum",i," = 0;"];
+ − ];
+ − Write[fp," int i;"];
+ − Write[fp," for (i = 0; i < nn; i++)"];
+ − Write[fp," {"];
+ − Write[fp," z = n1 * i - 1.L;"];
+ − list={"x[i]"};
+ − For[i=0, i<=p, i++,
+ − s1=" sum" <> ToString[i] <> " += ";
+ − wlist[s1, list, " * ", ";"];
+ − list=Append[list,"z"];
+ − ];
+ − Write[fp,"/* the above code is efficiently optimized by GCC4 */"];
+ − Write[fp,""];
+ − Write[fp," }"];
+ −
+ − s=sol[p];
+ − For[i=0, i<=p, i++,
+ − list={};
+ − For[j=0, j<=p, j++,
+ − If[!FreeQ[s[[i+1,j+1]],n],
+ − Write[fp, "temp", j, " = ", cform[s[[i+1,j+1]]], ";"];
+ − list=Append[list,"temp" <> ToString[j] <> " * sum" <> ToString[j]];
+ − ]
+ − ];
+ − s1 = "a" <> ToString[i] <> " = ";
+ − wlist[s1, list, " + ", ";"];
+ − ];
+ − Write[fp," for (i = 0; i < nn; i++)"];
+ − Write[fp," {"];
+ − Write[fp," z = n1 * i - 1.L;"];
+ − list={};
+ − For[i=0, i<=p, i++,
+ − list=Append[list, "a" <> ToString[i]];
+ − ];
+ −
+ − Write[fp," y[i] = x[i] - ", polyprint[list, "z"], ";"];
+ − Write[fp," }"];
+ − For[i=0, i<=p, i++,
+ − Write[fp," a[",i,"] = a",i,";"];
+ − ];
+ − Write[fp,"}"];
+ − Close[fp];
+ − ];
+ −
+ − For[i=0, i<=10, i++,
+ − out[i];
+ − Print["degree ",i," done."];
+ − ]
+ −