O(N*log(resolution)) complexity time, by yours truly. I believe that beats the world's best.
NOTE: this is the most up-to-date-revision posted. last revised 2013.07.05. also available at http://kbtrader.info/nbody.txt
* fixed a minor bug
* added expand and contract (for sparse quad-tree/oct-tree)
* added readWithLocalSubtract for concurrent "all particles minus one particle" calculations (see NOTE in code comments)
Code: Select all
//use unit square
import java.util.*;
import java.util.concurrent.*;
/*
WHAT THIS DOES:
this calculates the total force at a point x,y, where the force is a sum of forces from N points.
where the forces follow the inverse-square law. that is, the force is 1/(distance squared).
examples are gravity and the electro-magnetic force.
HOW TO USE:
*ADDING: to add a point of force, e.g. a charge carrier or a mass, at coordinates x,y, call add(x,y,vals),
where x,y is the coordinate, and vals is an array of force singularities at that point (e.g. electric charge, mass, etc.)
*SUBTRACTING: to subtract a point of force, call subtract(x,y,vals)
*MOVING/CHANGING: to move a point of force or change its force, first subtract the old values from the old spot, then add the new values to the new spot.
*CALCULATING: to find the total force on a point x,y, call readMaxContrib(x,y).
*NOTE: if you are calculating the force on a force point that you have already adding to the system,
* you want to caclulate it _without_ that point. to do that, call instead readWithLocalSubract(x,y,px,py,sub).
* where px and py are the coords of the point you want to subtract, and sub is that point's vals (it's array of force singularities).
* otherwise the calculated force will be singular at that point.
*
variations:
*1) fixed depth
**2) fixed depth w/threshold
***3) fixed depth w/threshold & bins
****4) dynamic depth w/threshold & bins (instead of binning, add or remove depth - bins are there so can be done)
*****5) dynamic compute space partitioning for non-uniform memory architecture; aka distributed n-body
***** to be designed
***** each bottom of a bin extension (thus holding x individual particles) is given a computation unit, with links to the up levels. make sure all levels are done before iterating.
for 4):
*make add and subtract also recurse to the replane
*make read also recurse to the replane
*make adjust x and y scale and translation for recursion.
for replaning, first add a distance scaling and shifting function (scale, xshift,yshift).
shown in the code below are the first three variations
1) use_bins = false, use read() function
2) use_bins = false, use readMaxContrib() function
3) use_bins = true , use readMaxContrib() function
4) all of the above, plus use_replane = true, fill in contract/expand thresholds
*/
public class Plane {
PlaneConfig planeConfig = null;
class PlaneConfig {
//params for variation 4.
ConcurrentLinkedQueue<Plane> planeReuseQueue = new ConcurrentLinkedQueue<Plane>();
boolean use_replane = false;
int contract_items_per_plane = 4; //OR
double contract_value_at_plane = 0.5;
int expand_items_per_plane = 8; //AND
double expand_value_at_plane = 1;
int num_values;
int cascade_depth;
boolean use_bins = false;
double max_contrib = 0.1;
boolean use_min_contrib_frac = false;
double min_contrib_frac = 0.0001;
}
double[][][][] vals = null;
Vector<double[]>[][] bins = null; //need bin members x and y coordinates!
Vector<double[]>[][] bin_coords = null; //need bin members x and y coordinates!
int item_count = 0;
Plane[][] re_plane = null; //need bin members x and y coordinates!
//int max_items_per_plane
double x_up_shift = 0;
double y_up_shift = 0;
double x_absolute_shift = 0;
double y_absolute_shift = 0;
double xy_down_scale = 1;
double xy_up_scale = 1;
double xy_absolute_scale = 1;
double translate_x_down(double x) { return (x - x_up_shift)*xy_down_scale; }
double translate_y_down(double y) { return (y - y_up_shift)*xy_down_scale; }
double translate_x_up(double x) { return x*xy_up_scale + x_up_shift; }
double translate_y_up(double y) { return y*xy_up_scale + y_up_shift; }
double translate_x_absolute(double x) { return x*xy_absolute_scale + x_absolute_shift; }
double translate_y_absolute(double y) { return y*xy_absolute_scale + y_absolute_shift; }
//note: replaning is relative to total value, not contribution.
//highest variation, each compute unit maintains its own map, depth is done relative to the compute unit.
//Queue<Plane> recycle_planes = new Queue<Plane>();
//cascade depth = log(resolution)
public void init(PlaneConfig planeConfig) {
this.planeConfig = planeConfig;
vals = new double[planeConfig.cascade_depth][][][];
int size = 1;
for( int i = 0; i < planeConfig.cascade_depth; i++) {
vals[i] = new double[size][size][planeConfig.num_values];
size *= 2;
}
if( planeConfig.use_bins) {
size /= 2;
bins = new Vector[size][size];
bin_coords = new Vector[size][size];
for( int j = 0; j < size; j++) {
for( int k = 0; k < size; k++) {
bins[j][k] = new Vector<double[]>();
bin_coords[j][k] = new Vector<double[]>();
}
}
}
}
public void expand(int xi, int yi) {
if( re_plane[yi][xi] != null) {
System.err.println("error! plane already expanded!");
return;
}
Plane p = planeConfig.planeReuseQueue.poll();
if( p == null) {
p = new Plane();
p.init(planeConfig);
} else {
p.reset();
}
re_plane[yi][xi] = p;
double scale = 1.0/(double)(0x01 << (planeConfig.cascade_depth-1));
p.xy_up_scale = scale;
p.xy_down_scale = 1.0/scale;
p.xy_absolute_scale = p.xy_up_scale*this.xy_absolute_scale;
p.x_absolute_shift = this.x_absolute_shift + p.xy_absolute_scale * (double)xi;
p.y_absolute_shift = this.y_absolute_shift + p.xy_absolute_scale * (double)yi;
p.x_up_shift = scale * (double)xi;
p.y_up_shift = scale * (double)yi;
Vector<double[]> bin = bins[yi][xi];
Vector<double[]> bin_coord = bin_coords[yi][xi];
for( int i = 0; i < bin.size(); i++) {
double[] xx = bin_coord.get(i);
double[] vv = bin.get(i);
p.add(xx[0], xx[1], vv);
}
}
public void contract(int ix, int iy) {
Plane p = re_plane[iy][ix];
re_plane[iy][ix] = null;
planeConfig.planeReuseQueue.offer(p);
}
public void reset() {
int size = 1;
for( int i = 0; i < planeConfig.cascade_depth; i++) {
for( int j = 0; j < size; j++) {
for( int k = 0; k < size; k++) {
for( int l = 0; l < planeConfig.num_values; l++) {
vals[i][j][k][l] = 0;
}
}
}
size *= 2;
}
if( planeConfig.use_bins) {
size /= 2;
for( int j = 0; j < size; j++) {
for( int k = 0; k < size; k++) {
bins[j][k] = new Vector<double[]>();
bin_coords[j][k] = new Vector<double[]>();
re_plane[j][k] = null;
}
}
}
}
public void subtract(double x, double y, double[] _vals) {
x = translate_x_down(x);
y = translate_y_down(y);
addScaled(x,y,_vals,-1);
if( planeConfig.use_bins) {
item_count--;
int depth = planeConfig.cascade_depth-1;
double scale = (double)(0x01 << (depth));
int xi0 = (int)(x*scale);
int yi0 = (int)(y*scale);
Vector<double[]> bin = bins[yi0][xi0];
Vector<double[]> bin_coord = bin_coords[yi0][xi0];
int ndx = bin.indexOf(_vals);
bin.remove(ndx);
bin_coord.remove(ndx);
Plane p = re_plane[yi0][xi0];
if( p != null) {
p.subtract(x,y,_vals);
boolean ok = true;
if( item_count > planeConfig.contract_items_per_plane) {
for( int i = 0; i < planeConfig.num_values; i++) {
if ( Math.abs(vals[depth][yi0][xi0][i]) > planeConfig.contract_value_at_plane) {
ok = false;
break;
}
}
}
if( ok)
contract(xi0,yi0);
}
}
}
public void add(double x, double y, double[] _vals) {
x = translate_x_down(x);
y = translate_y_down(y);
double xi = x;
double yi = y;
for( int i = 0; i < planeConfig.cascade_depth; i++) {
int xi0 = (int)Math.floor(xi);
int yi0 = (int)Math.floor(yi);
for( int j = 0; j < planeConfig.num_values; j++)
vals[i][yi0][xi0][j] += _vals[j];
xi *= 2;
yi *= 2;
}
if( planeConfig.use_bins) {
item_count++;
int depth = planeConfig.cascade_depth-1;
double scale = (double)(0x01 << (depth));
int xi0 = (int)(x*scale);
int yi0 = (int)(y*scale);
Vector<double[]> bin = bins[yi0][xi0];
bin.add(_vals);
Vector<double[]> bin_coord = bin_coords[yi0][xi0];
bin_coord.add(new double[]{x,y});
Plane p = re_plane[yi0][xi0];
if( p != null)
p.add(x,y,_vals);
if( p == null) {
boolean ok = false;
if( item_count > planeConfig.expand_items_per_plane) {
for( int i = 0; i < planeConfig.num_values; i++) {
if ( Math.abs(vals[depth][yi0][xi0][i]) > planeConfig.expand_value_at_plane) {
ok = true;
break;
}
}
}
if( ok)
expand(xi0,yi0);
}
}
}
public void addScaled(double x, double y, double[] _vals, double scale) {
double xi = x;
double yi = y;
if( scale != -1 && planeConfig.use_bins) {
System.err.println("don't use addScaled if you are using bins! instead subtract then add!");
return;
}
for( int i = 0; i < planeConfig.cascade_depth; i++) {
int xi0 = (int)Math.floor(xi);
int yi0 = (int)Math.floor(yi);
for( int j = 0; j < planeConfig.num_values; j++)
vals[i][yi0][xi0][j] += _vals[j]*scale;
xi *= 2;
yi *= 2;
}
}
public double[] read(double x, double y) {
double[] ret = new double[planeConfig.num_values];
read(x,y,0,0,0,ret);
return ret;
}
//enclosed volume must be less than squared distance at closest approach.
//otherwise, recurse in a level.
private void read(double x, double y, int depth, int xi, int yi, double[] _vals) {
double d = getSquaredDistance(x,y,depth,xi,yi);
if( d > getEnclosedVolume(depth) || depth == planeConfig.cascade_depth-1) {
for( int i = 0; i < planeConfig.num_values; i++) {
_vals[i] += vals[depth][yi][xi][i] / d;
}
} else {
xi *= 2;
yi *= 2;
depth++;
read(x,y,depth,xi+0,yi+0,_vals);
read(x,y,depth,xi+1,yi+0,_vals);
read(x,y,depth,xi+0,yi+1,_vals);
read(x,y,depth,xi+1,yi+1,_vals);
}
}
public double[] readMaxContrib(double x, double y) {
double[] ret = new double[planeConfig.num_values];
readMaxContrib(x,y,0,0,0,ret);
return ret;
}
public double[] readWithLocalSubtract(double x, double y, double px, double py, double[] sub) {
double[] ret = new double[planeConfig.num_values];
readWithLocalSubtract(x,y,0,0,0,ret,px,py,sub);
return ret;
}
private void readWithLocalSubtract(double x, double y, int depth, int xi, int yi, double[] _vals, double px, double py, double[] sub) {
double d = getSquaredDistance(x,y,depth,xi,yi);
double max = planeConfig.max_contrib * d;
double min_frac = planeConfig.min_contrib_frac * d;
boolean containsSubtractPoint = containsAbsolutePoint( depth, yi, xi, px, py);
//if it doesn't contain the subtract point, then we don't need to check that anymore going forward.
if( !containsSubtractPoint) {
readMaxContrib(x,y,depth,xi,yi,_vals);
return;
}
boolean ok = true;
//if using bins, then at max depth, add individual particles (becomes n^2 at max depth)
if( planeConfig.use_bins && depth == planeConfig.cascade_depth-1) {
if( re_plane[yi][xi] != null) {
re_plane[yi][xi].readWithLocalSubtract(x,y,0,0,0,_vals,px,py,sub);
return;
}
Vector<double[]> bin = bins[yi][xi];
Vector<double[]> bin_coord = bin_coords[yi][xi];
for( int i = 0; i < bin.size(); i++) {
double[] xx = bin_coord.get(i);
//skip if is the subtract point
if( xx[0] == px && xx[1] == py) {
continue;
}
double[] vv = bin.get(i);
double xd = xx[0]-x;
double yd = xx[1]-y;
double rsd = 1.0/(xd*xd+yd*yd);
for( int j = 0; j < _vals.length; j++) {
_vals[j] += vv[j]*rsd;
}
}
return;
}
if( depth != planeConfig.cascade_depth-1) {
for( int i = 0; i < planeConfig.num_values; i++) {
if( planeConfig.use_min_contrib_frac) {
if ( Math.abs((vals[depth][yi][xi][i]-sub[i])) < min_frac*_vals[i]) {
continue;
}
}
if ( Math.abs((vals[depth][yi][xi][i]-sub[i])) > max) {
ok = false;
break;
}
}
}
if( ok) {
for( int i = 0; i < planeConfig.num_values; i++) {
_vals[i] += (vals[depth][yi][xi][i]-sub[i]) / d;
}
return;
}
xi *= 2;
yi *= 2;
depth++;
readWithLocalSubtract(x,y,depth,xi+0,yi+0,_vals,px,py,sub);
readWithLocalSubtract(x,y,depth,xi+1,yi+0,_vals,px,py,sub);
readWithLocalSubtract(x,y,depth,xi+0,yi+1,_vals,px,py,sub);
readWithLocalSubtract(x,y,depth,xi+1,yi+1,_vals,px,py,sub);
}
public boolean containsAbsolutePoint(int depth, int yi, int xi, double px, double py) {
double mult = 1.0/(double)(0x01 << (depth));
double x0 = mult*(double)xi;
double y0 = mult*(double)yi;
double x1 = x0+mult;
double y1 = y0+mult;
x0 = translate_x_absolute(x0);
y0 = translate_y_absolute(y0);
x1 = translate_x_absolute(x1);
y1 = translate_y_absolute(y1);
if( px < x0 || px > x1 || py < y0 || py > y1)
return false;
return true;
}
//enclosed volume must be less than squared distance at closest approach.
//otherwise, recurse in a level.
//make read with local subtract function
private void readMaxContrib(double x, double y, int depth, int xi, int yi, double[] _vals) {
double d = getSquaredDistance(x,y,depth,xi,yi);
double max = planeConfig.max_contrib * d;
double min_frac = planeConfig.min_contrib_frac * d;
boolean ok = true;
//if using bins, then at max depth, add individual particles (becomes n^2 at max depth)
if( planeConfig.use_bins && depth == planeConfig.cascade_depth-1) {
if( re_plane[yi][xi] != null) {
re_plane[yi][xi].readMaxContrib(x,y,0,0,0,_vals);
return;
}
Vector<double[]> bin = bins[yi][xi];
Vector<double[]> bin_coord = bin_coords[yi][xi];
for( int i = 0; i < bin.size(); i++) {
double[] xx = bin_coord.get(i);
double[] vv = bin.get(i);
double xd = xx[0]-x;
double yd = xx[1]-y;
double rsd = 1.0/(xd*xd+yd*yd);
for( int j = 0; j < _vals.length; j++) {
_vals[j] += vv[j]*rsd;
}
}
return;
}
if( depth != planeConfig.cascade_depth-1) {
for( int i = 0; i < planeConfig.num_values; i++) {
if( planeConfig.use_min_contrib_frac) {
if ( Math.abs(vals[depth][yi][xi][i]) < min_frac*_vals[i]) {
continue;
}
}
if ( Math.abs(vals[depth][yi][xi][i]) > max) {
ok = false;
break;
}
}
}
if( ok) {
for( int i = 0; i < planeConfig.num_values; i++) {
_vals[i] += vals[depth][yi][xi][i] / d;
}
return;
}
xi *= 2;
yi *= 2;
depth++;
readMaxContrib(x,y,depth,xi+0,yi+0,_vals);
readMaxContrib(x,y,depth,xi+1,yi+0,_vals);
readMaxContrib(x,y,depth,xi+0,yi+1,_vals);
readMaxContrib(x,y,depth,xi+1,yi+1,_vals);
}
private double getSquaredDistance(double x, double y, int depth, int xi, int yi) {
double mult = 1.0/(double)(0x01 << (depth));
double x0 = mult*(double)xi + mult/2;
double y0 = mult*(double)yi + mult/2;
double xd = x - translate_x_absolute(x0);
double yd = y - translate_y_absolute(y0);
return xd*xd + yd*yd;
}
private double getEnclosedVolume(int depth) {
return xy_absolute_scale*xy_absolute_scale/(double)(0x01 << (depth*2));
}
}