Rcjp’s Weblog

February 1, 2007

globbing

Filed under: python, unix — rcjp @ 11:23 am

I hate globbing. Whoever thought that it was a good idea to let the shell snaffle the arguments you’ve given to a program needs to be pushed off somewhere high onto some burning cactii.

You can turn it off in the shell, but this makes most unix commands useless since they don’t do the file expansion themselves (like they should).

In ruby…

    for f in Dir.glob("/home/r/tmp/**/*")
        File.open(f) { |file|
            puts "filename #{f}"
        } if File.file?(f)
    end

python…

    import glob, os

    for root, dirs, files in os.walk('/home/r/tmp'):
        for f in files: print 'file -> ', os.path.join(root,f)

The Common Lisp pathname stuff can be a bit confusing, which is a pity because it does have some nice tricks up its sleeve.

You build the directory section of the pathname as a list taking keyword arguments :wild where you want to match at one directory level and :wild-inferiors when you want to recurse.

CL-USER> (make-pathname :name :wild
                         :type :wild
                         :directory (list :absolute "home" "r" "tmp" :wild-inferiors))
#P"/home/r/tmp/**/*.*"

so running (DIRECTORY #P"/home/r/tmp/**/*.*") should give you all the files in the tmp directory and below. Some compilers though take directories themselves to be files and others don’t, so you might need to use FILE-NAMESTRING to see if the filename part of the path is nil i.e. its a directory.

January 23, 2007

Symbol-function

Filed under: lisp — rcjp @ 9:36 am

You can make an adder function that changes its closure variable via an optional argument

    (defun make-variable-adder (n)
      (lambda (x &optional change)
        (if change
            (setq n x)
            (+ x n))))

    (setq add2 (make-variable-adder 2))
    (mapcar add2 '(1 2 3)) => (3 4 5)

But then whilst trying to change the closure variable I forgot that (add23 t) won’t work, you need (funcall add2 3 t) since (symbol-function 'add2) hasn’t been defined, only the symbol-value has. The only reason mapcar works is that as the hyperspec says… “If function is a symbol, it is coerced to a function as if by symbol-function.”

January 15, 2007

Swapping Spaces for Punctuation

Filed under: c, lisp, python — rcjp @ 10:17 am

Just a quicky. In Common Lisp…

    (defun spacify-punc (string)
      (substitute-if #\Space #'punc-p string))

    (defun punc-p (char) (find char "*,.:;-()"))

then use it…

    CL-USER> (spacify-punc "this, string. and: more")
    "this  string  and  more"

in python…

    def spacify_punc(char):
        if char not in "*,.:;-()":
            return char
        else:
            return ' '

    print ''.join(spacify_punc(c) for c in "this, string. and: more")

in C possibly something like…

    #include <stdio.h>
    #include <string.h>

    int main()
    {
      char words[] = "this, string. and: more";
      char *answer = strdup(words);
      unsigned int i;

      /* using sizeof (instead of strlen) so we include the NULL */
      for (i=0; i < sizeof(words); i++)
        answer[i] = strchr("*,.:;-()", words[i]) ? ' ': words[i];

      printf("'%s' becomes '%s'\n", words, answer);
      return 0;
    }

All three languages have functions to detect punctuation, but we are only looking for a subset in this case.

I think I like the CL solution best but the order of args to substitute-if always throws me, thank heavens for SLIME’s argument prompting.

January 11, 2007

Quick Graphviz Test

Filed under: python — rcjp @ 10:14 am

I occassionally need something to show a connection graph and graphviz is a handy tool especially using the python pydot interface, just define some node links…

    r@laptop:~/src/python/pydot-0.9.10$ sudo python setup.py install
    running install
    running build
    running build_py
    creating build
    creating build/lib
    copying pydot.py -> build/lib
    copying dot_parser.py -> build/lib
    running install_lib
    copying build/lib/pydot.py -> /usr/lib/python2.4/site-packages
    copying build/lib/dot_parser.py -> /usr/lib/python2.4/site-packages
    byte-compiling /usr/lib/python2.4/site-packages/pydot.py to pydot.pyc
    byte-compiling /usr/lib/python2.4/site-packages/dot_parser.py to dot_parser.pyc


then, running in ipython

    In [35]: import pydot

    In [36]: edges=[(1,2), (1,3), (1,4), (3,4)]

    In [37]: g=pydot.graph_from_edges(edges)

    In [38]: g.write_png('test-graph.png', prog='dot')
    Out[38]: True


test-graph

January 7, 2007

Renaming Groups of Files

Filed under: python, unix — rcjp @ 4:13 pm

I have done it in the times in the past, but I always forget the-right-way because different unix shells do loops differently etc, so I now always use ipython to e.g. lowercase all filenames in a directory…

    In [14]: import os

    In [15]: !! ls
    Out[15]:
    ['DEB1.GLE',
     'DEB1.eps',
     'DEBF1.RES',
     'DEBF2.RES',
     'DEBG1.RES',
     'DEBG2.RES',
     'DEBYE1.EXE',
     'DEBYE1.FOR',
     'DEBYE2.EXE',
     'DEBYE2.FOR']

    In [16]: for f in _15: os.rename(f, f.lower())


and for doing that recursively of course just do

    !! find . -type f

More nvi tricks

Filed under: unix — rcjp @ 2:31 pm

Came across Garrett Hildebrand’s nice collection of vi tricks including a cracking one to insert a previous occurring word – handy for long function names etc.

January 6, 2007

Twenty Questions Game in C

Filed under: c — rcjp @ 3:33 pm

Tree traversal is really really easy in Common Lisp so for a bit more of a challenge I wrote this in C. The game is for the user to think of an object and answer questions while the program tries to guess. The code is just a binary search tree, but you have to look ahead to see if you need to add in a piece of tree rather than just a straighforward leaf node.

    #include <stdio.h>
    #include <string.h>
    #include <stdbool.h>
    #include <stdlib.h>

    #define LINE_MAX 512

    typedef struct node NODE;
    struct node {
        char qu[LINE_MAX];
        NODE *yes, *no;
    };

    /*
     * print message and read stdin into LINE_MAX reply buffer
     * trimming newline
     */
    void getreply(char *message, char reply[LINE_MAX])
    {
        puts(message);
        fflush(stdout);
        fgets(reply, LINE_MAX, stdin);
        strtok(reply, "\n");        /* trim off \n */
    }

    bool ask(char *question, bool object)
    {
        char line[LINE_MAX], reply[LINE_MAX], qu[LINE_MAX];

        if (object)
            sprintf(qu, "Is it %s?", question);
        else
            strcpy(qu, question);
        getreply(qu, line);
        sscanf(line, "%s", reply);
        if (!strcmp(reply, "yes") || !strcmp(reply, "y"))
            return true;
        return false;
    }

    int main()
    {
        NODE info = { "Is it an animal?",
            &(NODE) {"Is it slimy?",
                     &(NODE) {"a fish", NULL, NULL},
                     &(NODE) {"a bear", NULL, NULL}},
            NULL
        };

        char object[LINE_MAX], clarify[LINE_MAX];
        bool yesno;
        do {
            NODE *curr = &info, *next, *newnode;
            while (curr) {
                yesno = ask(curr->qu, false);
                next = yesno ? curr->yes : curr->no;
                if (!next) {
                    /*  
                     * we need to add a leaf 
                     */
                    curr->no = malloc(sizeof(NODE));
                    getreply("What is it? ", curr->no->qu);
                    break;
                }
                if (!next->yes && !next->no) {
                    /* 
                     * we might need to add a new node in here 
                     */
                    if (ask(next->qu, true))
                        break;      /* the leaf was right.. new game */

                    newnode = malloc(sizeof(NODE));
                    getreply("What is it? ", object);
                    sprintf(clarify, "What question distiguishes %s from %s?",
                            object, next->qu);
                    getreply(clarify, newnode->qu);

                    /* check which way round this new question applies */
                    sprintf(clarify, "A %s, %s", object, newnode->qu);
                    if (ask(clarify, false)) {
                        newnode->yes = malloc(sizeof(NODE));
                        strcpy(newnode->yes->qu, object);
                        newnode->no = next;
                    } else {
                        newnode->no = malloc(sizeof(NODE));
                        strcpy(newnode->no->qu, object);
                        newnode->yes = next;
                    }
                    if (yesno)
                        curr->yes = newnode;
                    else
                        curr->no = newnode;
                    break;
                }
                curr = next;
            }
        }
        while (ask("New Game", false));
        return 0;
    }


    /*   Sample run...

        Is it an animal?
        y
        Is it slimy?
        n
        Is it a bear?
        y
        New Game
        y

        Is it an animal?
        n
        What is it?
        a box
        New Game
        y

        Is it an animal?
        n
        Is it a box?
        n
        What is it?
        a banana
        What question distiguishes a banana from a box?
        is it yellow?
        A a banana, is it yellow?
        y
        New Game
        y

        Is it an animal?
        n
        is it yellow?
        y
        Is it a banana?
        n
        What is it?
        the sun
        What question distiguishes the sun from a banana?
        is it cold?
        A the sun, is it cold?
        n
        New Game
        y

        Is it an animal?
        n
        is it yellow?
        y
        is it cold?
        n
        Is it the sun?
        y
        New Game
        n
        r@laptop:~/c/tmp$
    */

I started thinking about how this would look in a C++ object type program but, to be honest, I don’t enjoy trying to think out the right class structure; and what is annoying is that you can’t start coding and explore the problem as you code, because you’ll come to a grinding halt if you’ve got the structure wrong. I guess I just don’t think OOP is fun coding, at least not in C++.

December 30, 2006

Lyaponov Exponent

Filed under: c — rcjp @ 6:42 pm

Most introductions to chaos theory will talk about the logistic mapping x_{n+1}=r\, x_n(1-x_n) and then calculate the lyapunov exponent which is a measure of how rapidly small deviations in the initial conditions give exponential changes in results giving an indication how chaotic the system is. It is calculated from \lambda = \frac{1}{N} \sum^N_{n=1} \log(\frac{d x_{n+1}}{d x_n}) The program below calculates this for 2.95<r<4.0 after iterating for 1000 times to stabalise things then using N=30.

The code is slightly contrived; I wrote it to play with passing functions through C++ templates. The idea being that logis and deriv (its derivative) since they are templated, can get inlined by the compiler; though Id’ guess a much more straightforward C program would probably optimise just as well. In addition the std::bind1st(std::mem_fun(&IterMap::lyaponov_func),this) didn’t exactly trip off my fingers.

#include <iostream>
#include <fstream>
#include <algorithm>
#include <vector>
#include <numeric>
#include <functional>
#include <cmath>

typedef double (*CALC)(double, double);

class IterMap
{
    std::vector<double> data;
    const int nstabalise;
    const int nrecord;
    double r;
 public:
    IterMap(int rec, int stab=1000) : nstabalise(stab), nrecord(rec) {}
    template<CALC derivfunc>double lyaponov_func(double x);
    template<CALC mapfunc> void stabalise(double x);
    template<CALC mapfunc, CALC derivfunc> double lyapunov_exp(double r);
};

template<CALC mapfunc> void IterMap::stabalise(double x)
{
    for(int i = 0; i < nstabalise; ++i)
        x = mapfunc(r, x);

    for(int i = 0; i < nrecord; ++i) {
        x = mapfunc(r, x);
        data.push_back(x);
    }
}

template<CALC derivfunc> double IterMap::lyaponov_func(double x)
{
    return log(fabs( derivfunc(r,x) ));
}

template<CALC mapfunc, CALC derivfunc> double IterMap::lyapunov_exp(double ra)
{
    r = ra;
    stabalise<mapfunc>(0.1);

    // used mem_func, else have to declare func static

    transform(data.begin(), data.end(), data.begin(),
            std::bind1st(std::mem_fun(&IterMap::lyaponov_func<derivfunc>),this));
    double avg = accumulate(data.begin(), data.end(), 0.0) / data.size();
    data.clear();
    return avg;
}

////////////////////////////////////////////////////////////////////////////////

double logis(double r, double x)
{
    return r * x * (1.0 - x);
}

double deriv(double r, double x)
{
    return r - 2.0*r*x;
}

int main(int argc, char* argv[])
{
    std::ofstream out("ldata", std::ios::out);
    IterMap bir(30);
    for (double r = 2.95; r <= 4.0; r+=0.0005)
        out << r << " " << bir.lyapunov_exp<logis,deriv>(r) << std::endl;

    return 0;
}

December 29, 2006

Quick C++ Word Count Trick

Filed under: c — rcjp @ 2:19 pm

Probably not really useful in actual code but

#include <iostream>
#include <iterator>
#include <string>

typedef std::istream_iterator<std::string> in;

int main()
{
    std::cout << std::distance(in(std::cin), in());
}


will print the number of words piped in.

December 26, 2006

lotto numbers

Filed under: lisp, python — rcjp @ 9:26 am

I can’t remember how I came to be subscribed to a perl feed but I ended up reading http://www.oreillynet.com/onlamp/blog/2006/12/99_problems_in_perl_6.html
which was based on something similar in Lisp I think. The article showed some Scheme like Common Lisp solution to picking out the lottery numbers all cons’s etc.

My first somewhat crumby solution was

    (defun lotto (n m)
      "Draw N different random numbers from the set 1..M"
      (loop for num = (1+ (random m))
         for bag = nil then (adjoin num bag)
         while (< (length bag) n)
         finally (return bag)))

as I first thought of the adjoin function when I thought about adding numbers, but you shouldn’t be adjoining and then setting the result back to itself – pushnew would do that. Maybe its clearer just to use member

    (defun lotto1 (n m)
      (loop while (< (length bag) n)
         for num = (1+ (random m))
         unless (member num bag) collect num into bag
         finally (return bag)))

I vaguely remembered a thread on comp.lang.lisp on this subject and found a post from Alan Crowe did it prettier with ‘do’ – slapped wrists for me since I always dismiss ‘do’ as being harder to read than loop

    (defun lotto2 (n m)
      (do ((draw '()))
          ((= (length draw) n) draw)
        (pushnew (+ (random m) 1) draw)))

Of course these don’t work well if you want 49 out of 50 as there are too many non-unique numbers generated, Brian Downing posted a version using a nice idea of generating a list of numbers with an associated random number and sorting on that random number to shuffle the list

    (defun lotto3 (n m)
      (subseq (mapcar #'car (sort (loop for i from 1 upto m
                                     collect (cons i (random 1.0)))
                                  #'< :key #'cdr))
              0 n))

In python, as I was looking up the random function, the docs showed how you can do this sort of thing in one line (remember xrange, like range, generates numbers upto but not including the second arg)

    In [16]: random.sample(xrange(1,50),6)
    Out[16]: [29, 28, 32, 45, 42, 2]

December 22, 2006

String Searching/Replacing

Filed under: c, lisp, python — rcjp @ 9:04 am

just some quick notes on how various languages handle replacing/splicing a string…

In C++ chopping, searching and replacing in a string is fairly easy

    #include <iostream>
    #include <string>

    int main()
    {
        std::string s = "some one with more than one that ones.";

        s.erase(0,3);
        s.replace(s.find("one"), 3, "three");

        std::cout << s << std::endl;
    }

in straight C doing the same thing is a bit more fiddly

    #include <stdio.h>
    #include <string.h>

    int main()
    {
        char* s = "some one with more than one that ones.";

        char buf[255];
        strcpy(buf, s+3);  /* chop off the first 3 chars */

        char *p = strstr(buf, "one");
        if (p)
        {
            char tmp[255];
            *p = (char) 0;
            strcpy(tmp, buf);
            strcat(tmp, "three");
            strcat(tmp, p+3);  /* skip over length of "one" */
            strcpy(buf, tmp);
        }

        printf("%s\n", buf);
        return 0;
    }

you have to think about the size of the temporary buffer unless you malloc something based on the size of s – but then you have to know how much bigger your operations will make the string.

Perhaps suprisingly there isn’t any standard function to replace strings in Common Lisp – I guess its one of those things your are supposed to deal with yourself since if you know you your replacement word is the same size you can destructively alter the string (using setf on the subseq), otherwise you have to build a new string (since you may be replacing a word with a bigger word). So CL leaves things to you to figure out the best approach.

You can get the position of a string within another with

    * (search "one" "there one is one more than ones")
    6

and then build up the string a la C…

    (let* ((str "some one with more than one that ones.")
           (word "one")
           (p (search word str)))
      (if p
          (concatenate 'string (subseq str 3 p)
                               "three"
                               (subseq str (+ p (length word))))
          p))

In python

    In [3]: "there one is one more than ones".replace("one", "three", 1)
    Out[3]: 'there three is one more than ones'

where we are using a count=1 otherwise it would replace all instances. We can even chop off the first three chars and then replace all in one go

    In [5]: "there one is one more than ones"[3:].replace("one", "three")
    Out[5]: 're three is three more than threes'

Replacing all occurances in C++ is a bit more work

    std::string::size_type n;  // or   size_t n; 
    std::string word = "one";
    while((n=s.find(word)) != std::string::npos)
        s.replace(n, word.size(), "three");

Replacing all strings in C or Common Lisp is quite alot more work and probably better to write a utility function and keep it somewhere or use a library.
http://en.wikibooks.org/wiki/Programming:Common_Lisp/Strings has

    (defun replace-all (string part replacement &key (test #'char=))
      (with-output-to-string (out)
        (loop with part-length = (length part)
           for old-pos = 0 then (+ pos part-length)
           for pos = (search part string :start2 old-pos :test test)
           do
             (write-string string out :start old-pos
                                      :end (or pos (length string)))
           when pos do
             (write-string replacement out)
           while pos)))

or maybe use cl-ppcre

December 21, 2006

Compiler warnings GCC/SCC

Filed under: c — rcjp @ 8:51 am

Got Salford SCC running under WINE (except for the /break debugger option) which is impressive.

Just tried running it on a couple of test programs and it didn’t like

    typedef struct box {
        const char *name;
        const int count;
    } BOX;

    r@laptop:~/c/exercises$ scc bsearch.c /ansi_c/link
    [Salford SCC/WIN32 Ver 3.16 Copyright (c) Salford Software Ltd 2002]
       0008       const int count;
    *** A structure member may not be const or volatile
        1 ERRORS  [<BSEARCH> SCC/WIN32 Ver 3.16]
    *** Compilation failed

I’m not quite sure about this though. I was initialising with

    struct box ruler = { "ruler", 0 };

C++ has member initialisers which are required for const member objects as in

    typedef struct box {
        const char *name;
        const int count;
        box();
    } BOX;

    box::box () : name(""), count(5){}

as you can’t do

    box::box () {
        name = "";
        count = 5;
    }

because that would mean that the constructors for the member object are called for count and then altered in the box constructor and you can’t alter a const! But there aren’t the same concepts in C, hmmm.

I also realised that gcc’s -Wall doesn’t mean what I thought it did – warn about everything. Infact the recommended compiler options for finding problems are (see http://www.network-theory.co.uk/docs/gccintro/index.html)

gcc -ansi -pedantic -Wall -W -Wconversion -Wshadow -Wcast-qual -Wwrite-strings

using -std=c99 if using c99 features. I wonder why c99 isn’t the default yet?

December 19, 2006

Reading in Files

Filed under: c, lisp — rcjp @ 5:52 pm

I was catching up on Dr.Dobb’s columns and came across Walter Browns clean way of reading lines from files in Andrew Koenig’s c++ blog.

    #include <iostream>    // for cout 
    #include <fstream>     // for ifstream 
    #include <string>      // for string 

    int main()
    {
        std::ifstream file("t2.cpp");

        for(std::string s; getline(file, s);)
        {
            std::cout << s << '\n';
        }
    }


which is a nice alternative to the usual while loop as the string is scoped in the for loop that uses it.

Of course you can do something similar in C, provided you use -std=c99 compiler flag to gcc to allow declarations in for loops…

    #include <stdio.h>
    #define LINE_MAX 256

    int main()
    {
        FILE *fp=fopen("t1.c","r");

        for (char line[LINE_MAX]; fgets(line, LINE_MAX, fp);)
        {
            printf("line = %s\n", line);
        }
    }


Production code should always check the results of fopen, and fgets may return a NULL because of a problem reading the file – not because it reached the end, so you should check with feof. But in ‘quick and dirty’ mode I skip these and C will segmentation fault trying to read from a missing file. C++ will just quietly carry on, neither of which is ideal.

Common Lisp will throw you into the debugger when

    (defun test ()
      (with-open-file (stream "badfilename" :direction :input)
        (loop for line = (read-line stream nil nil)
           while line do
             (format t "line = ~A~%" line))))


fails to find the file. (I’ve never really liked those nil’s that follow read-line though, the first says the end of file isn’t an error and the second is what you want to signify the end of input. I think I’d prefer it if they were keyword arguments – more descriptive.)

If you are using a lisp compiler that supports restarts (unfortunately CMUCL/SBCL don’t as far as I can tell, but clisp and Lispworks do) you can quickly patch the running code if you don’t feel like going back to the source just yet…

    r@laptop:~/lisp/tmp$ clisp -q
    [1]> (load "test.lisp")
    ;; Loading file test.lisp ...
    ;; Loaded file test.lisp
    T
    [2]> (test)

    *** - OPEN: file #P"/home/r/tmp/badfilename" does not exist
    The following restarts are available:
    ABORT          :R1      ABORT
    Break 1 [3]>


first find out where you are with :w

    Break 1 [3]> :w
    <1> #<SYSTEM-FUNCTION OPEN>
    EVAL frame for form (OPEN "badfilename" :DIRECTION :INPUT)
    Break 1 [3]>


now you can see the actual open command the with-open-file macro got expanded into, and use :rt to return the correct filename by typing in the correct value of OPEN at the “Values:” prompt and carry on execution

    Break 1 [3]> :rt
    Values: (OPEN "correctfilename" :DIRECTION :INPUT)


There is an interesting document showing some more advanced restarting of code here using acl but is worth reading even if you don’t have that implementation.

November 16, 2006

Pattern Formation

Filed under: physics — rcjp @ 10:09 am

cccc
cccc program from 'Models of Biological Pattern Formation' Hans Meinhardt
cccc 

! gnuplot> set contour
! gnuplot> unset surface
! gnuplot> set style data lines
! gnuplot> set view 0,0
! gnuplot> set cntrparam levels 30
! gnuplot> splot 'ttt'

      parameter (igridx = 31, igridy = 31)
      real*8 g1(igridx,igridy), g2(igridx,igridy)
      real*8 s1(igridx,igridy), s2(igridx,igridy)
      real*8 vz(igridx,igridy)
      real*8 dg1(30), dg2(30), dr(30), ds1(30), ds2(30)
      kx = 30
      ky = 30
      qbb = 0.3
      
      ginit = sqrt(1.d0 - qbb)
      sinit = ginit
      
      ta = 0.1
      tb = 0.1
      dia = 0.005
      qd = 0.04 
      qe = 0.2
      qa = 0.001
      qaa = 0.001

cccc  initialise the grid

      do iy = 1, ky
         do ix = 1, kx
            g1(ix,iy) = ginit
            g2(ix,iy) = ginit
            s1(ix,iy) = sinit
            s2(ix,iy) = sinit
c            vz(ix,iy) = 1.d0 + rand()*0.0001
            vz(ix,iy) = 1.d0 + rand()*0.0001
         enddo
      enddo

      jy = (ky-1)/2 + 1

cccc  initial perterbation

      g1(kx,jy) = 1.1*ginit

      do iter = 1, 5
      dgc = 1.d0 - ta - 4.d0*dia
      dsc = 1.d0 - qd - 4.d0*qe
      drc = 1.d0 - tb

      do it = 1, 500
         do ix = 1, kx
            g1(ix, ky+1) = g1(ix, ky) ! lower border
            g2(ix, ky+1) = g2(ix, ky)    
            s1(ix, ky+1) = s1(ix, ky)    
            s2(ix, ky+1) = s2(ix, ky)
            
            dg1(ix) = g1(ix, 1) ! upper border
            dg2(ix) = g2(ix, 1)    
            ds1(ix) = s1(ix, 1)    
            ds2(ix) = s2(ix, 1)
         enddo
         do iy = 1, ky
            g1(kx+1, iy) = g1(kx, iy) ! right border
            g2(kx+1, iy) = g2(kx, iy)
            s1(kx+1, iy) = s1(kx, iy)
            s2(kx+1, iy) = s2(kx, iy)
         enddo
         do iy = 1, ky
            g1l = g1(1,iy)      ! left border
            g2l = g2(1,iy)
            s1l = s1(1,iy)
            s2l = s2(1,iy)

cccc  reaction and feedback

            do ix = 1, kx
               gf1 = g1(ix,iy)
               gf2 = g2(ix,iy)
               sf1 = s1(ix,iy)
               sf2 = s2(ix,iy)
               
               g1(ix,iy) = gf1*dgc+dia*
     &                     (dg1(ix)+g1(ix,iy+1)+g1l+g1(ix+1,iy))+
     &                     sf2*ta*vz(ix,iy)/(gf2**2+qbb)+qa
               g2(ix,iy) = gf2*dgc+dia*
     &                     (dg2(ix)+g2(ix,iy+1)+g2l+g2(ix+1,iy))+
     &                     ta*sf1/(gf1**2+qbb)+qaa
               s1(ix,iy) = sf1*dsc+qe*
     &                     (ds1(ix)+s1(ix,iy+1)+s1l+s1(ix+1,iy))+qd*gf1
               s2(ix,iy) = sf2*dsc+qe*
     &                     (ds2(ix)+s2(ix,iy+1)+s2l+s2(ix+1,iy))+qd*gf2

               g1l     = gf1
               g2l     = gf2
               dg1(ix) = gf1
               dg2(ix) = gf2
               s1l     = sf1
               s2l     = sf2
               ds1(ix) = sf1
               ds2(ix) = sf2
            enddo
         enddo
      enddo
      enddo

cccc  autocatalysis realised by inhibition of inhibition

      
C$$$      do iy = 1, igridy
C$$$         print '(31(g14.5,2x))', (s1(ix,iy), ix=1, igridx)
C$$$      enddo
      do iy = 1, ky
         do ix = 1, kx
            print *, s1(ix,iy)
         enddo
         print *
      enddo
      
      end

Pattern

November 8, 2006

Notes from TV Interview with Clive James

Filed under: misc — rcjp @ 4:37 pm

Its rare enough that I find anything worth watching on tv; even rarer to find something that I’d watch twice and positively nonexistent that I’d watch it a third time and scribble down some notes. Well, almost nonexistent, because here are some rough notes on a interview Clive James gave recently on TV on his life:

Born 1939, left Australia at 21, went to Pembroke Oxford and wrote a humorous tv column for the Observer in the 70’s, poems and a three volume autobiography.

Defining moment was his fathers death; he had survived the slave labour camps in Japan but died when, instead of waiting for the ship, he boarded an American plane laid on to take them back home to Australia which flew into a typhoon.

Saw Japanese ‘Endurance’ tv programme as a way to move the violent militaristic culture ‘Bushido’ onto television where it belongs.

Over-attached to his mother, but left for 16 years. “I didn’t get a very amazing degree from Sydney University and later on I didn’t get a very amazing degree from Cambridge. I knew some brilliant students and I wasn’t one of them – I spent too much time reading off the course, and I still do – its an inability to do what is set down.”

Terrible first jobs. Poems came back for two years. His old professor got him into Cambridge. “Very little happened because I willed it. On the whole I have to be asked to do things.”

Liked being a tv critic because it gave him a chance to talk about everything – since there is a programme about every subject. Critics these days get sent a dvd with a week’s tv, he used to watch the channels simultaneously (three) and scribble furiously as there was no way to wind back! Always thinks of pieces in terms of a running order – as he did when at Cambridge foot-lights. He’d take his notes from a week of watching tv and sit down at someone else’s typewriter at the Observer on a Friday morning and type.

“Whenever I thought of a joke, it was usually then, it was usually compressing an argument I already had – what’s the fastest way of saying this and I’ll think of a way that came out funny. And then sit back and laugh at it …”

This became legendary i.e. that he would sit there just laughing.

“…that was possibly a mistake – I was laughing because it had just hit me like a lightening bolt. All the other journalists which stretched into the distance, hated me because they already thought I only worked one day a week and there I was yodeling at my own stuff.”

Infamously described Arnold Schwargzenegar – Brown condom full of walnuts
and a tennis player as ‘a smile like a car crash’ (she had braces).

Some quotes (often abridged but hopefully representative) from the interview:

Ideally I’d like to make every answer an epigram. I was always trying to write the ideal paragraph; and sometimes the ideal paragraph is only a line long. I use a joke as a means of compression, its the fastest way of conveying an argument. But I always did love the tradition of the aphorism.

Some of my favourite things were said anonymously in Hollywood. Someone one said about a young actress “She’d be an nymphomaniac if only they could slow her down”

Rhyming slang I think is the dullest thing that ever happened to the world.

Get the joke in first – no one can say anything about my personal appearance because I got there first. Get your self deprecation in very early before the deprecation starts.

Interviews – Sit there mute, but put something on the desk.

It takes 10 years to build an image on Fleet St, I started finding a 1000 word column restrictive.

I will never stop writing poetry; thats a real test of seriousness because nobody really wants it. Hardly anybody reads poetry habitually, I do, and I like the people that do. But its a small market.

A woman’s world is what I look forward to.

The big lesson I’ve learned – A free society is simply bound to be full of things you don’t like – thats what freedom is.

He has a website http://clivejames.com

November 7, 2006

Series Solution Method in maxima

Filed under: physics — rcjp @ 5:56 pm

Just a small script to illustrate series solution

/***************************************************************************
  Solve the Hamiltonian

                       2  2    2
                    m q  w    p
          H(p, q) = ------- + ---
                       2      2 m

  by splitting to two coupled equations

                       d           p
                       -- (q(t)) = -
                       dt          m

                    d                  2
                    -- (p(t)) = - m q w
                    dt

  and using series method
****************************************************************************/
kill(p,q,t,m,w)$
eq1: diff(q(t), t)= p(t)/m;
eq2: diff(p(t), t)= -m*w^2*q(t);
assume(w>0)$
atvalue(p(t), t=0, 0)$
atvalue(q(t), t=0, 1)$
ans: desolve([eq1, eq2], [q(t), p(t)]);

/* And now do it just using power series */

kill(a,b,p,q,t,w,m)$
p:taylor(sum(a[i]*t^i,i,0,6),t,0,6)$
q:taylor(sum(b[i]*t^i,i,0,6),t,0,6)$

/* Note rearrange eqn=0 */

eq1:diff(q,t)-p/m$
eq2:diff(p,t)+m*w^2*q$
coeffs1:makelist(ratcoef(eq1,t,i),i,0,6)$
coeffs2:makelist(ratcoef(eq2,t,i),i,0,6)$
eqns:append(cons(subst(0,t,p) = 0,coeffs1),cons(subst(0,t,q) = 1,coeffs2))$
ab:append(makelist(a[i],i,0,6),makelist(b[i],i,0,6))$
ans:solve(eqns,ab)$
trunc(subst(ans[1],p));
trunc(subst(ans[1],q));

or inside an emaxima session:

% -*-EMaxima-*-
\documentclass{article}
\usepackage{a4wide}
\usepackage{emaxima}
\begin{document}
\section*{Series Solution of PDE}
$$
H(p,q) = {p^2\over 2 m} + {m\over 2} w^2 q^2
$$

$$
  \frac{dq}{dt}=\frac{\partial H}{\partial p} = \frac{p}{m} \qquad{\rm and}\qquad
  \frac{dp}{dt}=-\frac{\partial H}{\partial q} = -m\omega^2q
$$

\begin{maximasession}
  kill(p,q,t,m,w)$
  eq1: diff(q(t), t)= p(t)/m;
  eq2: diff(p(t), t)= -m*w^2*q(t);
  
  assume(w>0)$
  atvalue(p(t), t=0, 0)$
  atvalue(q(t), t=0, 1)$
  
  ans: desolve([eq1, eq2], [q(t), p(t)]);
\maximaoutput*
\end{maximasession}

To do that with series expansions only:

\begin{maximasession}
  p: taylor(sum(a[i]*t^i,i,0,inf),t,0,6)$ 
  q: taylor(sum(b[i]*t^i,i,0,inf),t,0,6)$ 

  eq1: diff(q, t) - p/m$
  eq2: diff(p, t) + m*w^2*q$

  coeffs1: makelist (ratcoeff(eq1, t, i), i , 0 , 6)$
  coeffs2: makelist (ratcoeff(eq2, t, i), i , 0 , 6)$

  /* add in boundary conditions p=0 and q=1 at t=0 */

  eqns: append(cons(subst(0,t,p)=0, coeffs1), cons(subst(0,t,q)=1, coeffs2))$
  coeffs:append(makelist(a[i],i,0,6), makelist(b[i],i,0,6))$
  ans: solve(eqns, coeffs)$

  trunc(subst(ans[1], p));
  trunc(subst(ans[1], q));

\maximaoutput*
\end{maximasession}
\end{document}

November 4, 2006

emaxima

Filed under: physics — rcjp @ 4:21 pm

Usually I run maxima just from the command line as rmaxima but occasionally it is useful to run it inside emacs in combination with Tex.

My .emacs has

;;;
;;; Maxima
;;;
(add-to-list 'load-path "/home/r/src/lisp/maxima/interfaces/emacs/emaxima")
(autoload 'maxima "maxima" "Running Maxima interactively" t)
(autoload 'maxima-mode "maxima" "Maxima editing mode" t)
(autoload 'emaxima-mode "emaxima" "eMaxima" t)

(add-hook 'inferior-maxima-mode-hook
          (lambda () (local-set-key [tab] 'inferior-maxima-complete)))
(add-hook 'inferior-maxima-mode-hook
          (lambda () (set (make-local-variable 'comint-prompt-read-only) t)))

and an example of a tex file that includes some maxima workings…

    % -*-EMaxima-*-
    \documentclass{article}
    \usepackage{a4wide}  % bit more room!
    \usepackage{amsmath} % used to get e.g. \boxed{...}
    \usepackage[breqn]{emaxima}
    \begin{document}
    \section*{Laplace Transform for Solving Integral Equation}

    % C-c +           to go to next cell.
    % C-u C-c C-u a   updates all the cells without asking - maxima output
    % C-u C-c C-u A   updates all the cells without asking - TeX output

    Solve $$\int_0^t cosh(a x) f(t-x)\,dx + b f(t)=t^3$$

    \begin{maximasession}
      'integrate(cosh(a*x)*f(t-x), x, 0, t) + b*f(t)=t**3;
      laplace(%,t,s);
      res: linsolve([%],['laplace(f(t),t,s)]);
    \end{maximasession}

    Use {\tt dpart} to put a box around which part you will get
    when using {\tt part}

    % use C-c C-u C  to get the tex output for just this one
    % C-c C-d        to delete the output

    \begin{maxima}[laplace1]
      dpart(res[1], 2);
      ilt(part(res[1], 2),s,t);
    \end{maxima}

    \end{document}

and a slightly more complex example…

    % -*-EMaxima-*-
    \documentclass{article}
    %\usepackage{a4wide}  % bit more room!
    \usepackage{amsmath} % used to get e.g. \boxed{...}
    \usepackage[breqn]{emaxima}
    \usepackage{vpage}
    \setpapersize{Afour}
    \setmarginsrb{1in}{1in}{1in}{1in}{0pt}{0mm}{0pt}{0mm}

    \begin{document}
    \section*{Perturbed DE Solution}
    Solve $u''+u=\epsilon u^2$ with $u(0)=1/r$ and $u'(0) = 0$ by looking
    for solutions of the form
    $$
    u(t)=u_0(t)+\epsilon u_1(t)+\epsilon^2u_2(t)+\ldots
    $$
    and initial conditions $u_0 = 1/r$, $u'_0=0$,
    $u_1=0$, $u'_1=0$, $u_2=0$, $u'_2=0$ at $t=0$

    \begin{maximasession}
      eqn_init: 'diff(u,t,2)+u=epsilon*u^2;
      v: u0 + epsilon*u1 + epsilon^2*u2;
      depends([u0,u1,u2],t);
      eq: subst(v,u,eqn_init);
      eq: ev(eq, diff);
      eq: ratexpand(eq);
      equ0: ev(eq, epsilon=0);
      solveu0: ode2(equ0, u0, t);
      solveu0: ic2(solveu0, t=0, u0=1/r, 'diff(u0,t)=0);
      uu0: rhs(%)$
      eq1: coeff(eq, epsilon, 1);
      eq1: ev(eq1, u0=uu0);
      solveu1: ode2(eq1, u1, t);
      solveu1: ic2(solveu1, t=0, u1=0, 'diff(u1,t)=0);
      uu1: rhs(%)$
      eq2: coeff(eq, epsilon, 2);
      eq2: ev(eq2, u0=uu0, u1=uu1);
      eq2: ratexpand(eq2);
      solveu2: ode2(eq2, u2, t);
      solveu2: ic2(solveu2, t=0, u2=0, 'diff(u2,t)=0);
      uu2: rhs(%)$
      u_final: ev(v, u0=uu0, u1=uu1, u2=uu2);
      trunc(%);
    \end{maximasession}
    \end{document}

October 21, 2006

Autogenerated Poetry

Filed under: python — rcjp @ 7:23 pm

#
# Sucks in a story like Milton's Paradise Lost and generates text
# using the word probabilty from the sample document
# translated from Common Lisp Graham p.140
#
import re
import random
import textwrap

words = {}

def read_text(filename):
    """Fill the words dictionary with unique words from 'filename',
    each word entry is itself a dictionary of the words that follow
    that word and their frequency"""

    sampletext = open(filename, 'r').read()
    previous = '.'  #  i.e. the start of a sentence
    for wordpunc in sampletext.lower().split():
        # keep the punctuation as part of the word e.g. 'bye.'
        # now split 'bye.' as separate words 'bye' and '.'  note (\W)
        # returns the splitting characters as well as the fields so we
        # get 'bye' '.' '' and we want to ignore the last seperator 

        # [Note normally we'd want to avoid punctuation so should do   
        #  words=re.compile(r'[\w'-]+') then  
        #  for word in words.finditer(line)
        #  do something to word.group(0) ]
        for word in [w for w in re.split(r'(\W)', wordpunc) if w != '']:
            if previous not in words:
                words[previous] = {word : 1}
            else:
                words[previous][word] = words[previous].get(word, 0) + 1
            previous = word
    print 'Using a vocabulary of', len(words), 'words'

def format_word(word, previous):
    if previous == '.':
        word = ' ' + word.capitalize()
    else:
        if word.isalpha(): word = ' ' + word
    return word

def generate_text(nwords, previous = '.'):
    """Prints 'nwords' random words chosen according to statistically
    likely order calculated in read_text"""

    text = ''
    for n in xrange(nwords):
        nextwords = words[previous]
        count = 0
        pick = random.randint(0, sum(nextwords.values()))
        for word, freq in nextwords.iteritems():
            count += freq
            if count >= pick:
                text += format_word(word, previous)
                break
        previous = word
    return text

if __name__ == '__main__':
    read_text('c:/tmp/testwords')
    print textwrap.fill(generate_text(100), 50)

"""
e.g.
In [210]: read_text('c:/tmp/ParadiseLost.txt')
Using a vocabulary of 8993 words

In [211]: print textwrap.fill(generate_text(100),50)
 Farewell, immutable; i give not: conviction to
torment me, and worse. Satan only disagree of
servants feet. To that i boast what was old night,
metals of belial came; lest the government well
thou for which we seek some glade, or who wrong,
but he more soft and evil ruin. So erroneous there
dwell permits, and inward faculties, and bliss),
or, to know. With design to spend, till then to
the fields, when he ended; of nature him hither
thrust
"""

October 11, 2006

Find classes in java jar files

Filed under: lisp, utils — rcjp @ 1:57 pm

;;;
;;; Find classes in jar files
;;; Usage (findjar classname directory)
;;; e.g (findjar "session" "c:/Novell/ndk/njclv2/lib")
;;;
(defun findjar (classname dir)
  ;;  (if (not (directory dir))
  ;;(format t "Invalid directory syntax - remember a trailing slash")
  (dolist (file (directory (concatenate 'string dir "/**/*")))
    (when (string-equal (pathname-type file) "JAR")
      (format t "processing ~S~%" file)
      (with-open-stream (str
        #+cmu
        (run-program (concatenate 'string "jar tf " (namestring file)) :output :stream)
        #+ clisp
        (ext:run-program "jar" :arguments (list "tf" (namestring file)) :output :stream)
        #+ lispworks
        (sys:open-pipe (concatenate 'string "jar tf " (namestring file)))
                          )
        (loop for line = (read-line str nil nil)
              while line
              do
              (when (search classname line :test 'string-equal)
                (format t "   found->~S in file (~S)~%" line (namestring file))
                (return)))))))

Nonograms in Common Lisp

Filed under: lisp — rcjp @ 1:47 pm

;;;
;;; Nonograms are puzzles from Sunday Telegraph: the length and order
;;; of blocks in a grid are given along both rows and columns the
;;; challenge is to work out how many spaces occur between the blocks
;;; and find the hidden picture (see example below)
;;;

;;; the strategy we use is to generate all the possible spacings
;;; of blocks and see if any parts of the grid, for all possible
;;; solutions, are always filled '1' or always empty '0'

(defun common-bits (bitvecs)
  "Returns a bit vector with '1' in positions where all supplied
bitvectors all have 1 or 0"
  (let ((firstvec (first bitvecs)))
    (reduce #'bit-and (loop for bv in bitvecs collect (bit-eqv firstvec bv)))))

(defun common-list-positions (lists)
  "Returns a list of (value . position) elements common to all lists"
  (let ((bitvecs (mapcar #'(lambda (list) (coerce list 'bit-vector)) lists)))
    (loop with common = (common-bits bitvecs)
       with vec = (first bitvecs)
       for start = 0 then (1+ pos)
       for pos = (position 1 common :start start)
       while pos collect (cons (sbit vec pos) pos))))

(defun make-section (lhs lspace blocksize rhs &optional not-end-section)
  "Make a section of a row (also used for columns), the end piece has
the rhs full of spaces else just pad a space ready for the next section"
  (append lhs
          (make-list lspace :initial-element 0)
          (make-list blocksize :initial-element 1)
          (if not-end-section
              '(0)
              (make-list rhs :initial-element 0))))

(defun blocks-spacings (blocklist nspaces &optional prevblocks)
  "Calculate all the arrangements of the list of blocks over nspaces"
  (loop
     with blocklen = (first blocklist)
     with moreblocks = (rest blocklist)
     with maxleft = (- nspaces (+ (apply #'+ moreblocks) (length moreblocks)))
     ;;
     ;; working left-to-right, work out all possible positions for the
     ;; current block recursively calling for blocks to the rhs
     ;;
     for blockpos upto (- maxleft blocklen)
     for section = (make-section prevblocks blockpos blocklen
                                 (- maxleft blocklen blockpos) moreblocks)
     ;;
     ;; build the row (column) in sections, appending sections to the
     ;; collect'ed right hand edge sections
     ;;
     when moreblocks append
       (blocks-spacings moreblocks (- nspaces blockpos blocklen 1) section)
     else collect section))

(defun find-spacings (blocklists nspaces)
  "Generate a vector from the blocklists where each element holds
all possible arrangements of each blocklist in nspaces"
  (let ((board (make-array (length blocklists) :fill-pointer 0)))
    (dolist (blocklist blocklists board)
      (vector-push (blocks-spacings blocklist nspaces) board))))

(defun multiple-solutions-p (boardx boardy)
  "If there is more than one possible arrangement anywhere - its not solved"
  (flet ((length1-p (seq) (= (length seq) 1)))
    (or (notevery #'length1-p boardx)
        (notevery #'length1-p boardy))))

(defun print-board (board)
  (loop for row across board
     do (format t "~{~&~{~[ ~;o~]~}~}" row)))

(defun delete-arrangement (array element test)
  "Delete from the list of arrangements in array element those
that fail the test"
  (setf (aref array element) (delete-if-not test (aref array element))))

(defun eltcheck (index val)
  "Returns a function which can be applied to a sequence to
check if the index value is val"
  #'(lambda (n) (= (elt n index) val)))

(defun nonogram (row-blocks col-blocks &optional (maxiterations 100))
  "Solve the nongram described by the list of row-blocks and
col-blocks, give up after maxiterations"
  (let ((possible-rows (find-spacings row-blocks (length col-blocks)))
        (possible-cols (find-spacings col-blocks (length row-blocks))))
    (flet ((find-apply-rules (board-find board-apply)
             (dotimes (line (length board-find))
               (dolist (rule (common-list-positions (aref board-find line)))
                 (delete-arrangement board-apply (cdr rule) (eltcheck line (car rule)))))))
      (loop while (multiple-solutions-p possible-rows possible-cols)
         repeat maxiterations
         do
           (find-apply-rules possible-rows possible-cols)
           (find-apply-rules possible-cols possible-rows))
      (print-board possible-rows))))

#|
e.g.
(nonogram '((5) (4) (3 3 3) (8 1) (7 3 3) (2 2 5 1) (1 6 3 1) (3 3 2) (3 2 1) 
            (2 2) (1 2) (1 2) (2) (2) (2) (2) (8) (11) (13) (4 4 5 3) (3 2 3 3) 
            (2 3) (20) (14 1 2) (4 1 12))
          '((6) (6) (2 3) (3 1 3) (3 1 2) (3 5 3 2) (3 4 4 3) (7 5 2) (4 4 3) 
            (4 2 3 3) (4 14 3) (20 3) (2 2 5 3) (1 3 5 3) (1 3 5 1 1) (1 3 3 1 1) 
            (1 1 3) (3 4 1) (1 6) (3 6)))
=>
          ooooo     
         oooo       
     ooo ooo   ooo  
    oooooooo     o  
   ooooooo ooo   ooo
   oo  oo ooooo    o
   o  oooooo ooo   o
     ooo ooo  oo    
     ooo  oo   o    
     oo   oo        
     o    oo        
     o    oo        
          oo        
          oo        
          oo        
          oo        
       oooooooo     
     ooooooooooo    
    ooooooooooooo   
oooo oooo  ooooo ooo
ooo   oo    ooo  ooo
oo               ooo
oooooooooooooooooooo
oooooooooooooo  o oo
oooo  o oooooooooooo

|#

« Newer PostsOlder Posts »

Blog at WordPress.com.