Archive for April 2010

IFS for the TI-89 Titanium   Leave a comment

Being bored this morning, I decided to implement the ifs algorithm on my TI-89 graphics calculator. It was not particularly difficult, and after 23 lines of TI-BASIC, I had this image on my screen.

Yes, rather low res, I know...

The program takes a n*7 matrix, like this

Where a, b, c, and d are the elements of a 2*2 transformation matrix, e and f are elements of an offset vector, and p is the probability of the rule occuring. Note that it runs extremely slowly, I waited about ten minutes for the picture shown to be generated. I’ve added it to the repo though; in case anyone else is also bored ;)

Prgm
rand()->a
rand()->b
0->i
1000000->iters
100->s
While ic
rand()->l
While pc
r[c,7]+p->p
EndWhile
r[c,1]*a+r[c,2]*b+r[c,5]->t
r[c,3]*a+r[c,4]*b+r[c,6]->b
t->a
If i>s
PtOn a,b
1+i->i
EndWhile
EndPrgm

Posted 30/04/2010 by Emmanuel Jacyna in Humor

Tagged with , , , ,

SICP – Section 2.2.4   Leave a comment

This section is a powerful demonstration of abstraction, using a “picture language”. For this section I used the sicp package from the PLaneT repository in DrScheme. Normally, I just use mzscheme and vim, but apparently you need DrScheme open to use that particular (graphical) package. Anyway, I found it pretty enjoyable, and even made some fractals with it.

;;Exercise 2.44

They could have made this one a little tougher ;)
Telling us to flip the arguments to right-split took any challenge out of it…

(define (up-split painter n)
  (if (= n 0)
      painter
      (let ((smaller (up-split painter (- n 1))))
        (below painter (beside smaller smaller)))))

;;Exercise 2.45

We follow the usual patter of making abstractions.

(define (split op1 op2)
  (define (rec painter n)
    (if (= n 0)
        painter
        (let ((smaller (rec painter (- n 1))))
          (op1 painter (op2 smaller smaller)))))
  rec)

(define up-split (split below beside))
(define right-split (split beside below))

;;Exercise 2.46

Not much to it, I used a list, but a cons would do fine too, however, I decided on list, because perhaps we’ll want to add a z coordinate for 3D vectors.

(define (make-myvect x y)
  (list x y))
(define (xcor-myvect vect)
  (car vect))
(define (ycor-myvect vect)
  (cadr vect))

(define (add-myvect a b)
  (make-myvect (+ (xcor-myvect a) (xcor-myvect b))
             (+ (ycor-myvect a) (ycor-myvect b))))

(define (sub-myvect a b)
  (make-myvect (- (xcor-myvect a) (xcor-myvect b))
             (- (ycor-myvect a) (ycor-myvect b))))

(define (scale-myvect s v)
  (make-myvect (* s (xcor-myvect v))
             (* s (ycor-myvect v))))

;;Exercise 2.47

First definition

(define (make-myframe origin edge1 edge2)
  (list origin edge1 edge2))

(define (origin-myframe f) (car f))

(define (edge1-myframe f) (cadr f))

(define (edge2-myframe f) (caddr f))

Second definition

(define (make-myframe origin edge1 edge2)
  (cons origin (cons edge1 edge2)))

(define (origin-myframe f) (car f))

(define (edge1-myframe f) (cadr f))

(define (edge2-myframe f) (cddr f))

;;Exercise 2.48

(define (make-mysegment start end) (cons start end))

(define (start-mysegment seg) (car seg))

(define (end-mysegment seg) (cdr seg))

;;Exercise 2.49


;;a.  The painter that draws the outline of the designated frame.
(define outline (segments->painter (list
  (make-segment (make-vect 0 0)
                (make-vect 0 1))
  (make-segment (make-vect 0 0)
                (make-vect 1 0))
  (make-segment (make-vect .99 0)
                (make-vect .99 1))
  (make-segment (make-vect 0 .99)
                (make-vect 1 .99)))))

;;b.  The painter that draws an ``X'' by connecting opposite corners of the frame.
(define X (segments->painter (list
  (make-segment (make-vect 0 0)
                (make-vect 1 1))
  (make-segment (make-vect 1 0)
                (make-vect 0 1)))))

