beta 2013.04.05 22:00

2 files changed, 296 insertions, 132 deletions

diff --git a/metapost/context/base/mp-grap.mpiv b/metapost/context/base/mp-grap.mpiv
index bc02e8610..cce63abe0 100644
--- a/metapost/context/base/mp-grap.mpiv
+++ b/metapost/context/base/mp-grap.mpiv
@@ -11,12 +11,13 @@
 %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 ;
+if known context_grap : endinput ; fi
 
 boolean context_grap ; context_grap := true ;
 
+% Instead we could include graph here and then clean it up as well as use private
+% variables in the grap_ namespace. After all, graph is frozen.
+
 input graph.mp ;
 
 vardef roundd(expr x, d) =
@@ -66,9 +67,23 @@ vardef dofmt_@#(expr f, x) =
   % textext(mfun_format_number(escaped_format(f),x))
 enddef ;
 
+% note that suffix @# is ignored above...
+
+vardef strfmt(expr f, x) =
+    "\MPgraphformat{" & escaped_format(f) & "}{" & (if string x : x else: decimal x fi) & "}"
+enddef ;
+
 % We redefine autogrid from graph.mp adding the possibility of differing X and Y
 % formats. Autoform is defined in graph.mp (by default "%g").
-%
+
+% graph.mp: string Autoform; Autoform = "%g";
+% graph.mp:
+% graph.mp: vardef autogrid(suffix tx, ty) text w =
+% graph.mp:   Gneedgr_:=false;
+% graph.mp:   if str tx<>"": for x=auto.x: tx(Autoform,x) w; endfor fi
+% graph.mp:   if str ty<>"": for y=auto.y: ty(Autoform,y) w; endfor fi
+% graph.mp: enddef;
+
 % string Autoform_X ; Autoform_X := "@.0e" ;
 % string Autoform_Y ; Autoform_Y := "@.0e" ;
 
@@ -86,9 +101,10 @@ vardef autogrid(suffix tx, ty) text w =
                 else :
                     Autoform
                 fi,
-                x) w ;
+                x
+            ) w ;
         endfor
-    fi ;
+    fi
     if str ty <> "" :
         for y=auto.y :
             ty (
@@ -101,16 +117,17 @@ vardef autogrid(suffix tx, ty) text w =
                 else :
                     Autoform
                 fi,
-                y) w ;
+                y
+            ) w ;
         endfor
-    fi ;
+    fi
 enddef ;
 
 % A couple of extensions:
 
-% Define a vector function sym returning a picture: 10 different shapes,
+% Define a function plotsymbol() returning a picture: 10 different shapes,
 % unfilled outline, interior filled with different shades of the background.
-% Thus, overlapping points on a plot are more clearly distinguishable.
+% This allows overlapping points on a plot to be more distinguishable.
 
 % grap_symsize := fontsize defaultfont ; % can be redefined
 %
@@ -124,11 +141,11 @@ enddef ;
 
 path grap_sym[] ; % (internal) symbol path
 
-grap_sym[0] := (0,0) ;                  % point
-grap_sym[1] := fullcircle ;             % circle
-grap_sym[2] := (up -- down) scaled .5 ; % vertical bar
+grap_sym[0] := (0,0) ;                      % point
+grap_sym[1] := fullcircle ;                 % circle
+grap_sym[2] := (up -- down) scaled .5 ;     % vertical bar
 
-for i = 3 upto 9 :                  % polygons
+for i = 3 upto 9 :                          % polygons
     grap_sym[i] :=
         for j = 0 upto i-1 :
             (up scaled .5) rotated (360j/i) --
@@ -179,9 +196,9 @@ endfor cycle ;
 
 grap_sym[34] := grap_sym[24] rotated 45 ;
 
-    % usage: gdraw p plot plotsymbol(1,red,1) ;       % a filled red circle
-    % usage: gdraw p plot plotsymbol(4,blue,0) ;      % a blue square
-    % usage: gdraw p plot plotsymbol(14,green,0.5) ;  % a 50% filled green diamond
+% usage: gdraw p plot plotsymbol(1,red,1) ;       % a filled red circle
+% usage: gdraw p plot plotsymbol(4,blue,0) ;      % a blue square
+% usage: gdraw p plot plotsymbol(14,green,0.5) ;  % a 50% filled green diamond
 
 def plotsymbol(expr n,c,f) = % (number,color,color|number)
     if known grap_sym[n] :
