公式樹開源庫分析

1.公式樹模塊

公式樹模塊的作用是，從訓練集X和function_set中進行隨機采樣，生成一棵公式樹，同時提供子樹變異、 crossover、hoist變異和點變異的方法。

1.1基本屬性

1.2方法

1.3 build_tree算法原理

用到的數據結構:

terminal_stack: 存儲是幾元運算的一個棧

symbol_tree: lisp_tree 列表樹, Lisp列表是基于廣義表的結構，所以很容易將一個列表表達成樹結構。 S-表達式可能以其在Lisp家族的編程語言中的使用而為人所知，約翰·麥卡錫發明LISP于1958年，首次由史蒂夫·拉塞爾實施在IBM704計算機上，它特別適合用于人工智能方案，因為它有效地處理的符號信息。

在前綴表示法，運算符在自己操作數前寫。例如，表達式

a * ( b + c ) / d

被寫成

(/ (* a (+ b c) ) d)

例如公式：

它也可以寫作：

寫成S表達式就變成了這個

對應的二叉樹

也就是說s表達式對應于符號樹的先序遍歷

算法輸入:function_set[‘add’, ‘sub’, ‘mul’] , arities{2:[‘add’, ‘sub’, ‘mul’]}, method: grow , max_depth:3(2-4內的一個隨機數)

method: grow n_features: 10 max_depth:2 function_index:0 function地址: function_name:add program:[add] terminal_stack:[2] 循環部分 ############LOOP########################## 第1次 depth: 1 # program的長度也就是符號樹的list表長度 等于 len(terminal_stack) choice: 13 choice_index: 1 #depth < max_depth或者choice_index <= len(self.function_set)就會從function_set里面選擇一個比如這里選擇add function_index: 0 function: add 1_program: [add, add] 1_terminal_stack: [2, 2] 第2次 depth: 2 choice: 13 choice_index: 11 2_terminal: 10 # 這里terminal是10和n_features相等，所以會生成一個隨機的float數值 3_terminal: 0.8395650516882855 2_program: [add, add, 0.8395650516882855] 2_terminal_stack: [2, 1]# 加入了常數 第二個值就減去1 第3次 depth: 2 choice: 13 choice_index: 0 2_terminal: 8 # 這里terminal是10和n_features如果不相等的話就加入一個特征列 3_terminal: 8 2_program: [add, ddd, 0.8395650516882855, 8] 2_terminal_stack: [2, 0] 2_terminal_stack.pop(): 0 # 等于0的時候就被pop掉了 然后stack減去1 第4次 depth: 1 choice: 13 choice_index: 5 2_terminal: 0 3_terminal: 0 2_program: [add, add, 0.8395650516882855, 8, 0] 2_terminal_stack: [0] 2_terminal_stack.pop(): 0 最后樹的形狀變成這個 tree1: add(add(0.840, X8), X0) 生成第2顆樹的過程 method: grow max_depth:4 function_index:0 function的地址: function_name:add program:[add] terminal_stack:[2] #############LOOP######################### 第1次 depth: 1 choice: 13 choice_index: 4 2_terminal: 3 3_terminal: 3 2_program: [add, 3] 2_terminal_stack: [1] 第2次 depth: 1 choice: 13 choice_index: 4 2_terminal: 6 3_terminal: 6 2_program: [add, 3, 6] 2_terminal_stack: [0] 2_terminal_stack.pop(): 0 最后樹的形狀變成這個 tree2: [add, 3, 6]

畫成流程圖就是

1.4 get_subtree 獲取隨機子樹

給symbol_tree里面的元素賦予權重如果是算子就是0.9 不是算子就是0.1

比如tree是

tree1: mul(X6, sub(X3, X1))

那么久賦予權重

probs: [0.9 0.1 0.9 0.1 0.1]

然后進行歸一化就是除以sum

_probs: [0.42857143 0.47619048 0.9047619 0.95238095 1. ]