;;c.  The painter that draws a diamond shape by connecting the midpoints of the sides of the frame.
(define diamond (segments->painter (list
  (make-segment (make-vect 0 .5)
                (make-vect .5 0))
  (make-segment (make-vect .5 0)
                (make-vect 1 .5))
  (make-segment (make-vect 1 .5)
                (make-vect .5 1))
  (make-segment (make-vect .5 1)
                (make-vect 0 .5)))))

;;d.  The wave painter.
(define wave (segments->painter (list
  (make-segment (make-vect 0 .84) (make-vect .13 .6))
  (make-segment (make-vect .13 .6) (make-vect .3 .66))
  (make-segment (make-vect .3 .66) (make-vect .4 .66))
  (make-segment (make-vect .4 .66) (make-vect .36 .84))
  (make-segment (make-vect .36 .84) (make-vect .4 1))
  (make-segment (make-vect .6 1) (make-vect .64 .84))
  (make-segment (make-vect .64 .84) (make-vect .6 .66))
  (make-segment (make-vect .6 .66) (make-vect .71 .66))
  (make-segment (make-vect .71 .66) (make-vect 1 .34))
  (make-segment (make-vect 1 .16) (make-vect .6 .44))
  (make-segment (make-vect .6 .44) (make-vect .74 0))
  (make-segment (make-vect .6 0) (make-vect .5 .3))
  (make-segment (make-vect .5 .3) (make-vect .4 0))
  (make-segment (make-vect .24 0) (make-vect .34 .5))
  (make-segment (make-vect .34 .5) (make-vect .3 .6))
  (make-segment (make-vect .3 .6) (make-vect .14 .4))
  (make-segment (make-vect .14 .4) (make-vect 0 .64)))))

;;Exercise 2.50
The key to figuring these out is to draw a diagram of the unit square, draw the original picture, and then draw the rotated one, and note the coordinates of the origin, edge1, and edge2.

(define (myflip-horiz painter)
  ((transform-painter (make-vect 1 0)
                      (make-vect 0 0)
                      (make-vect 1 1)) painter))

(define (rotate-90 painter)
  ((transform-painter (make-vect 1 0)
                     (make-vect 1 1)
                     (make-vect 0 0)) painter))

(define (rotate-180 painter)
  ((transform-painter (make-vect 1 1)
                      (make-vect 0 1)
                      (make-vect 1 0)) painter))

(define (rotate-270 painter)
  ((transform-painter (make-vect 0 1)
                      (make-vect 0 0)
                      (make-vect 1 1)) painter))

;;Exercise 2.51

(define (my-below painter1 painter2)
  (let ((split-point (make-vect 0 .5)))
    (let ((paint-top
           ((transform-painter split-point
                               (make-vect 1 .5)
                               (make-vect 0 1)) painter1))
          (paint-bottom
           ((transform-painter (make-vect 0 0)
                               (make-vect 1 0)
                               split-point) painter2)))
      (lambda (frame)
        (paint-top frame)
        (paint-bottom frame)))))

(define (my-below2 painter1 painter2)
  (rotate-90 (beside (rotate-270 painter1) (rotate-270 painter2))))

Okay, this has been a pretty boring article… lack of pictures, &c. However, I do have a couple pictures of some nice fractals that I generated with some functions I wrote in this picture language.

The sierpinski painter generates a sierpinski triangle with n iterations using an MRCM (Multiple Reduction Copy Machine) process. This is based on the idea from the book, Computational Beauty of Nature.

(define (sierpinski painter n)
  (if (= n 0)
      painter
            (let ((top ((transform-painter (make-vect .25 0)
                      (make-vect .75 0)
                      (make-vect .25 1))
                      painter)))
        (sierpinski (below (beside painter painter) top) (- n 1))))) 

And here it is generating it with the einstein painter.

There’s also another neat fractal, which is defined as

(define (crystal painter n)
  (if (= n 0)
      painter
            (let ((top ((transform-painter (make-vect .25 0)
                      (make-vect .75 0)
                      (make-vect .25 1))
                      painter)))
        (crystal (below (beside (rotate-90 painter)
                                   (rotate-270 painter))
                           top)
                    (- n 1))))) 

And looks like

As you can see, it’s pretty easy to define MRCM fractals using this picture language.

