محاسبات تکاملی چیست؟ – به زبان ساده

1function varargout = tsp_ga(varargin)
2    
3    % Initialize default configuration
4    defaultConfig.xy          = 10*rand(50,2);
5    defaultConfig.dmat        = [];
6    defaultConfig.popSize     = 100;
7    defaultConfig.numIter     = 1e4;
8    defaultConfig.showProg    = true;
9    defaultConfig.showResult  = true;
10    defaultConfig.showWaitbar = false;
11    
12    % Interpret user configuration inputs
13    if ~nargin
14        userConfig = struct();
15    elseif isstruct(varargin{1})
16        userConfig = varargin{1};
17    else
18        try
19            userConfig = struct(varargin{:});
20        catch
21            error('Expected inputs are either a structure or parameter/value pairs');
22        end
23    end
24    
25    % Override default configuration with user inputs
26    configStruct = get_config(defaultConfig,userConfig);
27    
28    % Extract configuration
29    xy          = configStruct.xy;
30    dmat        = configStruct.dmat;
31    popSize     = configStruct.popSize;
32    numIter     = configStruct.numIter;
33    showProg    = configStruct.showProg;
34    showResult  = configStruct.showResult;
35    showWaitbar = configStruct.showWaitbar;
36    if isempty(dmat)
37        nPoints = size(xy,1);
38        a = meshgrid(1:nPoints);
39        dmat = reshape(sqrt(sum((xy(a,:)-xy(a',:)).^2,2)),nPoints,nPoints);
40    end
41    
42    % Verify Inputs
43    [N,dims] = size(xy);
44    [nr,nc] = size(dmat);
45    if N ~= nr || N ~= nc
46        error('Invalid XY or DMAT inputs!')
47    end
48    n = N;
49    
50    % Sanity Checks
51    popSize     = 4*ceil(popSize/4);
52    numIter     = max(1,round(real(numIter(1))));
53    showProg    = logical(showProg(1));
54    showResult  = logical(showResult(1));
55    showWaitbar = logical(showWaitbar(1));
56    
57    % Initialize the Population
58    pop = zeros(popSize,n);
59    pop(1,:) = (1:n);
60    for k = 2:popSize
61        pop(k,:) = randperm(n);
62    end
63    
64    % Run the GA
65    globalMin = Inf;
66    totalDist = zeros(1,popSize);
67    distHistory = zeros(1,numIter);
68    tmpPop = zeros(4,n);
69    newPop = zeros(popSize,n);
70    if showProg
71        figure('Name','TSP_GA | Current Best Solution','Numbertitle','off');
72        hAx = gca;
73    end
74    if showWaitbar
75        hWait = waitbar(0,'Searching for near-optimal solution ...');
76    end
77    for iter = 1:numIter
78        % Evaluate Each Population Member (Calculate Total Distance)
79        for p = 1:popSize
80            d = dmat(pop(p,n),pop(p,1)); % Closed Path
81            for k = 2:n
82                d = d + dmat(pop(p,k-1),pop(p,k));
83            end
84            totalDist(p) = d;
85        end
86        
87        % Find the Best Route in the Population
88        [minDist,index] = min(totalDist);
89        distHistory(iter) = minDist;
90        if minDist < globalMin
91            globalMin = minDist;
92            optRoute = pop(index,:);
93            if showProg
94                % Plot the Best Route
95                rte = optRoute([1:n 1]);
96                if dims > 2, plot3(hAx,xy(rte,1),xy(rte,2),xy(rte,3),'r.-');
97                else plot(hAx,xy(rte,1),xy(rte,2),'r.-'); end
98                title(hAx,sprintf('Total Distance = %1.4f, Iteration = %d',minDist,iter));
99                drawnow;
100            end
101        end
102        
103        % Genetic Algorithm Operators
104        randomOrder = randperm(popSize);
105        for p = 4:4:popSize
106            rtes = pop(randomOrder(p-3:p),:);
107            dists = totalDist(randomOrder(p-3:p));
108            [ignore,idx] = min(dists); %#ok
109            bestOf4Route = rtes(idx,:);
110            routeInsertionPoints = sort(ceil(n*rand(1,2)));
111            I = routeInsertionPoints(1);
112            J = routeInsertionPoints(2);
113            for k = 1:4 % Mutate the Best to get Three New Routes
114                tmpPop(k,:) = bestOf4Route;
115                switch k
116                    case 2 % Flip
117                        tmpPop(k,I:J) = tmpPop(k,J:-1:I);
118                    case 3 % Swap
119                        tmpPop(k,[I J]) = tmpPop(k,[J I]);
120                    case 4 % Slide
121                        tmpPop(k,I:J) = tmpPop(k,[I+1:J I]);
122                    otherwise % Do Nothing
123                end
124            end
125            newPop(p-3:p,:) = tmpPop;
126        end
127        pop = newPop;
128        
129        % Update the waitbar
130        if showWaitbar && ~mod(iter,ceil(numIter/325))
131            waitbar(iter/numIter,hWait);
132        end
133        
134    end
135    if showWaitbar
136        close(hWait);
137    end
138    
139    if showResult
140        % Plots the GA Results
141        figure('Name','TSP_GA | Results','Numbertitle','off');
142        subplot(2,2,1);
143        pclr = ~get(0,'DefaultAxesColor');
144        if dims > 2, plot3(xy(:,1),xy(:,2),xy(:,3),'.','Color',pclr);
145        else plot(xy(:,1),xy(:,2),'.','Color',pclr); end
146        title('City Locations');
147        subplot(2,2,2);
148        imagesc(dmat(optRoute,optRoute));
149        title('Distance Matrix');
150        subplot(2,2,3);
151        rte = optRoute([1:n 1]);
152        if dims > 2, plot3(xy(rte,1),xy(rte,2),xy(rte,3),'r.-');
153        else plot(xy(rte,1),xy(rte,2),'r.-'); end
154        title(sprintf('Total Distance = %1.4f',minDist));
155        subplot(2,2,4);
156        plot(distHistory,'b','LineWidth',2);
157        title('Best Solution History');
158        set(gca,'XLim',[0 numIter+1],'YLim',[0 1.1*max([1 distHistory])]);
159    end
160    
161    % Return Output
162    if nargout
163        resultStruct = struct( ...
164            'xy',          xy, ...
165            'dmat',        dmat, ...
166            'popSize',     popSize, ...
167            'numIter',     numIter, ...
168            'showProg',    showProg, ...
169            'showResult',  showResult, ...
170            'showWaitbar', showWaitbar, ...
171            'optRoute',    optRoute, ...
172            'minDist',     minDist);
173        
174        varargout = {resultStruct};
175    end
176    
177end
178
179% Subfunction to override the default configuration with user inputs
180function config = get_config(defaultConfig,userConfig)
181    
182    % Initialize the configuration structure as the default
183    config = defaultConfig;
184    
185    % Extract the field names of the default configuration structure
186    defaultFields = fieldnames(defaultConfig);
187    
188    % Extract the field names of the user configuration structure
189    userFields = fieldnames(userConfig);
190    nUserFields = length(userFields);
191    
192    % Override any default configuration fields with user values
193    for i = 1:nUserFields
194        userField = userFields{i};
195        isField = strcmpi(defaultFields,userField);
196        if nnz(isField) == 1
197            thisField = defaultFields{isField};
198            config.(thisField) = userConfig.(userField);
199        end
200    end
201    
202end

1import numpy as np, random, operator, pandas as pd, matplotlib.pyplot as plt
2
3class City:
4    def __init__(self, x, y):
5        self.x = x
6        self.y = y
7    
8    def distance(self, city):
9        xDis = abs(self.x - city.x)
10        yDis = abs(self.y - city.y)
11        distance = np.sqrt((xDis ** 2) + (yDis ** 2))
12        return distance
13    
14    def __repr__(self):
15        return "(" + str(self.x) + "," + str(self.y) + ")"
16
17class Fitness:
18    def __init__(self, route):
19        self.route = route
20        self.distance = 0
21        self.fitness= 0.0
22    
23    def routeDistance(self):
24        if self.distance ==0:
25            pathDistance = 0
26            for i in range(0, len(self.route)):
27                fromCity = self.route[i]
28                toCity = None
29                if i + 1 < len(self.route):
30                    toCity = self.route[i + 1]
31                else:
32                    toCity = self.route[0]
33                pathDistance += fromCity.distance(toCity)
34            self.distance = pathDistance
35        return self.distance
36    
37    def routeFitness(self):
38        if self.fitness == 0:
39            self.fitness = 1 / float(self.routeDistance())
40        return self.fitness
41
42def createRoute(cityList):
43    route = random.sample(cityList, len(cityList))
44    return route
45
46def initialPopulation(popSize, cityList):
47    population = []
48
49    for i in range(0, popSize):
50        population.append(createRoute(cityList))
51    return population
52
53def rankRoutes(population):
54    fitnessResults = {}
55    for i in range(0,len(population)):
56        fitnessResults[i] = Fitness(population[i]).routeFitness()
57    return sorted(fitnessResults.items(), key = operator.itemgetter(1), reverse = True)
58
59def selection(popRanked, eliteSize):
60    selectionResults = []
61    df = pd.DataFrame(np.array(popRanked), columns=["Index","Fitness"])
62    df['cum_sum'] = df.Fitness.cumsum()
63    df['cum_perc'] = 100*df.cum_sum/df.Fitness.sum()
64    
65    for i in range(0, eliteSize):
66        selectionResults.append(popRanked[i][0])
67    for i in range(0, len(popRanked) - eliteSize):
68        pick = 100*random.random()
69        for i in range(0, len(popRanked)):
70            if pick <= df.iat[i,3]:
71                selectionResults.append(popRanked[i][0])
72                break
73    return selectionResults
74
75def matingPool(population, selectionResults):
76    matingpool = []
77    for i in range(0, len(selectionResults)):
78        index = selectionResults[i]
79        matingpool.append(population[index])
80    return matingpool
81
82def breed(parent1, parent2):
83    child = []
84    childP1 = []
85    childP2 = []
86    
87    geneA = int(random.random() * len(parent1))
88    geneB = int(random.random() * len(parent1))
89    
90    startGene = min(geneA, geneB)
91    endGene = max(geneA, geneB)
92
93    for i in range(startGene, endGene):
94        childP1.append(parent1[i])
95        
96    childP2 = [item for item in parent2 if item not in childP1]
97
98    child = childP1 + childP2
99    return child
100
101def breedPopulation(matingpool, eliteSize):
102    children = []
103    length = len(matingpool) - eliteSize
104    pool = random.sample(matingpool, len(matingpool))
105
106    for i in range(0,eliteSize):
107        children.append(matingpool[i])
108    
109    for i in range(0, length):
110        child = breed(pool[i], pool[len(matingpool)-i-1])
111        children.append(child)
112    return children
113
114def mutate(individual, mutationRate):
115    for swapped in range(len(individual)):
116        if(random.random() < mutationRate):
117            swapWith = int(random.random() * len(individual))
118            
119            city1 = individual[swapped]
120            city2 = individual[swapWith]
121            
122            individual[swapped] = city2
123            individual[swapWith] = city1
124    return individual
125
126def mutatePopulation(population, mutationRate):
127    mutatedPop = []
128    
129    for ind in range(0, len(population)):
130        mutatedInd = mutate(population[ind], mutationRate)
131        mutatedPop.append(mutatedInd)
132    return mutatedPop
133
134def nextGeneration(currentGen, eliteSize, mutationRate):
135    popRanked = rankRoutes(currentGen)
136    selectionResults = selection(popRanked, eliteSize)
137    matingpool = matingPool(currentGen, selectionResults)
138    children = breedPopulation(matingpool, eliteSize)
139    nextGeneration = mutatePopulation(children, mutationRate)
140    return nextGeneration
141
142def geneticAlgorithm(population, popSize, eliteSize, mutationRate, generations):
143    pop = initialPopulation(popSize, population)
144    print("Initial distance: " + str(1 / rankRoutes(pop)[0][1]))
145    
146    for i in range(0, generations):
147        pop = nextGeneration(pop, eliteSize, mutationRate)
148    
149    print("Final distance: " + str(1 / rankRoutes(pop)[0][1]))
150    bestRouteIndex = rankRoutes(pop)[0][0]
151    bestRoute = pop[bestRouteIndex]
152    return bestRoute
153
154cityList = []
155
156for i in range(0,25):
157    cityList.append(City(x=int(random.random() * 200), y=int(random.random() * 200)))
158
159geneticAlgorithm(population=cityList, popSize=100, eliteSize=20, mutationRate=0.01, generations=500)
160
161def geneticAlgorithmPlot(population, popSize, eliteSize, mutationRate, generations):
162    pop = initialPopulation(popSize, population)
163    progress = []
164    progress.append(1 / rankRoutes(pop)[0][1])
165    
166    for i in range(0, generations):
167        pop = nextGeneration(pop, eliteSize, mutationRate)
168        progress.append(1 / rankRoutes(pop)[0][1])
169    
170    plt.plot(progress)
171    plt.ylabel('Distance')
172    plt.xlabel('Generation')
173    plt.show()
174
175geneticAlgorithmPlot(population=cityList, popSize=100, eliteSize=20, mutationRate=0.01, generations=500)

1%This class impliments a string that mutates to a target
2classdef EvolutionaryAlgorithm
3 
4    properties
5 
6        target;
7        parent;
8        children = {};
9        validAlphabet;
10 
11        %Constants
12        numChildrenPerIteration;
13        maxIterations;
14        mutationRate;
15 
16    end
17 
18    methods
19 
20        %Class constructor
21        function family = EvolutionaryAlgorithm(target,mutationRate,numChildren,maxIterations)
22 
23            family.validAlphabet = char([32 (65:90)]); %Space char and A-Z
24            family.target = target;
25            family.children = cell(numChildren,1);
26            family.numChildrenPerIteration = numChildren;
27            family.maxIterations = maxIterations;
28            family.mutationRate = mutationRate;
29            initialize(family);
30 
31        end %class constructor
32 
33        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
34        %Helper functions and class get/set functions
35        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
36 
37        %setAlphabet() - sets the valid alphabet for the current instance
38        %of the EvolutionaryAlgorithm class.
39        function setAlphabet(family,alphabet)
40 
41            if(ischar(alphabet))
42                family.validAlphabet = alphabet;
43 
44                %Makes change permanent
45                assignin('caller',inputname(1),family); 
46            else
47                error 'New alphabet must be a string or character array';
48            end            
49 
50        end
51 
52        %setTarget() - sets the target for the current instance
53        %of the EvolutionaryAlgorithm class.
54        function setTarget(family,target)
55 
56            if(ischar(target))
57                family.target = target;
58 
59                %Makes change permanent
60                assignin('caller',inputname(1),family); 
61            else
62                error 'New target must be a string or character array';
63            end            
64 
65        end
66 
67        %setMutationRate() - sets the mutation rate for the current instance
68        %of the EvolutionaryAlgorithm class.
69        function setMutationRate(family,mutationRate)
70 
71            if(isnumeric(mutationRate))
72                family.mutationRate = mutationRate;
73 
74                %Makes change permanent
75                assignin('caller',inputname(1),family); 
76            else
77                error 'New mutation rate must be a double precision number';
78            end            
79 
80        end
81 
82        %setMaxIterations() - sets the maximum number of iterations during
83        %evolution for the current instance of the EvolutionaryAlgorithm class.
84        function setMaxIterations(family,maxIterations)
85 
86            if(isnumeric(maxIterations))
87                family.maxIterations = maxIterations;
88 
89                %Makes change permanent
90                assignin('caller',inputname(1),family); 
91            else
92                error 'New maximum amount of iterations must be a double precision number';
93            end            
94 
95        end
96 
97        %display() - overrides the built-in MATLAB display() function, to
98        %display the important class variables
99        function display(family)
100            disp([sprintf('Target: %s\n',family.target)...
101                  sprintf('Parent: %s\n',family.parent)...
102                  sprintf('Valid Alphabet: %s\n',family.validAlphabet)...
103                  sprintf('Number of Children: %d\n',family.numChildrenPerIteration)...
104                  sprintf('Mutation Rate [0,1]: %d\n',family.mutationRate)...
105                  sprintf('Maximum Iterations: %d\n',family.maxIterations)]);                 
106        end
107 
108        %disp() - overrides the built-in MATLAB disp() function, to
109        %display the important class variables
110        function disp(family)
111            display(family);
112        end
113 
114        %randAlphabetElement() - Generates a random character from the
115        %valid alphabet for the current instance of the class.
116        function elements = randAlphabetElements(family,numChars)
117 
118            %Sample the valid alphabet randomly from the uniform
119            %distribution
120            N = length(family.validAlphabet);
121            choices = ceil(N*rand(1,numChars));
122 
123            elements = family.validAlphabet(choices);
124 
125        end
126 
127        %initialize() - Sets the parent to a random string of length equal
128        %to the length of the target
129        function parent = initialize(family)
130 
131            family.parent = randAlphabetElements(family,length(family.target));
132            parent = family.parent;
133 
134            %Makes changes to the instance of EvolutionaryAlgorithm permanent
135            assignin('caller',inputname(1),family); 
136 
137        end %initialize
138 
139        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
140        %Functions required by task specification
141        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
142 
143        %mutate() - generates children from the parent and mutates them
144        function mutate(family)
145 
146            sizeParent = length(family.parent);
147 
148            %Generate mutatant children sequentially
149            for child = (1:family.numChildrenPerIteration)
150 
151                parentCopy = family.parent;
152 
153                for charIndex = (1:sizeParent) 
154                    if (rand(1) < family.mutationRate)
155                        parentCopy(charIndex) = randAlphabetElements(family,1);
156                    end
157                end
158 
159                family.children{child} = parentCopy;
160 
161            end
162 
163            %Makes changes to the instance of EvolutionaryAlgorithm permanent
164            assignin('caller',inputname(1),family);  
165 
166        end %mutate
167 
168        %fitness() - Computes the Hamming distance between the target
169        %string and the string input as the familyMember argument
170        function theFitness = fitness(family,familyMember)
171 
172            if not(ischar(familyMember))
173                error 'The second argument must be a string';
174            end
175 
176            theFitness = sum(family.target == familyMember);
177        end
178 
179        %evolve() - evolves the family until the target is reached or it 
180        %exceeds the maximum amount of iterations        
181        function [iteration,mostFitFitness] = evolve(family)
182 
183            iteration = 0;
184            mostFitFitness = 0;
185            targetFitness = fitness(family,family.target);
186 
187            disp(['Target fitness is ' num2str(targetFitness)]);
188 
189            while (mostFitFitness < targetFitness) && (iteration < family.maxIterations)
190 
191                iteration = iteration + 1;
192 
193                mutate(family);
194 
195                parentFitness = fitness(family,family.parent);                
196                mostFit = family.parent;
197                mostFitFitness = parentFitness;
198 
199                for child = (1:family.numChildrenPerIteration)
200 
201                    childFitness = fitness(family,family.children{child});
202                    if childFitness > mostFitFitness
203                        mostFit = family.children{child};
204                        mostFitFitness = childFitness;
205                    end
206 
207                end               
208 
209                family.parent = mostFit;
210                disp([num2str(iteration) ': ' mostFit ' - Fitness: ' num2str(mostFitFitness)]);
211 
212            end
213 
214            %Makes changes to the instance of EvolutionaryAlgorithm permanent
215            assignin('caller',inputname(1),family);
216 
217        end %evolve
218 
219    end %methods
220end %classdef

1>> instance = EvolutionaryAlgorithm('METHINKS IT IS LIKE A WEASEL',.08,50,1000)
2Target: METHINKS IT IS LIKE A WEASEL
3Parent: UVEOCXXFBGDCSFNMJQNWTPJ PCVA
4Valid Alphabet:  ABCDEFGHIJKLMNOPQRSTUVWXYZ
5Number of Children: 50
6Mutation Rate [0,1]: 8.000000e-002
7Maximum Iterations: 1000
8 
9>> evolve(instance);

1from string import letters
2from random import choice, random
3 
4target  = list("METHINKS IT IS LIKE A WEASEL")
5charset = letters + ' '
6parent  = [choice(charset) for _ in range(len(target))]
7minmutaterate  = .09
8C = range(100)
9 
10perfectfitness = float(len(target))
11 
12def fitness(trial):
13    'Sum of matching chars by position'
14    return sum(t==h for t,h in zip(trial, target))
15 
16def mutaterate():
17    'Less mutation the closer the fit of the parent'
18    return 1-((perfectfitness - fitness(parent)) / perfectfitness * (1 - minmutaterate))
19 
20def mutate(parent, rate):
21    return [(ch if random() <= rate else choice(charset)) for ch in parent]
22 
23def que():
24    '(from the favourite saying of Manuel in Fawlty Towers)'
25    print ("#%-4i, fitness: %4.1f%%, '%s'" %
26           (iterations, fitness(parent)*100./perfectfitness, ''.join(parent)))
27 
28def mate(a, b):
29    place = 0
30    if choice(xrange(10)) < 7:
31        place = choice(xrange(len(target)))
32    else:
33        return a, b
34 
35    return a, b, a[:place] + b[place:], b[:place] + a[place:]
36 
37iterations = 0
38center = len(C)/2
39while parent != target:
40    rate = mutaterate()
41    iterations += 1
42    if iterations % 100 == 0: que()
43    copies = [ mutate(parent, rate) for _ in C ]  + [parent]
44    parent1 = max(copies[:center], key=fitness)
45    parent2 = max(copies[center:], key=fitness)
46    parent = max(mate(parent1, parent2), key=fitness)
47que()

1from random import random, randint, seed
2from statistics import mean
3from copy import deepcopy
4import matplotlib.pyplot as plt
5from IPython.display import Image, display
6from graphviz import Digraph, Source 
7
8POP_SIZE        = 60    # population size
9MIN_DEPTH       = 2     # minimal initial random tree depth
10MAX_DEPTH       = 5     # maximal initial random tree depth
11GENERATIONS     = 250   # maximal number of generations to run evolution
12TOURNAMENT_SIZE = 5     # size of tournament for tournament selection
13XO_RATE         = 0.8   # crossover rate 
14PROB_MUTATION   = 0.2   # per-node mutation probability 
15BLOAT_CONTROL   = False # True adds bloat control to fitness function
16
17def add(x, y): return x + y
18def sub(x, y): return x - y
19def mul(x, y): return x * y
20FUNCTIONS = [add, sub, mul]
21TERMINALS = ['x', -2, -1, 0, 1, 2] 
22
23def target_func(x): # evolution's target
24    return x*x*x*x + x*x*x + x*x + x + 1
25
26def generate_dataset(): # generate 101 data points from target_func
27    dataset = []
28    for x in range(-100,101,2): 
29        x /= 100
30        dataset.append([x, target_func(x)])
31    return dataset
32
33class GPTree:
34    def __init__(self, data = None, left = None, right = None):
35        self.data  = data
36        self.left  = left
37        self.right = right
38        
39    def node_label(self): # return string label
40        if (self.data in FUNCTIONS):
41            return self.data.__name__
42        else: 
43            return str(self.data)
44    
45    def draw(self, dot, count): # dot & count are lists in order to pass "by reference" 
46        node_name = str(count[0])
47        dot[0].node(node_name, self.node_label())
48        if self.left:
49            count[0] += 1
50            dot[0].edge(node_name, str(count[0]))
51            self.left.draw(dot, count)
52        if self.right:
53            count[0] += 1
54            dot[0].edge(node_name, str(count[0]))
55            self.right.draw(dot, count)
56        
57    def draw_tree(self, fname, footer):
58        dot = [Digraph()]
59        dot[0].attr(kw='graph', label = footer)
60        count = [0]
61        self.draw(dot, count)
62        Source(dot[0], filename = fname + ".gv", format="png").render()
63        display(Image(filename = fname + ".gv.png"))
64
65    def compute_tree(self, x): 
66        if (self.data in FUNCTIONS): 
67            return self.data(self.left.compute_tree(x), self.right.compute_tree(x))
68        elif self.data == 'x': return x
69        else: return self.data
70            
71    def random_tree(self, grow, max_depth, depth = 0): # create random tree using either grow or full method
72        if depth < MIN_DEPTH or (depth < max_depth and not grow): 
73            self.data = FUNCTIONS[randint(0, len(FUNCTIONS)-1)]
74        elif depth >= max_depth:   
75            self.data = TERMINALS[randint(0, len(TERMINALS)-1)]
76        else: # intermediate depth, grow
77            if random () > 0.5: 
78                self.data = TERMINALS[randint(0, len(TERMINALS)-1)]
79            else:
80                self.data = FUNCTIONS[randint(0, len(FUNCTIONS)-1)]
81        if self.data in FUNCTIONS:
82            self.left = GPTree()          
83            self.left.random_tree(grow, max_depth, depth = depth + 1)            
84            self.right = GPTree()
85            self.right.random_tree(grow, max_depth, depth = depth + 1)
86
87    def mutation(self):
88        if random() < PROB_MUTATION: # mutate at this node
89            self.random_tree(grow = True, max_depth = 2)
90        elif self.left: self.left.mutation()
91        elif self.right: self.right.mutation() 
92        
93    def size(self): # tree size in nodes
94        if self.data in TERMINALS: return 1
95        l = self.left.size()  if self.left  else 0
96        r = self.right.size() if self.right else 0
97        return 1 + l + r
98
99    def build_subtree(self): # count is list in order to pass "by reference"
100        t = GPTree()
101        t.data = self.data
102        if self.left:  t.left  = self.left.build_subtree()
103        if self.right: t.right = self.right.build_subtree()
104        return t
105                        
106    def scan_tree(self, count, second): # note: count is list, so it's passed "by reference"
107        count[0] -= 1            
108        if count[0] <= 1: 
109            if not second: # return subtree rooted here
110                return self.build_subtree()
111            else: # glue subtree here
112                self.data  = second.data
113                self.left  = second.left
114                self.right = second.right
115        else:  
116            ret = None              
117            if self.left  and count[0] > 1: ret = self.left.scan_tree(count, second)  
118            if self.right and count[0] > 1: ret = self.right.scan_tree(count, second)  
119            return ret
120
121    def crossover(self, other): # xo 2 trees at random nodes
122        if random() < XO_RATE:
123            second = other.scan_tree([randint(1, other.size())], None) # 2nd random subtree
124            self.scan_tree([randint(1, self.size())], second) # 2nd subtree "glued" inside 1st tree
125# end class GPTree
126                   
127def init_population(): # ramped half-and-half
128    pop = []
129    for md in range(3, MAX_DEPTH + 1):
130        for i in range(int(POP_SIZE/6)):
131            t = GPTree()
132            t.random_tree(grow = True, max_depth = md) # grow
133            pop.append(t) 
134        for i in range(int(POP_SIZE/6)):
135            t = GPTree()
136            t.random_tree(grow = False, max_depth = md) # full
137            pop.append(t) 
138    return pop
139
140def error(individual, dataset):
141    return mean([abs(individual.compute_tree(ds[0]) - ds[1]) for ds in dataset])
142
143def fitness(individual, dataset): 
144    if BLOAT_CONTROL:
145        return 1 / (1 + error(individual, dataset) + 0.01*individual.size())
146    else:
147        return 1 / (1 + error(individual, dataset))
148                
149def selection(population, fitnesses): # select one individual using tournament selection
150    tournament = [randint(0, len(population)-1) for i in range(TOURNAMENT_SIZE)] # select tournament contenders
151    tournament_fitnesses = [fitnesses[tournament[i]] for i in range(TOURNAMENT_SIZE)]
152    return deepcopy(population[tournament[tournament_fitnesses.index(max(tournament_fitnesses))]]) 
153            
154def prepare_plots():
155    fig, axarr = plt.subplots(2, sharex=True)
156    fig.canvas.set_window_title('EVOLUTIONARY PROGRESS')
157    fig.subplots_adjust(hspace = 0.5)
158    axarr[0].set_title('error', fontsize=14)
159    axarr[1].set_title('mean size', fontsize=14)
160    plt.xlabel('generation', fontsize=18)
161    plt.ion() # interactive mode for plot
162    axarr[0].set_xlim(0, GENERATIONS)
163    axarr[0].set_ylim(0, 1) # fitness range
164    xdata = []
165    ydata = [ [], [] ]
166    line = [None, None]
167    line[0], = axarr[0].plot(xdata, ydata[0], 'b-') # 'b-' = blue line    
168    line[1], = axarr[1].plot(xdata, ydata[1], 'r-') # 'r-' = red line
169    return axarr, line, xdata, ydata
170
171def plot(axarr, line, xdata, ydata, gen, pop, errors, max_mean_size):
172    xdata.append(gen)
173    ydata[0].append(min(errors))
174    line[0].set_xdata(xdata)
175    line[0].set_ydata(ydata[0])
176    sizes = [ind.size() for ind in pop]
177    if mean(sizes) > max_mean_size[0]:
178        max_mean_size[0] = mean(sizes)
179        axarr[1].set_ylim(0, max_mean_size[0])
180    ydata[1].append(mean(sizes))
181    line[1].set_xdata(xdata)
182    line[1].set_ydata(ydata[1])
183    plt.draw()  
184    plt.pause(0.01)
185
186def main():      
187    # init stuff
188    seed() # init internal state of random number generator
189    dataset = generate_dataset()
190    population= init_population() 
191    best_of_run = None
192    best_of_run_error = 1e20 
193    best_of_run_gen = 0
194    fitnesses = [fitness(ind, dataset) for ind in population]
195    max_mean_size = [0] # track maximal mean size for plotting
196    axarr, line, xdata, ydata = prepare_plots()
197
198    # go evolution!
199    for gen in range(GENERATIONS):        
200        nextgen_population=[]
201        for i in range(POP_SIZE):
202            parent1 = selection(population, fitnesses)
203            parent2 = selection(population, fitnesses)
204            parent1.crossover(parent2)
205            parent1.mutation()
206            nextgen_population.append(parent1)
207        population=nextgen_population
208        fitnesses = [fitness(ind, dataset) for ind in population]
209        errors = [error(ind, dataset) for ind in population]
210        if min(errors) < best_of_run_error:
211            best_of_run_error = min(errors)
212            best_of_run_gen = gen
213            best_of_run = deepcopy(population[errors.index(min(errors))])
214            print("________________________")
215            best_of_run.draw_tree("best_of_run",\
216                                  "gen: " + str(gen) + ", error: " + str(round(best_of_run_error,3)))
217        plot(axarr, line, xdata, ydata, gen, population, errors, max_mean_size)
218        if best_of_run_error <= 1e-5: break
219    
220    endrun = "_________________________________________________\nEND OF RUN (bloat control was "
221    endrun += "ON)" if BLOAT_CONTROL else "OFF)"
222    print(endrun)
223    s = "\n\nbest_of_run attained at gen " + str(best_of_run_gen) + " and has error=" + str(round(best_of_run_error,3))
224    best_of_run.draw_tree("best_of_run",s)
225    
226if __name__== "__main__":
227  main()