Tuesday, December 22, 2015

Maxima in Common Lisp via ASDF

I have previously used Lisp code that ran in Maxima. This is fairly straightforward if you know Common Lisp. Basic instructions can be found in Chapter 37.1 of the Maxima help files (look for Contents link at the top of any of the help pages). But, how about the other way around—Maxima from Lisp side?

I recently looked into running Maxima inside Common Lisp using ASDF and found this. That gave me hope that it was possible and a little insight into what it was trying to do. I have sometimes wanted to create a user interface in some other environment (such as LispWorks) and call upon Maxima to do some computations for me.

To get yourself set up to access Maxima from Common Lisp you need:
  1. The source, including the maxima.asd file.
  2. Some idea where to put stuff and, perhaps, a few tweaks to the maxima.asd file.
  3. Some file search variable set up to allow Maxima load statements in your common lisp code so you can access functionality found in the share packages or your own packages.

Tweaks

I put the source under my common-lisp directory so that a simple call to (asdf:load-system :maxima) would work.

I found that I needed to edit the maxima.asd file. I don't want to tell ASDF where to put the result or interfere with its "normal" operation. "ASDF, you're in charge!" is basically how I look at it. Accordingly, I commented out the output-files :around method. I suppose you could change the *binary-output-dir* here as well, if you would like to—perhaps to suit whichever Lisp you are compiling to.

Fig. 1. I'm pretty sure I don't need this code for my purposes, so I have it commented out.

Accessing Shared Packages

You might not realize how much stuff is in the share packages if you have been using the wxMaxima interface only. It appears that a number of widely used share packages are automatically loaded in wxMaxima, but maxima.asd only loads the core components. This is not a serious impediment—you can determine which packages to load from the Maxima help files.

But before you can use the Maxima load command, you need to set the file_search_maxima and file_search_lisp Maxima variables. I combined a call to asdf:load-system with the setting of these variables in a file under my common-lisp directory. I called this file load-max.lisp and it looks a bit like this (note how hairy it gets if you scroll right):



So, if I want to use Maxima in my Lisp session, I just need to call (load "~/common-lisp/load-max.lisp"). I can move to the Maxima package by calling (in-package :maxima).

