Last modified on 3 December 2014, at 15:01

Fractals/Computer graphic techniques/2D/optimisation

"premature optimization is the root of all evil" - Donald Knuth in paper "Structured Programming With Go To Statements"

Software optimisationEdit

Inner loopEdit

  • inner loop in wikipedia [1]

Optimisation of numerical computationsEdit

DistanceEdit

Distance approximation

SymmetryEdit

Mirror symmetryEdit

Parameter plane and the Mandelbrot set is symmetric with respect to the x-axis in the plane :[2]

colour(x,y) = colour(x,-y)

Here is a C code that shows how to divide symmetric image (with axis in the middle of its height) into 3 parts :

  • axis of symmetry
  • 2 symmetric parts : above and below axis


// fill array using mirror symmetry of image 
// uses global var :  ...
int FillArray(unsigned char data[] )
{
 
 unsigned char Color; // gray from 0 to 255 
 
// draw axis of symmetry
iy = iyAxisOfSymmetry; 
for(ix=ixMin;ix<=ixMax;++ix) PlotPoint(ix, iy, GiveColor(ix, iy)); 
 
 
// draw symmetric parts : above and below axis 
for(iyAbove = iyAboveMin; iyAbove<=iyAboveMax; ++iyAbove) 
  for(ix=ixMin; ix<=ixMax; ++ix) 
 
  { // above axis compute color and save it to the array
    iy = iyAxisOfSymmetry + iyAbove;
    Color = GiveColor(ix, iy);
    PlotPoint(ix, iy, Color ); 
    // below the axis only copy Color the same as above without computing it 
    PlotPoint(ix, iyAxisOfSymmetry - iyAbove , Color ); 
   }   
 return 0;
}


See also Pascal code ( Delphi) for general case [3]

Rotational symmetryEdit

Dynamical plane and the Julia set have rotational symmetry. It can be used to speed up computing :

// fill array using symmetry of image 
// uses global var :  ...
int FillArraySymmetric(unsigned char data[] )
{
 
 unsigned char Color; // gray from 0 to 255 
 
printf("axis of symmetry \n"); 
iy = iyAxisOfSymmetry; 
 
for(ix=ixMin;ix<=ixMax;++ix) PlotPoint(ix, iy, GiveColor(ix, iy)); 
 
// above and below axis 
for(iyAbove = iyAboveMin; iyAbove<=iyAboveMax; ++iyAbove) 
  {printf(" %d from %d\n", iyAbove, iyAboveMax); //info 
  for(ix=ixMin; ix<=ixMax; ++ix) 
 
  { // above axis compute color and save it to the array
    iy = iyAxisOfSymmetry + iyAbove;
    Color = GiveColor(ix, iy);
    PlotPoint(ix, iy, Color ); 
    // below the axis only copy Color the same as above without computing it 
    PlotPoint(ixMax-ix, iyAxisOfSymmetry - iyAbove , Color ); 
   } 
}  
 return 0;
}

Complex numbersEdit

Computations without explicit definition of complex numbers or without complex type are usually faster.

Parallel computingEdit

In Mandelbrot and Julia sets each point/pixel is calculated independently, so it can be easly divided into an smaller tasks and carried out simultaneously.[4]

Types of parallel computing : Wikipedia description[5]

  • Multi-core (computing)
  • General-purpose computing on graphics processing units (GPGPU):

VectorisationEdit

CEdit

REdit

This image is made with use of vectorised code. It's creation time is very short.

Vectorisation [22][23]


Here is R code of above image with comments: [24]

## based on the R code by Jarek Tuszynski
## http://commons.wikimedia.org/wiki/File:Mandelbrot_Creation_Animation_%28800x600%29.gif

iMax=20; # number of frames and iterations

## world coordinate ( float )
cxmin = -2.2;
cxmax = 1.0;
cymax = 1.2;
cymin = -1.2;

## screen coordinate ( integer)

dx=800; dy=600 ;           # define grid size

## sequences of  values = rows and columns
cxseq = seq(cxmin, cxmax, length.out=dx);
cyseq = seq(cymin, cymax, length.out=dy);
## sequences of  values = rows * columns
cxrseq = rep(cxseq, each = dy);
cyrseq = rep(cyseq,  dx);
C = complex( real=cxrseq, imag = cyrseq); # complex vector
C = matrix(C,dy,dx)       # convert from vector to matrix
# Now C is a matrix of c points coordinates (cx,cy)

# allocate memory for all the frames
F = array(0, dim = c(dy,dx,iMax)) # dy*dx*iMax array filled with zeros

Z = 0                     # initialize Z to zero
for (i in 1:iMax)         # perform iMax iterations
{
  Z = Z^2+C               # iteration; uses n point to compute n+1 point
  
  F[,,i]  = exp(-abs(Z))              # an omitted index is used to represent an entire column or row
     # store magnitude of the complex number, scale it using an exponential function to emphasize fine detail, 
} 


# color and save F array to the file 
library(caTools)          # load library with write.gif function
# mapped to color palette jetColors . Dark red is a very low number, dark blue is a very high number
jetColors = colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))
write.gif(F, "Mandelbrot.gif", col=jetColors, delay=100, transparent=0)  

