As many (or very few in the real life haha) people know, today is the Pi Approximation Day ! So it’s time to make a contribution to celebrate this funny day =)

My contribution is to use Python and Pyevolve to approximate Pi number using Genetic Programming approach. I’ve created the functions gp_add(+), gp_sub(-), gp_div(/), gp_mul(*) and gp_sqrt (square root) to use as non-terminals of the GP. The fitness function is very simple too, it simple returns the absolute difference between the Python math.pi and the evaluated individual. I’ve used also a population size of 1k individuals with max tree depth of 8 and the random ephemeral constants as random integers. The best approximation I’ve got while running the GP for about 8 minutes (40 generations) was 3.1416185511, best for 3 digits, you can improve it and let it run for more time to get better approximations.

Here is the formulae I’ve got with the GP (click to enlarge):

And here is the output of the script:

Best (0): 3.1577998365 Error: 0.0162071829 Best (10): 3.1417973679 Error: 0.0002047143 Best (20): 3.1417973679 Error: 0.0002047143 Best (30): 3.1417973679 Error: 0.0002047143 Best (40): 3.1416185511 Error: 0.0000258975 - GenomeBase Score: 0.000026 Fitness: 15751.020831 Params: {'max_depth': 8, 'method': 'ramped'} Slot [Evaluator] (Count: 1) Slot [Initializator] (Count: 1) Name: GTreeGPInitializator - Weight: 0.50 Doc: This initializator accepts the follow parameters: *max_depth* The max depth of the tree *method* The method, accepts "grow" or "full" .. versionadded:: 0.6 The *GTreeGPInitializator* function. Slot [Mutator] (Count: 1) Name: GTreeGPMutatorSubtree - Weight: 0.50 Doc: The mutator of GTreeGP, Subtree Mutator .. versionadded:: 0.6 The *GTreeGPMutatorSubtree* function Slot [Crossover] (Count: 1) Name: GTreeGPCrossoverSinglePoint - Weight: 0.50 - GTree Height: 8 Nodes: 21 GTreeNodeBase [Childs=1] - [gp_sqrt] GTreeNodeBase [Childs=2] - [gp_div] GTreeNodeBase [Childs=2] - [gp_add] GTreeNodeBase [Childs=0] - [26] GTreeNodeBase [Childs=2] - [gp_div] GTreeNodeBase [Childs=2] - [gp_mul] GTreeNodeBase [Childs=2] - [gp_add] GTreeNodeBase [Childs=2] - [gp_sub] GTreeNodeBase [Childs=0] - [34] GTreeNodeBase [Childs=2] - [gp_sub] GTreeNodeBase [Childs=0] - [44] GTreeNodeBase [Childs=0] - [1] GTreeNodeBase [Childs=2] - [gp_mul] GTreeNodeBase [Childs=0] - [49] GTreeNodeBase [Childs=0] - [43] GTreeNodeBase [Childs=1] - [gp_sqrt] GTreeNodeBase [Childs=0] - [18] GTreeNodeBase [Childs=0] - [16] GTreeNodeBase [Childs=2] - [gp_add] GTreeNodeBase [Childs=0] - [24] GTreeNodeBase [Childs=0] - [35] - GTreeGP Expression: gp_sqrt(gp_div(gp_add(26, gp_div(gp_mul(gp_add(gp_sub(34, gp_sub(44, 1)), gp_mul(49, 43)), gp_sqrt(18)), 16)), gp_add(24, 35)))

And finally, here is the source code:

from __future__ import division from pyevolve import * import math def gp_add(a, b): return a+b def gp_sub(a, b): return a-b def gp_div(a, b): return 1 if b==0 else a/b def gp_mul(a, b): return a*b def gp_sqrt(a): return math.sqrt(abs(a)) def eval_func(chromosome): code_comp = chromosome.getCompiledCode() ret = eval(code_comp) return abs(math.pi - ret) def step_callback(engine): gen = engine.getCurrentGeneration() if gen % 10 == 0: best = engine.bestIndividual() best_pi = eval(best.getCompiledCode()) print "Best (%d): %.10f" % (gen, best_pi) print "\tError: %.10f" % (abs(math.pi - best_pi)) return False def main_run(): genome = GTree.GTreeGP() genome.setParams(max_depth=8, method="ramped") genome.evaluator += eval_func ga = GSimpleGA.GSimpleGA(genome) ga.setParams(gp_terminals = ['ephemeral:random.randint(1, 50)'], gp_function_prefix = "gp") ga.setMinimax(Consts.minimaxType["minimize"]) ga.setGenerations(50000) ga.setCrossoverRate(1.0) ga.setMutationRate(0.09) ga.setPopulationSize(1000) ga.stepCallback.set(step_callback) ga.evolve() best = ga.bestIndividual() best.writeDotImage("tree_pi.png") print best if __name__ == "__main__": main_run()

If you are interested why today is the Pi Approximation day, see some resources:

Little Cartoon

Some Background History

Some Pi Approximations