summaryrefslogtreecommitdiff
path: root/metapost
diff options
context:
space:
mode:
Diffstat (limited to 'metapost')
-rw-r--r--metapost/context/base/mp-grap.mpiv173
1 files changed, 172 insertions, 1 deletions
diff --git a/metapost/context/base/mp-grap.mpiv b/metapost/context/base/mp-grap.mpiv
index 98f537315..bc02e8610 100644
--- a/metapost/context/base/mp-grap.mpiv
+++ b/metapost/context/base/mp-grap.mpiv
@@ -11,6 +11,8 @@
%C therefore copyrighted by \PRAGMA. See licen-en.pdf for
%C details.
+% maybe we should have another Gerr ... something grph_error_message
+
if known context_grap : endinput ; fi ;
boolean context_grap ; context_grap := true ;
@@ -130,7 +132,7 @@ for i = 3 upto 9 : % polygons
grap_sym[i] :=
for j = 0 upto i-1 :
(up scaled .5) rotated (360j/i) --
- endfor cycle ;
+ endfor cycle ;
endfor
grap_sym[12] := grap_sym[2] rotated +90 ; % horizontal line
@@ -203,3 +205,172 @@ def plotsymbol(expr n,c,f) = % (number,color,color|number)
nullpicture
fi
enddef ;
+
+% Here starts a section with some extensions that come in handy when drawing
+% polynomials. We assume that metapost is run in double number mode.
+
+% Least-squares "fit" to a polynomial
+%
+% Example of use:
+%
+% path p[] ;
+% numeric a[] ; a0 := 1 ; a1 := .1 ; a2 := .01 ; a3 := .001 ; a4 := 0.0001 ;
+% for i=0 upto 10:
+% x1 := 5i/10 ;
+% y1 := poly.a(4,x1) ;
+% augment.p0(z1) ;
+% augment.p1((x1,y1+.005normaldeviate)) ;
+% endfor
+% gdraw p0 ;
+% gdraw p1 plot plotsymbol(1,black,.5) ;
+%
+% numeric chisq, b[] ;
+% polyfit.p1(chisq, b, 4) ;
+% for i=0 upto length p1 :
+% x1 := xpart(point i of p1) ;
+% y1 := poly.b(4,x1) ;
+% augment.p2(z1) ;
+% endfor
+% gdraw p2 ;
+%
+% numeric c[] ;
+% linefit.p1(chisq, c) ;
+% for i=0 upto length p1 :
+% x1 := xpart(point i of p1) ;
+% y1 := line.c(x1) ;
+% augment.p3(z1) ;
+% endfor
+% gdraw p3 dashed evenly ;
+
+vardef det@# (expr n) = % find the determinant of a (n+1)*(n+1) matrix
+ % indices run from 0 to n.
+ % first, we make a copy so as not to corrupt the matrix.
+ save copy ; numeric copy[][] ;
+ for k=0 upto n :
+ for j=0 upto n :
+ copy[k][j] := @#[k][j] ;
+ endfor
+ endfor
+ numeric determinant, jj ; determinant := 1 ;
+ boolean zero ; zero := false ;
+ for k=0 upto n :
+ if copy[k][k] = 0 :
+ for 0 = k upto n :
+ if copy[k][j]=0 :
+ zero := true ;
+ else :
+ jj := j ;
+ fi
+ exitunless zero ;
+ endfor
+ if zero :
+ determinant := 0 ;
+ fi
+ exitif zero ;
+ for j = k upto n : % interchange the columns
+ temp := copy[j][jj] ;
+ copy[j][jj] := copy[j][k] ;
+ copy[j][k] := temp ;
+ endfor
+ determinant = -determinant ;
+ fi
+ exitif zero ;
+ determinant := determinant * copy[k][k] ;
+ if k < n : % subtract row k from lower rows to get a diagonal matrix
+ for j = k + 1 upto n :
+ for i = k + 1 upto n :
+ copy[j][i] := copy[j][i] - copy[j][k] * copy[k][i] / copy[k][k] ;
+ endfor
+ endfor
+ fi
+ endfor ;
+ determinant % no ;
+enddef ;
+
+% least-squares fit of a polynomial $ of order n to a path @#
+
+vardef polyfit@# (suffix chisq, $) (expr n) =
+ if not path begingroup @# endgroup :
+ Gerr(begingroup @# endgroup, "Cannot fit--not a path") ;
+ elseif length @# < n :
+ Gerr(begingroup @# endgroup, "Cannot fit--not enough points") ;
+ else:
+ chisq := 0 ;
+ % calculate sums of the data
+ save sumx, sumy ; numeric sumx[], sumy[] ;
+ save nmax ; numeric nmax ; nmax := 2*n ;
+ for i = 0 upto nmax :
+ sumx[i] := 0 ;
+ endfor
+ for i = 0 upto n :
+ sumy[i] := 0 ;
+ endfor
+ save xp, yp ; numeric xp, yp ;
+ save zi ; pair zi ;
+ for i = 0 upto length @# :
+ zi := point i of @# ;
+ x0 := xpart zi ; y0 := ypart zi ;
+ x1 := 1 ;
+ for j = 0 upto nmax :
+ sumx[j] := sumx[j] + x1 ;
+ x1 := x1 * x0 ;
+ endfor
+ y1 := y0 ;
+ for j = 0 upto n :
+ sumy[j] := sumy[j] + y1 ;
+ y1 := y1 * x0 ;
+ endfor
+ chisq := chisq + y0*y0 ;
+ endfor
+ % construct matrices and calculate the polynomial coefficients
+ save m ; numeric m[][] ;
+ for j = 0 upto n :
+ for k = 0 upto n :
+ i := j + k ;
+ m[j][k] := sumx[i] ;
+ endfor
+ endfor
+ save delta ; numeric delta ;
+ delta := det.m(n) ;
+ if delta = 0 :
+ chisq := 0 ;
+ for j=0 upto n :
+ $[j] := 0 ;
+ endfor
+ else :
+ for l = 0 upto n :
+ for j = 0 upto n :
+ for k = 0 upto n :
+ i := j + k ;
+ m[j][k] := sumx[i] ;
+ endfor
+ m[j][l] := sumy[j] ;
+ endfor
+ $[l] := det.m(n) / delta ;
+ endfor
+ for j = 0 upto n :
+ chisq := chisq - 2*$[j]*sumy[j] ;
+ for k = 0 upto n :
+ i := j + k ;
+ chisq := chisq + $[j]*$[k]*sumx[i] ;
+ endfor
+ endfor
+ % normalize by the number of degrees of freedom
+ chisq := chisq / (length(@#) - n) ;
+ fi
+ fi
+enddef ;
+
+vardef poly@#(expr n, x) =
+ for j = 0 upto n :
+ + @#[j]*(x**j)
+ endfor % no ;
+enddef ;
+
+vardef line@#(expr x) =
+ poly@# (1,x)
+enddef ;
+
+vardef linefit@#(suffix chisq, $) =
+ polyfit@#(chisq, $, 1) ;
+enddef ;