Tag Archive: programming


OOP in LUA

Currently I’m writing a game for Android using Corona SDK. And its main scripting language is LUA. My first impression of LUA
was that it is extremely easy to learn. However, being a lightweight scripting language, LUA doesnt support OOP in the direct way like Java/C++. Instead you have to go the round way using its metatable feature. Detail can be found via http://www.lua.org/pil/16.html.

So what if I want to implement my classes in separate files. It brought me surprise that there’s not much information about this on the Internet. So, feeling a bit frustrated, I decided to do it on my own.

First let’s start with the standard way. I.e. declare everything in your main.lua

    Account = {balance = 0}
    
    function Account:new (o)
      o = o or {}
      setmetatable(o, self)
      self.__index = self
      return o
    end
    
    function Account:deposit (v)
      self.balance = self.balance + v
    end
    
    function Account:withdraw (v)
      if v > self.balance then error"insufficient funds" end
      self.balance = self.balance - v
    end

    SpecialAccount = Account:new()
    s = SpecialAccount:new{limit=1000.00}
    function SpecialAccount:withdraw (v)
      if v - self.balance >= self:getLimit() then
        error"insufficient funds"
      end
      self.balance = self.balance - v
    end
    
    function SpecialAccount:getLimit ()
      return self.limit or 0
    end

    function s:getLimit ()
      return self.balance * 0.10
    end

Now let’s do the partitioning
Account class (account.lua)

module  (...,package.seeall)

function new (ret)
	local Account = {balance=0}

	function Account:new (ret) 
		ret = ret or {}
		setmetatable (ret,self)
		self.__index = self
		return ret
	end

	function Account:deposit (amount)
		self.balance = self.balance + amount
	end

	function Account:withdraw (amount)
		if (amount>self.balance) then
			print ("not sufficient money")
		else
			self.balance = self.balance - amount
		end
	end

	return Account:new (ret)
end

SpecialAccount (specialaccount.lua)

module (...,package.seeall)
Account = require ("account")

function new (arg)
	local SpecialAccount = Account.new (arg)
	
        -- method override
	function SpecialAccount:withdraw (amount)
		print ("asdf")
		if (amount>self.balance+self.getLimit()) then
			print ("not surficient money& limit")
		else
			self.balance = self.balance - amount
		end
	end

	function SpecialAccount:getLimit ()
		return self.limit
	end

	return SpecialAccount
end

To use these 2 classes, first we need to import them

SpecialAccount = require ("specialaccount")
Account = require ("account")

The rest is just to create their instances (notice the use of “.” not “:”) and calling corresponding methods

local s =  SpecialAccount.new ({balance=200,limit=500})
s:deposit (50)
s:withdraw (700)
print (s.balance)
local a = Account.new ({balance=100})
a:withdraw (101)
print (a.balance)
Advertisements

2D n-body problem

The n-body problem is a classic example used in parallel programming. The straight forward solution takes O(N^2) time complexity, whereas the Barnes-Hut reduce it to O(NlogN) based on the assumption that, a group of bodies that is far enough from a body can be clustered in to a BIG body.

Based on that assumption, we can divide the space into sub 4 equal rectangles (2D) or 8 cubes (3D). Then we build a quad tree (2D) or oct tree (3D) whose root is the bounding box of all n bodies and divide each node into 4(2D)/8(3D) until each leaf contains no more than 1 body. Each node will need to store its bounding box, its mass (total mass of the children) and center of gravity, which is calculated by the formula:

COM = (C1*M1+C2*M2+C3*M3+C4*M4)/M (for quad tree)
COM = (C1*M1+…+C8*M8)/M (for quad tree)
where Ci is the position of Child ith, Mi is the respective mass and M is the total mass of all children

To compute the total force acting on body X, we just need to traverse down the tree, at each node Y, we calculate the distant R of the node to X. if R >= d*theta, where d is the largest size of the bounding box enclosing Y, theta is a constant (in my code I choose theta = 1.2), we treat Y as a body and calculate its force on X. Otherwise, we go down to Y’s children.

I have written 2 versions. One is Quadtree, and the other uses OpenMp. The quadtre version is no doubt faster but openmp is much shorter and easy to implement. And it considerably faster than naive O(N^2) version.

