diff --git a/README.md b/README.md index 9acebc4c..6e2cfc7b 100644 --- a/README.md +++ b/README.md @@ -59,11 +59,11 @@ which gives: ### API What follows is the API reference for running the numpy interface. -Note that most parameters here -have been tuned with ~1000 trials over several example -equations, so you likely don't need to tune them yourself. +You likely don't need to tune the hyperparameters yourself, +but if you would like, you can use `hyperopt.py` as an example. However, you should adjust `threads`, `niterations`, -`binary_operators`, `unary_operators`, and `maxsize` to your requirements. +`binary_operators`, `unary_operators`, and `maxsize` +to your requirements. The program will output a pandas DataFrame containing the equations, mean square error, and complexity. It will also dump to a csv @@ -148,7 +148,7 @@ pd.DataFrame, Results dataframe, giving complexity, MSE, and equations # TODO - [ ] Make scaling of changes to constant a hyperparameter -- [ ] Update hall of fame every iteration +- [ ] Update hall of fame every iteration? - [ ] Calculate feature importances of future mutations, by looking at correlation between residual of model, and the features. - Store feature importances of future, and periodically update it. - [ ] Implement more parts of the original Eureqa algorithms: https://www.creativemachineslab.com/eureqa.html @@ -162,6 +162,8 @@ pd.DataFrame, Results dataframe, giving complexity, MSE, and equations - Current most expensive operations: - [ ] Calculating the loss function - there is duplicate calculations happening. - [x] Declaration of the weights array every iteration +- [x] Add a node at the top of a tree +- [x] Insert a node at the top of a subtree - [x] Record very best individual in each population, and return at end. - [x] Write our own tree copy operation; deepcopy() is the slowest operation by far. - [x] Hyperparameter tune diff --git a/eureqa.jl b/eureqa.jl index 88ecc26d..9c06fddd 100644 --- a/eureqa.jl +++ b/eureqa.jl @@ -3,6 +3,22 @@ import Optim const maxdegree = 2 const actualMaxsize = maxsize + maxdegree + +# Sum of square error between two arrays +function SSE(x::Array{Float32}, y::Array{Float32})::Float32 + diff = (x - y) + return sum(diff .* diff) +end + +# Mean of square error between two arrays +function MSE(x::Array{Float32}, y::Array{Float32})::Float32 + return SSE(x, y)/size(x)[1] +end + +const len = size(X)[1] +const avgy = sum(y)/len +const baselineSSE = SSE(y, convert(Array{Float32, 1}, ones(len) .* avgy)) + id = (x,) -> x const nuna = size(unaops)[1] const nbin = size(binops)[1] @@ -211,7 +227,6 @@ end # Evaluate an equation over an array of datapoints function evalTreeArray(tree::Node)::Array{Float32, 1} - len = size(X)[1] if tree.degree == 0 if tree.constant return ones(Float32, len) .* tree.val @@ -225,21 +240,10 @@ function evalTreeArray(tree::Node)::Array{Float32, 1} end end -# Sum of square error between two arrays -function SSE(x::Array{Float32}, y::Array{Float32})::Float32 - diff = (x - y) - return sum(diff .* diff) -end - -# Mean of square error between two arrays -function MSE(x::Array{Float32}, y::Array{Float32})::Float32 - return SSE(x, y)/size(x)[1] -end - # Score an equation function scoreFunc(tree::Node)::Float32 try - return MSE(evalTreeArray(tree), y) + countNodes(tree)*parsimony + return SSE(evalTreeArray(tree), y)/baselineSSE + countNodes(tree)*parsimony catch error if isa(error, DomainError) return 1f9 @@ -290,17 +294,55 @@ function appendRandomOp(tree::Node)::Node return tree end -# Select a random node, and replace it an the subtree -# with a variable or constant -function deleteRandomOp(tree::Node)::Node +# Add random node to the top of a tree +function popRandomOp(tree::Node)::Node + node = tree + choice = rand() + makeNewBinOp = choice < nbin/nops + left = tree + + if makeNewBinOp + right = randomConstantNode() + newnode = Node( + binops[rand(1:length(binops))], + left, + right + ) + else + newnode = Node( + unaops[rand(1:length(unaops))], + left + ) + end + node.l = newnode.l + node.r = newnode.r + node.op = newnode.op + node.degree = newnode.degree + node.val = newnode.val + node.constant = newnode.constant + return node +end + +# Insert random node +function insertRandomOp(tree::Node)::Node node = randomNode(tree) - # Can "delete" variable or constant too - if rand() > 0.5 - val = Float32(randn()) + choice = rand() + makeNewBinOp = choice < nbin/nops + left = copyNode(node) + + if makeNewBinOp + right = randomConstantNode() + newnode = Node( + binops[rand(1:length(binops))], + left, + right + ) else - val = rand(1:nvar) + newnode = Node( + unaops[rand(1:length(unaops))], + left + ) end - newnode = Node(val) node.l = newnode.l node.r = newnode.r node.op = newnode.op @@ -310,6 +352,120 @@ function deleteRandomOp(tree::Node)::Node return tree end +function randomConstantNode()::Node + if rand() > 0.5 + val = Float32(randn()) + else + val = rand(1:nvar) + end + newnode = Node(val) + return newnode +end + +# Return a random node from the tree with parent +function randomNodeAndParent(tree::Node, parent::Union{Node, Nothing})::Tuple{Node, Union{Node, Nothing}} + if tree.degree == 0 + return tree, parent + end + a = countNodes(tree) + b = 0 + c = 0 + if tree.degree >= 1 + b = countNodes(tree.l) + end + if tree.degree == 2 + c = countNodes(tree.r) + end + + i = rand(1:1+b+c) + if i <= b + return randomNodeAndParent(tree.l, tree) + elseif i == b + 1 + return tree, parent + end + + return randomNodeAndParent(tree.r, tree) +end + +# Select a random node, and replace it an the subtree +# with a variable or constant +function deleteRandomOp(tree::Node)::Node + node, parent = randomNodeAndParent(tree, nothing) + isroot = (parent == nothing) + + if node.degree == 0 + # Replace with new constant + newnode = randomConstantNode() + node.l = newnode.l + node.r = newnode.r + node.op = newnode.op + node.degree = newnode.degree + node.val = newnode.val + node.constant = newnode.constant + elseif node.degree == 1 + # Join one of the children with the parent + if isroot + return node.l + elseif parent.l == node + parent.l = node.l + else + parent.r = node.l + end + else + # Join one of the children with the parent + if rand() < 0.5 + if isroot + return node.l + elseif parent.l == node + parent.l = node.l + else + parent.r = node.l + end + else + if isroot + return node.r + elseif parent.l == node + parent.l = node.r + else + parent.r = node.r + end + end + end + return tree +end + +# Simplify tree +function combineOperators(tree::Node)::Node + # (const (+*) const) already accounted for + # ((const + var) + const) => (const + var) + # ((const * var) * const) => (const * var) + # (anything commutative!) + if tree.degree == 2 && (tree.op == plus || tree.op == mult) + op = tree.op + if tree.l.constant || tree.r.constant + # Put the constant in r + if tree.l.constant + tmp = tree.r + tree.r = tree.l + tree.l = tmp + end + topconstant = tree.r.val + # Simplify down first + tree.l = combineOperators(tree.l) + below = tree.l + if below.degree == 2 && below.op == op + if below.l.constant + tree = below + tree.l.val = op(tree.l.val, topconstant) + elseif below.r.constant + tree = below + tree.r.val = op(tree.r.val, topconstant) + end + end + end + end + return tree +end # Simplify tree function simplifyTree(tree::Node)::Node @@ -355,12 +511,15 @@ function iterate( tree = mutateOperator(tree) elseif mutationChoice < cweights[3] && n < maxsize tree = appendRandomOp(tree) - elseif mutationChoice < cweights[4] - tree = deleteRandomOp(tree) + elseif mutationChoice < cweights[4] && n < maxsize + tree = insertRandomOp(tree) elseif mutationChoice < cweights[5] + tree = deleteRandomOp(tree) + elseif mutationChoice < cweights[6] tree = simplifyTree(tree) # Sometimes we simplify tree + tree = combineOperators(tree) # See if repeated constants at outer levels return tree - elseif mutationChoice < cweights[6] + elseif mutationChoice < cweights[7] tree = genRandomTree(5) # Sometimes we simplify tree else return tree @@ -608,14 +767,15 @@ function fullRun(niterations::Integer; debug(verbosity, "-----------------------------------------") debug(verbosity, "Complexity \t MSE \t Equation") println(io,"Complexity|MSE|Equation") - for size=1:maxsize + for size=1:actualMaxsize if hallOfFame.exists[size] member = hallOfFame.members[size] - numberSmallerAndBetter = sum([member.score > hallOfFame.members[i].score for i=1:(size-1)]) + curMSE = MSE(evalTreeArray(member.tree), y) + numberSmallerAndBetter = sum([curMSE > MSE(evalTreeArray(hallOfFame.members[i].tree), y) for i=1:(size-1)]) betterThanAllSmaller = (numberSmallerAndBetter == 0) if betterThanAllSmaller - debug(verbosity, "$size \t $(member.score-parsimony*size) \t $(stringTree(member.tree))") - println(io, "$size|$(member.score-parsimony*size)|$(stringTree(member.tree))") + debug(verbosity, "$size \t $(curMSE) \t $(stringTree(member.tree))") + println(io, "$size|$(curMSE)|$(stringTree(member.tree))") push!(dominating, member) end end diff --git a/eureqa.py b/eureqa.py index 85f1d789..0f933605 100644 --- a/eureqa.py +++ b/eureqa.py @@ -6,27 +6,26 @@ import pandas as pd # Dumped from hyperparam optimization -default_alpha = 5.429178 -default_annealing = 1.000000 -default_fractionReplaced = 0.415076 -default_fractionReplacedHof = 0.031713 -default_ncyclesperiteration = 1614.000000 -default_niterations = 25.000000 -default_npop = 300.000000 -default_parsimony = 0.001861 -default_topn = 9.000000 -default_weightAddNode = 3.788920 -default_weightDeleteNode = 2.553399 -default_weightDoNothing = 0.498528 -default_weightMutateConstant = 2.081372 -default_weightMutateOperator = 2.003413 -default_weightRandomize = 4.679883 -default_weightSimplify = 0.009620 -default_result = -1.183938 +default_alpha = 5 +default_fractionReplaced = 0.30 +default_fractionReplacedHof = 0.05 +default_npop = 200 +default_weightAddNode = 1 +default_weightInsertNode = 1 +default_weightDeleteNode = 1 +default_weightMutateConstant = 10 +default_weightMutateOperator = 1 +default_weightRandomize = 1 +default_weightSimplify = 0.1 +default_weightDoNothing = 1 +default_result = 1 +default_topn = 10 +default_parsimony = 0.0 + def eureqa(X=None, y=None, threads=4, niterations=20, - ncyclesperiteration=int(default_ncyclesperiteration), + ncyclesperiteration=10000, binary_operators=["plus", "mult"], unary_operators=["cos", "exp", "sin"], alpha=default_alpha, @@ -40,6 +39,7 @@ def eureqa(X=None, y=None, threads=4, shouldOptimizeConstants=True, topn=int(default_topn), weightAddNode=default_weightAddNode, + weightInsertNode=default_weightInsertNode, weightDeleteNode=default_weightDeleteNode, weightDoNothing=default_weightDoNothing, weightMutateConstant=default_weightMutateConstant, @@ -84,6 +84,7 @@ def eureqa(X=None, y=None, threads=4, constants (Nelder-Mead/Newton) at the end of each iteration. :param topn: int, How many top individuals migrate from each population. :param weightAddNode: float, Relative likelihood for mutation to add a node + :param weightInsertNode: float, Relative likelihood for mutation to insert a node :param weightDeleteNode: float, Relative likelihood for mutation to delete a node :param weightDoNothing: float, Relative likelihood for mutation to leave the individual :param weightMutateConstant: float, Relative likelihood for mutation to change @@ -141,6 +142,7 @@ def eureqa(X=None, y=None, threads=4, {weightMutateConstant:f}, {weightMutateOperator:f}, {weightAddNode:f}, + {weightInsertNode:f}, {weightDeleteNode:f}, {weightSimplify:f}, {weightRandomize:f}, @@ -199,10 +201,18 @@ def eureqa(X=None, y=None, threads=4, parser.add_argument("--maxsize", type=int, default=20, help="Max size of equation") parser.add_argument("--niterations", type=int, default=20, help="Number of total migration periods") parser.add_argument("--npop", type=int, default=int(default_npop), help="Number of members per population") - parser.add_argument("--ncyclesperiteration", type=int, default=int(default_ncyclesperiteration), help="Number of evolutionary cycles per migration") + parser.add_argument("--ncyclesperiteration", type=int, default=10000, help="Number of evolutionary cycles per migration") parser.add_argument("--topn", type=int, default=int(default_topn), help="How many best species to distribute from each population") parser.add_argument("--fractionReplacedHof", type=float, default=default_fractionReplacedHof, help="Fraction of population to replace with hall of fame") parser.add_argument("--fractionReplaced", type=float, default=default_fractionReplaced, help="Fraction of population to replace with best from other populations") + parser.add_argument("--weightAddNode", type=float, default=default_weightAddNode) + parser.add_argument("--weightInsertNode", type=float, default=default_weightInsertNode) + parser.add_argument("--weightDeleteNode", type=float, default=default_weightDeleteNode) + parser.add_argument("--weightMutateConstant", type=float, default=default_weightMutateConstant) + parser.add_argument("--weightMutateOperator", type=float, default=default_weightMutateOperator) + parser.add_argument("--weightRandomize", type=float, default=default_weightRandomize) + parser.add_argument("--weightSimplify", type=float, default=default_weightSimplify) + parser.add_argument("--weightDoNothing", type=float, default=default_weightDoNothing) parser.add_argument("--migration", type=bool, default=True, help="Whether to migrate") parser.add_argument("--hofMigration", type=bool, default=True, help="Whether to have hall of fame migration") parser.add_argument("--shouldOptimizeConstants", type=bool, default=True, help="Whether to use classical optimization on constants before every migration (doesn't impact performance that much)") diff --git a/hyperparamopt.py b/hyperparamopt.py index 3fb4741f..e2716d49 100644 --- a/hyperparamopt.py +++ b/hyperparamopt.py @@ -37,9 +37,14 @@ def run_trial(args): for key in 'niterations npop'.split(' '): args[key] = int(args[key]) - total_steps = 10*100*5000 + + total_steps = 20*100*5000 niterations = args['niterations'] npop = args['npop'] + if niterations == 0 or npop == 0: + print("Bad parameters") + return {'status': 'ok', 'loss': np.inf} + args['ncyclesperiteration'] = int(total_steps / (niterations * npop)) args['topn'] = 10 args['parsimony'] = 1e-3 @@ -52,7 +57,7 @@ def run_trial(args): args['weightDoNothing'] = 1.0 - maxTime = 2*60 + maxTime = 3*60 ntrials = 2 equation_file = f'.hall_of_fame_{np.random.rand():f}.csv' @@ -112,6 +117,7 @@ def run_trial(args): 'weightMutateConstant': hp.lognormal('weightMutateConstant', np.log(4.0), 1.0), 'weightMutateOperator': hp.lognormal('weightMutateOperator', np.log(0.5), 1.0), 'weightAddNode': hp.lognormal('weightAddNode', np.log(0.5), 1.0), + 'weightInsertNode': hp.lognormal('weightInsertNode', np.log(0.5), 1.0), 'weightDeleteNode': hp.lognormal('weightDeleteNode', np.log(0.5), 1.0), 'weightSimplify': hp.lognormal('weightSimplify', np.log(0.05), 1.0), 'weightRandomize': hp.lognormal('weightRandomize', np.log(0.25), 1.0), diff --git a/operators.jl b/operators.jl index 727170b4..db635b49 100644 --- a/operators.jl +++ b/operators.jl @@ -1,6 +1,6 @@ # Define allowed operators. Any julia operator can also be used. -plus(x::Float32, y::Float32)::Float32 = x+y -mult(x::Float32, y::Float32)::Float32 = x*y +plus(x::Float32, y::Float32)::Float32 = x+y #Do not change the name of this operator. +mult(x::Float32, y::Float32)::Float32 = x*y #Do not change the name of this operator. pow(x::Float32, y::Float32)::Float32 = sign(x)*abs(x)^y div(x::Float32, y::Float32)::Float32 = x/y loga(x::Float32)::Float32 = log(abs(x) + 1)