## Sunday, June 24, 2007

### Creating fractal images using C# 3.0 features

In this post I'm going to try to implement the basic Escape time algorithm for creating fractal images. My main goal is to use as many C# 3.0 features as possible .

This post is inspired by the very nice blog post A Ray Tracer in C#3.0.

The Escape time algorithm is very simple, but depending on its input it can generate very complex images.

Although not mandatory, the first thing to create is a class from complex numbers. The complex number class contains the following members:

``public class Complex{   public Complex(double real, double img)   public Complex Multiply(Complex c)   public Complex Add(Complex c)   public double Norm()   public double NormSquared()   public static Complex operator *(Complex c1,Complex c2)    public static Complex operator +(Complex c1, Complex c2)   public double Real   public double Img }``

Now we need something to convert from the image coordinate system to the real coordinate system. To do this the following function was created:

``public static Func<double, double> InterpFunc(double x1,double y1,double x2,double y2){   double m = (y2 - y1) / (x2 - x1);   double b = y1 - (m * x1);   return (x => m * x + b);}``

The `InterpFunc` method creates a function that maps values from one coordinate system to the other.

The code that creates the image is the following:

``void CreateFractal(double minX, double minY,double maxX, double maxY,int imageWidth, int imageHeight){   Func<double, double> xF = MathUtils.InterpFunc(0, minX, imageWidth, maxX);   Func<double, double> yF = MathUtils.InterpFunc(0, minY, imageHeight, maxY);   foreach (var p in from yi in Enumerable.Range(0, imageHeight)                     from xi in Enumerable.Range(0, imageWidth)                     select new                       {                          x = xi,                          y = yi,                          xD = xF(xi),                          yD = yF(yi)                       })   {      Complex p0 = new Complex(p.xD, p.yD);      Func<Complex, Complex> function = functionConstructor(p0);      int i = ApplyFunction(function, p0)              .TakeWhile(                (x, j) => j < maxIteration && x.NormSquared() < 4.0)              .Count();      HandlePixel(p.x, p.y, i);   }}``

Basically `select ` expression generates the coordinates for all the points in the image . The `functionConstructor` function is used to create the function that will be applied to generate the fractal. One example of this functions is the one used for the Mandelbrot fractal `f(x) = x^2 + p0` . The escape time is calculated by counting the number of elements in the generated sequence of function applications before the value escaped.

The `ApplyFunction` method is interesting since is the one that creates the sequence of recursive function applications required for this algorithm. The method looks like this:

``IEnumerable<Complex> ApplyFunction(Func<Complex, Complex> function,Complex initial){   Complex last = initial;   while (true)   {      last = function(last);      yield return last;   }}``

By running this algorithm using the Mandelbrot formula:

``Func<Complex,Func<Complex,Complex>> fc =       p0 => (c => c.Multiply(c).Add(p0));EscapeTimeFractal p = new EscapeTimeFractal(300, 300, fc);p.MinX = -0.03;p.MinY = 0.68;p.MaxX = 0.03;p.MaxY = 0.62;p.Run();``

The generated image is the following: By running it using the Szegedi Butterfly 1 formula:

``Func<Complex,Func<Complex,Complex>> fc =       p0 =>        (c => (new Complex(                  ((c.Img * c.Img) - Math.Sqrt(Math.Abs(c.Real))),                  ((c.Real * c.Real) - Math.Sqrt(Math.Abs(c.Img)))))                + p0);EscapeTimeFractal p = new EscapeTimeFractal(300, 300,fc);p.MinX = -2;p.MinY = 2;p.MaxX =2;p.MaxY = -2;p.MaxIteration = 127;p.OutputFileName = @"c:\temp\output.bmp";p.Run();``

The generated image is the following: The code for this experiment can be found here.

## Thursday, June 21, 2007

### Why AutoCAD uses Lisp as extension/scripting language

The following links provide historic information on AutoLISP:

1. Why Lisp? from The Autodesk File

2. The AutoLISP Wikipedia entry

3. Introduction to AutoLISP