這里是采用輪盤賭法來選擇切割點的

步驟

1）隨機生成一個（0,1）內的隨機值比如生成

s_index:0.8299421213898753

2）找到隨機值在probs中應該在的位置這里這個位置就是start的位置

start: 2

3）初始化 end = start=2

stack = 1

如果end - start < stack那么

node = program[end]

如果node是算子的話那么stack要加上元數

stack += node.arity

end 自身一直加一直到program最后

畫成流程圖就是

1.5 crossover算法模塊原理

crossover的目的是進行子樹交叉

第一步從符號樹模塊獲得隨機子樹

start, end = self.get_subtree(random_state)

第二步從其他符號樹個體中獲得隨機子樹

donor_start, donor_end = self.get_subtree(random_state, donor)

第三步獲得交叉后的符號樹

self.symbol_tree[: start] + donor[donor_start : donor_end] + self.symbol_tree[end : ]

tree1: add(mul(X9, X4), X6) start, end: 1, 4 removed:range(1, 4) donor: [mul, 6, 0.6656811846792283] donor_start, donor_end:(0, 3) donor_removed:[] 結合生成的子樹 new_tree[add, mul, 6, 0.6656811846792283, 6] 最后得到的結果 add mul x6 0.6656811846792283 x6 self.symbol_tree[: start] + donor[donor_start : donor_end] + self.symbol_tree[end : ]

1.6 subtree_mutation子樹變異

由p_subtree_mutation參數控制。這是一種更激進的變異策略：優勝者的一棵子樹將被另一棵完全隨機的全新子樹代替。

chicken = self.build_tree(random_state) 直接生產一棵新的子樹 self.crossover(chicken, random_state) # 然后用crossover算法直接進行交叉生成子樹

1.7 hoist_mutation hoist變異

hoist變異是一種對抗公式樹膨脹（bloating，即過于復雜）的方法：從優勝者公式樹內隨機選擇一個子樹A，再從A里隨機選擇一個子樹B，然后把B提升到A原來的位置，用B替代A。hoist的含義即「升高、提起」。

第一步獲取一個隨機子樹A

start, end = self.get_subtree(random_state) subtree = self.symbol_tree[start:end]

第二步從子樹A中獲取一個子樹B

# 獲取隨機子樹B sub_start, sub_end = self.get_subtree(random_state, subtree) hoist = subtree[sub_start:sub_end]

第三步 把B提升到A原來的位置，用B替代A

self.symbol_tree[:start] + hoist + self.symbol_tree[end:]

tree1: add(X6, sub(add(X0, X2), X6)) start, end: 0, 7 subtree : [add, x6, sub, add, x0, x2, x6] mutation: sub_start, sub_end: 3, 6 mutation_hoist : [add, x0, x2] removed: [0, 1, 2, 6] new_tree: [add, x0, x2] 第二次 tree1: mul(X8, X0) mutation_start, end: 0, 3 mutation_subtree : [mul, x8, x0] mutation: sub_start, sub_end: 1, 2 mutation_hoist : [x8] removed: [0, 2] new_tree: [8]

1.8 point_mutation 點變異

由p_point_replace參數控制。一個隨機的節點將會被改變，比如加法可以被替換成除法，變量X0可以被替換成常數-2.5。點變異可以重新加入一些先前被淘汰的函數和變量，從而促進公式的多樣性。

第一步復制符號樹，并獲取一個隨機的點

program = copy(self.symbol_tree) # 自己復制一份 # 隨機生成符號樹長度個點，然后找到其中小于點變異概率的點組成一個list mutate = np.where(random_state.uniform(size = len(program)) < self.p_point_replace)[0] # 獲取一個隨機的點

第二步遍歷mutate的node節點，如果node是一個Function就替換，不是的話就加入常量或者feature

