gmane.comp.mathematics.maxima.general

Re: lisp-stat for Maxima

Robert Dodier — 2010-03-19T21:55:31

On Fri, Mar 19, 2010 at 9:42 AM, Richard Fateman
<fateman< at >cs.berkeley.edu> wrote:


I looked at the stuff at github and it seems to be be pretty
limited in functionality. I don't think I'd want to import the
whole thing into Maxima, and I don't see any discrete chunks
that could be imported. Maxima already has a lot of that
stuff scattered in different corners (statistical stuff, linear
algebra, plotting).

From what I can tell Lisp-stat is much more limited than R.


Yeah, that's worth considering. Another possibility is
to write a parser for R scripts in Lisp such that they can
be parsed and loaded into Maxima. Well, to make that
really workable one would also need FFI to link the C
code which is part of many R packages.

FWIW

Robert Dodier

lisp-stat for Maxima

Richard Fateman — 2010-03-19T15:42:10

In looking for the graphing code (see previous message),
I came across a directory I had of a large body of statistical
computing code. This code, written by experts, is a competitor to "S"
and was written in  lisp, a number of years ago. In fact, a whole
lisp system was built on it: Lisp-stat.   It lost out in competition
to "S", for reasons that have to do more with market pressure rather
than technical prowess.  (there is a history on this..)

anyway

http://github.com/blindglobe/common-lisp-stat

is a common lisp port of the system. 

I have an older version from Luke Tierney, which includes explicit
ports with pieces of C code to optimize KCL (now GCL), as well as
a number of other lisps.  I have not looked at the github code.

But numerical statistical computing in lisp is available, and
perhaps sufficiently portable as to become useful to people.

(if you take this up seriously, you should also consider whether
linking to "S" is feasible or even preferable.)


RJF

Re: graphical interface to commands?

Raymond Toy — 2010-03-19T15:41:22

Quite a while ago, I played around with maxima and McCLIM (a gui
system).  I got as far as having maxima command line interface showing
up in a mcclim window.  I wanted to be able to do more than just that,
but never got around to it. 

Since McCLIM is in Lisp, it's reasonably portable and probably works
with any Lisp supported by maxima, except, perhaps gcl.  McCLIM also
comes with SciGraph, a 2D graphing package.  I think it can provide the
information you want from graphs.  I know that after plotting a graph,
you can use the mouse to drag the plotted points around and get
coordinates of the new point.

But if you don't know Lisp, this is probably not a good place to start.

Ray

Re: graphical interface to commands?

Richard Fateman — 2010-03-19T15:28:10

Here is a paper describing the program:

http://www.cs.berkeley.edu/~fateman/papers/graphing7.pdf

the paper claims the program is 217 lines long, but I can't find it 
right now. However,
Ken Cheatham at Franz Inc, the vendor for Allegro Common Lisp, improved 
upon it and
set it up as a tutorial for learning about the user interface.

The version from Ken, in 2002, is stored in 
http://www.cs.berkeley.edu/~fateman/lisp/graphcode.lisp
and  may still work, but if you want to experiment, I suggest you get 
the current ACL version
and look at the tutorial there.  One of the features of GUI-code is that 
it tends to be a lot of
boiler plate, and so programs that support GUI-building also tend to 
have programs that automatically
generate the boilerplate code.  That is, a graphical interface where you 
place pictures of widgets here and there,
and fill in the blanks, and out comes a program.  ACL does that for you 
in the GUI "IDE" (Interactive
Development Environment). 

I understand that at least some other lisps have IDEs, but they probably 
differ in important details.
The nice thing about GUI standards is that there are so many to choose from.

Linking from Lisp to gplot is a way to "make things work" but is 
probably not the way to do interaction.
For example, with my graphing program, I can click on something on the 
screen and trigger an arbitrary
lisp function.  I can drag things around, change their color, delete 
them, etc.  And when I am happy with
the way things appear, I can take the data structure that is the 
internal form of the display, and run more
lisp or maxima programs on it.

The particular task discussed in the paper above has to do with graph 
layout.  e.g. take the organization chart
of some company and automatically (or with human assistance) arrange it 
so that it is easy to understand.

There are lots and lots of algorithms for this, and as mentioned in the 
paper, competitions between programs
at annual conferences.

RJF





Sheldon Newhouse wrote:

Re: graphical interface to commands?

Sheldon Newhouse — 2010-03-19T14:25:11

This question about the

]1,1.2000 -:00, y]

has been resolved.  It involved restarting maxima.  Some bug in the 
functioning, I suppose.

-sen

Re: graphical interface to commands?

