PyGEP的使用说明
本帖最后由 Rainyboy 于 2010-12-12 21:14 编辑PyGEP说用说明
简介
PyGEP包中的文件如下：
里面最核心的程序就是population.py和chromosome.py，只要读懂这两个程序，整个GEP的思想也在其中了
还是先看一下该方法的框架吧!
从流程图可以看出，基本上的步骤就是
1、生成初始种群（由染色体组成chromosomes），每个染色体表达了问题的一个解。这本质上就是一个测试空间。我们所有的工作就是在该测试空间下通过一系列的操作以寻找出符合我们要求的解。
2、对每个染色体进行解码（Execute），得到我们的一个评估结果。每个染色体对应于一个适应性函数（fitness），和遗传算法（GA）类似。
3、对当前种群的改造，生成新的种群，重复2-3直到得到我们需要的解。
PyGEP结构介绍
一、生成种群
初始化种群时需要如下参数的设置，这些参数将决定如何对现在种群进行改造或者叫遗传。
- exclusion_level: selection pressure (typically 1.5)
- mutation_rate: mutation rate (2 per chromosome)
- inversion_rate: inversion probability (0.1)
- is_transposition_rate: non-root transposition (0.1)
- is_transposition_length: IS transposition lengths (1,2,3)
- ris_transposition_rate: root transposition rate (0.1)
- ris_transposition_length: RIS transposition lengths (1,2,3)
- gene_transposition_rate: gene transposition rate (0.1)
- crossover_one_point_rate: 1-point crossover rate (0.3)
- crossover_two_point_rate: 2-point crossover rate (0.3)
- crossover_gene_rate: full gene crossover (0.1)__init__()函数：def __init__(self, cls, size, head, genes=1, linker=default_linker)
cls表示种群对应的染色体类型，
size表示种群大小，
head种群头部长度，
genes 基因组数，
linker基因链接的方式
这些就是我们调用population的接口。整个改造过程是有cycle()函数完成的。
类中还定义了一个best属性用以标示当前的最好染色体。
best = property(
# Gives preference to later individuals tied for best
lambda self: max(reversed(self.population)),
doc='The best Chromosome of the current generation'
)
二、染色体
问题的描述主要是在染色体中进行的。所以实质上染色体的构造是我们对模型的构造。而种群是对算法的构造。
之所以说GEP算法在并行上有很大的优势是因为她天生就是并行的。因为对于每个染色体的解码是与其他染色体无关的。
分析chromosome类的外部作用，我们需要知道的是如何生成chromosome类，并且需要定义如何评估（fitness及solved）
chromosome类提供了一个generate函数：
def generate(cls, head, genes=1, linker=default_linker)
head是头部长度，genes是基因组数，linker是指链接方式。
同样还有：
def __init__(self, genes, head, linker=default_linker)同上。
关于fitness和solved，PyGEP提供了两个属性来连接。
fitness = property(lambda self: self._fitness(), doc='Fitness value')
solved= property(lambda self: self._solved(),doc='Problem solved')
而所有具体问题的描述都将在fitness和solved中进行。
当然在进行如上操作之前，我们必须要设置符号空间 functions和 terminals ，用以表示表达式。
三、一个例子
这个例子是包中的例子（majority），理解一下就清楚如何使用这个程序包了。
首先定义一个类：
from pygep import *
from pygep.functions.linkers import or_linker
from pygep.functions.logical import LOGIC_ALL
# Demo for the boolean majority function
class Majority(object):
SAMPLE = []
@staticmethod
def populate():
for a in 0, 1:
for b in 0, 1:
for c in 0, 1:
Majority.SAMPLE.append(Majority(a, b, c))
def __init__(self, a, b, c):
self.a = a
self.b = b
self.c = c
self.majority = (a and b) or (a and c) or (b and c)
该方程实质上是给了Majority这个bool运算函数的一个运算法则。我们想让gep通过符号运算来得到这法则。
在这个地方我所知道的其实就是不同的参数对应的Majority的值。
class BooleanFunction(Chromosome):
functions = LOGIC_ALL
terminals = 'a', 'b', 'c'
def _fitness(self):
# Fitness function: number of hits
hits = 0
for i in Majority.SAMPLE:
if i.majority == bool(self(i)):
hits += 1
return hits
def _solved(self):
return self.fitness >= len(Majority.SAMPLE)
这部分继承chromosome类，初始化了functions和terminals的算符表，并重写了_fitness()和_solved()函数。
if __name__ == '__main__':
Majority.populate() #生成数据
# Search for a solution
p = Population(BooleanFunction, 10, 3, 3, or_linker)#建立种群
print p
for _ in xrange(100):
if p.best.solved:#（检查是否满足要求）
break
p.cycle() #不满足就继续
print p
if p.best.solved:
print 'SOLVED:', p.best #输出最终计算结果
这部分所做的，调用了整个过程。由这个例子可以看出，如果我们要使用该包去计算一个问题时，我们所要做的是
1、建立一个变量列表所对应的数据表（如Majority中给出的是a，b，c，及对应的majority）
2、建立一个chromosome的类的继承，给定算符空间（functions and terminals）在继承中重写_fitness()和_solved()函数，
以满足我们具体问题的需要。
3、生成初始种群P0，并调用_solved()函数检查是否满足要求，调用cycle()函数进行种群改造。
这些就是我们所要做的事，如果你对比如选择算子，变异算子，等的概率设置不满意，同样可以到population类代码中进行修改。
同样对其他不满意也可以自己去修改，因为代码就在你手中。
又一个例子：
寻找一个多项式y=x**4+x**3+x**2+x
数据存在data.txt文件里。data.txt文件内容：
9.500366 9103.54918799079
-6.130432 1213.47815867463
3.252685 160.18146840485
7.88797 4432.23529247904
9.090484 7671.79392962917
1.485199 11.8327154232663
-3.950531 193.570365560718
10.003326 11124.3786278048
-.607453 -.326443121119579
5.469299 1093.78835980998
源代码：from pygep import *
from pygep.functions.mathematical.arithmetic import add_op, subtract_op, multiply_op, divide_op
from pygep.functions.linkers import sum_linker
# example for y=x**4+x**3+x**2+x
HH=add_op, subtract_op, multiply_op, divide_op
class mp(object):
SAMPLE = []
@staticmethod
def populate():
data=open("data.txt").read().split()
for i in range(0,5):
mp.SAMPLE.append(mp(float(data),float(data)))
def __init__(self, x,y):
self.x =x
self.value = y
class mpf(Chromosome):
functions = HH
terminals = tuple('x')
def _fitness(self):
# Fitness function: number of hits
hits = 0
for i in mp.SAMPLE:
if abs(i.value -float(self(i)))/abs(i.value)<0.01:
hits += 1
return hits
def _solved(self):
return self.fitness >= len(mp.SAMPLE)
if __name__ == '__main__':
mp.populate()
# Search for a solution
p = Population(mpf, 10, 10, 1,sum_linker)
for _ in xrange(1000):
if p.best.solved:
break
p.cycle()
if p.best.solved:
print p
print 'SOLVED:', p.bestresult:
+*x++x*x/*xxxxxxxxxxx : 4
+x*x++x***xxxxxxxxxxx : 5
+*x++x*x/*xxxxxxxxxxx : 4
+*+++x*x/*xxxxxxxxxxx : 0
+*/++x*x/*xxxxxxxxxxx : 0
+*x++x*x/*xxxxxxxxxxx : 4
+*x++x*x/*xxxxxxxxxxx : 4
+*x++x*-+*xxxxxxxxxxx : 0
+*x++x*x/*xxxxxxxxxxx : 4
+*x++x*x/*xxxxxxxxxxx : 4
SOLVED: +x*x++x***xxxxxxxxxxx
so we get the result:x**4+x**3+x**2+x.
这里需要重新修改一下乘法和除法运算，在pygep.functions.mathematical.arithmetic 这个文件里。做如下修改：
#divide_op = symbol('/')(lambda i, j: float(i) / j)
#modulus_op= symbol('%')(lambda i, j: i % j)
def divfun(i,j):
if abs(j)<10**(-12):
return float(i)*10**12
else:
return float(i)/j
def modfun(i,j):
if abs(j)<10**(-12):
return float(i)%10**(-12)
else:
return float(i)/j
divide_op = symbol('/')(divfun)
modulus_op= symbol('%')(modfun)
代码分析与应用
继续我们的讨论哈!
#--encoding=UTF-8--
#分析PyGEP
一. Chromosome类分析
1.MetaChromosome
这个metachromosome类只重写了type的new方法,并在new方法中为meta实例
添加了symbols字段,用以表示符号集
arity字段,用以识别符号集中的最大操作数目.比如'+'法函数又两个参数'a+b',则arity=2
_fitness字段用以得到适应性函数
2.Chromosome
__metaclass__ = MetaChromosome #得到类型名
__next_id = 1
gene_type = KarvaGene #表示gene类型是KarvaGene的
#下面重点分析一下代码
@classmethod
def generate(cls, head, genes=1, linker=default_linker):
tail = head * (cls.arity - 1) + 1 #得到尾部长度
while True:
new_genes = * genes #新建一个列表
for i in xrange(genes):
#调用KarvaGene,新建染色体序列
new_genes = cls.gene_type(
+ \
, head
)
yield cls(new_genes, head, linker) #返回该染色体,也就是Chromosome的一个实例
#注意这里使用yield生成的是一个代码序列.
#以下__init__代码:
def __init__(self, genes, head, linker=default_linker):
''
if head < 0:
raise ValueError('Head length must be at least 0')
if not genes:
raise ValueError('Must have at least 1 gene')
self.genes= genes
self.head = head
self.linker = linker
# Unique number of the organism
self.__id = type(self).__next_id #每一个实例对应一个id
type(self).__next_id += 1 #这里设置了一个迭代,这个迭代是在上面的while True:..中实现.
#__call__函数
def __call__(self, obj):
return self.linker(*)
#该函数只有一句,在return中执行,却是得到的我们染色体序列的值.具体需要分析karvaGene.这里不做解释.
#其他是一些操作.暂时不做解释.这些不影响我们的使用.
#现在有一个问题,如果我希望将一些染色体本身作为染色体序列中的一个gene来进行处理,那么该如何去做呢?
#这个问题的背景可以设想为,比如我知道当前的符号表达式是一些三角多项式,我要固定三角函数的形式,比如
#sin(x),那么我希望在生成的序列中要么是一些变量比如x,要么是sin(x)的形式.或者更复杂一些,我希望是
#sin(f(a)*x)的形式,其中f(a)是一些常数的组合,这些常数是需要我们在基因生成过程中来选择的.
#那么就需要考虑刚才的问题,这让我想到了<盗梦空间>中的多层梦境.
#和rainboy约定sin(x**2+x)*x来演示sin(f(x))*h(x),其中f,h需要在GEP中模拟.
#发现分别用种群来模拟,很难加入到种群遗传序列里.
#最终不得不用了一个最为简单的办法.重写__call__()函数.
#这里是两个例子.可见该方法还是有效的,目前看是最简单的,显然也是有很大限制的,这个一方面是由于我用的这个
#勉强使用的办法的缺陷,另一方面,由于需要在整体上进行进化,这也是不得已的办法.希望rainboy和各位朋友能想
#出更好的办法.
例子:sin(f(x))*h(x)
1.f(x)=x**2+x,h(x)=x
代码:(有很多与以上例子重复的地方,但是为了清晰起见还是把他放在这里.
第三个例子
HH=add_op, subtract_op, multiply_op, divide_op
class mp(object):
SAMPLE = []
@staticmethod
def populate():
for i in range(0,5):
x=(i+1)/math.e*math.pi
mp.SAMPLE.append(mp(x,math.sin(x**2+x)*(x**3+x**2)))
def __init__(self, x,y):
self.x =x
self.value = y
class mpf(Chromosome):
functions = HH
terminals = tuple('x')
def _fitness(self):
# Fitness function: number of hits
hits = 0
for i in mp.SAMPLE:
if abs(i.value -float(self(i)))/abs(i.value)<0.01:
hits += 1
return hits
def _solved(self):
return self.fitness >= len(mp.SAMPLE)
def __call__(self, obj):
'''
Evaluates a given GEP chromosome against some instance.The
terminals in the chromosome are assumed to be attributes on
the object instance provided (unless they are numeric constants).
@param obj: an object instance with terminal attributes set
@return: result of evaluating the chromosome
'''
return math.sin(float(self.genes(obj)))*float(self.genes(obj))
if __name__ == '__main__':
mp.populate()
# Search for a solution
p = Population(mpf, 10, 10, 2,sum_linker)
for _ in xrange(1000):
if p.best.solved:
break
p.cycle()
if p.best.solved:
print p
print 'SOLVED:', p.best
运行结果:
*++xxxxx+xxxxx : 1
*++xxxxx/xxxxx : 1
*++xxxxxx*xxxx : 1
x+xxxxxx//xxxx : 0
*++xxxxx/xxxxx : 1
*+*xxxxxxxxxxx : 0
-++xxxxxx-xxxx : 0
*++xxxxx++xxxx : 1
*++xxxxx+xxxxx : 1
+x*xxxxx+xxxxx : 5
SOLVED: +x*xxxxx+xxxxx
这个代码就是表示x**2+x 和x
另外也算了一个h(x)=x**3+x**2的例子.
只需要在以上代码里将mp类中的初始数据改一下就可以了.
运行结果:
012345678901234567890123456789012345678901
------------------------------------------
+-*--/+*//xxxxxxxxxxx*-*/+/*-+/xxxxxxxxxxx : 1
+-*--/+*//xxxxxxxxxxx*-*/+/*-+/xxxxxxxxxxx : 1
+-*+-/+*//xxxxxxxxxxx*-*/+-*-+/xxxxxxxxxxx : 0
+-**-/+*//xxxxxxxxxxx*-*/+//-+/xxxxxxxxxxx : 0
+-*--/+*//xxxxxxxxxxx*-*/+/+-+/xxxxxxxxxxx : 0
+-*--/--+*xxxxxxxxxxx*-*/+/*-+/xxxxxxxxxxx : 5
+-*-//+*/-xxxxxxxxxxx--*-*/+**-xxxxxxxxxxx : 0
+-*--/+*//xxxxxxxxxxx*--/+/*-+/xxxxxxxxxxx : 0
+-*--/+*//xxxxxxxxxxx*-*/+/*-+/xxxxxxxxxxx : 1
*-*/+/*-+/xxxxxxxxxxx+-*--/+*//xxxxxxxxxxx : 0
SOLVED: +-*--/--+*xxxxxxxxxxx*-*/+/*-+/xxxxxxxxxxx这个符号序列解码出来就是x**2+x和x**3+x**2.
希望有更好的办法,开始的时候一直希望能通过设置算符,来满足要求,结果却没有实现哈.d
本帖最后由 Rainyboy 于 2010-12-11 00:46 编辑
我来给smtmobly的第二个例子写个注释，呵呵。
========================代码分割线===============================
# -*- coding: cp936 -*-
from pygep import *
from pygep.functions.mathematical.arithmetic import add_op, subtract_op, multiply_op, divide_op
from pygep.functions.linkers import sum_linker
# example for y=x**4+x**3+x**2+x
HH=(add_op, subtract_op, multiply_op, divide_op)
#此处要准备欲拟合的差值基函数
class mp(object):
#这个类的作用主要是构造欲拟合的数据
#这相当于是一个public static变量，类型为list
SAMPLE = []
@staticmethod
def populate():
#从文件读入两列数据，并充满SAMPLE
data=open("data.txt").read().split()
for i in range(0,5):
mp.SAMPLE.append(mp(float(data),float(data)))
def __init__(self, x,y):
#这里要指明的变量名与后文密切相关
self.x =x
self.value = y
class mpf(Chromosome):
functions = HH
#由于有Self.x,所以这里是'x'，不信你改成'z'试试
terminals = ('x',)
def _fitness(self):
# Fitness function: number of hits
hits = 0
for i in mp.SAMPLE:
#比较个体（其实是一个函数）在拟合点的值与原始值的差异
#由于有self.value,所以这里有i.value
#至于self(i)手法，由于Chromosome定义了如下函数:
#|__call__(self, obj)
#| Evaluates a given GEP chromosome against some instance.The
#| terminals in the chromosome are assumed to be attributes on
#| the object instance provided (unless they are numeric constants).
#|
#| @param obj: an object instance with terminal attributes set
#| @return: result of evaluating the chromosome
#相当于这个类的对象同时具有函数的性质，以下两种写法等效
if abs(i.value -float(self.__call__(i)))/abs(i.value)<0.01:
#if abs(i.value -float(self(i)))/abs(i.value)<0.01:
hits += 1
return hits
def _solved(self):
#与所有拟合点的相差都满足条件才算满足
return self.fitness >= len(mp.SAMPLE)
if __name__ == '__main__':
mp.populate()
# Search for a solution
#下面的这行代码所设置的是：
#Chromosome的派生类
#代规模：population size
#个体编码的长度：chromosome head length
#一代中个体的个数：number of genes
#多个个体适应值累计方法：multigenic results linker function
p = Population(mpf, 10, 10, 1,sum_linker)
#后面的代码几乎是定式
for _ in xrange(1000):
if p.best.solved:
break
p.cycle()
if p.best.solved:
print p
print 'SOLVED:', p.best
========================================================
不知道个人这样理解对不对，还望各位及smtmobly指正！
本帖最后由 Rainyboy 于 2010-12-12 21:15 编辑
第三个例子也可以这样完成，采用定义新的linker的方式，不用重写__Call()__:
====================代码分割线================================
# -*- coding: cp936 -*-
from pygep import *
from pygep.functions.mathematical.arithmetic import add_op, subtract_op, multiply_op, divide_op
#from pygep.functions.linkers import sum_linker
import numpy as math
HH=add_op, subtract_op, multiply_op, divide_op
#自定义的linker
def udf_linker(*args):
return math.sin(float(args))*float(args);
class mp(object):
SAMPLE = []
@staticmethod
def populate():
for i in range(0,5):
x=(i+1)/math.e*math.pi
mp.SAMPLE.append(mp(x,math.sin(x**2+x)*(x**3+x**2)))
def __init__(self, x,y):
self.x =x
self.value = y
class mpf(Chromosome):
functions = HH
terminals = tuple('x')
def _fitness(self):
# Fitness function: number of hits
hits = 0
for i in mp.SAMPLE:
if abs(i.value -float(self(i)))/abs(i.value)<0.01:
hits += 1
return hits
def _solved(self):
return self.fitness >= len(mp.SAMPLE)
if __name__ == '__main__':
mp.populate()
# 创建两个种群
# 用自定义的linker表示两个种群之间的适应值关系
p = Population(mpf, 10, 10, 2,udf_linker)
#以下代码为定式
for _ in xrange(2000):
if p.best.solved:
break
p.cycle()
if p.best.solved:
print p
print 'SOLVED:', p.best
==============================================================
计算结果：
012345678901234567890123456789012345678901
------------------------------------------
--+-++++*xxxxxxxxxxxx*-x++x-**xxxxxxxxxxxx : 1
-x+-++++*xxxxxxxxxxxx*-x++x-+*xxxxxxxxxxxx : 5
+--+-++++*xxxxxxxxxxx*-x++x-**xxxxxxxxxxxx : 0
--+-++++*xxxxxxxxxxxx*-*++x-**xxxxxxxxxxxx : 0
--x-++++*xxxxxxxxxxxx*-x++x-**xxxxxxxxxxxx : 0
++*--+-+++xxxxxxxxxxx*-x++x-**xxxxxxxxxxxx : 0
--+-++++*xxxxxxxxxxxx*-x++x-**xxxxxxxxxxxx : 1
-++-++++*xxxxxxxxxxxx*xx++**-xxxxxxxxxxxxx : 0
--+/++++*xxxxxxxxxxxx*-x++x-**xxxxxxxxxxxx : 0
--+-+/++*xxxxxxxxxxxx**x**x++x*xxxxxxxxxxx : 0
SOLVED: -x+-++++*xxxxxxxxxxxx*-x++x-+*xxxxxxxxxxxx
请各位及smtmobly指正！
奇怪了，输入法没反应 睡觉前看了一下代码，其实你这样写和我的重写__call__()没有本质的区别，差别就在于你是直接修改了linker，而我是将修改的代码写到调用里去。
linker函数有一个基本要求就是需要有对称性，比如or,a or b =b or a ,+号也是。
之前我们在qq里讨论的是，使用f（G1）*G2的形式，G1作为一个单独的基因，G2作为单独的基因。但是那天弄了一下午还是不行，后来发现，其实这些基因调用全部都是call函数来指导的，因为call函数本质上excute了那个genes。先用这个办法先做做看看呵呵。想到好的再说。现在觉得这个GEP真的是不错，有意思 回复 5 # smtmobly 的帖子
恩，我也试了下，分别作为单独的基因我也没弄成……不过linker函数可以不对称啊，没有地方说他非得对称啊？ 回复 6 # Rainyboy 的帖子
是啊 ！没有必要。这个不是硬性规定。呵呵，这个东西想怎么玩就怎么玩 哈哈!就是哈!google第一个是pycode(google的)下面就是我们的帖子!
问下斑竹，在chromosome.py中_call_方法中的返回值return self.linker(*)是什么意思，谢谢 回复 9 # dahai66001 的帖子
就是将所有种群按照给定的linker组合起来，以便计算它们组合产生的SAMPLE在_fitness函数计算下的的适应值。 还有，我想做一个多输入单输出的系统辨识，我改了下代码，但是怎么弄都没有结果，不知道怎么回事，请斑竹指点下，源码：
from pygep import *
from pygep.functions.mathematical.arithmetic import add_op, subtract_op, multiply_op, divide_op
from pygep.functions.linkers import sum_linker
# example for y1=x1**2+x2**2+x1*x2+x2*x3
HH=add_op, subtract_op, multiply_op, divide_op
class mp(object):
SAMPLE = []
@staticmethod
def populate():
data=open("mimo.dat").read().split()
for i in range(0,5):
mp.SAMPLE.append(mp(float(data),\
float(data),\
float(data),\
float(data)))
def __init__(self, x1,x2,x3,y1,y2):
self.x1 =x1
self.x2 =x2
self.x3 =x3
self.value1 = y1
class mpf(Chromosome):
functions = HH
terminals = 'x1','x2','x3'
def _fitness(self):
# Fitness function: number of hits
hits = 0
for i in mp.SAMPLE:
if abs(i.value -float(self(i)))/abs(i.value)<0.001:
hits += 1
return hits
def _solved(self):
return self.fitness >= len(mp.SAMPLE)
if __name__ == '__main__':
mp.populate()
# Search for a solution
p = Population(mpf, 10, 15, 1,sum_linker)
for _ in xrange(1000):
if p.best.solved:
break
p.cycle()
if p.best.solved:
print p
print 'SOLVED:', p.best 那在调用的时候是这样调用的：if abs(i.value -float(self(i)))/abs(i.value)<0.001:
但是如果我的返回值是两个，该怎么做？就像多输入多输出的辨识，你在smtmobly的第二个例子写个注释，是个但输入单输出的例子，把他改成多输入多输出的时候得返回值该怎么处理？ 回复 12 # dahai66001 的帖子
将Value设置成元组应该就行，就是value = (y1 , y2) 俺去试一下 回复 11 # dahai66001 的帖子
我看了一下，其实你的问题适用第三个例子，也可参见我在3楼的做法，将你的“多输入”处理为“多种群”，然后通过定义linker表示它们的组合“输入”，在你的问题中，可以这样定义：
#自定义的linker
def udf_linker(*args):
(x1,x2,x3) = (float(args),float(args),float(args));
return x1**2+x2**2+x1*x2+x2*x3
其它地方的修改类似。
页:
[1]
2