MulticoreEdit

Image and src code using OpenMP and symmetry

"I always create as many worker threads as I have cores. More than that and your system spends too much time task switching" - Duncan C[25][26]


Using OpenMP :

// fill array using symmetry of image 
// uses global var :  ...
int FillArraySymmetric(unsigned char data[] )
{
 
 unsigned char Color; // gray from 0 to 255 
 
printf("axis of symmetry \n"); 
iy = iyAxisOfSymmetry; 
#pragma omp parallel for schedule(dynamic) private(ix,Color) shared(ixMin,ixMax, iyAxisOfSymmetry)
for(ix=ixMin;ix<=ixMax;++ix) PlotPoint(ix, iy, GiveColor(ix, iy)); 
 
 
/*
The use of ‘shared(variable, variable2) specifies that these variables should be shared among all the threads.
The use of ‘private(variable, variable2)’ specifies that these variables should have a seperate instance in each thread.
*/
 
#pragma omp parallel for schedule(dynamic) private(iyAbove,ix,iy,Color) shared(iyAboveMin, iyAboveMax,ixMin,ixMax, iyAxisOfSymmetry)
 
// above and below axis 
for(iyAbove = iyAboveMin; iyAbove<=iyAboveMax; ++iyAbove) 
  {printf(" %d from %d\n", iyAbove, iyAboveMax); //info 
  for(ix=ixMin; ix<=ixMax; ++ix) 
 
  { // above axis compute color and save it to the array
    iy = iyAxisOfSymmetry + iyAbove;
    Color = GiveColor(ix, iy);
    PlotPoint(ix, iy, Color ); 
    // below the axis only copy Color the same as above without computing it 
    PlotPoint(ixMax-ix, iyAxisOfSymmetry - iyAbove , Color ); 
   } 
}  
 return 0;
}


CellEdit

"The Sony PlayStation 3 is one of the cheapest parallel computers available on the consumer market."[27][28]

GPGPUEdit

Canvas : list of supporting browsers

EffectsEdit

Effects on 2 core CPU ( with 1 thread per core )

Method time [minutes] relative speed
CPU & no optimisation 80 1
CPU & symmetry 39 2
CPU & symmetry & OpenMP 24 3
GPU  ? ( to do )

GPU should 62 time faster than CPU [31]


Using WinXP Intel 2.9GHz CPU (1 CPU used) with a GTX 480 GPU I get the following using 1000x1000 plot with 1000 max iterations : [32]

type time description
gpu 0.07s gpu is a pure CUDA solution on the GPU
gpuarray 3.4s gpuarray uses a numpy-like CUDA wrapper in Python on the GPU
numpy 43.4s numpy is a pure Numpy (C-based) solution on the CPU
python (serial) 1605.6s python is a pure Python solution on the CPU with numpy arrays

ReferencesEdit

  1. Inner loop in wikipedia
  2. wikipedia : Reflection symmetry
  3. Mirror symmetry around x axis at fraktal.republika.pl
  4. embarrassingly parallel problem
  5. wikipedia : Parallel_computers - Classes_of_parallel_computers
  6. Mandelbrot using OpenCl and Parallella
  7. wikipedia : Intel MIC
  8. parallel mandelbrot set (C Code with Message Passing Interface (MPI) library) by Omar U. Florez
  9. Guide into OpenMP: Easy multithreading programming for C++ by Joel Yliluoma
  10. Parallel Mandelbrot with OpenMP by dcfrogle
  11. claudiusmaximus  : exponential mapping and openmp
  12. Part 1: OpenCL™ – Portable Parallelism By manythreads
  13. A Mandelbrot Set on the GPU in Matlab by Loren Shure
  14. progressive-julia-fractal using webgl by Felix Woitzel
  15. glsl sandbox at heroku.com
  16. gcc  : Vector Extensions
  17. First code example using gcc vector support by bert hubert
  18. Mandelbrot calculation using SIMD by EJRH Edmund Horner
  19. The Computer Language Benchmarks Game
  20. Intel Advanced Vector Extensions
  21. ICC and Mandelbrot by Tim Horton
  22. Vectorization (parallel computing)
  23. Matlab - Code Vectorization Guide
  24. mandelbrot set in r by J.D. Long
  25. /programming/what-is-the-general-approach-to-threading-for-2d-plotting/ Fractal forums discussion : What is the general approach to threading for 2d plotting
  26. fractalforums discussion : most-powerful-computer-possible-for-a-reasonable-price/?PHPSESSID=7db41da33f499eb25ef85f46866438f2
  27. Algorithms for Specific Hardware -Mathematical Image Analysis Group Saarland University, Germany
  28. Jenks Parallel Programming Blog
  29. Programming high-performance applications on the Cell BE processor, Part 1: An introduction to Linux on the PLAYSTATION 3
  30. [par4all]
  31. MandelCPU vs MandelGPU Written by David Bucciarelli
  32. [http://wiki.tiker.net/PyCuda/Examples/Mandelbrot | Mandelbrot calculate using GPU, Serial numpy and faster numpy Use to show the speed difference between CPU and GPU calculations by Ian Ozsvald ianozsvald.com July 2010]