diff options
Diffstat (limited to 'metapost')
| -rw-r--r-- | metapost/context/base/mp-grap.mpiv | 173 |
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 ; |