@@ -206,171 +223,303 @@ def plotsymbol(expr n,c,f) = % (number,color,color|number)
     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.
+% convert a polygon path to a smooth path (useful, e.g. as a guide to the eye)
+
+def smoothpath (suffix $) =
+    if path $ :
+    (for i=0 upto length $ :
+        if i>0 : .. fi
+        (point i of $)
+     endfor )
+    fi
+enddef ;
+
+% return a path of a function func(x) with abcissa running from f to t over n intervals
 
-% Least-squares "fit" to a polynomial
+def makefunctionpath (expr f, t, n) (text func) =
+    (for x=f step ((t-f)/n) until t :
+        if x<>f : .. fi
+        (x, func)
+    endfor )
+enddef ;
+
+% shift a path, point by point
 %
-% Example of use:
+% example:
+%
+%   p1 := addnoisetopath(p0,(.1normaldeviate,.1normaldeviate)) ;
+
+vardef addnoisetopath (suffix p) (text t) =
+    if path p :
+    hide(pair pi)
+    (for i=0 upto length p :
+        if i>0 : -- fi
+        hide(pi := point i of p; x := xpart pi; y := ypart pi)z shifted t
+    endfor)
+    fi
+enddef ;
+
+% return a new path of a function func(x) using the same abcissa as an existing path
+
+vardef functionpath (suffix p) (text t) =
+    (for i=0 upto length p :
+        if i>0 : .. fi
+        (hide(x := xpart(point i of p))x,t)
+     endfor )
+enddef ;
+
+% least-squares "fit" to a polynomial
+%
+% example:
 %
 %   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
+%   p0 := makefunctionpath(0,5,10,poly(a,4,x)) ;
+%   p1 := addnoisetopath(p0,(0,.001normaldeviate)) ;
 %   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 b[] ;
+%   polyfit(p1, b, 4, 1) ;
+%   gdraw functionpath(p1,poly(b,4,x)) ;
 %
 %   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 ;
+%   linefit(p1, c, 1) ;
+%   gdraw functionpath(p1,line(c,x)) dashed evenly ;
+
+% a polynomial function:
+%
+% y = a0 + a1 * x + a2 * x^2 + ... + a[n] * x^n
+
+vardef poly (suffix $) (expr n, x) =
+    (for j=0 upto n : + $[j]*(x**j) endfor) % no ;
+enddef ;
+
+% find the determinant of a (n+1)*(n+1) matrix; indices run from 0 to n
+
+vardef det (suffix $) (expr n) =
+    hide(
+        numeric determinant ; determinant := 1 ;
+        save jj ; numeric jj ;
+        for k=0 upto n :
+            if $[k][k]=0 :
+                jj := -1 ;
+                for j=0 upto n :
+                    if $[k][j]<>0 :
+                        jj := j ;
+                        exitif true ;
+                    fi
+                endfor
+                if jj<0 :
+                    determinant := 0 ;
+                    exitif true ;
                 fi
-                exitunless zero ;
-            endfor
-            if zero :
-                determinant := 0 ;
+                for j=k upto n : % interchange the columns
+                    temp := $[j][jj] ;
+                    $[j][jj] := $[j][k] ;
+                    $[j][k] := temp ;
+                endfor
+                determinant = -determinant ;
             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] ;
+            exitif determinant=0 ;
+            determinant := determinant * $[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:
+                        $[j][i] := $[j][i]-$[j][k]*$[k][i]/$[k][k] ;
+                    endfor
                 endfor
-            endfor
-        fi
-    endfor ;
+            fi
+        endfor ;
+    )
     determinant % no ;
 enddef ;
 
-% least-squares fit of a polynomial $ of order n to a path @#
+numeric Gchisq_ ; % global - use graph.mp naming convention (until we clean things up).
+% why not grap_pf_temp? answer: because we want this to be available to the user.
 