Sheldon Newhouse — 2010-03-19T14:13:53

Where is your program?
  Thanks,
  -sen

Re: graphical interface to commands?

Richard Fateman — 2010-03-19T14:01:30

Unfortunately, the interface capabilities of different common lisp 
systems tend to be different, because
the common lisp standard has not been updated to include this kind of 
stuff.  There are some lisps
with very good graphical interfaces on some or all host operating 
systems.  I wrote a system that
manipulates directed graphs interactively and/or by command from lisp.  
It works on Windows and
Unixes.  It does not work on CLISP or GCL, just Allegro. (commercial or 
free download).

I assume some other systems have some other particular GUI setups, and 
so if you wanted to write
a GUI, you could pick one such system and write one.

RJF

Re: graphical interface to commands?

Sheldon Newhouse — 2010-03-19T13:09:36

Hi Mario,
  Thanks for the code.  It gives me a beginning.

Here is a question about the code.

When I do the plot, instead of seeing the coordinates on the screen, I 
only see the second coordinate.
The shown coordinates look like
   ]1,1.2000 -:00, y]  where y is the second coordinate.
Why is this? Also, how do I correct it so that I see both coordinates?

TIA,
  -sen

Re: graphical interface to commands?

Mario Rodriguez — 2010-03-19T11:12:42




Hello,

I don't know if it is possible to send information generated by Gnuplot
back to Maxima through the pipe. That would be very nice. But it is
possible to write some information into a file from Gnuplot and then let
Maxima read it and make something with it.

The script below is an example. First, you middle clic on the window to
select the points; after every clic, press the 'x' key to save the
coordinates into file 'data_for_maxima'. Once you have finished your
selections, press the 'd' key and the polygonal line should be plotted.

Once you press 'd', Gnuplot writes into file 'message_for_maxima.mac'
the code 'message:1$'. Meanwhile, Maxima is repeating an endless
while-loop; during every iteration, Maxima reads if variable 'message'
has changed from 0 to 1. When it is equal to 1, Maxima reads the
contents of 'data_for_maxima' and plots the segments.

See also embedded comments.

/**************** begin Maxima code *****************/
/*

An example of interacting plot.

The script creates two files:
  data_for_maxima, where points are saved after clicking
                   with middle button and then pressing the
                   'x' key.
  message_for_maxima.mac, where gnuplot saves messages for
                   Maxima; message:0 means don't plot, and
                   message:1 means read data from data_for_maxima
                   and plot the line segments.

When interactive_draw is called, both files are initialized, so that
no interferences occur between two consecutive calls.

Once you have selected the point, press the 'd' key.

Communication channels are:
    Maxima  --> Gnuplot, via pipes
    Gnuplot --> Maxima,  via files

*/


load ("draw.lisp") $