I always wanted to know why Lisp was used as the scripting/extension language of AutoCAD. Mainly because Lisp may be perceived as a Computer Science only language.

## Wednesday, June 20, 2007

### Calling the Boo compiler from a Boo program

In this post I'm going to show a little demo of creating a piece of a Boo program on the fly, compile it and execute it as part of the main program.

The Boo language implementation comes with two very useful and interesting libraries: `Boo.Language.Parser.dll` and `Boo.Language.Compiler.dll`. These libraries give the developer access to Boo's parser and compiler for use within any Boo(or .NET) program. Also the AST classes(` Boo.Lang.Compiler.Ast`) and visitors(`Boo.Lang.Compiler.Ast.Visitors`) are available to create and manipulate pieces of programs .

The `ast-to-string.boo,ast-to-xml.boo,run-ast.boo` and ` run-ast-without-compiler.boo` examples included with the Boo distributions shows how to use both the AST and the compiler.

The little demo I worked on, is a program that shows the content of an AST tree using a Winforms TreeView control. In order create this little experiment I wanted to create a class that inherits from `Boo.Lang.Compiler.Ast.DepthFirstVisitor` that will create the TreeNode instances used to visualize the tree.

Creating a class that inherits from `Boo.Lang.Compiler.Ast.DepthFirstVisitor` is a tedious work since you have to create a method for each kind of AST node with the same body. Because of this the Ast classes and the compiler to create the class and compile it on the fly. Here's the code.

