一文詳細分析公式樹開源庫

語言: CN / TW / HK
摘要:公式樹模塊的作用是,從訓練集X和function_set中進行隨機採樣,生成一棵公式樹,同時提供子樹變異、 crossover、hoist變異和點變異的方法。

本文分享自華為雲社區《公式樹開源庫分析》,作者:鯉魚君 。

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地址:<functions._Function object at 0x000001E538356EB0>
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的地址:<functions._Function object at 0x000001E538356EB0>
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初始化模塊

4.1.1基本屬性

4.1.2方法

_verbose_reporter:控制日誌輸出

4.1.3 fit模塊

接口

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

入參:

第1步:對數據進行校驗

校驗:檢查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進行類型檢查和範圍檢查

第2步:取出參數並對其進行賦值

 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()# 輸出日誌

第3步:並行的運行程序

循環層(從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[<symboltree._SymbolTree object at 0x00000118ABB89E20>, <symboltree._SymbolTree object at 0x00000118ABB89E80>, <symboltree._SymbolTree object at 0x00000118ABB89F10>, <symboltree._SymbolTree object at 0x00000118ABB89FD0>]
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

第4步進行運行的信息

對應代碼

# 記錄運行信息
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

到這裏循環結束,得到所有的世代。

第5步如果是變換的話

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]]

 

點擊關注,第一時間瞭解華為雲新鮮技術~