/* [x1, x2] and [y1, y2] are the x and y ranges*/
interactive_draw(x1, x2, y1, y2) :=

  block([plotted: false, message, data],

    /* create file with messages for Maxima */
    /* initial message is "no" */
    with_stdout (
      "message_for_maxima.mac",
      print("message:0$")),

    /* initial data file is empty */
    with_stdout (
      "data_for_maxima"),

    /* set draw defaults */
    set_draw_defaults(
      /*terminal = wxt,*/
      xrange = [x1, x2],
      yrange = [y1, y2],
      grid   = true),

    /* open empty window */
    draw2d(
      xy_file = "data_for_maxima",
      user_preamble = 
       "bind 'd' 'system \"echo message:1$ > message_for_maxima.mac\" ",
      points([[x1,y1]])),

    /* while-loop until a '1' messsage is read */
    while (not plotted) do (
      load("message_for_maxima.mac"),
      if message = 1
        then (plotted: true,
              data: read_matrix("data_for_maxima"),
              draw2d(
                points_joined = true,
                points(data))   )),
    'done ) $



interactive_draw(1,10,1,10);

/**************** end Maxima code *****************/



I have tested this in Linux, both with x11 and wxt terminals.

In Windows, the user_preamble option should be changed according to
msdos commands. Since in Windows the pipe concept doesn't exist, the
global behaviour may be different.

--
Mario

Re: numerical integration

dlakelan — 2010-03-19T02:59:45

Yes, I am convinced that it is right, but this kind of problem is 
related to (for example) boundary layer theories in fluid mechanics, and 
it would be interesting to see whether some alternative type of method 
could produce a similarly accurate answer.

Re: Gamma Incomplete integrator

Richard Hennessy — 2010-03-19T02:44:31

I have put this on my site (gm.mac) and also added the special case below.

intgamma(gamma_incomplete(b, a*x)/x, x);
->  -%e^-(a*x)*(a^b*gamma(1-b)*x^(2*b)*log(a*x)-a^b*psi[0](1-b)*gamma(1-b)*x^(2*b)
                                           +(-a*log(a)*expintegral_e(b,a*x)-a*hypergeometric_regularized([1-b, 
1-b],[2-b, 2-b],-a*x)*gamma(1-b)^2)*x^(b+1))

Unfortunately I don't know how to eliminate the hypergeometric_regularized() function from the answer assuming that is 
possible.  This is a symbolic solution. Maxima as far as I know cannot evaluate the hypergeometic_regularized function. 
You also cannot diff the answer to see if you get the original problem back.  Anyway I guess it is right.  I took a 
limit of an expression involving expintegral_e() to get this answer (which Maxima can do).  So this is the output from 
the limit() function.

Rich


From: Richard Hennessy
Sent: Tuesday, March 16, 2010 2:21 AM
To: Maxima List
Subject: [Maxima] Dammar Incomplete integrator


Forgot the year.

Rich




_______________________________________________
Maxima mailing list
Maxima< at >math.utexas.edu
http://www.math.utexas.edu/mailman/listinfo/maxima

Re: numerical integration

Richard Hennessy — 2010-03-19T01:57:56

Hi,

I increased the size of the interval by 66% in these two cases and it only affects the answer in the 15th digit.  I 
think that is a sound reason to believe it is very close to being right.

Rich


--------------------------------------------------
From: "dlakelan" <dlakelan< at >street-artists.org>
Sent: Thursday, March 18, 2010 9:36 PM
To: <maxima< at >math.utexas.edu>
Subject: Re: [Maxima] numerical integration

Re: numerical integration

dlakelan — 2010-03-19T01:36:34

Yes, this is similar to Raymond Toy's solution. I get a similar answer 
with trapezoid rule and an estimate of the truncation error that assumes 
similitude between the truncated part and the calculated part.

Is there a solution in terms of an asymptotic series? we assume 
integral(f(x),x,5,6) ~ integral(f(x),x,6-eps,6), expand the 
approximation in terms of the small parameter eps, and then truncate the 
expansion at some small number of terms (perhaps 2 or 3) and get some 
kind of decent approximate result?

I tried to do this via integration by parts taking

x^x^x = x * x^(x^x-1)

take x dx = dv and x^(x^x-1) = u, plug into uv - int(v du), iterate 
again, and then chop off the remaining integral... but it got me 
nowhere, the resulting series explodes.

I still feel like there should be a "clever" way to get good accuracy 
using some kind of asymptotic method. This type of extremely fast 
growing integral is a very special sort of thing.

Re: numerical integration

Richard Hennessy — 2010-03-18T22:15:28

I got it this way.


bromberg(x^x^x, x, 599975/100000b0, 6b0),brombergit=15,fpprec:2000;
->  1.1026649993619407788421378964b36300
bromberg(x^x^x, x, 599985/100000b0, 6b0),brombergit=15,fpprec:2000;

When you are this close to 6 the area under the curve between 5 and 5.99985 is irrelevant since it does not affect the 
answer (at least not the first 8-10 digits.)

Rich


--------------------------------------------------
From: "dlakelan" <dlakelan< at >street-artists.org>
Sent: Thursday, March 18, 2010 5:46 PM
To: <maxima< at >math.utexas.edu>
Subject: Re: [Maxima] numerical integration

Re: numerical integration

dlakelan — 2010-03-18T21:46:53


(to remind everyone we're trying to integrate x^x^x from x=5 to x=6)

I thought about this issue as well. It's similar but even more extreme 
than calculation of a hypervolume in a very high dimensional space. 
That's done in statistical mechanics with ~10^23 dimensions, and the 
integral is approximated as simply the value of the integrand at the 
outer radius. That is for something which grows simply exponentially. 
here we have even faster growth.

Your calculation is very sound with the only thing I'd like to check 
being the accuracy of the bfromberg step. I tried a similar thing where 
I calculate a point where g(x)/g(6) < epsilon for a small epsilon, and 
then use

integrate(g(x),x,5,6) = integrate(g(x),x,5,6-epsilon) + 
integrate(g(x),x,6-epsilon,6).

We can approximate the first term as 1*g(6-epsilon) for an upper bound 
on the truncation error, and then spend some time getting different 
estimates of the second term. I did that as follows (see below) and I 
get an estimate that is 1.09b36300 + O(1b36298) using a 15th order 
taylor series, which is consistent with your 1.1b36300 estimate.

I am trying to figure out how to use renormalization group theory to 
solve this problem, assuming that the left side integral is self-similar 
to the right side integral.

For example, when I require that f(x-eps) is the average value of the 
integral as calculated from a partition [5,6-eps] and [6-eps,6] with the 
right side integral calculated using simpson's rule, and the left side 
integral calculated using self-similarity with the right side integral, 
I get the estimate: 2.58b36300 which is also relatively consistent.

------------------ output of my maxima script --------

(%i212) kill(all);

(%o0) done
(%i1) display2d:false;

(%o1) false
(%i2) fpprec:60;

(%o2) 60
(%i3) fpprintprec:8;

(%o3) 8
(%i4) ratepsilon:1e-15;

(%o4) 1.e-15
(%i5) ratprint:false;

(%o5) false
(%i6) f(x) := x^x^x;

(%o6) f(x):=x^x^x
(%i7) g(x) := f(x)/f(6.0b0);

(%o7) g(x):=f(x)/f(6.0b0)
(%i8) 
draw2d(xrange=[5.9999,6.00001],explicit(lambda([x],float(g(bfloat(x)))),x,5.9999,6));

(%o8) [gr2d(?explicit)]
(%i9) epslower:5b-8;

(%o9) 5.0b-8
(%i10) 
epsb3:bfloat(6-find_root(lambda([x],float(g(bfloat(x))-epslower)),x,5.9995,5.99999));

(%o10) 6.9717956b-5
(%i11) matchdeclare(a,atom);

(%o11) done
(%i12) matchdeclare(x,atom);

(%o12) done
(%i13) tellsimpafter(a*O(x),O(a*x));

(%o13) [?\*rule8,?simptimes]
(%i14) trapest:bfloat((g(6-eps) + g(6))/2*eps + O(g(6-eps)));

(%o14) O(3.7606429b-36306*(6.0b0-1.0b0*eps)^(6-eps)^(6-eps))
         +5.0b-1*eps
                *(3.7606429b-36306*(6.0b0-1.0b0*eps)^(6-eps)^(6-eps)
                 +9.9999999b-1)
(%i15) trapest,eps=epsb3;

(%o15) O(5.0b-8)+3.4858979b-5
(%i16) expand(bfloat(%*f(6.0b0)));

(%o16) O(1.3295598b36298)+9.2694202b36300
(%i17) simpest:expand(bfloat((g(6-eps) + g(6) + 4*g(6-eps/2))/6*eps + 
O(g(6-eps)))*f(6.0b0));

(%o17) 2.6591197b36305*O(3.7606429b-36306*(6.0b0-1.0b0*eps)^(6-eps)^(6-eps))
         +6.6666666b-1*eps*(6.0b0-5.0b-1*eps)^(6-eps/2)^(6-eps/2)
         +1.6666666b-1*eps*(6.0b0-1.0b0*eps)^(6-eps)^(6-eps)
         +4.4318662b36304*eps
(%i18) simpest,eps=epsb3;

(%o18) O(1.3295598b36298)+3.0925691b36300
(%i19) tayg: bfloat(taylor(g(x),x,6-epsb3/2,15));

(%o19) 9.2721644b64*(x-5.9999651b0)^15+5.7669108b60*(x-5.9999651b0)^14
                                       +3.3477067b56*(x-5.9999651b0)^13
                                       +1.8045635b52*(x-5.9999651b0)^12
                                       +8.9792514b47*(x-5.9999651b0)^11
                                       +4.0956685b43*(x-5.9999651b0)^10
                                       +1.6983299b39*(x-5.9999651b0)^9
                                       +6.3382174b34*(x-5.9999651b0)^8
                                       +2.1026401b30*(x-5.9999651b0)^7
                                       +6.1034602b25*(x-5.9999651b0)^6
                                       +1.5186083b21*(x-5.9999651b0)^5
                                       +3.1487593b16*(x-5.9999651b0)^4
                                       +5.2231011b11*(x-5.9999651b0)^3
                                       +6.4980634b6*(x-5.9999651b0)^2
                                       +5.3895621b1*(x-5.9999651b0)
                                       +2.2351069b-4
(%i20) 
draw2d(yrange=[-1,1],explicit(lambda([x],float(g(bfloat(x)))),x,6-epsb3,6+epsb3/1000),color='blue,
   explicit(float(subst('x=bfloat(x),tayg)),x,6-epsb3,6));

(%o20) [gr2d(?explicit,?explicit)]
(%i21) taygs:subst(s,bfloat(x-(6-epsb3/2)),tayg);

(%o21) 
9.2721644b64*s^15+5.7669108b60*s^14+3.3477067b56*s^13+1.8045635b52*s^12
 
+8.9792514b47*s^11+4.0956685b43*s^10+1.6983299b39*s^9
                         +6.3382174b34*s^8+2.1026401b30*s^7+6.1034602b25*s^6
                         +1.5186083b21*s^5+3.1487593b16*s^4+5.2231011b11*s^3
                         +6.4980634b6*s^2+5.3895621b1*s+2.2351069b-4
(%i22) taygsest:bfloat(integrate(taygs,s,-epsb3/2,epsb3/2));

(%o22) 4.1127528b-6
(%i23) tayest:expand((% + O(epslower))*f(6.0b0));

(%o23) O(1.3295598b36298)+1.0936302b36300
(%i24) simpint2:(g(6-eps) + g(6)+4*g(6-eps/2))*eps/6;

(%o24) (1.5042571b-36305*(6-eps/2)^(6-eps/2)^(6-eps/2)
         +3.7606429b-36306*(6-eps)^(6-eps)^(6-eps)+9.9999999b-1)
         *eps
         /6
(%i25) epseqn:g(6-eps) * 1 = simpint2 + (1-eps)/eps*simpint2*g(6-eps);

(%o25) 3.7606429b-36306*(6-eps)^(6-eps)^(6-eps)
          = (1.5042571b-36305*(6-eps/2)^(6-eps/2)^(6-eps/2)
          +3.7606429b-36306*(6-eps)^(6-eps)^(6-eps)+9.9999999b-1)
          *eps
          /6
          +6.2677382b-36307*(1.5042571b-36305*(6-eps/2)^(6-eps/2)^(6-eps/2)
                            +3.7606429b-36306*(6-eps)^(6-eps)^(6-eps)
                            +9.9999999b-1)*(1-eps)*(6-eps)^(6-eps)^(6-eps)
(%i26) epsexpr: lhs(epseqn) - rhs(epseqn);

(%o26) -(1.5042571b-36305*(6-eps/2)^(6-eps/2)^(6-eps/2)
         +3.7606429b-36306*(6-eps)^(6-eps)^(6-eps)+9.9999999b-1)
         *eps
         /6
         -6.2677382b-36307*(1.5042571b-36305*(6-eps/2)^(6-eps/2)^(6-eps/2)
                           +3.7606429b-36306*(6-eps)^(6-eps)^(6-eps)
                           +9.9999999b-1)*(1-eps)*(6-eps)^(6-eps)^(6-eps)
         +3.7606429b-36306*(6-eps)^(6-eps)^(6-eps)
(%i27) 
foundeps:find_root(lambda([x],float(ev(epsexpr,'eps=bfloat(x)))),x,1e-8,1e-3);

(%o27) 4.78631498e-5
(%i28) g(6.0b0-foundeps)*f(6.0b0);

(%o28) 2.5836545b36300

Re: choicesin in sin.lisp

Robert Dodier — 2010-03-18T21:27:25


If you can verify that it's equivalent to some formulation
in terms of REMOVE or whatever, then it's OK by me to
replace it.

Robert Dodier

Re: discontinuos function

Richard Hennessy — 2010-03-18T20:34:27

_______________________________________________
Maxima mailing list
Maxima< at >math.utexas.edu
http://www.math.utexas.edu/mailman/listinfo/maxima

Re: choicesin in sin.lisp

Raymond Toy — 2010-03-18T18:36:49

Given the experience Dieter had with one file where changes shouldn't
have had any affect but did, perhaps we can leave choicesin but replace
the body with remove.

But I note that choicesin and remove are not exactly the same. 
(choicesin 1 (list 1 2 3)) returns a list which is EQ to the CDR of the
second arg.  REMOVE doesn't, at least on cmucl.  Don't know if that
matters or not.

Ray

Re: choicesin in sin.lisp

Richard Fateman — 2010-03-18T18:34:44

I haven't looked at the definition, but a reason why it is there may be 
that in Maclisp in 1964, there
was no built-in remove with :count.

choicesin in sin.lisp

Andreas Eder — 2010-03-18T17:40:40

Hi,

does anyone know why there is choicesin in sin.lisp?
It is used only there, and apart from the apparent stack overflow
(choicesin  a b) seems to be equivalent to  (remove a b :count 1).

Is anyone opposed to replacing it?

Andreas

Re: discontinuos function

Jaime Villate — 2010-03-18T17:30:37

Do you mean something like this?:
(%i1) f(a,x):= if x > a then x-a else 0$
(%i2) plot2d(f(2,x), [x,-4,4], [y,-1,5])$

Regards,
Jaime