``import Systemimport System.Reflectionimport System.IOimport Boo.Lang.Compilerimport Boo.Lang.Compiler.Astimport Boo.Lang.Parserimport System.Windows.Forms from System.Windows.Forms// Get the type instance for DepthFirstVisitora = Assembly.Load("Boo.Lang.Compiler")t = a.GetType("Boo.Lang.Compiler.Ast.DepthFirstVisitor")// Create a class definition for a custom visitormyVisitor = ClassDefinition(Name:"DynamicVisitor")myVisitor.BaseTypes.Add(SimpleTypeReference("Boo.Lang.Compiler.Ast.DepthFirstVisitor"))// Add a new field for our node stackmyVisitor.Members.Add(Field(Name:"stck",                            Type:SimpleTypeReference("System.Collections.Stack"),                            Modifiers:TypeMemberModifiers.Public))// Create all 'Enter' methodsfor m as MethodInfo in [m for m in t.GetMethods() if m.Name.StartsWith("Enter")]:   nm = Method(Name: m.Name,  Modifiers: TypeMemberModifiers.Override,                            Body:Block(),ReturnType:SimpleTypeReference("System.Boolean"))   nm.Parameters.Add(ParameterDeclaration(Name:"p",                                         Type:SimpleTypeReference(                                            Name:m.GetParameters().ParameterType.FullName)))   // Node for TreeNode('')   tnCreation = MethodInvocationExpression(ReferenceExpression("TreeNode"))                                     tnCreation.Arguments.Add(            StringLiteralExpression(m.GetParameters().ParameterType.Name))      // Create node for 't = TreeNode('')'    decl = DeclarationStatement(           Declaration: Declaration(Name:"t",                      Type: SimpleTypeReference("System.Windows.Forms.TreeNode")),           Initializer: tnCreation)      nm.Body.Add(decl)      // Create nodes for '(stck.Peek() as TreeNode).Nodes.Add(t)'   at =  MethodInvocationExpression(           MemberReferenceExpression(              Target:MemberReferenceExpression(                 Target: TryCastExpression(                          Target:MethodInvocationExpression(                                MemberReferenceExpression(                                     Target:ReferenceExpression("stck"),                                     Name:"Peek")),                          Type:SimpleTypeReference("System.Windows.Forms.TreeNode")),                 Name:"Nodes"),              Name:"Add"))   at.Arguments.Add(ReferenceExpression("t"))      nm.Body.Add(ExpressionStatement(at))      // Create node for 'stck.Push(t)'   mie = ExpressionStatement(              MethodInvocationExpression(                         MemberReferenceExpression(                               Target:ReferenceExpression("stck"),Name:"Push")))                                                (mie.Expression as MethodInvocationExpression).Arguments.Add(ReferenceExpression("t"))         nm.Body.Add(mie)      // Add 'return true' statement   nm.Body.Add(ReturnStatement(Expression:BoolLiteralExpression(Value:true)))      myVisitor.Members.Add(nm)// Create all 'Leave' methodsfor m as MethodInfo in [m for m in t.GetMethods() if m.Name.StartsWith("Leave")]:   // Create 'stck.Pop()' node   mie = ExpressionStatement(              MethodInvocationExpression(                         MemberReferenceExpression(                               Target:ReferenceExpression("stck"),Name:"Pop")))   // CreateNode   nm = Method(Name: m.Name,  Modifiers: TypeMemberModifiers.Override,                           Body:Block())   nm.Body.Add(mie)   nm.Parameters.Add(ParameterDeclaration(Name:"p",                                          Type:SimpleTypeReference(                                             Name:m.GetParameters().ParameterType.FullName)))   myVisitor.Members.Add(nm)   // Compile unit and module creationcu = CompileUnit()mod = Boo.Lang.Compiler.Ast.Module(Name:"Module")cu.Modules.Add(mod)mod.Imports.Add(Import(Namespace:"System.Windows.Forms"))mod.Members.Add(myVisitor)   // Initialize compiler preferencespipeline = Pipelines.CompileToMemory()ctxt = CompilerContext(cu)// Run the compilerpipeline.Run(ctxt)// Check the resultsif (ctxt.GeneratedAssembly != null):   // Create an instance of the new visitor and initialize it   i as object = ctxt.GeneratedAssembly.CreateInstance("DynamicVisitor");   t = ctxt.GeneratedAssembly.GetType("DynamicVisitor")   tStack = System.Collections.Stack()   tStack.Push(System.Windows.Forms.TreeNode("root"))   t.GetField("stck").SetValue(i,tStack );   v as DepthFirstVisitor = i   // Parse a file an run the visitor   //tast = BooParser.ParseString("myUnit","x = w.foo(1)")   tast = BooParser.ParseFile(argv)   tast.Accept(i)      //Initialize a form an show the results   frm = Form(Text: "AST content",Width: 300,Height: 300)   tv = TreeView(Dock: DockStyle.Fill)   tv.Nodes.Add(tStack.Pop() as TreeNode)   frm.Controls.Add(tv)   Application.Run(frm)else:   for e in ctxt.Errors:      print e``

Here's a sample of the output. ## Friday, June 8, 2007

### Mandelbrot Set Fractal in Fortress

One of my favorite embarrassingly parallel problems is the rendering of the Mandelbrot set fractal. It is a simple program that produces very interesting images.

There are many ways to optimize this program, however I'm going to avoid these optimizations and try to create a naive implementation of the problem, just to see to the program looks like in Fortress. Also I'm going to try to take advantage of the parallel for loops that Fortress provides.

The first thing to define is something to represent complex numbers (remember this is a non-optimize version of the problem) . It seems that in the current reference implementation of there's no native support complex number (that I could find). So this gave me the opportunity to implement it.

``value object Complex( real:RR64, img:RR64)   opr+(self,other:Complex) = Complex(real+other.real,img+other.img)   opr juxtaposition(self,other:Complex) =           Complex((real other.real) - (img other.img),                  (real other.img) + (other.real img))   toString() = ("Complex(" real ", " img ")")endopr |x:Complex| = SQRT(x.real^2 + x.img^2)complexExp (o : Complex,n : Integral) = do   if (n = 1)       then o      else (o complexExp(o,n - 1))   endendopr^(o:Complex,n:Integral) =    if (n = 2)      then (o o)     else complexExp(o,n)   end ``

Note that this implementation is not complete, only the required elements for the Mandelbrot set are implemented. It is interesting that multiplication in Fortress doesn't use the '*' operation but it uses juxtaposition of elements to be multiplied. That is, instead of saying (x*y) you say (x y).

