はじめに

Cartesian Genetic Programmingを実装してみたので，それについて．
今回実装するのは，回帰問題(単目的)用です．全体のソースコードはここにあげてます．
修正点・ミス・バグあったら教えてください．

Cartesian Genetic Programming(CGP)の実装

一般的にCGPは，多層のフィードフォワード型の合成関数で表現されます．

しかし，今回は実装が簡単なので，一層当たり1ノードのフィードフォワードネットワークにします．CGPは自分より前の層ならどこからでも接続できるので，ノード総数が同じならば，一層当たり1ノードのフィードフォワードネットワークでもそんなに問題ないはずです．

個体表現

個体は，固定長の配列で表現します．要素数がノードの総数です．このとき，すべてのノードが使われることはないので，表現型は可変長ネットワークになります．

また，各ノードは

functionノードの種類
子ノードの接続関係
定数(constノードのみ使用)

をパラメータとして持つ構造体です．さらに，Ephemeral Constant(ERC)と呼ばれるconstantノードを一定数持つことにします．ERCとは，個体の生成時と同時にランダムに決まる定数です．

class Gene:
	def __init__(self, func_name, inputs_arg, arity, const):
		self.func_name = func_name
		self.inputs_arg = inputs_arg
		self.arity = arity
		self.const = const

class Individual(list):
	def __init__(self, content=[]):
		list.__init__(self, content)
		self.fitness = None

	def create(self, gene_num, output_num, nodeSet, leafSet, ERC_func=lambda:random.random()*2-1):
		self.output_num = output_num
		variable_leaf = [(name,leaf_info) for name,leaf_info in leafSet.items() if leaf_info["const"]!=None]
		variable_leaf.sort(key=lambda obj:obj[1]["const"])
		const_leaf = [(name,leaf_info) for name,leaf_info in leafSet.items() if leaf_info["const"]==None]
		const_leaf.sort(key=lambda obj:obj[0])

		for n, item in enumerate(variable_leaf):
			name, leaf_info = item
			self.append(Gene(name, [], 0, leaf_info["const"]))
		for n, item  in enumerate(const_leaf):
			name, leaf_info = item
			self.append(Gene(name, [], 0, ERC_func()))

		function_node = [(name,node_info) for name,node_info in nodeSet.items() if name!="out"]
		max_arg_num = max([node_info["arity"] for node_info in nodeSet.values()])
		for n in range(gene_num):
			name,node_info = random.choice(function_node)
			inputs_arg = [random.randrange(len(self)) for x in range(max_arg_num)]
			self.append(Gene(name, inputs_arg, node_info["arity"], None))

		function_num = len(self)
		self.function_num = function_num
		for n in range(output_num):
			inputs_arg = [random.randrange(function_num) for x in range(max_arg_num)]
			self.append(Gene("out", inputs_arg, 1, None))

あとで，交叉しやすいように，子ノードの数は揃えてます．ただし，子ノードは先頭から使われていき，そのノードの関数の規定まで使われます．また，関数のノードセット(nodeSet)と入力ノードのノードセット(leafSet)は，各要素に関数名，入力数，定数が定義された辞書in辞書となります．

個体関数の実行

実行ですが，本実装の場合は，配列の先頭から実行しすれば大丈夫です．
しかし，実行しなくてもよいノードも存在するため，実行するノードだけを先に計算します．これは配列の最後の要素であるoutputノードから子ノードを辿ることで計算できます．

	def calc_path(self):
		path = [False for n in range(len(self))]
		for n in range(len(self))[::-1]:
			gene = self[n]
			if gene.func_name == "out":
				path[gene.inputs_arg[0]] = True
				path[n] = True
			else:
				if path[n]:
					for i in range(gene.arity):
						path[gene.inputs_arg[i]] = True
		return path

そして，実行するノードだけを，合成関数として計算します．

	def run(self, X, nodeSet, leafSet):
		path = self.calc_path()
		calc_rets = [None for _ in range(len(self))]
		for n, use in enumerate(path):
			if use:
				gene = self[n]
				name = gene.func_name
				if gene.arity == 0:
					if "x" in name:
						ret = leafSet[name]["function"](X, gene.const)
					else:
						ret = leafSet[name]["function"](gene.const, len(X))
				else:
					ret = nodeSet[name]["function"](calc_rets, gene.inputs_arg)
				calc_rets[n] = ret
		for n in range(1,self.output_num+1):
			ret = calc_rets[-1*n]
			mask = np.isfinite(ret)
			m = np.mean(ret[mask])
			ret[mask==False] = m
			calc_rets[-1*n] = ret

		return calc_rets[-1*self.output_num:]