Posted 24/04/2010 by Emmanuel Jacyna in Code, Scheme, SICP

Tagged with , , , ,

SICP – Section 2.2.3   Leave a comment

This section had just enough challenging problems, along with just enough hints to guide me along, to make it quite fun. I did this section about two weeks ago, but, due to school starting, and some other things in my life, it’s taken a bit longer to put it up. Also, I’ve put up my solutions in a github repo for people to check out, critique, and possibly correct. Anyway, without further ado, here are my solutions.

;;Exercise 2.34

The horner-eval exercise was fairly straightforward, you just have to remember that
(...(anx + an-1) x + ... + a1) x + a0
is equivalent to
a0 + x (a1 + x( ... an-1 + x (an))...)


(define (horner-eval x coefficient-sequence)
  (accumulate (lambda (this-coeff higher-terms)
                (+ this-coeff (* x higher-terms)))
              0
              coefficient-sequence))

;;Exercise 2.35

I felt a bit sneaky with my solution for this, my map was pretty much an identity function, it had no affect on the end result…

(define (count-leaves t)
  (accumulate (lambda (x y) (+ (if (not (pair? x)) 1 (count-leaves x)) y))
              0
              (map (lambda (x) x) t)))

;;Exercise 2.36

I had to define a couple auxiliary functions for this one, car-tree and cdr-tree take the car and cdr of the lists in a list… that is,
> (car-tree ((1 2 3) (2 3 4) (3 4 5)))
> (1 2 3)

(define (car-tree t)
  (map (lambda (l) (car l)) t))
(define (cdr-tree t)
  (map (lambda (l) (cdr l)) t))

(define (accumulate-n op init seqs)
  (if (null? (car seqs))
      '()
      (cons (accumulate op init  (map (lambda (l) (car l)) seqs))
            (accumulate-n op init (map (lambda (l) (cdr l)) seqs)))))

;;Exercise 2.37
I haven’t done much matrix and vector math, so some of my definitions may be incorrect, however, I’m pretty sure they work like they’re meant to. I’m also pretty impressed with the way these functions work so nicely. I also added a bit of error/dimensions checking, so you can’t multiply a 2×3 by a 1×4 matrix, etc.

(define (dot-product v w)
  (accumulate + 0 (map * v w)))

(define (element r c)
  (accumulate + 0 (accumulate-n * 1 (list r c))))

(define (matrix-*-vector m v)
  (map (lambda (r) (element r v)) m))

(define (dimensions m)
  (cons (length m) (length (car m))))