Try out the diff function:

 MAXIMA> ($diff #$x^2$ #$x$)
((MTIMES SIMP) 2 $X)

 Note the syntax for referring to the diff function ($diff) and the syntax for an entire maxima expression (#$expression$). There's also some lower to uppercase switching to watch out for, but see chapter 37 of the Maxima help documentation for more about this.

 Suppose I want to use something in a share package or one of my own custom Maxima packages. I can use the Maxima version of load, namely, $load.

 MAXIMA> ($load "stringproc")

MAXIMA> ($parse_string "x^2")
((MEXPT) $X 2)

One more example, to show how to deal with Maxima-specific features that are beyond the scope of Common Lisp. In general, use the meval function as in:

 (setf n (* 3 4 5 6 7))
(maxima::meval `((maxima::$divisors) ,n))

 ((MAXIMA::$SET MAXIMA::SIMP) 1 2 3 4 5 6 7 8 9 10 12 14 15 18 20 21 24 28 30 35 36 40 42 45 56 60 63 70 72 84 90 105 120 126 140 168 180 210 252 280 315 360 420 504 630 840 1260 2520)

 (If you're interested, see Maxima-discuss archives. HT Robert Dodier.)
    Posted by Darren Irvine at 1:26 PM 0 comments
    Labels: , ,

    Sunday, September 6, 2015

    Plotting a LWPOLYLINE in Maxima

    In this post, I want to show how to produce a representation of LWPOLYLINEs in Maxima that is suitable for plotting to the plot2d function. The direct usefulness of such an operation is dubious. Think of this as a worked example to help build understanding of a few mathematical and computing tools.

    Recap
    The LWPOLYLINE is a 2d representation of polylines that represents each vertex as having x-value, y-value, and bulge value. The bulge value indicates the path of travel from the current point to the next one. Either it will go straight (bulge = 0), veer to the left and back as a CW arc (bulge < 0), or it will veer to the right and back as a CCW arc (bulge > 0).

    So, it will
    • bulge left (< 0),
    • go straight (= 0), or
    • bulge right (> 0).
    More details on the mathematics of the bulge factor can be found here and here.

    The Pieces
    First of all, we define our problem and break its solution up into pieces. A polyline is of such a nature that is does not make a good "\(f(x)\)" type function (plotting \(y\) versus \(x\)). However, it lends itself to parametric expression well enough. So, we will have an \(x(s)\) function and a \(y(s)\) function. Ideally, I would like to specify a length along the polyline (measured from the start) and get a point \((x(s), y(s))\). (If you have a surveying background, think \(s = \) station number.)

    So, I need to have parametric equations for lines and for arcs that are parameterized for arc length and initial value of the parameter. I then need to stitch these together to make both an \(x\) and a \(y\) piecewise function.

    Linear

    Parametric equations for lines involve a direction and a starting point. Simply normalize the direction vector to parameterize for length. This looks like
    \[(x(s), y(s)) = A + \frac{(B-A)}{||B-A||}s,\]
    where \(A\) is the starting point and \(B\) is the final point. We will need to adjust the initial value to account for all of the length of the polyline prior to the initial point on this segment, say, \(s_i\) which changes our equations to
    \[(x(s), y(s)) = A + \frac{(B-A)}{||B-A||}(s-s_i).\]

    Arc

    Parametric equations for arcs use a center, radius, initial angle, direction of rotation, and total angle subtended (some variation of feasible input parameters admitted). Parameterizing for arc length isn't as obvious as it was for lines, so let's start with a basic circle:
    \[x(s) = c_x + r \cos{s}\]
    \[y(s) = c_y + r \sin{s},\]
    where \((c_x, c_y)\) is the center and \(r\) is the radius of the circle. For our reparameterization to work we require that radian measure be used (or else require degrees and do the conversion in the function, but why make work for ourselves and introduce inefficiency?). Recall that the definition of radians is the ratio of arc length to radius (\(\theta = s/r\)). So, our parameter becomes \(s/r\) and so we have

    \[x(s) = c_x + r \cos{\left(\frac{s}{r}\right)}\]
    \[y(s) = c_y + r \sin{\left(\frac{s}{r}\right)}.\]

    To get the right starting point for our function we need the parameter inside the trig functions to be equal to the correct initial angle when the input value of \(s_i\) (length of polyline so up to the start of the arc) is given. We need, then, the initial angle, say \(\theta_i\), plus the angle past initial that the given \(s\) corresponds to. Thus,

    \[x(s) = c_x + r \cos{\left(\theta_i + \sigma \frac{s-s_i}{r}\right)}\]
    \[y(s) = c_y + r \sin{\left(\theta_i + \sigma \frac{s-s_i}{r}\right)},\]
    where \(\sigma\) is 1 or -1 as per the sign of the bulge factor.

    Piecewise

    Making a piecewise function parameterized for length requires that we keep track of the length of the polyline up to the segment we are about to address. Implementing it in Maxima also requires a bit of acquaintance with the methods available to us for doing so. It is possible to use if-then-else logic to achieve this and obtain something that works fine for plotting purposes. As an example of this, consider a really simple example:

    f(x) := if x < 1 then x^2 else x^3;
    plot2d(f(x),[x,-1,2]);

    Fig 1. If-then-else logic in functions plots fine.
    There's another, somewhat handier approach to this problem, namely, the characteristic function. The characteristic function function is a way of expressing conditions as multiplication by 1 or 0 in order to include and exclude the appropriate expression. charfun2(x,a,b) returns 1 if \(x\in [a,b)\), else 0. So, we can modify the above function to be as follows and obtain the same result. 

    load(interpol);
    f(x) := charfun2(x,-inf,1)*x^2 + charfun2(x,1,inf)*x^3;
    plot2d(f(x),[x,-1,2]);

    One nice thing about these expressions is that it is possible to create them separately and add them together later and get a relatively normal Maxima expression. (There are some hiccups if you try to differentiate or integrate them.) I can't comment on the performance of the characteristic function approach versus the if-then-else approach, but the characteristic function approach may be more manageable to produce.

    The Program
    There are a number of details that don't specifically relate to this problem that I will not go into detail about and other details I have given in the two previously mentioned posts on bulge factors. Here is the complete code for producing the parametric equations and plotting them followed by the result for a specific example polyline.



    Posted by Darren Irvine at 1:26 PM 0 comments
    Labels: , , , ,

    Circle Fitting in Maxima

    In previous posts, I've considered the problem of fitting a circle to 2D data points. Fitting problems such as this require:
    1. An equation form to fit to. 
    2. A function to quantify the error, which is to be minimized.
    For a circle, the form of equation we want to fit to is
    \[r^2=(x-a)^2 + (y-b)^2.\]
    Maxima's lsquares package removes the burden of item number 2 for us. But for interest's sake, the error function we want to minimize can be given as
    \[\mathrm{r}\left( a,b\right) =\frac{\sum_{i=1}^{n}\sqrt{{\left( {x}_{i}-a\right) }^{2}+{\left( {y}_{i}-b\right) }^{2}}}{n}\]
    \[\mathrm{SSE}\left( a,b\right) =\sum_{i=1}^{n}{\left( \mathrm{r}\left( a,b\right) -\sqrt{{\left( {x}_{i}-a\right) }^{2}+{\left( {y}_{i}-b\right) }^{2}}\right) }^{2}.\]
    Recalling that \(\sum_{i=1}^{n}(\bar{x}-x_i)^2 = \sum_{i=1}^{n}{x_i}^2 - (\sum_{i=1}^{n}{x_i})/n\), we see that
    \[\mathrm{SSE}\left( a,b\right) =\sum_{i=1}^{n}\left({\left( {x}_{i}-a\right) }^{2}+{\left( {y}_{i}-b\right) }^{2}\right)-\frac{{\left(\sum_{i=1}^{n}\sqrt{{\left({x}_{i}-a\right) }^{2}+{\left( {y}_{i}-b\right) }^{2}}\right)}^{2}}{n}.\]

    lsquares_estimates()

    Given a set of points in matrix form (m), we need only to use the equation of a circle and indicate which are the data point variables and which are the "solve for" variables. In our problem,

    lsquares_estimates(m, [x,y], (x-a)^2 + (y-b)^2 = r^2, [a,b,r]);

    will suffice. The following function can be used to check the SSE.

    sse(a, b, r, pts) := sum((r-sqrt((pts[i][1]-a)^2+(pts[i][2]-b)^2))^2,i,1,length(pts));

    Testing for Reasonableness

    It is reasonable to ask whether least squares (or other data fitting methods) will reliably produce good results with a particular quality of data for a particular application. This is true for any problem. What if the math works fine, but the likelihood of being able to "math filter" the  likely errors in the data to reach an accurate answer is not so high? In other words, how much inaccuracy in the collected data can you tolerate and still have a reasonable expectation of getting to the right answer? I'm not going to attempt to deal with probabilities (per se) or confidence intervals, but look at a simulation approach that may be able to give you an idea of the reasonableness of the circle fit based on an (intuitively) estimated amount of random error in your data.

    Simulation of Circular Arc Data

    We will consider the problem of data collected, using a total station, of a circular arc, with an "as-built" scenario particularly in view. As such, auto-generated sample points should simulate the expected behaviors of a rod/prism person stepping out regular sample locations from start to end of an ostensibly circular arc. This can be implemented as random normal perturbations of a circular arc.

    One of the first pieces to a realistic perturbed circle is to be able to apply randomization that follows a normal distribution given an expected value and a standard deviation. Maxima has a function that behaves this way. random_normal() can be called with an expected value and standard deviation and optionally with a number of values to generate. Here is a sample use and output of the function, pictured using the histogram function:

    load(distrib)$
    load(draw)$
    histogram(
        random_normal(50.0, 1.0, 1000),
        nclasses = 10,
        xrange = [40,60]);

    Fig. 1. Normally distributed data
    The standard deviation should be chosen such that roughly two-thirds of the time (≈68%), the sample value will be within one standard deviation of the mean; that is, on the interval \((\bar{x}-\sigma, \bar{x}+\sigma).\) 

    Where the radius is concerned, we have the relatively predictable variations of the prism person's actions and the relatively unpredictable variations of the ostensible circular arc, which could turn out to not be an arc at all. In the absence of any sort of meaningful solution to this latter variation, we could possibly just bump up the standard deviation. But I'm not clear on whether this is a meaningful thing to do in assessing the acceptability of the data collection precision. The error should be "small" if the thing we measured conforms well to a circle. In order to check for whether this is "small" you need to produce a standard deviation from the SSE: \[\sigma_{SSE} = \sqrt{\frac{SSE}{n-1}}\] In the scenario we're dealing with here, someone would probably draw the computed arc in a CAD program and compare it to the measured data points and in an "arm-waving" sort of way decide whether it looks good or not.

    Where the point sampling distances are concerned, we will assume that the rod person can step out distances (or angle sweeps) within some percentage of that intended. We will also make the pragmatic assumption that the prism person has decided a definite number of samples per given arc and will try to make them evenly spaced. The prism person will try to subdivide the remaining space at each sample point, and two thirds of the time, they will be within some percentage of the distance (or angle sweep) they intended.

    Here is the perturbedArc() function with a sample usage:


    Fig. 2 - This is a perturbed circular arc. Doesn't look very perturbed does it? If we bumped up the standard deviation of the radius enough, we'd get something more clearly perturbed.
    To produce the matrix form of the data I use a simple function I created (see here) called MakeMatrix(). Beyond that, it is a simple as

    m: MakeMatrix(pts)$
    lsquares_estimates(m, [x,y], (x-a)^2 + (y-b)^2 = r^2, [a,b,r]);
    subst(%[1],sse(a,b,abs(r),pts));
    Posted by Darren Irvine at 1:26 PM 0 comments
    Labels: ,

    Thursday, August 6, 2015

    Polylines: Radius-Bulge Turnaround

    The dxf representation of LWPOLYLINEs represents arcs using the end points of the arc and a bulge factor which is the tangent of a quarter of the central angle of the arc (see this previous post). We want to interchange between the radius and the bulge factor (both ways). Let's pull up the earlier diagram.

    Looking at the diagram, we can produce several useful formulas. Taking the bulge factor as \(b\), observe that
    \[b = \tan {\theta \over 4} = {i \over u/2} = {2 i \over u}\tag{1},\]
    where we are temporarily ignoring the issue of sign.
    The distance \(u\) is simply the distance from \(A\) to \(B\) and these point are always given to us in the context of a polyline whether we are given the bulge factor or the radius. We write this as 
    \[u = ||B-A||.\tag{2}\]
    Using the Pythagorean theorem we have
    \[r^2 = {u^2 \over 4} + a^2 .\tag{3}\]
    Using the substitution \(a=r-i\), we obtain
    \[r^2= {u^2 \over 4} + (r-i)^2 = {u^2 \over 4} + r^2 - 2 r i + i^2\]
    and so
    \[2 r i = {u^2 \over 4} + i^2.\tag{4}\]

    Radius/bulge direction => Bulge factor

    Using equation (3), we can solve for \(a\)
    \[a = \sqrt{r^2 - \frac{u^2}{4}}.\]
    From this we can obtain \(i = r - a\) which allows us to calculate \(b\) according to equation 1. Additionally, set the sign of \(b\) to be positive for a CCW arc and negative for a CW arc. Alternatively stated, \(b\) is positive if the peak of the arc is a right offset from ray \(\overrightarrow{AB}\) and negative if the peak is a left offset from ray \(\overrightarrow{AB}\). In our diagram above, if point A comes before point B in the order they occur in the polyline, then we have a CW arc (P offset to the left of \(\overrightarrow{AB}\)) and so \(b\) should be negative. A negative bulge value indicates bulging to the left and a positive value indicates bulging to the right.

    Bulge factor => Radius/bulge direction

    We calculate \(u\) from \(A\) and \(B\) and then substitute equation 1 (solved for \(i\)) into equation 4 to get
    \[2 r \frac{b u}{2} = \frac{u^2}{4} + \frac{b^2 u^2}{4}\]
    \[4 r b =  u + b^2 u\]
    \[r = u \left(\frac{1 + b^2}{4 b}\right).\]
    Note that we would actually use the absolute value of \(b\) to obtain a positive radius. We can use the same logic as above for determining CW versus CCW (it's "iff" logic).

    Locating the Center Point

    Whether we are given the radius or bulge factor, we may want to find the center point. From above, we know we have \(u\), \(i\), \(b\), \(r\), and \(a\) all easily available to us if we need them.

    We take our unit vector for the direction of the line segment taken from \(A\) to \(B\) as \(d = (B-A)/u\). Then the normal direction (90° right turn for offsetting) is given by
    \[N = (d_y, -d_x).\]
    Taking \(\sigma = {\rm signum}(b)\),
    \[C = \frac{A+B}{2} - \sigma a N\]
    \[P = \frac{A+B}{2} + \sigma i N\].

    In my next post, I show a way to implement this in Maxima.
    Posted by Darren Irvine at 1:25 PM 0 comments
    Labels: , ,

    Wednesday, May 6, 2015

    >> Imitation Forward Pipe for Common Lisp

    Since learning a little F# about a year ago, I took a liking to its forward pipe operator (|>). The semantics of the forward pipe work great with F#'s sequence data type (with so-called, "lazy" computation) and functional style programming. It allows you to express a progression of consecutive steps in the order they will be done, while allowing you to directly feed the results of the last function into the input of the next function, avoiding redundant variable declarations. I suppose you could call these types of consecutive actions "transitive" combinations of operations. Imperative style programming would preserve the apparent sequence of events (meaning, the order you type the operations in, is the order the operations happen in), but it also involves declaring intermediate variables, which perhaps you will only use once.

    I'm not an expert on functional style programming and don't intend, right now, to explain or justify its merits, except to say that sometimes it produces cleaner, clearer code. But while I have a current preference for Common Lisp above most other (programming) languages at the present, I really like that forward pipe operator in F#. This morning, I thought it was time to cross that bridge—try to make an imitation forward pipe in Common Lisp.

    In F#, if you wanted to turn a into b into c, you might do something like

    let a2c x =
        a2b x
        |> b2c

    where b2c also takes one parameter. In Lisp, you end up with something more like

    (defun a->c (x)
      (b->c (a->b x)))

    This is fine for being concise but it doesn't preserve the (apparent) logical ordering of \(a\to b\to c\). We could do that with a change to something more imperative, such as

    (defun a->c (x)
      (let ((b (a->b x)))
        (b->c b)))

    and this works okay. However, I would like to write

    (defun a->c (x)
      (>> 
        (a->b x)
        (b->c)))

    Or,

    (defun a->e (x)
      (>> :pipe-in
        (a->b x)
        (b->c)
        (c->d :pipe-in 0.0625d0)
        (d->e)))

     where we want this last one to expand into

    (D->E (C->D (B->C (A->B X)) 0.0625))

    I've introduced a keyword so that we can dictate which parameter of the next function the previous result goes into. You can specify any keyword to use as a pipe-in parameter, but make sure you don't use a keyword of one of the functions being called in the sequence—it'll really mess with your mind. Here's the macro definition along with an error condition:



    To test that the code works, we can try this at the REPL,

    CL-USER> (macroexpand-1 '(>> :pipe-in
     (a->b x)
     (b->c) 
     (c->d :pipe-in 0.0625) 
     (d->e)))
    (D->E (C->D (B->C (A->B X)) 0.0625))

     This macro does not attempt to handle multiple return values—though such a development might lead to an interesting result. It might be as easy as introducing more keywords.
    Posted by Darren Irvine at 1:23 PM 0 comments
    Labels: , , ,

    Monday, April 6, 2015

    Dynamic Zoom in LispWorks' CAPI

    Using the LispWorks' CAPI, I have implemented a simple program to demonstrate the basic elements in the implementation of a 2D dynamic zoom feature. The program is basically uninteresting except for the dynamic zoom implementation. You can start with the code immediately, or skip to the other side for a few comments and refer back.



    Every time the user clicks the left button (:button-1), the mouse location is trapped and turned into a "model" location for a line starting position. If the user presses the right button (:button-3) the common terminal point for all of the lines is changed—the same principle of storing a "model" point is followed. Basically it is a spider with as many legs as you want originating from whatever point you select.

    When the program commences, the top-left corner is (0, 0) and view scale is 1:1. The rolling of the scroll wheel of the mouse will change both of these parameters by scaling relative to the current location of the mouse pointer. Since the scroll event does not appear to provide a means to get the current mouse pointer location, I am tracking the mouse :motion "gesture" (basically an "event") and storing the pointer location, in screen coordinates, as is moves.

    The change-start (I guess I should have said "add-start"—ah, scope creep) and change-center handler functions do basically the same thing. Using the current mouse position (these two handlers do receive mouse position information), determine the model point corresponding to the screen position of the mouse. The top-left corner is an actual ("model") location. So the model location that corresponds to the given screen coordinates is the top-left corner plus the model distance between the mouse pointer and the top left corner of the screen. Thus, \(M' = (M/s + TL)\), where \(M'\) is the model location. We manipulate the same formula to obtain screen coordinates (\(M\)) from the model coordinates when we display to the screen.

    The part that is a little trickier is to adjust the top-left corner when dynamic zoom is applied so that you can maintain the mapping between model and screen coordinates. This is an important part of change-scale's job. Our basic formula is
    $$TL_i + A_i = TL_f + A_f,$$ where \(TL_i, TL_f\) are the initial and final top-left corners (i.e., before and after zoom) and \(A_i, A_f\) are the initial and final distances between the model position of the mouse and the top-left corner. In other words, zooming does not change the mouse location in either screen or model coordinates. Hence,
    $$TL_f = TL_i + A_i - A_f = TL_i + M {{s_f - s_i}\over{s_f s_i}},$$where \(s_i, s_f\) are the initial and final scale factors.

    LispWorks Personal Edition for non-commercial use is available here: http://www.lispworks.com/downloads/.
    Posted by Darren Irvine at 1:22 PM 0 comments
    Labels: , , ,
    Subscribe to: Posts (Atom)