遺伝的プログラミング(GP)を用いて、特定の点を通る関数近似をしたい

実現したいこと

少し数学よりの質問になってしまうのですが、特定の点を通る関数を遺伝的プログラミング(以下、GP)で生成したいです。
事前にあるデータ(xi, yi)が与えられているときに、データの中の複数個を選択し、それらの点を必ず通るという制約を設けたうえでGPによる回帰を行いたいです。
GPの流れを簡単に説明すると、
①初期解(遺伝子)生成。（今回は関数近似が目的なので、木構造で数式を表したもの）
②適合度関数で遺伝子を評価(適合度関数は絶対誤差の総和)
③遺伝子を選択。トーナメント選択を採用。ランダムに指定した個数分遺伝子を取り出し、その中から最も最適な遺伝子を次の世代に残す。
④遺伝子交叉。遺伝子同士を組み合わせて新たな解を作成し、次世代に残す。
⑤遺伝子突然変異。
⑥②～⑤の操作を繰り返し、最終的に準最適解(数式)を出力
というものです。GPについてはこちらのサイトがわかりやすく説明されています。

該当のソースコード

Python
1# 遺伝子(木の生成)
2def randomTreeMake(nodeSet,leafSet,_min,_max):
3
4    height = random.randint(_min,_max)　# 生成する木の深さ
5    items = []
6    stack = [0]
7    while len(stack) != 0:
8        depth = stack.pop()
9        if depth == height:
10            node = random.choice(leafSet)    # leafSetは木の終端記号(定数、変数xなど．内容については後述)
11        else:
12            node = random.choice(nodeSet) # nodeSetは木の非終端記号(四則演算記号や三角関数など)
13            stack.extend([depth+1]*node[2])
14        items.append(node) 
15    tree = Tree(items)
16    return tree
17
18# 遺伝子(数式)にデータxを代入し、計算
19def get_y(tree,nodeSet,leafSet,x):
20    stack = []
21    for node in tree[::-1]:
22        name,func,arity = node
23        if name == "x":
24            output = x
25        else:
26            if arity == 0:
27                output = func*np.ones(x.shape[0])
28            if arity == 1:
29                outout = func(stack.pop())
30            if arity == 2:
31                output = func(stack.pop(),stack.pop())
32        stack.append(output)
33    return stack.pop()
34
35# 遺伝子の評価
36def evaluate(population,nodeSet,leafSet,x,y):
37
38    for tree in population:
39        if tree.fitness == None:
40            _y = get_y(tree,nodeSet,leafSet,x)   # _y：GPで生成したランダムな数式にデータのx座標を代入して求めた値
41            fitness = np.sum(np.abs(_y - y))　# 適合度関数．データxにおけるGPで生成した(_y)とデータyの絶対誤差の総和
42            tree.fitness = fitness
43
44def main():
45    """main function
46
47    Args:
48        N: 個体数
49        M: トーナメント選択数
50        maxG: 最大世代数
51        Pcx: 交叉率
52        Pmut: 突然変異率
53        nodeSet: 非終端記号
54        leafSet: 終端記号
55        pupulation: 個体群
56    """
57    N = 500
58    M = 8
59    maxG = 200
60    Pcx = 0.5
61    Pmut = 0.1
62    nodeSet,leafSet = makeNodeSet()
63    x, y = generate_dataset(input_data.dataset)  # データ(x, y)を取得
64
65    population = initial_population(N,nodeSet,leafSet, MIN_DEPTH, MAX_DEPTH)
66    g = 0
67    while True:
68        evaluate(population,nodeSet,leafSet,x,y)
69        printLog(population,g)
70        if g == maxG:
71            break
72        population = select(population,M)
73        population = crossover(population,Pcx)
74        population = mutation(population,Pmut,nodeSet,leafSet)
75        g+=1
76    best = min(population,key=lambda tree:tree.fitness)

また、leafSet, nodeSetは木構造の数式を生成するときに使用する終端記号、非終端記号を代入したものです。

Python
1def makeNodeSet():
2    nodeSet = []
3    leafSet = []
4    nodeSet.append(("add",np.add,2))   # 3つ目の引数の数字は子ノードの数を表す。
5    nodeSet.append(("sub",np.subtract,2))
6    nodeSet.append(("mul",np.multiply,2))
7    nodeSet.append(("div",protectedDiv,2))
8    # nodeSet.append(("pow",np.power,2))
9    nodeSet.append(("sin",np.sin,1))
10    nodeSet.append(("cos",np.cos,1))
11    nodeSet.append(("tan",np.tan,1))
12    nodeSet.append(("log",np.log,1))
13    nodeSet.append(("exp",np.exp,1))
14    nodeSet.append(("sqrt",np.sqrt,1))
15
16    leafSet.append(("-1",-1,0))
17    leafSet.append(("0",0,0))
18    leafSet.append(("1",1,0))
19    leafSet.append(("0.5",0.5,0))
20    leafSet.append(("0.1",0.1,0))
21    leafSet.append(("pi",np.pi,0))
22    leafSet.append(("e",np.e,0))
23    leafSet.append(("x",None,0))

現状と修正したい点

＜現状＞
ランダムにleafSet, nodeSetを用いて木構造(数式)を生成。生成した数式にデータxiにおける_yを求める。
求めた_yと元のデータyiの絶対誤差の総和を取り、この値が最小となるような関数を探索。

<修正したい点>
使用するデータ(xi, yi)のうち、複数個のデータを指定し、それらは必ず通るという制約を設けたうえで関数近似をしたいがこの部分の実装をどのようにすればよいのかわからない。
いじるとすれば、randomTreeMake、get_y，evaluate関数あたりかなと考えています。
回答するうえでわからない点があればすぐに補足いたしますのでおっしゃってください。