(define (transpose mat)
  (accumulate-n (lambda (x y) (cons x y)) '() mat))

(define (matrix-*-matrix m n)
  (if (not (= (cdr (dimensions m)) (car (dimensions n))))
      (error "Dimension mismatch")
      (let ((cols (transpose n)))
      (map (lambda (row)
             (map (lambda (col)
                    (element col row))
                  cols))
           m))))

;;Exercise 2.38

> (fold-right / 1 (list 1 2 3))
> 3/2
> (fold-left / 1 (list 1 2 3))
> 1/6

As we can see, we get different answers, the property is fairly obvious, op must have the associative property.

;;Exercise 2.39

(define (right-reverse sequence)
  (fold-right (lambda (x y) (append y (list x))) '() sequence))

(define (left-reverse sequence)
  (fold-left (lambda (x y) (append (list y) x)) '() sequence))

Just have to flip the order in either one.

;;Exercise 2.40

(define (unique-pairs n)
  (filter (lambda (i) (not (eqv? i '())))
          (flatmap (lambda (i)
                   (map (lambda (j) (if (< j i) (list i j) '()))
                        (enumerate-interval 1 (- n 1))))
                   (enumerate-interval 1 n))))

(define (prime-sum-pairs n)
  (map make-pair-sum
       (filter prime-sum? (unique-pairs n))))

This does indeed greatly simplify the definition of prime-sum-pairs, from 8 to 3 lines.

;;Exercise 2.41

I found this one a little a little tricky, but after studying the chapter again, I got it. The key is the nested mapping, and the use of flatmap, also note the abstraction in in-order?

(define (triples n)
  (flatmap (lambda (i)
             (flatmap (lambda (j)
                    (map (lambda (k)
                           (list i j k))
                         (enumerate-interval 1 n)))
                  (enumerate-interval 1 n)))
            (enumerate-interval 1 n)))

(define (in-order? cmp l)
  (cmp (car l) (cadr l) (caddr l)))

(define (gt? l)
  (in-order? > l))

(define (ordered-triples n)
  (filter gt? (triples n)))

(define (sum-to? triple s)
  (= (in-order? + triple) s))

(define (ordered-triples-which-sum-to-s n s)
  (filter (lambda (t) (sum-to? t s)) (ordered-triples n)))

;;Exercise 2.42

I haven’t completed this one yet, I think I have the general idea, but there are some issues with my implementation.

Posted 24/04/2010 by Emmanuel Jacyna in Code, Scheme, SICP

Tagged with , ,

Using python generators   Leave a comment

I just made a generator in python.

Wow. Just, Wow.

Generators are really awesome.
I’ve been having a couple issues with my IFS algorithm, because it’s not quite as fast as I would like, and I can’t seem to do everything I want just by passing callback functions. Here’s the original algorithm:

def run_IFS(iters, skip, rules, draw_method, speak_method=None, interval = 10000):
    """Calculate the points in an iterated fractal system, drawing them with draw_method,
    and announcing results with speak_method every interval iterations."""
    probs = [r[PROB_I] for r in rules]
    rules = [r[RULE_I] for r in rules]

    x = random_float()
    y = random_float()

    while i < iters + skip
        chosen = choose_rule(probs)
        x, y = calculate_transform(x, y, rules[chosen])
        if i > skip:
            draw_method(x, y, chosen)
        #Announce the results
        if not speak_method: pass
        elif i%interval == 0:
            speak_method(x, y, i)
        i += 1

This isn’t bad code, and it’s not particularly difficult to use. I know, because it was very easy for me to add a tkinter interface in less than ten minutes. However, when I wanted to write a function to calculate the best size settings for a fractal, I found it to be rather messy:


def old_calculate_best_settings(rules, scale, iters = 10000, skip = 100):
    """Return a tuple of width and height and x_off and y_off settings
((w,h),(x_off, y_off))"""
    #note, I know this is ugly, I don't know how else to do it though
    global h_x, h_y, l_x, l_y
    h_x = h_y = l_x = l_y = 0

    def draw_method(x, y, chosen):
        global h_x, h_y, l_x, l_y
        if x*scale > h_x: h_x = x*scale
        elif x*scale < l_x: l_x = x*scale
        if y*scale > h_y: h_y = y*scale
        elif y*scale < l_y: l_y = y*scale

    run_IFS(iters, skip, rules, draw_method)
    width = int(abs(l_x)+abs(h_x))
    height = int(abs(l_y)+abs(h_y))
    x_off = int(abs(l_x))
    y_off = int(abs(l_y))

    return ((width, height), (x_off, y_off))

Yes, that’s right.

    Global Variables


Because of an interesting feature in python, it’s not possible for draw_method to change the values of h_x (and friends) without resorting to globals. This is because you can access variables from outside scope, but assigning to them makes a new local variable with the same name, which is not really what we want (thanks to bob2 from #python on freenode.net for the explanation). Of course, I could have made these variables classes, and then called some __set__ method to change their values. Purge your mind friend. The very fact that we’re thinking of such roundabout, convoluted methods is a sure sign that something is not being done right.

Enter generators

Python, of course, has a nice solution to this problem. We can see that this function is basically iterating and producing points. A possible solution would be to store the points in a list, and then iterate over these to see what the low/high points are. This is halfway there, however, when we’re talking millions of points, memory quickly becomes an issue (at least on my 486 :). Instead, we can use a generator, which behaves like a list, and makes for more “pythonic” code. Here’s a simple example:


def generate_x():
    i = 0
    while i < 10:
        print i
        yield i
        i += 1

for x in generate_x():
    print x

When you run that code (or on codepad) you can see that the prints are interleaved in the generator and the for loop, demonstrating how it works.
To change a function to a generator, we just rearrange the function a bit and make it yield rather than return the points.
Here’s the run_IFS function as a generator:


def IFS_generator(iters, skip, rules):

    probs = [r[PROB_I] for r in rules]
    rules = [r[RULE_I] for r in rules]

    x = random_float()
    y = random_float()
    i = 0

    while i < iters + skip:         chosen = choose_rule(probs)         x, y, = calculate_transform(x, y, rules[chosen])         if i > skip:
            yield (x, y, chosen)
        i += 1

That’s a lot simpler, isn’t it? We can now modify the calculate_best_settings function to use the generator:


def calculate_best_settings(rules, scale, iters = 10000, skip = 100):
    """Return a tuple of width and height and x_off and y_off settings
((w,h),(x_off, y_off))"""
    h_x = l_x = h_y = l_y = 0

    for points in IFS_generator(iters, skip, rules):
        x, y, _ = points
        if (x*scale) > h_x: h_x = x*scale
        elif (x*scale) < l_x: l_x = x*scale         if (y*scale) > h_y: h_y = y*scale
        if (y*scale) < l_y: l_y = y*scale

    width = int(abs(l_x)+abs(h_x))
    height = int(abs(l_y)+abs(h_y))
    x_off = int(abs(l_x))
    y_off = int(abs(l_y))

    return ((width, height), (x_off, y_off))

See how this is now nice, pythonic code? No more messing with global variables and callbacks. Of course, there’s also a nice speed bonus in getting rid of those callbacks.

After timing both functions, I got these results:
Old:
[4.5442321300506592, 4.4281911849975586, 4.3002359867095947] min: 4.30023598671
New:
[5.3141071796417236, 4.0443418025970459, 3.9916889667510986] min: 3.99168896675

This is a pretty significant speed up for such a small change! Note that the calculate_best_settings function only runs ~10000 iterations. When you generate a fractal with 1000000 iterations this difference becomes even greater. The code gets a boost once we’re rid of draw_method and speak_method. Of course, the next step would be to implement this function as a C extension.

Further reading:
Pep 255 – Simple Generators

Posted 19/04/2010 by Emmanuel Jacyna in Code, Python

Tagged with , , ,

pyafg – transform editor   Leave a comment

I was at school on thursday, and had IT Software Development in the afternoon. Surprisingly, it turned out to be quite enjoyable. I tested out pyafg on the school computers running windows xp, and it worked quite nicely. I can now say that pyafg is cross platform! (w00t)

Since we were being “taught” how to use tkinter, I decided to hack up a little gui for pyafg using it. It’s a rather simple affair, but is both a testament to the power of tkinter, and the success of my abstraction in libifs. I plan to turn run_IFS into a generator, and generalize it even further. While I was working on the (simplistic) gui, I had the idea to make a transform editor for pyafg. I think it should be pretty easy, and for those who’ve seen transform_viewer.py in the git repo, you’ll know that I already have the means to do so :)

I’ve also updated the pyafg page into something looking a little more like a project page, and plan to add a TODO list to that. Transform Editor will be first on the list :)

Posted 17/04/2010 by Emmanuel Jacyna in pyafg

Tagged with , , , , ,

On Driving & Coding   4 comments

Recently I’ve had the opportunity to try out driving a manual car. It was a pretty big change at first, I spent an afternoon figuring out how to start it properly and change gears. It reminded me a bit of switching from something like python to C. First it feels very cumbersome and difficult, but once you’ve gotten used to it, it just feels like a different way of driving. You’ll still get caught up every now and then, like hitting an invisible clutch when you’re in an automatic, but once you learn it, you don’t forget it. In the same way, I think hacking in C is just a different way of accomplishing the same thing. I’m careful not to stretch the analogy too far, but I think it’s a nice one. :)

Posted 14/04/2010 by Emmanuel Jacyna in General Life

Tagged with , ,

A new project   Leave a comment

I’ve started work on a new project, pyafg (pyafg is an affine fractal generator). It’s a little python program that generates affine fractals using the iterated function system algorithm. I was inspired to start it whilst reading the computational beauty of nature; I’ll write a post up about my experiences hacking it later. For now, download it, play around with it, and enjoy!

(Note, I did only hack this up in 2 days. It’s nowhere near complete, but I’ve made a very good start to it. Expect more features!)

Computer Generated Fern

A Fern I generated

Posted 08/04/2010 by Emmanuel Jacyna in Code, pyafg, Python

Tagged with , , , ,