-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 ;
+% least-squares fit of a polynomial $ of order n to a path p (unweighted)
+%
+% reference: P. R. Bevington, "Data Reduction and Error Analysis for the Physical
+% Sciences", McGraw-Hill, New York 1969.
+
+vardef polyfit (suffix p, $) (expr n) (text t) =
+    if not path p :
+        Gerr(p, "Cannot fit--not a path") ;
+    elseif length p < n :
+        Gerr(p, "Cannot fit--not enough points") ;
+    else :
+        Gchisq_ := 0 ;
         % calculate sums of the data
         save sumx, sumy ; numeric sumx[], sumy[] ;
-        save nmax ; numeric nmax ; nmax := 2*n ;
-        for i = 0 upto nmax :
+        save w ; numeric w ;
+        for i=0 upto 2n :
             sumx[i] := 0 ;
         endfor
-        for i = 0 upto n :
+        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 :
+        for i=0 upto length p :
+            clearxy; z = point i of p ;
+            w := if length(t) > 0 : t else : 1 fi ; % weight
+            x1 := w ;
+            for j=0 upto 2n :
                 sumx[j] := sumx[j] + x1 ;
-                x1 := x1 * x0 ;
+                x1 := x1 * x ;
             endfor
-            y1 := y0 ;
-            for j = 0 upto n :
+            y1 := y * w ;
+            for j=0 upto n :
                sumy[j] := sumy[j] + y1 ;
-               y1 := y1 * x0 ;
+               y1 := y1 * x ;
             endfor
-            chisq := chisq + y0*y0 ;
+            Gchisq_ := Gchisq_ + y*y*w ;
         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] ;
+        for j=0 upto n :
+            for k=0 upto n :
+                m[j][k] := sumx[j+k] ;
             endfor
         endfor
         save delta ; numeric delta ;
-        delta := det.m(n) ;
+        delta := det(m,n) ; % this destroys the matrix m[][], which is OK
         if delta = 0 :
-            chisq := 0 ;
+            Gchisq_ := 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] ;
+            for i=0 upto n :
+                for j=0 upto n :
+                    for k=0 upto n :
+                        m[j][k] := sumx[j+k] ;
                     endfor
-                    m[j][l] := sumy[j] ;
+                    m[j][i] := sumy[j] ;
                 endfor
-                $[l] := det.m(n) / delta ;
+                $[i] := det(m,n) / delta ; % matrix m[][] gets destroyed...
             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] ;
+            for j=0 upto n :
+                Gchisq_ := Gchisq_ - 2sumy[j]*$[j] ;
+                for k=0 upto n :
+                    Gchisq_ := Gchisq_ + $[j]*$[k]*sumx[j+k] ;
                 endfor
             endfor
             % normalize by the number of degrees of freedom