if isinstance(program[node], _Function): arity = program[node].arity # 算子元數 replacement_len = len(self.arities[arity]) # 找到和arity元數一樣的算子有多少個 replacement_index = random_state.randint(replacement_len) # 從里面隨機選擇一個 replacement = self.arities[arity][replacement_index] # 找到index對應的算子 program[node] = replacement # 進行替換

如果不是function

第一種情況

# 不是算子的話就是常量或者端點 加入一個常量 if self.const_range is not None: terminal = random_state.randint(self.n_features + 1) else: terminal = random_state.randint(self.n_features) # 隨機生成一個（0，n_features）內的一個數字terminal if terminal == self.n_features: # 如果terminal和n_features相等的話就替換為一個（0,1）內float的數字 terminal = random_state.uniform(*self.const_range) if self.const_range is None: raise ValueError('A constant was produced with const_range=None.') program[node] = terminal

2. fitness 模塊獲得符號樹的適應性

2.1 get_all_indices 接口獲得所有數據的index

第一步：進行參數校驗

if self._indices_state is None and random_state is None: # 如果_indices_state和random_state都是None raise ValueError('The program has not been evaluated for fitness\n yet, indices not available.') if n_samples is not None and self._n_samples is None: #如果n_samples不為None self._n_samples = n_samples if max_samples is not None and self._max_samples is None: # n_samples代表數據集的行數 self._max_samples = max_samples # max_samples最大采樣樹 if random_state is not None and self._indices_state is None: self._indices_state = random_state.get_state()

第二步 獲得隨機種子，然后獲得袋外數據和袋內數據的index

indices_state = check_random_state(None) indices_state.set_state(self._indices_state) # 得到random_state not_indices = sample_without_replacement( self._n_samples, self._n_samples - self._max_samples, random_state=indices_state) # 袋外數據 這里是從[0,self._n_samples]中選出self._n_samples - self._max_samples個數據 sample_counts = np.bincount(not_indices, minlength=self._n_samples) # 找到每個index出現的次數了 indices = np.where(sample_counts == 0)[0] # 出現次數是零的index就是被留下的數據，在袋內的數據了

其他函數

sample_without_replacement(n_population, n_samples, random_state,method): 采樣函數，隨機獲取袋外數據，從集合[0，n_population]中選擇n_samples個數據，有放回的抽樣

參數介紹

2.2 raw_fitness

接口

raw_fitness(self, X, y, sample_weight)

先執行X的算法得到y_pred，然后根據y,y_pred以及權重計算誤差

# 根據x，y 評估符號樹的適用性 返回fitness y_pred = self.execute(X) raw_fitness = self.metric(y, y_pred, sample_weight)

2.3 fitness模塊

接口

fitness(self, parsimony_coefficient=None)

先執行X的算法得到y_pred，然后根據y,y_pred以及權重計算誤差

if parsimony_coefficient is None: parsimony_coefficient = self.parsimony_coefficient penalty = parsimony_coefficient * len(self.symbol_tree) * self.metric.sign # 這里是節約系數乘以公式樹的長度如果越大越好sign是1否則是-1 return self.raw_fitness_ - penalty # fitness進行了約減，這里懲罰了樹過于膨脹的公式

3.并行模塊

3.1并行parallel_evolve

接口

_parallel_evolve(n_programs, parents, X, y, sample_weight, seeds, params)

入參

屬性：

3.2內置接口_tournament

目的：找到表現最好的符號樹

contenders = random_state.randint(0, len(parents), tournament_size) # 生成tournament_size個（0，len(parents)）的數字相當于從父類中隨機選多少個 fitness = [parents[p].fitness_ for p in contenders] # 得到這些被選中的符公式樹的評分 if metric.greater_is_better: # 判斷指標是不是越大越好還是越小越好 parent_index = contenders[np.argmax(fitness)] # np.argmax找到最大的那個值的index else: parent_index = contenders[np.argmin(fitness)] # 越小越好 return parents[parent_index], parent_index # 返回最大的對象和他的index