To compile the sourcode under Linux, you need to add “-lX11 -lm” as the program will rely on XWindows to draw
To compile OpenMp version, add -openmp

QuadTree version

#include <X11/Xlib.h>
#include <X11/Xutil.h>
#include <X11/Xos.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <iostream>
#include <sys/time.h>
#include <vector>

#define N 1000000
#define G 6.673e-11
#define T_MAX 100
#define dt 1
#define THETA 1.2
#define  X_RESN	800       /* x resolution */
#define  Y_RESN	800       /* y resolution */

using namespace std;


typedef struct complextype 
{
    float real, imag;
} Compl;

struct  QTreeNode {
	float x0,y0,x1,y1,mass;
	float COM[2]; // center of mass
	bool is_leaf;
	vector<int> id;	
	struct QTreeNode** children;
};



Window        win;  /* initialization for a window */
GC black_gc; // to draw the body
GC white_gc; // to clear the body
GC gc;
Display       *display;
unsigned int  screen; /* which screen */

float m[N]; // mass 
float f[N][2],v[N][2],p[N][2]; // forces, velocity, position
int n,t_max;

void init_display () {
    unsigned 
    int           width, height,  /* window size */
                  x, y, /* window position */
                  border_width, /*border width in pixels */
                  display_width, display_height;
                  /* size of screen */
                  char *window_name = "Mandelbrot Set", 
                  *display_name = NULL;
    unsigned
    long          valuemask = 0;
    XGCValues     values;
    XSizeHints    size_hints;
    Pixmap        bitmap;
    XPoint        points[800];
    FILE          *fp, *fopen ();
    char          str[100];

	XSetWindowAttributes attr[1];

    /* Mandlebrot variables */
    int           i, j, k;
    
    /* connect to Xserver */

    if ( (display = XOpenDisplay (display_name)) == NULL ) {
	  fprintf (stderr, 
                  "Cannot connect to X server %s\n",
                  XDisplayName (display_name) );
	  exit (-1);
    }
	
    /* get screen size */

    screen = DefaultScreen(display);
    display_width = DisplayWidth(display, screen);
    display_height = DisplayHeight(display, screen);

    /* set window size */

    width = X_RESN;
    height = Y_RESN;

    /* set window position */

    x = 0;
    y = 0;

    /* create opaque window */
	border_width = 4;
    win = XCreateSimpleWindow(display, 
                  RootWindow(display, screen),
                  x, y, width, height, border_width, 
                  BlackPixel(display, screen), 
                  WhitePixel(display, screen));

    size_hints.flags = USPosition|USSize;
    size_hints.x = x;
    size_hints.y = y;
    size_hints.width = width;
    size_hints.height = height;
    size_hints.min_width = 300;
    size_hints.min_height = 300;

    XSetNormalHints(display, win, &size_hints);
    XStoreName(display, win, window_name);
    /* create graphics context */
	/*
	values.foreground = WhitePixel (display, screen);
    white_gc = XCreateGC(display, win, valuemask, &values);
	values.foreground = BlackPixel (display, screen);
	black_gc = XCreateGC (display, win, valuemask, &values);
	*/
	gc = XCreateGC (display, win, valuemask, &values);
    XSetBackground(display, gc, 
                  WhitePixel (display, screen));
    XSetForeground(display, gc, 
                  BlackPixel (display, screen));
    XSetLineAttributes(display, gc, 1, 
                  LineSolid, CapRound, JoinRound);

    attr[0].backing_store = Always;
    attr[0].backing_planes = 1;
    attr[0].backing_pixel = BlackPixel(display, screen);

    XChangeWindowAttributes(display, win, 
                  CWBackingStore | CWBackingPlanes | CWBackingPixel, 
                  attr);

    XMapWindow(display, win);
    XSync(display, 0);

}


float abs (float x) {
	return (x>=0)?x:-x;
}