一応，対応できるように書いてあります．単出力の場合は，返り値の最初の要素が出力列となります．

進化アルゴリズム

実装した進化アルゴリズムはGenetic Algorithm(GA)とEvolutionary Strategy(ES)の二つです．
また，この二つのアルゴリズムのパラメータは後述するsklearn_objectが持つパラメータを使います．

Genetic Algorithm

GAの場合は交叉・突然変異・選択の3つの操作を行います．トーナメント選択のみしか実装していないです．
交叉は，1点交叉と一様交叉の2つが選べます．突然変異は，使用ノードのみを変異させる方法と，すべてのノードを変異させる方法がとれます．
以後，増やしていきたいです．
また，各手法はエリート戦略として，エリートを一個体保存するようにしました．

def GA_evolution(sklearn_object):
	ind_num = sklearn_object.ind_num
	nodeSet = sklearn_object.nodeSet
	leafSet = sklearn_object.leafSet
	gene_num = sklearn_object.gene_num
	output_num = sklearn_object.output_num
	generation_num = sklearn_object.generation_num
	tour_size = sklearn_object.tour_size

	CXPB = sklearn_object.CXPB
	MUTPB = sklearn_object.MUTPB
	CXMode = sklearn_object.CXMode
	MUTMode = sklearn_object.MUTMode

	X = sklearn_object.train_X
	y = sklearn_object.train_y
	evaluator_args = (X,y,nodeSet,leafSet)

	pop = op.makeInitialPopulation(ind_num, nodeSet, leafSet, gene_num, output_num)
	elite = None

	for g in range(generation_num):
		pop = op.evaluation(pop, evaluator, evaluator_args)

		new_elite = min([ind for ind in pop if ind.fitness!=None], key=lambda ind:ind.fitness)
		if elite == None or elite.fitness == None or new_elite.fitness < elite.fitness:
			elite = deepcopy(new_elite)
		pop = op.tournament_selection(pop, tour_size, elite)

		#print(g,*getGenerationInfo(pop))
		if CXMode == "OnePoint":
			pop = op.one_point_crossover(pop, CXPB)
		if CXMode == "Uniform":
			pop = op.uniform_crossover(pop, CXPB)
		if MUTMode == "Point":
			pop = op.point_mutate(pop, MUTPB, nodeSet, leafSet)
		if MUTMode == "UsePoint":
			pop = op.use_point_mutate(pop, MUTPB, nodeSet, leafSet)


	pop = op.evaluation(pop, evaluator, evaluator_args)
	new_elite = min([ind for ind in pop if ind.fitness!=None], key=lambda ind:ind.fitness)
	if new_elite.fitness < elite.fitness:
		elite = deepcopy(new_elite)
	return elite

Evolutionary Strategy

ESの場合は，(1+λ)ESのみ実装しました．論文見てると(1+λ)ESが強いようなので．
ESの場合は，交叉ができないようにしました．突然変異はGAと同じく，2つの方法が選べます．

def ES_evolution(sklearn_object):
	ind_num = sklearn_object.ind_num
	nodeSet = sklearn_object.nodeSet
	leafSet = sklearn_object.leafSet
	gene_num = sklearn_object.gene_num
	output_num = sklearn_object.output_num
	generation_num = sklearn_object.generation_num

	MUTPB = sklearn_object.MUTPB
	MUTMode = sklearn_object.MUTMode
	stop_eval = sklearn_object.stop_eval

	X = sklearn_object.train_X
	y = sklearn_object.train_y
	evaluator_args = (X,y,nodeSet,leafSet)

	pop = op.makeInitialPopulation(ind_num, nodeSet, leafSet, gene_num, output_num)
	elite = None

	eval_num = 0
	for g in range(generation_num):
		pop = op.evaluation(pop, evaluator, evaluator_args)
		eval_num += len(pop)
		if eval_num > stop_eval:
			break
		new_elite = min([ind for ind in pop if ind.fitness!=None], key=lambda ind:ind.fitness)
		if elite == None or elite.fitness == None or  new_elite.fitness < elite.fitness:
			elite = deepcopy(new_elite)

		pop = [deepcopy(elite) for _ in range(ind_num)]
		if MUTMode == "Point":
			pop = op.point_mutate(pop, MUTPB, nodeSet, leafSet)
		if MUTMode == "UsePoint":
			pop = op.use_point_mutate(pop, MUTPB, nodeSet, leafSet)
		#print(g,*getGenerationInfo(pop))

	new_elite = min([ind for ind in pop if ind.fitness!=None], key=lambda ind:ind.fitness)
	if elite == None or elite.fitness == None or  new_elite.fitness < elite.fitness:
		elite = deepcopy(new_elite)
	return elite