3.3運行流程：

第一步：n_programs表示種群里的一個組內有多少顆樹，這里要循環n_programs次

初始化

method = random_state.uniform() # method從crossover subtree hoist point變異里選中的概率 parent, parent_index = _tournament() # 找到了表現好的公式樹

第二步根據method的概率選擇突變的類型

method_probs定義

self._method_probs = np.array([self.p_crossover, self.p_subtree_mutation, self.p_hoist_mutation, self.p_point_mutation]) self._method_probs = np.cumsum(self._method_probs)

if method < method_probs[0]: # 如果method小于crossover概率的話 # 走crossover方法 donor, donor_index = _tournament() # 次好的公式樹作為子樹 program, removed, remains = parent.crossover(donor.symbol_tree, random_state) # 兩者交叉 genome = {'method':'Crossover', 'parent_idx': parent_index, 'parent_nodes':removed, 'donor_idx':donor_index, 'donor_nodes':remains } elif method < method_probs[1]:# 如果method小于crossover概率的話 # subtree突變 program, removed, _ = parent.subtree_mutation(random_state) genome = {'method':'Subtree Mutation', 'parent_idx':parent_index, 'parent_nodes':removed } elif method < method_probs[2]: # hoist突變 program, removed = parent.hoist_mutation(random_state) genome = {'method':'Hoist Mutation', 'parent_idx': parent_index, 'parent_node': removed } elif method < method_probs[3]: # point突變 program, mutated = parent.point_mutation(random_state) genome = {'mehtod':'Point Mutation', 'parent_idx':parent_index, 'parent_nodes':mutated } else: # 自身拷貝 program = parent.reproduce() genome = {'mehtod': 'Reproduction', 'parent_idx':parent_index, 'parent_nodes':[] }

第三步 根據參數和第二步得到的program生成公式樹

program = _SymbolTree(function_set=function_set, arities=arities, init_depth=init_depth, init_method=init_method, n_features=n_features, metric=metric, const_range=const_range, p_point_replace=p_point_replace, parsimony_coefficient=parsimony_coefficient, feature_names=feature_names, random_state=random_state, symbol_tree = program)

然后

program.parents = genome

這里的genome存儲的是之前生成子樹的過程中刪掉的信息把他賦值給parents

第四步 根據sample_weight中的權重信息對特征列賦予權重。

if sample_weight is None: # 計算袋內數據權重 curr_sample_weight = np.ones((n_samples,)) # 如果沒有權重信息那么所有樣本權重都是1 else: curr_sample_weight = sample_weight.copy() oob_sample_weight = curr_sample_weight.copy() # 袋外數據

計算袋外數據的fitness

indices, not_indices = program.get_all_indices(n_samples, max_samples, random_state) # 得到選擇的在袋外的數據index curr_sample_weight[not_indices] = 0 # 原數據里面屬于袋外數據的其index對應的權重置為零 oob_sample_weight[indices] = 0 # 袋外數據里面不在原數據里的其index對應的權重置為零 program.raw_fitness_ = program.raw_fitness(X, y, curr_sample_weight) # 計算袋內數據的fitness

計算袋外數據的fitness

if max_samples < n_samples: # 計算袋外數據的適用性 program.oob_fitness_ = program.raw_fitness(X, y, oob_sample_weight) # 計算袋內數據的fitness

最后循環n次就得到了n顆變異后的子樹programs，它里面有兩個私有屬性raw_fitness_，oob_fitness_分別存儲了袋內袋外數據的適用性

4.SymbolicTransformer模塊

4.1初始化模塊

_verbose_reporter：控制日志輸出

接口

fit(self, X, y, sample_weight = None)

入參：

校驗：檢查X和y的長度是否一致、hall_of_fame、function_set、_arities是不是正常以及metric是不是Fitness類型 自身是否繼承于TransformerMixin抽象類

然后把概率放到list里面，逐步加