Then we need something that helps us to convert from screen coordinates to real coordinates. To do this I created the following function:

``lFunc(x1 : RR64,y1 : RR64, x2 : RR64, y2 :RR64) = do   m = (y2 - y1) / (x2 - x1)   b = y1 - (m x1)   fn(x) => (m x) + bend``

Note that this function returns another function that preforms the conversion.

Now the following code shows the implementation of the Mandelbrot set algorithm for one row in the image.

``lineMandelbrot[\W\](startIndex : ZZ32, endIndex : ZZ32,                     lineData :Array[\ZZ32,ZZ32\],                     f :RR64 -> RR64, y : RR64) = do   for i <- startIndex:endIndex do      p0 = Complex( f(i) , y)      p : Complex := p0       j : ZZ32 := 0      maxIteration = 255      while ( |p| < 2.0  AND j < maxIteration ) do         p := p^2 + p0         j := j+1      end            lineData[i] := j   endend``

Note that I use a parallel for loop for the external for statement in order to say that each point can be calculated independently in a separate thread.

This code looks very nice when formated with Fortify: Finally the following code show the main program. For the graphic generation, a PPM file was used since is the easiest image format to generate!.

``run(a:String...) =   do     imageWidth = 100     imageHeight = 100     startY = -1.0     endY = 1.0     startX = -1.0     endX = 1.0     image =  array[\ZZ32\](imageWidth)     fout: BufferedWriter = outFileOpen "output.ppm"     outFileWrite(fout,"P3\n")     outFileWrite(fout,"" imageWidth " " imageHeight "\n")     outFileWrite(fout,"255\n")          startIY = 1     endIY = imageHeight     fY = lFunc(startIY,startY,endIY,endY)     fX = lFunc(0,startX,imageWidth - 1,endX)     for iY <- seq(startIY#endIY) do                  print "Line: " iY "\n"         lineMandelbrot(0,imageWidth - 1,image,fX,fY(iY))             for i <- seq(0#imageWidth) do            outFileWrite (fout,"0 0 " image[i]  "   ")         end         outFileWrite(fout,"\n")     end          outFileClose(fout)     ()end``

Since the reference implementation focus is not performance, I think is not fair to publish benchmarks or something like that. However I noticed that by playing around with the `NumFortressThreads` environment variable (found by looking at the source code) I got different execution times in my dual core machine .

## Sunday, June 3, 2007

### Parallel For Loops in Fortress

On interesting and promising feature of Fortress is that `for` loops are parallel by default. This is documented in the 2.8 For Loops Are Parallel by Default section of the Fortress Language Specification document.

It seems that the reference implementation already supports this feature. For example when running this code:

``component ForLoopsexport Executablerun(args:String...):() = do   for i <- 1:10 do      print("Iteration " i ""\n")   endendend``

The output shows:

``Parsing /home/ldfallas/blog/fortress/for/forloopexperiment.fss with the Rats! parser: 283 millisecondsRead /home/ldfallas/fortress/FortressLibrary.tfs: 566 millisecondsIteration 10Iteration 9Iteration 6Iteration 4Iteration 8Iteration 7Iteration 2Iteration 5Iteration 1Iteration 3finish runProgram1711 milliseconds``

According to the specification the generator section of the `for` loop controls this behavior (section 13.15 and 13.14) . If the `sequencial` generator is used, the loop is executed in the classic order. For example:

``component ForLoopsexport Executablerun(args:String...):() = do   for i <- sequential(1:10) do      print("Iteration " i ""\n")   endendend``

The output is:

``Parsing /home/ldfallas/blog/fortress/for/forloopexperimentSeq.fss with the Rats! parser: 282 millisecondsRead /home/ldfallas/fortress/FortressLibrary.tfs: 546 millisecondsIteration 1Iteration 2Iteration 3Iteration 4Iteration 5Iteration 6Iteration 7Iteration 8Iteration 9Iteration 10finish runProgram1523 milliseconds``