void build_qtree (struct QTreeNode* aNode) {
	if (aNode->id.size () < 2) {
		aNode->is_leaf = true;
		return;
	}
	aNode->children = new QTreeNode*[4];
	float x0,y0,x1,y1;
	x0 = aNode->x0; 
	x1 = aNode->x1;
	y0 = aNode->y0;
	y1 = aNode->y1;
	for (int i=0;i<4;++i) { 
		aNode->children[i] = new QTreeNode;
		aNode->is_leaf = false;
		if (i==0||i==2) {
			aNode->children[i]->x0 = x0;
			aNode->children[i]->x1 = (x0+x1)*0.5;
		}
		else {
			aNode->children[i]->x0 = (x0+x1)*0.5;
			aNode->children[i]->x1 = x1;
		}
		if (i<2) {
			aNode->children[i]->y0 = y0;
			aNode->children[i]->y1 = (y0+y1)*0.5;
		}
		else {
			aNode->children[i]->y0 = (y0+y1)*0.5;
			aNode->children[i]->y1 = y1;
		}
	}
	for (int i = 0;i<aNode->id.size ();++i) {
		// decide which child the body belongs to 
		int pid = aNode->id[i];
		float x = p[pid][0];
		float y = p[pid][1];
		for (int j=0;j<4;++j) {
			if (x>=aNode->children[j]->x0 && x<=aNode->children[j]->x1 && y>=aNode->children[j]->y0 && y<=aNode->children[j]->y1) {
				aNode->children[j]->id.push_back (pid);
				break;
			}
		}
	}
	for (int i=0;i<4;++i)
		build_qtree (aNode->children[i]);
}

void compute_mass (struct QTreeNode* aNode) {
	if (aNode->id.size () == 0) return;
	aNode->mass = 0;
	aNode->COM[0] = aNode->COM[1] = 0;
	if (aNode->is_leaf) {
		aNode->mass = m[aNode->id[0]];
		aNode->COM[0] = p[aNode->id[0]][0];
		aNode->COM[1] = p[aNode->id[0]][1];
	}
	else {
		for (int i=0;i<4;++i) {
			struct QTreeNode* child = aNode->children[i];
			compute_mass (child);
			aNode->mass+=child->mass;
			aNode->COM[0]+=child->mass*child->COM[0];
			aNode->COM[1]+=child->mass*child->COM[1];
		}
		aNode->COM[0]/= aNode->mass;
		aNode->COM[1]/= aNode->mass;
	}
}