self._method_probs = np.array([self.p_crossover, self.p_subtree_mutation, self.p_hoist_mutation, self.p_point_mutation]) self._method_probs = np.cumsum(self._method_probs)

然后校驗_method_probs、init_method、const_range、init_depth、feature_names進行類型檢查和范圍檢查

params = self.get_params() print(f'params: {params}') params['_metric'] = self._metric params['function_set'] = self._function_set params['arities'] = self._arities params['method_probs'] = self._method_probs

如果不是warm_start模式就準備好_programs和 run_details_字典

if not self.warm_start or not hasattr(self, '_programs'): # warm_start為false時就重新開始，釋放分配的內存 self._programs = [] # _programs里面存的是每一代的優勝者 list[list1,list2... listn] self.run_details_ = { 'generation':[], 'average_length':[], 'average_fitness':[], 'best_length':[], 'best_fitness': [], 'best_oob_fitness':[], 'generation_time':[] }

prior_generations = len(self._programs) # _programs里面存的是每一代的優勝者 list[list1,list2... listn] 其長度代表迭代了多少代 n_more_generations = self.generations - prior_generations # 還有多少代需要迭代

然后是對n_more_generations進行校驗

population_size代表每一個世代中種群的數目

if self.warm_start: # 丟棄之前用過的隨機種子 for i in range(len(self._programs)): _ = random_state.randint(MAX_INT, size = self.population_size)

if self.verbose: self._verbose_reporter()# 輸出日志

循環層（從prior_generations到generations）

1）記錄時間找到父類

start_time = time() if gen == 0: # 如果是第一代的話parents是None parents = None else: parents = self._programs[gen - 1] # _programs里面的最后一代

2）并行運行程序

n_jobs, n_programs, starts = _partition_estimators(self.population_size, self.n_jobs) # 把種群分成n_job份n_programs表示第幾份中有多少個數據 starts記錄了每組數據的起點 seeds = random_state.randint(MAX_INT, size = self.population_size) # 產生population_size個隨機種子 population = Parallel(n_jobs=n_jobs, verbose=int(self.verbose > 1))( delayed(_parallel_evolve)(n_programs[i],parents, X, y,sample_weight, seeds[starts[i]:starts[i + 1]], params) for i in range(n_jobs))

對數據進行合并，得到下一代變異的種群

population[, , , ]

population = list(itertools.chain.from_iterable(population))

得到種群的所有的個體的fitness和length是一個list

fitness = [program.raw_fitness_ for program in population] length = [program.length_ for program in population]

3）懲罰系數對fitness進行約束

parsimony_coefficient = None if self.parsimony_coefficient == 'auto': parsimony_coefficient = (np.cov(length, fitness)[1, 0] / np.var(length)) # 取出（length, fitness）協方差矩陣的第2行1列除以方差 for program in population: program.fitness_ = program.fitness(parsimony_coefficient) # 收縮后的適應度 self._programs.append(population) #新生成的這一代的信息放入_programs

4）刪除被淘汰的個體

if not self.low_memory: for old_gen in np.arange(gen, 0, -1): # 把到gen的世代數倒序排成list類似[5 4 3 2 1] indices = [] for program in self._programs[old_gen]: # 找到上一代的種群每一個符號樹 if program is not None:# 不是None的話 for idx in program.parents: # 找到他的parents_idx parents_idx里面存的是其表現最好的父類 if 'idx' in idx:# 找到其中的parent_idx indices.append(program.parents[idx]) indices = set(indices) # 去重復 for idx in range(self.population_size):# 種群內每一個個體 if idx not in indices: # 如果該個體不在最優集合里面就把他置為None self._programs[old_gen - 1][idx] = None elif gen > 0: self._programs[gen - 1] = None #不然就把上一代置為None

對應代碼

