Wednesday, December 29, 2010

Notespad 1

I've turned my dynamic compilation code towards audio. I wrote a program recently that is a "lab" for creating audio effects: you type in a few lines of code and press Go, and the code is compiled and run. These processing snippets can be saved, along with parameter settings. I jokingly call it Windows Notespad.

With this tool I can quickly prototype audio effects. As a source I will use Cathy Berberian's reading from Joyce's Ulysses. (This is an excerpt from Luciano Berio's piece Thema, Omaggio a Joyce).

Some interesting changes can be made in the frequency domain. I divide the audio into many overlapping pieces. If I take the dft of each piece, change the phases every 256 samples, and ifft back, this imparts a high tone to the sound (frequency = 44100hz / 256 = 689hz), turning the speaking into singing. If I make each piece 0.1 seconds, ramp the pieces for a smooth transition, and set the phases to random values, I can 'smear' the sound, keeping the frequency content within each piece but losing the envelopes. This resembles adding a lot of reverb.

Hear:


The original, for comparison:

I wrote a script to plot amplitude (average signal power) over time. Plotted on a log scale, it looks like this:

One thought was to turn this curve into pitch. What if one played a sine wave corresponding to these pitches? I experimented with code until I found an effect like a bird's song.

Another non-linear transformation is to change the pitch based on power. Imagine a poor vocalist who goes sharp on loud notes and sings flat on soft notes. What if the whole band played this way? The effect is too comical to include.

We can find the loudness at a certain point in the audio clip. The idea struck me: what does it sound like if all of the quiet parts are made loud, and the loud parts made quiet? If find the effect to be rather terrifying:

There is much left to explore.

Saturday, November 6, 2010

drawdraw


drawdraw is a web app for recursive drawing. Try it!

Instructions: Click the buttons to make things. Drag the green dots to move things. Experiment until it looks cool. Send links.



Now that you have spent time playing with the program:
The light-blue arrow in the background represents the initial scale and rotation. At every step, the orange arrows are replaced with a scaled and rotated copy of the figure. If an orange arrow is longer than the light-blue arrow, the figure "explodes" because the size increases with every step.

In the end, if enough steps are done (the pencil+ button adds more steps), look at the result - the location where each orange arrow will contain the entire figure (picture a tree, where each branch can have the form of the entire tree). It is self-similar, a property held by fractals. Because this program represents a set of affine mappings, it can be used to draw all of the traditional IFS fractals in a new manner.

Have fun!

(I wonder how one would make a non-linear version...)

Saturday, October 23, 2010

plotmydata

Need to plot some data quickly? Try plotmydata.blogspot.com.



Your data is plotted in real time as you type, or you can copy/paste in your data in csv format.

This is for the times when you don't want to wait for Matlab to open, and it's not worth it to write a python matplotlib script/use a command line tool.

What features do you want - uploading text files? scatter plots? zooming? comparing multiple trends? plotting an arbitrary function? anything else? I'm all ears.

Just one of those projects. It makes use of Stanford's protovis library.

Friday, September 3, 2010

pixels

I've taken a break from writing music this summer and have done this instead.















































Wednesday, August 18, 2010

Thread Unsafety

#!/usr/bin/env python
import time,itertools
from threading import Thread

class FizzBuzz():
  g_place = ''
  def th(self, t, tminc, fn):
    start = time.clock()
    while True:
      t +=tminc
      while time.clock()<start+t: time.sleep(0.01)
      self.g_place = fn(self.g_place)
    
  def go(self, counter, printer):
    thds = ((0.7, 1.0, lambda c: printer(c)), 
      (0.0, 1.0, lambda c: counter.next()),
      (0.1, 3.0, lambda c: 'Fizz'), 
      (0.2, 5.0, lambda c: 'FizzBuzz' if c=='Fizz' else 'Buzz'))
    for thd in thds: Thread(target=self.th, args=thd).start()

def _print(s): print s
FizzBuzz().go(itertools.count(1), _print)
As far as I know, this is a unique way to code FizzBuzz. It's also a pretty bad way, for a variety of reasons. But at least it's multithreaded!

Saturday, August 14, 2010

Physics into audio

I've been toying with the idea of using physics simulations to generate audio.

First, imagine a audio in a circular loop, as on a record player. As the needle moves clockwise, the music plays. If the needle moves counter-clockwise, the music is played backwards. The faster the needle moves, the faster (and higher pitched) the sound.

Now, model the chaotic motion of a double pendulum, and use that as the needle that is driving the music:

It can also be interesting to use the data directly. Take a physics model of a double pendulum, and record the height (y) of the longer leg over time. Then, interpret that data as raw audio samples and play it!

Double pendulum as raw audio by jamonben

I inadvertently created a mashup while working on a related project, but it's not very good.

Monday, August 2, 2010

LnzScript 0.4

The LnzScript 0.4 release is here! This scripting tool can automate repetitve tasks in Windows, work with files and folders, and easily perform text manipulation.

The editor makes it easy to write quick one-off scripts. Standard output appears in a pane to the right. Let's say you have some text that you want shuffled randomly. If one were to write a Python script, you'd have to save the text in a temporary file, or load some library to access the Windows clipboard. With LnzScript, the clipboard is directly accessible without importing anything. One could just write Clipboard.set(Math.shuffleArray(Clipboard.get().split('\r\n'))) Furthermore, the LnzScript Editor can run a script before it has even been saved to a file, a unique feature, so you can run a script immediately.

LnzScript 0.4 includes more methods for working with files. File.dirListFiles(@'c:\myfolder', '*.txt','size') returns all of the text files in the directory, sorted by size. File.copyMany can copy full directories. File.createShortcut and File.getShortcutTarget can work with shortcuts. The Rename.renameFn will call a function you provide to rename files. To rename .jpeg to .jpg, say, all you need to do is:
File.cd(@'c:\myfolder');
Rename.rename('.jpeg','.jpg');

Rename.renamePreview can be used to see the results before renaming anything.

And, it's still good at automating tasks. I used it to access Webmail (for which there was not a 'remember me' option); much faster than typing in the fields every time. An old demo of LnzScript opening Firefox tabs is here.

The installer associates script files so that you can run them by double-clicking. Just by using argv in a script, you can now drag and drop a file into the script and it will run with that argument in argv[1], a 'droplet.' Lately I've been using Clavier+ to assign scripts to hotkeys. Other programs such as Launchy could also be used.

Why JavaScript? I find it to be much cleaner than Basic-derived semi-languages. Lambdas and the good parts can be very useful. I've added @'string' literals to the language, that can contain single backslashes, so that one can type @'c:\foo' instead of 'c:\\foo'. There are other additions like include(), alert(), print(), and argv.

Parsing .url shortcuts in a directory and printing out URLs:
File.cd(@'C:\folder');
var arFiles = File.dirListFiles('.', '*.url')
for(var i=0; i<arFiles.length;i++)
{
  var pageName = arFiles[i].slice(0,-4)
  var url = File.readFile(arFiles[i]).split('URL=')[1].split('\n')[0]
  print(pageName)
  print(url + '\n')  
}



An example using JavaScript's closures to rename files in the pattern 00, 01, 02, and so on, by last modification time:
include('<std>')
File.cd(@'C:\folder')
var i=0;
function renameNumbered(s)
{
  i++;
  return (i<10)?'0'+i+s : i+s;
}
//first, preview to see consequences
Rename.renameFnPreview(renameNumbered, 'time')
//then,
//Rename.renameFn(renameNumbered, 'time')

Wednesday, July 14, 2010

Audiovisual

I few months ago I was struck with an idea. What would it look like if traditional audio effects were applied to videos?

I could treat every pixel as a separate signal in time, with its brightness being the value. These become sequences of values, just like streams of audio data. I could put these signals through various audio filters, like high-pass, low-pass, and so on. (I solely apply effects in time, not space.)

I wrote code in C that worked on individual frame images at the pixel level in order to have complete control.

Applying a low pass filter to all pixels created a motion-blur like effect, because it was essentially blurring together consecutive frames. A high pass filter looked more interesting - the video became all dark except for the areas with the fastest motion. Things really started to get interesting, though, when I used other audio effects. Watch this:



A flange effect can be created with is the sum of two signals: one is the input signal, unchanged, and the other is the input signal alternating between being played slightly faster, and slightly slower. Because of the alternation, they stay in sync, but the effect is heard as a moving, shimmering sound. (Used in guitars in psychedelic rock).

I applied an exaggerated version of this effect. I combine the video signal with an altered version of itself that speeds up and slows down, just like a flange effect. This means that sometimes the altered signal is a ghostly signal sometimes ahead in the future, and sometimes lagging behind.

The effect with rainbow colors is related to a high pass filter - the areas that are changing the most are the ones that are allowed to overflow into vivid colors.

In audio, there are methods of looping a non-periodic sound. One way is to play many copies of the sound offset in time, but giving each a envelope, so that each piece on its own fades in and out, and isn't very noticeable. I used this in the last effect of the video, creating a loop.

Tuesday, July 6, 2010

Launchorz scripts

My Launchorz project makes it easy to automate repetitive tasks in Windows. Here are some of the scripts I've found to be the most useful. (Note that they are intended for Windows 7 and might not be compatible with something earlier.)

Open command-line to the current Explorer directory

include('getExplorerDirectory.js');
Window.activate( {'class':'CabinetWClass'});
// activate the most-recently-used Explorer window.
Time.sleep(100);
var strDir = getCurrentExplorerDirectory();
if (strDir)
  Process.openFile('cmd.exe', strDir);

You can download getExplorerDirectory.js.

Creating an index of all files in a directory

This script creates a text file index of all files in a directory, including subdirectories. The output is formatted nicely by indentation level. It uses the path of the currently open Explorer window and creates tree_files.txt.
include('getExplorerDirectory.js');

Window.activate( {'class':'CabinetWClass'})
Time.sleep(100);
var strDir = getCurrentExplorerDirectory();
if (strDir)
  Process.open('cmd.exe /c tree /f /a > tree_files.txt', strDir);

Expand repeated code

When working on a quick project, or writing in a language like Verilog, sometimes the same line of code is repeated a few times, with a different index. The following script takes a code fragment like
a[#] = b[#] & c[#];
and expands it to
a[0] = b[0] & c[0];
a[1] = b[1] & c[1];
a[2] = b[2] & c[2];
...
(This isn't a useful example but it depicts what the script does). All it does is replace the '#' character with a number.
include('<std>');
var sClipboard = Clipboard.getText();
if (sClipboard  && sClipboard.contains('#'))
{  
  var sRep = Dialog.input('Replace', 'How many repetitions? n=0 to ?', '3');
  if (isNumber(sRep))
  {    
    var sResult = '';
    for (var i=0; i<sRep; i++)
      sResult += sClipboard.replace(/#/g, i) + '\r\n';
    Clipboard.setText(sResult);
  }
}

Sort lines in the clipboard

var strClipboard = Clipboard.getText();
if (strClipboard) 
{
  var arr = strClipboard.split('\r\n');
  arr.sort();
  var strResult = arr.join('\r\n');
  Clipboard.setText(strResult);
}
(You could also normalize the line endings with strClipboard.replace(/\r\n/g, '\n') and split on \n).

Displaying an ANSI table

//opens an ascii table in notepad.
var s ='';
for (var i=32; i<128; i++)
{
  var si = (i<100)?' '+i : i.toString();
  s+= si + '\t' + String.fromCharCode(i);
  s+= '\r\n';
}
File.writeFile( File.getPathTemp()+'\\tmpascii.txt', s);
Process.open('notepad.exe '+File.getPathTemp()+'\\tmpascii.txt');
Showing a simple reference of ansi characters.

I also have a clipboard replacement script that gives you 26 different clipboards, each recalled by a letter. However, there exist many other clipboard-handling tools.

To learn more about Launchorz, you can watch some screencasts here. (It is similar to AutoIt, but has a full JavaScript language and namespaces).

Wednesday, June 2, 2010

Chaotic Maps!

This is the first release of "fastmapmap", a program I wrote that plots 2D chaotic maps. A discrete chaotic map will take the point (x,y) to a new (x_, y_) based on a function. For example, the Henon map takes the point (x,y) to (y+1-a*x*x, b*x), where a and b are constants. One can draw a "portrait" of a map by choosing an (x0, y0) and repeatedly evaluating the function, shading in each point that is landed on. Depending on the values of a and b, sometimes the point moves off to become very large, or sometimes the point alternates jumping between two or more values. Other times, the points shade in a region in a complicated way, never repeating, which is a chaotic orbit.

The program fastmapmap can draw these portraits in real time. (I put some effort into optimization). The parameters a and b can be quickly changed with the arrow keys. The plot can be easily navigated; hold alt and drag the mouse to zoom in on a region.

The program is called fastmapmap because an additional plot can be drawn that is a "map" of the behavior of the map. The x axis represents the constant a and the y axis is b. In this plot, black areas are periodic, red areas escape to infinity, and colored areas are chaotic. So, the plot can depict which parameter values may be interesting to observe.
It's cross-platform, tested in Windows and Ubuntu. Try it! Download Windows | Source (GPL)

A short tutorial:
  • Open fastmap.exe
  • Tap the right arrow a bit to change the value of a.
  • Press Ctrl-s, type 'test', and press Enter to save as test.cfg
  • Ctrl-click on the shape to zoom in, Shift-click to zoom out
  • Press Ctrl-o, type 'test' and press Enter to open our test.cfg
  • Press Alt-b and watch the shapes move
  • Press Alt-b again to stop movement
  • Click the cross-hairs icon to show the diagram. Click in the diagram.
  • Click the "..." icon to read what else can be done.
The program can easily create animations. Hit Alt-B and watch the shapes start to move. Use Ctrl-F1 through Ctrl-F9 to set a keyframe. Refer to readme.txt for all of the other features.

 

 

 

I came up with other drawing modes. In the following, x and y represent x0 and y0, and the color is chosen according to the value after 20 iterations.




Assigning colors based on quadrant:





I invented this drawing method originally as way of depicting the action of the Henon map.

 


 

Remember to download: Windows | Source

SDL Text and Dialogs

For one of my recent projects, I wanted dialogs in an sdl program. This small library draws simple "dialogs" to ask the user for text or a number.



Usage:
// where SDL_Surface* pSurface is the screen surface.
char * s; 
s= Dialog_GetText("Enter your name:", "Fred", pSurface);
Dialog_Messagef(pSurface, "Welcome, %s", s);
free(s);

double val = 1.5;
Dialog_GetDouble("Enter a value", pSurface, &val);
int n = 1;
Dialog_GetInt("Choose a number.", pSurface, &n);
if (Dialog_GetBool("Double the number?", pSurface))
  n*=2;
Dialog_Message("Goodbye", pSurface);
freeFonts();


To draw text, I used the code at this page. (I added the ability to have newlines and tabs). It supports lower ascii chars, and code could be improved, but currently suits my purposes. There's also SFont and others.

Source (includes a bitmap font).

Sunday, April 25, 2010

SSE instructions

For the past few weeks I have been experimenting with Streaming SIMD Extension instructions. It takes a while to write at this lower level, but I've been able to get my code to run 3 times faster!

Most processors since about 1999 contain eight 128 bit registers that can be used to do multiple computations in one instruction. 128 bits can hold 4 32-bit floating point numbers, and so there could be a theoretical 4x improvement in speed, as long as the same operation is done to each. For example, if you have two arrays a and b , and want to sum them, instead of:
float a[4]; float b[4]; float c[4];
//...
c[0] = a[0]+b[0];
c[1] = a[1]+b[1];
c[2] = a[3]+b[2];
c[3] = a[3]+b[3];

, this can be expressed as:
__m128 a, b, c;
// ...
c = _mm_add_ps (a, b);
The _mm_add_ps "instrinsic" becomes 1 instruction on the cpu (addps) instead of 4. Similar code can be written with 32-bit integers, and a few other data types. (Technology like CUDA is faster, though). In case you are interested, here is an example of optimizing with SSE instructions.

The following code will draw a chaotic map (c1 and c2 are parameters). For this particular map, we know that initial points or "seeds" close to (0.0,0.0) are good starting points. Let's say, for fun, that we want this to run as fast as possible.
const int nseeds=4;
float xseeds[nseeds] = {0.0f, 0.0f, 0.0f, 0.0f};
float yseeds[nseeds] = {0.000001f, 0.000002f,-0.0000011f,-0.0000011f};
float x, y, x_, y_;
for (int seed=0; seed<nseeds; seed++)
{
  x=xseeds[seed]; y=yseeds[seed];
  // ... settling code ...
  
  for (int i=0; i<PixelsToDrawPerSeed; i++)
  {
    x_ = c1*x - y*y; 
    y_ = c2*y + x*y;
    x=x_;
    y=y_;
    if (ISTOOBIG(x)||ISTOOBIG(y)) break; 
    //if escaping to infinity, there won't be more points drawn, so quit early.
    //quitting early is important, far slower without this.

    //scale to pixel coordinates, and plot point if it is within bounds. X1 is greatest x, and so on.
    int px = (int)(SIZE * ((x - X0) / (X1 - X0)));
    int py = (int)(SIZE * ((y - Y0) / (Y1 - Y0)));
    if (py >= 0 && py < SIZE && px>=0 && px<SIZE)
      arr[px+py*SIZE] = 1; //plot point at px, py.
  }
}
Runtime for 1,310,720 different parameters is Test 1: 14,747ms, Test 2: 25,026ms. First, I added threading to to take advantage of the two cpu cores. I made sure to balance the work so that each thread is kept busy. This improved time to Test 1:7,452ms Test 2:12,617ms.

At first glance, it's not obvious how this could be optimized. I was a bit concerned about the (int) casts becoming a slow _ftol call, but looking at the generated assembly showed that this wasn't the case. (Read here,here for why casts can be slow, although with modern compilers and instruction sets this doesn't seem to be an issue).

In visual studio, I noticed a setting for Streaming SIMD Extensions 2 (/arch:SSE2). With this enabled, the compiler was able to vectorize some of the code (I verified in the asm that it was using the xmm registers), which resulted in an ok speed increase. Test 1:5,981ms Test 2:10,987ms

Can I do better?

Here is my first attempt at a vectorized version, moving all 4 points at once.
__m128 x128 = _mm_setr_ps( 0.0f, 0.0f, 0.0f, 0.0f);
__m128 y128 = _mm_setr_ps( 0.000001f, 0.000002f,-0.0000011f,-0.0000011f);
__m128 mConst_a = _mm_set1_ps(c1); //sets all 4 fields to c1
__m128 mConst_b = _mm_set1_ps(c2);
__m128 nextx128, tempx, tempy;

__m128i const_size = _mm_set1_epi32(SIZE), const_zero = _mm_set1_epi32(0);
// ... settling code ...
for (int i=0; i<PixelsToDrawPerSeed; i++)
{
  nextx128 = _mm_sub_ps(_mm_mul_ps(mConst_a, x128), _mm_mul_ps(y128,y128));
  y128 = _mm_mul_ps(y128, _mm_add_ps(x128, mConst_b)); 
  x128 = nextx128;

  if (ISTOOBIG(x128.m128_f32[0]) ||ISTOOBIG(x128.m128_f32[1]) ||
    ISTOOBIG(x128.m128_f32[2]) ||ISTOOBIG(x128.m128_f32[3]) ||
    ISTOOBIG(y128.m128_f32[0]) ||ISTOOBIG(y128.m128_f32[1]) ||
    ISTOOBIG(y128.m128_f32[2]) ||ISTOOBIG(y128.m128_f32[3] ))
  {
    break;
  } //assume that if one is too big, all four will be too big, nearly always valid for this map

  //scale to pixel coordinates. assume that pixels are square, that X1-X0==Y1-Y0.
  //the following computes px = (int)(SIZE * ((x - X0) / (X1 - X0)));
  __m128 mConst_scale = _mm_set1_ps(SIZE / (X1 - X0));
  tempx = _mm_sub_ps(x128, _mm_set1_ps(X0));
  tempx = _mm_mul_ps(tempx, mConst_scale); 
  tempy = _mm_sub_ps(y128, _mm_set1_ps(Y0));
  tempy = _mm_mul_ps(tempy, mConst_scale); 
  __m128i xPt = _mm_cvttps_epi32(tempx); //cast to int32
  __m128i yPt = _mm_cvttps_epi32(tempy); 
  __m128i xInBounds =  _mm_and_si128(_mm_cmpgt_epi32(xPt, const_zero), _mm_cmplt_epi32(xPt, const_size));
  __m128i yInBounds =  _mm_and_si128(_mm_cmpgt_epi32(yPt, const_zero), _mm_cmplt_epi32(yPt, const_size));
  __m128i inBounds = _mm_and_si128(xInBounds, yInBounds);
  
  if (inBounds.m128i_i32[3] != 0)
    arr[xPt.m128i_i32[3] +yPt.m128i_i32[3]*SIZE] = 1; //plot point
  if (inBounds.m128i_i32[2] != 0)
    arr[xPt.m128i_i32[2] +yPt.m128i_i32[2]*SIZE] = 1; //plot point
  if (inBounds.m128i_i32[1] != 0)
    arr[xPt.m128i_i32[1] +yPt.m128i_i32[1]*SIZE] = 1; //plot point
  if (inBounds.m128i_i32[0] != 0)
    arr[xPt.m128i_i32[0] +yPt.m128i_i32[0]*SIZE] = 1; //plot point
}
It is convinient that conditions can be expressed in SSE with the and instruction and comparison instruction. If a comparison is true, the field is filled with 0xffffffff, and if false, 0x00000000. This version was faster. Test 1: 4,650ms Test 2: 7,293ms

It was clear from the generated assembly that the ISTOOBIG(x128.m128_f32[0]) ||ISTOOBIG(x128.m128_f32[1] was inefficient, as it had to move data out of the xmm registers before comparing. I then learned about _mm_movemask_ps, which does exactly what I want. It condenses the 128 bits into 4 bits of an integer in one step, much faster than moving out data four times, and so I could write
//instead of ISTOOBIG(x128.m128_f32[0])||ISTOOBIG(x128.m128_f32[1]) and so on
__m128 estimateMag = _mm_add_ps(x128,y128); //estimate magnitude as (x+y)(x+y)
estimateMag = _mm_mul_ps(estimateMag,estimateMag);
__m128 isTooBig = _mm_cmpgt_ps(estimateMag, constThreshold);
if (_mm_movemask_ps(isTooBig) != 0) //if any of the four fields are true
  return 0;

This increased speed dramatically! The score is now Test 1: 1,956ms, Test 2: 3,323ms! I later saw on this site assembly for a quick absolute value hack. I used the idea and rewrote as
__m128 estimateMag = _mm_add_ps(x128,y128); //estimate magnitude as abs(x+y)
estimateMag = _mm_andnot_ps(_mm_set1_ps(-0.0f), estimateMag); //credit: Concatenative language wiki
__m128 isTooBig = _mm_cmpgt_ps(estimateMag, constToobig);
if (_mm_movemask_ps(isTooBig) != 0) //if any of the four fields are true
  return 0;

A further optimization was to vectorize finding the array pixel location. Instead of "if (inBounds.m128i_i32[3] != 0)...", I could write:
//Make size a power of 2 so multiplication is a shift left. Example: SIZE=128,LOG2SIZE=7
__m128i pixel = _mm_add_epi32(xPt, _mm_slli_epi32(yPt, LOG2SIZE)); //pixel=xPt+yPt*size 
pixel = _mm_and_si128(inBounds, pixel); //draw out-of-bounds pixels at 0,0
arr[pixel.m128i_i32[3]] = 1; //plot point
arr[pixel.m128i_i32[2]] = 1; //plot point
arr[pixel.m128i_i32[1]] = 1; //plot point
arr[pixel.m128i_i32[0]] = 1; //plot point

Note that I've eliminated some conditionals. The trick is to mask out the point if inBounds is false with the line "pixel = _mm_and_si128(inBounds, pixel);". False is represented by 0x00000000, and so if inBounds is false, pixel is set to 0. All out-of-bounds pixels are drawn onto (0,0), but this pixel can be cleared later. Speed is improved to Test 1: 1,726ms, Test 2: 2,991ms.

If it were likely that a pixel would be out-of-bounds, I could change the order and check inbounds much earlier, before casting to int, and do things like checking if x is in bounds before calculating y. Because this case not occur often, it is not worth attention.

I notice some redundancy in "int py = (int)(SIZE * ((y - Y0) / (Y1 - Y0)))" and "pixel=x+y*size." In both, y is multiplied by a constant. I could eliminate a multiply with something like:
int py_times_size = (int)(SIZE*SIZE* ((y - Y0) / (Y1 - Y0)))
if py_times_size>0 and py_times_size<SIZE*SIZE and x>0...
pixel=x+py_times_size.

This didn't work in sse vectors because of overflow with _mm_cvttps_epi32.

My final step is the least elegant: loop unrolling. Repeating the loop body 4 times gave the fastest results. (Faster because less overhead and maybe better scheduling, there is a trade off because of cache misses as size increases). There are additional savings because I only need one "isTooBig test" per every four iterations, and so I deleted all but one of the tests from the loop. It doesn't matter if the point is too big for a few iterations as long as we realize it within a few iterations. This final step also made a significant difference.

Our final times are Test 1: 1,346ms, Test 2: 2,285ms. This amount of optimization, by about a factor of 4, can be the difference in a real-time program. The code in this example, by the way, makes it possible to quickly draw this rather striking figure:

Monday, March 8, 2010

Bifurcation

I wrote a program to draw bifurcation diagrams of any 1-dimensional map.


Recently, I studied 1-dimensional maps in my nonlinear dynamics course. I recommend reading about this topic: it is fascinating that an expression as simple as p_next = r*p*(1-p) can lead to what is known as "chaos." r is a constant, p is a number between 0 and 1, and the expression is repeatedly evaluated. The resulting sequence of p can either settle to a fixed number, periodically cycle between values, or continue aperiodically and seemingly in a random fashion. A bifurcation diagram concisely depicts the behavior as the parameter r changes. Read more about the logistic map on Wikipedia or elsewhere.

I used Matlab for most numerical work pertaining to the course. In my spare time, though, I wrote this program to experiment with different maps and parameters. Also, because the program runs so much more quickly, I could create better images without waiting. (This program builds upon my earlier project that dynamically compiles code, and so it runs efficiently.)

Download (Unzip and run CsBifurcation.exe).
Source (GPLv3 license)

To use the program, type in an expression, and click Go. By default, the famous logistic map is shown. Use the mouse to zoom in. By holding the Alt key and dragging, you can zoom in on a non-square rectangle of the plot. Refer to the readme.txt file for more information. I encourage you to run this program and explore the intricate details of the logistic map, and to invent your own 1D maps. The program can also create animations by slowly changing one of the parameters over time. I will post some of the animations I've made at a later date.

As I wrote new expressions and drew bifurcation diagrams, I began to search for interesting shapes and details. I began to compose images based on appearance rather than mathematical properties. A selection of this work:



Tuesday, February 9, 2010

An Excessively Fast Graphing Calculator

Lately I've been exploring the code-generation features of the .NET framework. It's cool to be able to compile and run new code at runtime. The CSharpCodeProvider and System.CodeDom.Compiler classes provide this ability.

Evaluating code at runtime is especially useful when accepting user scripts, or evaluating a new mathematical expression. In Python, one can use eval() or exec to run code, but in C# the process is more complicated. So, I wrote a wrapping class with a convenient interface, that can be reused in other projects.

For example, public double simpleMathEval(string strExp, out string strErr) can easily evaluate an expression: double ret = simpleMathEval("1+3+Math.sin(4.6)", out strErr); The more advanced method mathEval accepts any amount of code and inner functions.

As a proof of concept, I wrote a graphing calculator. It's really fast - compiling your expression. When you click Plot, a new program is generated that is essentially a loop containing your expression. (Compilation or runtime errors are caught and relayed back to you). On the graph, zoom in easily by clicking and dragging to draw a rectangle.


(Note that I try to be intelligent when drawing asymptotes, making this plot of cosecant look cleaner than if it'd been plotted in Matlab). I based the code on the efforts of Mike Gold's "CodeDom Calculator", but made a new interface and extended it to allow arrays. As the title of this blog suggests, though, I haven't polished or thoroughly tested it.


A more advanced example - you can put all sorts of code in here.

To see how to use CodedomEvaluator, CodedomTest.cs contains some examples. CodedomEvaluator Usage:
double simpleMathEval(string strExp, out string strErr)
-evaluates a string, not a full statement. example: "3.0*5.1"returns 15.3.
double simpleMathEval(string strExp, string strVarname, double varValue, out string strErr)
-same as above, but provide a variable. 
example: strVarname="x", varValue=3, "x*x" returns 9.0.
double mathEval(string strExp, Dictionary<string, double> vars, out string strErr)
-evaluates a string that is a full statement. must assign to "ans". 
the vars dictionary can be used to provide variables.
example: vars["x"] = 3.0, "ans=x+4;" returns 7.0.
double[] mathEvalArray(string strExp, Dictionary<string, double> vars, int arrayLen, out string strErr)
-like above, but will return an array of values. The array arrAns[] will be created of length arrayLen.
example: arrayLen=40, "for(int i=0;i<40;i++) arrAns[i]=i*i;"
Because the inner loop is compiled, 
this is a fast way to evaluate an expression many times.


Windows binary
C# src, GPL v3.