void force (const float& x1, const float& y1, const float& x2, const float& y2, const float& m1, const float& m2, float& Fx, float& Fy) {
	float r = sqrt ((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
	if (r>0) {
		r*=r*r;
		Fx +=G*m1*m2*(x2-x1) / r;
		Fy +=G*m1*m2*(y2-y1) / r;
	}
}
void cal_force (struct QTreeNode* aNode, const float& mm, const float& x, const float& y, float& Fx, float& Fy) {
	if (aNode->is_leaf) {
		if (aNode->id.size () == 0) return;
		int pid = aNode->id[0];
		force (x,y,p[pid][0],p[pid][1],mm,m[pid], Fx, Fy);
		return;
	} 
	float r = sqrt (pow (aNode->COM[0]-x,2) + pow (aNode->COM[1]-y,2));	
	float d = max (abs (aNode->x1-aNode->x0), abs (aNode->y1-aNode->y0));

	// if only cluster of bodies is far enough
	if (r>d*THETA) {
		force (x,y,aNode->COM[0],aNode->COM[1],mm,aNode->mass,Fx,Fy);
		return;
	}
	for (int i=0;i<4;++i) 
		cal_force (aNode->children[i], mm, x,y,Fx,Fy);
}

void compute_force (struct QTreeNode* aNode) {
	for (int i=0;i<n;++i) {
		cal_force (aNode, m[i], p[i][0], p[i][1], f[i][0],f[i][1]);	
	}
}
/*
void traverse (struct QTreeNode* aNode) {
	if (aNode->id.size ()==0) cout<<"nil"<<endl;
	else if (aNode->is_leaf) cout<<"L: "<<aNode->id[0]<<endl;
	else {
		for (int i = 0;i<aNode->id.size();++i) 
			cout<<aNode->id[i]<<" ";
		cout<<endl;
		for (int i=0;i<4;++i)
			traverse (aNode->children[i]);
	}
}
*/

main (int argc, char* argv[])
{
	init_display ();     	 
	FILE* fin = fopen ("nbody.txt","r");
    /* Calculate and draw points */
    Compl         z, c;
    float         lengthsq, temp;
	fscanf (fin,"%d %d",&n,&t_max);

		
	for (int i = 0;i<n;++i) {
		m[i] = 1;
		f[i][0] = f[i][1] = v[i][0] = v[i][1] = 0;
		fscanf (fin, "%f %f",&p[i][0],&p[i][1]);
	}

	float p_new[n][2],v_new[n][2];
	cout<<n<<endl;

	struct timeval tv1,tv2;
	struct timezone tz1,tz2;
	gettimeofday (&tv1,&tz1);
	for (int t=0;t<t_max;++t) {
		// build quad tree
		struct QTreeNode* root = new QTreeNode;
		root->x0 = root->y0 = 0;
		root->x1 = root->y1 = 1;
		root->is_leaf = false;
		for (int i=0;i<n;++i) 
			root->id.push_back (i);
		build_qtree (root);
		//traverse (root);
		compute_mass (root);
		compute_force (root);
		
		for (int i=0;i<n;++i) {
			v_new[i][0] = v[i][0] + f[i][0]*dt/m[i];
			v_new[i][1] = v[i][1] + f[i][1]*dt/m[i];
			p_new[i][0] = p[i][0]+v_new[i][0]*dt/m[i];
			p_new[i][1] = p[i][1]+v_new[i][1]*dt/m[i];
		}

		for (int i=0;i<n;++i) {
			int x,y;
			x = floor (p[i][0]*X_RESN);
			y = floor (p[i][1]*Y_RESN);
			//XSetForeground (display, gc, WhitePixel (display,screen));
			//XDrawPoint (display, win, gc, x,y);
			v[i][0] = v_new[i][0];
			v[i][1] = v_new[i][1];
			p[i][0] = p_new[i][0];
			p[i][1] = p_new[i][1];
			x = floor (p[i][0]*X_RESN);
			y = floor (p[i][1]*Y_RESN);
			//cout<<x<<" "<<y<<endl;
			//XSetForeground (display, gc, BlackPixel (display,screen));
			//XDrawPoint (display, win, gc, x, y);
		}
		//sleep (1);
		
	}
	gettimeofday (&tv2,&tz2);
	printf ("Eslaped time = %ld\n", (tv2.tv_sec-tv1.tv_sec)*1000000 + tv2.tv_usec-tv1.tv_usec);
    XFlush(display);

    /* Program Finished */

}

OpenPM version

#include <X11/Xlib.h>
#include <X11/Xutil.h>
#include <X11/Xos.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <iostream>
#include <sys/time.h>
#include <omp.h>

#define N 10   				// number of bodies
#define G 6.673e-11
#define T_MAX 100
#define N_THREADS 7			// default number of threads
#define dt 1
#define  X_RESN	800       /* x resolution */
#define  Y_RESN	800       /* y resolution */

using namespace std;


typedef struct complextype 
{
    float real, imag;
} Compl;


Window        win;  /* initialization for a window */
GC black_gc; // to draw the body
GC white_gc; // to clear the body
GC gc;
Display       *display;
unsigned int  screen; /* which screen */

void init_display () {
    unsigned 
    int           width, height,  /* window size */
                  x, y, /* window position */
                  border_width, /*border width in pixels */
                  display_width, display_height;
                  /* size of screen */
                  char *window_name = "Mandelbrot Set", 
                  *display_name = NULL;
    unsigned
    long          valuemask = 0;
    XGCValues     values;
    XSizeHints    size_hints;
    Pixmap        bitmap;
    XPoint        points[800];
    FILE          *fp, *fopen ();
    char          str[100];

	XSetWindowAttributes attr[1];

    /* Mandlebrot variables */
    int           i, j, k;
    
    /* connect to Xserver */

    if ( (display = XOpenDisplay (display_name)) == NULL ) {
	  fprintf (stderr, 
                  "Cannot connect to X server %s\n",
                  XDisplayName (display_name) );
	  exit (-1);
    }
	
    /* get screen size */

    screen = DefaultScreen(display);
    display_width = DisplayWidth(display, screen);
    display_height = DisplayHeight(display, screen);

    /* set window size */

    width = X_RESN;
    height = Y_RESN;

    /* set window position */

    x = 0;
    y = 0;

    /* create opaque window */

    border_width = 4;
    win = XCreateSimpleWindow(display, 
                  RootWindow(display, screen),
                  x, y, width, height, border_width, 
                  BlackPixel(display, screen), 
                  WhitePixel(display, screen));

    size_hints.flags = USPosition|USSize;
    size_hints.x = x;
    size_hints.y = y;
    size_hints.width = width;
    size_hints.height = height;
    size_hints.min_width = 300;
    size_hints.min_height = 300;

    XSetNormalHints(display, win, &size_hints);
    XStoreName(display, win, window_name);
    /* create graphics context */
	/*
	values.foreground = WhitePixel (display, screen);
    white_gc = XCreateGC(display, win, valuemask, &values);
	values.foreground = BlackPixel (display, screen);
	black_gc = XCreateGC (display, win, valuemask, &values);
	*/
	gc = XCreateGC (display, win, valuemask, &values);
    XSetBackground(display, gc, 
                  WhitePixel (display, screen));
    XSetForeground(display, gc, 
                  BlackPixel (display, screen));
    XSetLineAttributes(display, gc, 1, 
                  LineSolid, CapRound, JoinRound);

	//cout<<"here"<<endl;
    attr[0].backing_store = Always;
    attr[0].backing_planes = 1;
    attr[0].backing_pixel = BlackPixel(display, screen);

    XChangeWindowAttributes(display, win, 
                  CWBackingStore | CWBackingPlanes | CWBackingPixel, 
                  attr);

    XMapWindow(display, win);
    XSync(display, 0);
}

main (int argc, char* argv[])
{
	init_display ();     	 
    /* Calculate and draw points */
    Compl         z, c;
    float         lengthsq, temp;
	int n,t_max;
	FILE* fin = fopen ("nbody.txt","r");
	fscanf (fin, "%d %d",&n,&t_max);
	int n_threads = (argc==2)?atoi(argv[3]):N_THREADS;
	double m[n];
	double f[n][2],v[n][2],p[n][2];
	
	#pragma omp parallel for 
	for (int i = 0;i<n;++i) {
		m[i] = 1;
		f[i][0] = f[i][1] = v[i][0] = v[i][1] = 0;
		fscanf (fin,"%lf %lf",&p[i][0],&p[i][1]);
	}
	double p_new[n][2],v_new[n][2];
	omp_set_num_threads (n_threads);
	struct timeval tv1,tv2;
	struct timezone tz1,tz2;
	gettimeofday (&tv1,&tz1);
	for (int t=0;t<t_max;++t) {
		//cout<<t<<endl;
		#pragma omp parallel for 
		for (int i=0;i<n;++i) {
			for (int j=0;j<n;++j) {
				double r = sqrt (pow (p[j][0]-p[i][0],2) + pow(p[j][1]-p[i][1],2));
				r*=r*r;
				if (r>0) {
					f[i][0]+=G*m[i]*m[j]*(p[j][0]-p[i][0]) / r;
					f[i][1]+=G*m[i]*m[j]*(p[j][1]-p[i][1]) / r;
				}
			}
			v_new[i][0] = v[i][0] + f[i][0]*dt/m[i];
			v_new[i][1] = v[i][1] + f[i][1]*dt/m[i];
			p_new[i][0] = p[i][0]+v_new[i][0]*dt/m[i];
			p_new[i][1] = p[i][1]+v_new[i][1]*dt/m[i];

		}
		#pragma omp barrier

		#pragma omp parallel for 
		for (int i=0;i<n;++i) {
			int x,y;
			x = floor (p[i][0]*X_RESN);
			y = floor (p[i][1]*Y_RESN);
			#pragma omp critical 
			{
				XSetForeground (display, gc, WhitePixel (display,screen));
				XDrawPoint (display, win, gc, x,y);
			}
			v[i][0] = v_new[i][0];
			v[i][1] = v_new[i][1];
			p[i][0] = p_new[i][0];
			p[i][1] = p_new[i][1];
			x = floor (p[i][0]*X_RESN);
			y = floor (p[i][1]*Y_RESN);
			#pragma omp critical 
			{
				XSetForeground (display, gc, BlackPixel (display,screen));
				XDrawPoint (display, win, gc, x, y);
			}
		}
		//sleep(1);
	}
	gettimeofday (&tv2,&tz2);
	printf ("Eslaped time = %ld\n", (tv2.tv_sec-tv1.tv_sec)*1000000 + tv2.tv_usec-tv1.tv_usec);
    XFlush(display);

    /* Program Finished */

}