-            chisq := chisq / (length(@#) - n) ;
+            Gchisq_ := Gchisq_ / (length(p) - n) ;
+        fi
+    fi
+enddef ;
+
+% y = a0 + a1 * x
+%
+% of course a line is just a polynomial of order 1
+
+vardef line    (suffix $)    (expr x) = poly($,1,x) enddef ; % beware: potential name clash
+vardef linefit (suffix p, $) (text t) = polyfit(p, $, 1, t) ; enddef ;
+
+% and a constant is polynomial of order 0
+
+vardef const    (suffix $)    (expr x) = poly($,0,x) enddef ;
+vardef constfit (suffix p, $) (text t) = polyfit(p, $, 0, t) ; enddef ;
+
+% y = a1 * exp(a0*x)
+%
+% exp and ln defined in metafun
+
+vardef exponential (suffix $) (expr x) = $1*exp($0*x) enddef ;
+
+% since we take a log, this only works for positive ordinates
+
+vardef exponentialfit (suffix p, $) (text t) =
+    save a ; numeric a[] ;
+    save q ; path q ; % fit to the log of the ordinate
+    for i=0 upto length p :
+        if ypart(point i of p)>0 :
+            augment.q(xpart(point i of p),ln(ypart(point i of p))) ;
+        fi
+    endfor
+    linefit(q,a,t) ;
+    $0 := a1 ;
+    $1 := exp(a0) ;
+enddef ;
+
+% y = a1 * x**a0
+
+vardef power (suffix $) (expr x) = $1*(x**$0) enddef ;
+
+% since we take logs, this only works for positive abcissae and ordinates
+
+vardef powerfit (suffix p, $) (text t) =
+    save a ; numeric a[] ;
+    save q ; path q ; % fit to the logs of the abcissae and ordinates
+    for i=0 upto length p :
+        if (xpart(point i of p)>0) and (ypart(point i of p)>0) :
+            augment.q(ln(xpart(point i of p)),ln(ypart(point i of p))) ;
         fi
+    endfor
+    linefit(q,a,t) ;
+    $0 := a1 ;
+    $1 := exp(a0) ;
+enddef ;
+
+% gaussian: y = a2 * exp(-ln(2)*((x-a0)/a1)^2)
+%
+% a1 is the hwhm; sigma := a1/sqrt(2ln(2)) or a1/1.17741
+
+numeric lntwo ; lntwo := ln(2) ; % brrr, why not inline it
+
+vardef gaussian (suffix $) (expr x) =
+    if $1 = 0 :
+        if x = $0 : $2 else : 0 fi
+    else :
+        $2 * exp(-lntwo*(((x-$0)/$1)**2))
+    fi
+    if known $3 :
+        + $3
     fi
 enddef ;
 
-vardef poly@#(expr n, x) =
-    for j = 0 upto n :
-        + @#[j]*(x**j)
-    endfor % no ;
+% since we take a log, this only works for positive ordinates
+
+vardef gaussianfit (suffix p, $) (text t) =
+    save a ; numeric a[] ;
+    save q ; path q ; % fit to the log of the ordinate
+    for i=0 upto length p :
+        if ypart(point i of p)>0 :
+            augment.q(xpart(point i of p), ln(ypart(point i of p))) ;
+        fi
+    endfor
+    polyfit(q,a,2,if t > 0 : ln(t) else : 0 fi) ;
+    $1 := sqrt(-lntwo/a2) ;
+    $0 := -.5a1/a2 ;
+    $2 := exp(a0-.25*a1*a1/a2) ;
+    $3 := 0 ; % polyfit will NOT work with a non-zero background!
 enddef ;
 
-vardef line@#(expr x) =
-    poly@# (1,x)
+% lorentzian: y = a2 / (1 + ((x - a0)/a1)^2)
+
+vardef lorentzian (suffix $) (expr x) =
+    if $1 = 0 :
+        if x = $0 : $2 else : 0 fi
+    else :
+        $2 / (1 + ((x - $0)/$1)**2)
+    fi
+    if known $3 :
+        + $3
+    fi
 enddef ;
 
-vardef linefit@#(suffix chisq, $) =
-    polyfit@#(chisq, $, 1) ;
+vardef lorentzianfit (suffix p, $) (text t) =
+    save a ; numeric a[] ;
+    save q ; path q ; % fit to the inverse of the ordinate
+    for i=0 upto length p :
+        if ypart(point i of p)<>0 :
+            augment.q(xpart(point i of p), 1/ypart(point i of p)) ;
+        fi
+    endfor
+    polyfit(q,a,2,if t <> 0 : 1/(t) else : 0 fi) ;
+    $0 := -.5a1/a2 ;
+    $2 := 1/(a0-.25a1*a1/a2) ;
+    $1 := sqrt((a0-.25a1*a1/a2)/a2) ;
+    $3 := 0 ; % polyfit will NOT work with a non-zero background!
 enddef ;
diff --git a/metapost/context/base/mp-tool.mpiv b/metapost/context/base/mp-tool.mpiv
index a947143c4..5b53dcdef 100644
--- a/metapost/context/base/mp-tool.mpiv
+++ b/metapost/context/base/mp-tool.mpiv
@@ -2307,3 +2307,18 @@ enddef ;
 %         setbounds currentpicture to mfun_bleed_box ;
 %     )
 % enddef ;
+
+%D Dimensions have bever been an issue as traditional MP can't make that large
+%D pictures, but with double mode we need a catch:
+
+newinternal maxdimensions ; maxdimensions := 14000 ;
+
+def mfun_apply_max_dimensions = % not a generic helper, we want to protect this one
+    if bbwidth currentpicture > maxdimensions :
+        currentpicture := currentpicture if bbheight currentpicture  > bbwidth currentpicture : ysized else : xsized fi maxdimensions ;
+    elseif bbheight currentpicture  > maxdimensions :
+        currentpicture := currentpicture ysized maxdimensions ;
+    fi ;
+enddef;
+
+extra_endfig := extra_endfig & "mfun_apply_max_dimensions ;" ;