scikit-learnからの呼び出し

sklearnのcross validationやgrid searchから呼び出せるように，sklearnのbaseの要件を満たします．
baseオブジェクトを継承をしない場合は，

fit関数
predict関数
get_params関数
set_params関数
score関数

を定義すればいいはずです(2019/07/12現在)．もしかしたらpredict関数はいらないかもしれません．
また，grid searchなどでscore関数を参照する場合には，戻り値のスコアを最大化します．今回は回帰問題を対象とするため，スコアに-1を掛けて返します．

class CGPRegressor:
	def __init__(self,ind_num=4, nodeSet=None,
		const_num=0, gene_num=100, generation_num=100, tour_size=7,
		CXPB=0.7, CXMode="OnePoint", MUTPB=0.05, MUTMode="UsePoint",
		ea_type="ES",stop_eval=100000):
		self.ind_num = ind_num
		if nodeSet == None:
			self.nodeSet = cgp_base.CreateNodeSet({"add","sub","mul","div"})
		else:
			self.nodeSet = nodeSet
		self.const_num = const_num
		self.gene_num = gene_num
		self.generation_num = generation_num
		self.tour_size = tour_size

		self.CXPB = CXPB
		if ea_type == "GA":
			self.MUTPB = 1.0 - CXPB
		else:
			self.MUTPB = MUTPB
		self.CXMode = CXMode
		self.MUTMode = MUTMode
		self.ea_type = ea_type
		self.output_num = 1
		self.stop_eval = stop_eval

	def fit(self, X, y):
		self.train_X = X
		self.train_y = y
		variable_num= X.shape[1]
		output_num = 1
		self.leafSet = cgp_base.CreateLeafSet(variable_num, self.const_num)
		if self.ea_type == "ES":
			self.best_ind = ea.ES_evolution(self)
		if self.ea_type == "GA":
			self.best_ind = ea.GA_evolution(self)
		return self

	def predict(self, X):
		_y = self.best_ind.run(X, self.nodeSet, self.leafSet)[0]
		return _y

	def get_params(self, deep=True):
		return {
		"nodeSet":self.nodeSet,
		"ind_num":self.ind_num,
		"const_num":self.const_num,
		"gene_num":self.gene_num,
		"generation_num":self.generation_num,
		"tour_size":self.tour_size,
		"CXPB":self.CXPB,
		"MUTPB":self.MUTPB,
		"CXMode":self.CXMode,
		"MUTMode":self.MUTMode,
		"ea_type":self.ea_type
		}

	def set_params(self, **parameters):
		for parameter, value in parameters.items():
			setattr(self,parameter, value)
		return self

	def score(self,X,y,sample_weight=None):
		_y = self.predict(X)
		s =  np.sum((y-_y)**2)/len(y)
		return -1*s

実行結果

実行してみます．対象の回帰問題ですが，データセットを用意するのが面倒なので，乱数列とします．
そのうちちゃんとした問題に適用してみます．

if __name__ == "__main__":# test code
	seed = 6
	inputs_num = 3
	const_num = 1
	instances_num = 5

	random.seed(seed)
	np.random.seed(seed)
	X = np.random.random((instances_num, inputs_num))
	y = np.random.random(instances_num)

	model = CGPRegressor(ind_num=4, nodeSet=None,
		const_num=const_num, gene_num=10, generation_num=100,
		MUTPB=0.05, MUTMode="UsePoint",
		ea_type="ES")
	model.fit(X,y)
	s = model.score(X,y)
	print(s)

実行結果

-0.08112019089484938

とりあえず，動いているようです．

おわりに

今回は，CGPを実装してみました．ソースコードはこちら
次は，CGPにはGAがよいのか，それともESがよいのかを検証しようと思っています．
そのときは，実際の回帰問題で試してみることにします．