# 記錄運行信息 if self._metric.greater_is_better: # 如果是越大越好的話 best_program = population[np.argmax(fitness)] else: best_program = population[np.argmin(fitness)] self.run_details_['generation'].append(gen) self.run_details_['average_length'].append(np.mean(length)) self.run_details_['average_fitness'].append(np.mean(fitness)) self.run_details_['best_length'].append(best_program.length_) self.run_details_['best_fitness'].append(best_program.raw_fitness_) oob_fitness = np.nan if self.max_samples < 1.0: oob_fitness = best_program.oob_fitness_ self.run_details_['best_oob_fitness'].append(oob_fitness) # 袋外數據準確性 generation_time = time() - start_time self.run_details_['generation_time'].append(generation_time) # 運行時間

處理early stopping

if self._metric.greater_is_better: best_finess = fitness[np.argmax(fitness)] if best_finess >= self.stopping_criteria: # 達到停止標準的時候 break else: best_finess = fitness[np.argmix(fitness)] if best_finess <= self.stopping_criteria: # 達到停止標準的時候 break

到這里循環結束，得到所有的世代。

a）首先得到hall_of_fame個索引

fitness = np.array(fitness) # 找到這一代種群的fitness if self._metric.greater_is_better: # 越大越好的好就倒序選擇 hall_of_fame = fitness.argsort()[::-1][:self.hall_of_fame] #得到hall_of_fame個fitness的索引 else: hall_of_fame = fitness.argsort()[:self.hall_of_fame] # 越小越好就正序選擇

對最后一代的種群里所有的個體（其中屬于hall_of_fame的）進行計算得到預測的值

evaluation = np.array([gp.execute(X) for gp in [self._programs[-1][i] for i in hall_of_fame]])

如果指標是spearman系數的話，計算得到evaluation每一組數據的排序值

if self.metric == 'spearman': evaluation = np.apply_along_axis(rankdata, 1, evaluation)

from scipy.stats import rankdata evaluation = np.array([[1,2,3,4], [6,5,7,8], [9,10,11,12]]) print(np.apply_along_axis(rankdata, 1, evaluation)) #輸出 [[1. 2. 3. 4.] [2. 1. 3. 4.] [1. 2. 3. 4.]]

然后計算相關系數矩陣

with np.errstate(divide = 'ignore', invalid = 'ignore'): # 去掉除0 無效值等警告 correlations = np.abs(np.corrcoef(evaluation)) # 得到相關系數矩陣 如果是spearman系數這里就是spearman相關系數

[[1. 2. 3. 4.] [2. 1. 3. 4.] [1. 2. 3. 4.]] [[1. 0.8 1. ] [0.8 1. 0.8] [1. 0.8 1. ]]

np.fill_diagonal(correlations, 0.) # 對角線元素填充0

components = list(range(self.hall_of_fame)) # hall_of_frame個0到hall_of_frame的數字 indices = list(range(self.hall_of_fame)) # X_shape(354, 13) # evaluation:(50, 354)

# 迭代刪除最小相關的的元素 while len(components) > self.n_components: # 如果最小相關個體大于要留的個數 # correlations_shape(50, 50) most_correlated = np.unravel_index(np.argmax(correlations), correlations.shape) # 得到最大值的索引 相關性越大越不好 # 通過適應度對相關矩陣進行排序，確定最不合適的索引 worst = max(most_correlated) # worst就是索引大的那一列 components.pop(worst) indices.remove(worst) # 從序列號中刪除 correlations = correlations[:, indices][indices, :] indices = list(range(len(components))) self._best_programs = [self._programs[-1][i] for i in hall_of_fame[components]]

AI EI企業智能 機器學習 運籌優化

版權聲明：本文內容由網絡用戶投稿，版權歸原作者所有，本站不擁有其著作權，亦不承擔相應法律責任。如果您發現本站中有涉嫌抄襲或描述失實的內容，請聯系我們jiasou666@gmail.com 處理，核實后本網站將在24小時內刪除侵權內容。