💻

遺伝的アルゴリズムで条件を満たすいい感じの部分集合をとってくる

に公開
遺伝的アルゴリズム
Bioinformatics
NSGA-II
生命科学
化合物群
tech

多目的最適化というとNSGAIIが有名.
遺伝的アルゴリズムのパッケージDEAP, pymoo辺りはもちろんのこと, 最近だとOptunaにも実装がある.
これらの例題については情報が豊富で, google先生に聞けば多くの解説サイトがある.
一方, 生命科学への具体的な事例へどう適用するか？についてはまだ少なく感じたので, その点について簡単にまとめる.

遺伝的アルゴリズムはどのような時に使うか？

そもそも遺伝的アルゴリズムとは？的な話は多くの良い記事があるので割愛.
使い時をわかっているとは言い難いが, 個人的には 「厳密解は存在しないものの8割くらい条件を満たすものがあれば事足りる」 時に使っている.
特に多目的最適化の場合のアウトカムはパレート最適解であり, そもそも厳密解は存在しない.
パレート最適解とはざっくり言えば「他により優越する解がない解」の集合のこと, 詳細はgoogle先生に譲る(こちらのサイトさんなど).

遺伝的アルゴリズムの生物・数学の対応

集団(population), 個体(individual), 遺伝子(gene)の三つの階層を考える.
他に世代(generation), 操作(operation)が重要なターミノロジー.
訳語などがGA界隈で厳密なのかはわからないものの, 生物情報学者としてはgenomeとgeneが一緒くたになっているのは見過ごせないのでこのように表記 (これはこれで気になるところも).

生物的な面と数学的な面と対応させるとこんな感じか:

  • 個体
    • (生物) 集団の構成要素, genomeで定義される
    • (数学) 多変量ないし部分集合, 最適化すべき解の実体
  • 遺伝子
    • (生物) genomeの構成要素
    • (数学) 多変量の各変数, 部分集合の要素
  • 集団
    • (生物) 個体の集合, 集団の中で個体が交叉したり変異したりして集団に属する個体が進化する
    • (数学) 解候補の集合, 解候補同士で変数・要素を交換するなど操作してより良い解を探索する
  • 世代
    • (生物) ある時期の集団を表すラベル
    • (数学) イテレーション
  • 操作
    • (生物) 変異, 交叉, など
    • (数学) 各変数・要素に対する乱数置換など確率的な操作, 多変量・部分集合同士の変数・要素の入れ替え, など

ざっくりした流れ

  1. 初期集団を準備する
  2. 各個体の良し悪しをスコアリングする
  3. 個体を操作して新しい個体を作る
  4. 新しい個体とこれまでの集団の個体から新しい集団を作る
  5. 特定の基準まで2-4を繰り返す

NSGAII

NSGA-Ⅱ(Elitist Non-dominated Sorting Genetic Algorithm)は多目的遺伝的アルゴリズム(multi-objective GA)の一種で, 2002年に開発されてから今でも第一線で使われている.
それまでの多目的GAに対して以下の工夫を加え, 速度面などを改善している.

  1. 高速優越ソート
  2. シェアリングの代わる混雑度指標の導入
  3. エリート主義の導入
    色々なところで紹介されているので詳細は割愛, こちらなどがGA系の歴史もわかってわかりやすい.

NSGAIIの生命科学っぽい課題への適用：部分集合の選定

特定の化合物群や遺伝子群といった「全体集合から複数の条件をいい感じに満たす部分集合抽出したい」ときに使いやすい印象.
この辺りで述べられているよう特徴量選択に関する紹介はちょこちょこある.

定式化

集合Vの部分集合S=\{s_1,\cdots,s_{k} \}を入力とする互いに競合するp個の目的関数f_{i} \colon S \to \mathbb{R} (i = 1, \cdots ,p)を, 最小化(または最大化)するS^*を求める問題として定式化できる.
肝は具体的にどうやって表現するか. これ次第で変わるんだと思う.
ナップサック問題的な問題を例によくある定式化と, 集合を使った定式化を紹介する.

よくある？定式化

普通ナップサック問題というと, 「詰め込める最大の重量だか容積が決まったナップサックに対し, 異なる価値・重さのモノを何個ずつ入れるとナップサック全体の価値を最大化できるか」といった問題.
この問題を, 「有限な集合からいくらかの条件を満たす指定した要素数の部分集合を選定する問題」に見立てるために, 選べる数が1つのみというケースを考える.

全体の集合(ex. ナップサックの中身候補, データベース中の全化合物)は有限なので, 各要素にインデックスを振り, 固定長の二値ベクトルで選択を表現する.
普通のナップサック問題ではここが二値ベクトルではなく, 個数を表現するベクトルとなっている(こちらのサイトさんとか).
全体の集合をV_n = \{v_1, \cdots, v_n \}, その部分集合(ex. 条件を満たすナップサックの中身達, 化合物群)をSとする.

Individual:=\boldsymbol{x}=[x_1, \cdots, x_n]

, where x_i = \begin{cases} 1 & (v_i \in S)\\ 0 & (v_i \notin S) \end{cases} and |S| = k

以下の問題となる.

{\underset{x}{\argmin}} f_{i} (\boldsymbol{x}) (i = 1, \cdots ,p)

ちなみにOptunaだとsuggestの候補にベクトルが存在しないので現状(220929)ではできないっぽい (suggest_intなどでindex指定する形はできるため重複ありなら組めるが, 抽出先を弄ると動的はダメと怒られる).
遺伝的アルゴリズムと言えばDEAPがやはり良い (もちろんpymooも良さそう, 単に私が使っていないだけ).

集合を使った定式化

さてよくあると書いたものの, 上記のタイプだと実は固定長が扱いづらい.
なぜなら交叉や変異といった操作に際し, 要素数を保存するようにgenomeを作っていく必要があるため.
この点固定長でも使いやすいのが集合を使った定式化, 単に集合を見つけてくるだけになる.

Individual:=S

, where S \subset V and |S| = k

DEAPの公式さんが普通のナップサック問題を例に, 集合を使ったバージョンでの実装わかりやすく書いてくれている.
交叉と変異の関数であるcxSetとmutSetの部分で固定長になるように工夫すれば, 目的とする固定長の部分集合をとってくる最適化アルゴリズムができる.

ただDEAP自体癖があるので説明しながら書くのに時間がかかりそう…なのでとりあえずCLIコードだけぶん投げて説明は追記する.
DEAPの癖についてはこちらがわかりやすい.

import argparse
import random
import numpy as np
import pandas as pd

from deap import algorithms
from deap import base
from deap import creator
from deap import tools

# input preparation, necessary to be modified
def prep_input(n_total):
    """ prepare items for defining knapsack problem """
    items = dict()
    for i in range(n_total):
        items[i] = (random.randint(1, 10), random.uniform(0, 100))
    return items

# objective function, necessary to be modified
def eval(individual, items, max_weight):
    weight = 0.0
    value = 0.0
    for item in individual:
        weight += items[item][0]
        value += items[item][1]
    if weight > max_weight: # penalty for overweight
        return 10000, 0
    return weight, value

# functions for handling fixed length set in GA
def prepInd(INSTANCE, n_total, n_set):
    """
    prepare individual inheriting a creator instance

    """
    whole = list(range(n_total))
    np.random.shuffle(whole)
    return INSTANCE(whole[:n_set])

def cxSet(ind1, ind2):
    """
    crossover keeping length of individual
    note ind1 and ind2 should be creator instances

    """
    symdif = list(ind1 ^ ind2)
    if len(symdif)==0:
        pass
    else:
        np.random.shuffle(symdif)
        point = int(len(symdif)/2)
        l, r = symdif[:point], symdif[point:]
        tmp = ind1.copy()
        ind1 &= ind2
        ind1 |= set(l)
        ind2 &= tmp
        ind2 |= set(r)
    return ind1, ind2

def mutSet(individual, n_total):
    """
    mutation keeping length of individual
    note individual should be a creator instance
    
    """
    whole = set(range(n_total)) - individual
    individual.remove(random.choice(list(individual)))
    whole = list(whole)
    np.random.shuffle(whole)
    individual.add(whole[0])
    return individual,

# main runner
def main():
    # prepare arguments
    parser = argparse.ArgumentParser()
    # IO
    parser.add_argument("-i", "--input", type=int, required=True)
    parser.add_argument("-o", "--output", type=str, required=True)
    # conditions
    parser.add_argument("-s", "--nset", type=int, required=True)
    parser.add_argument("-t", "--ntotal", type=int, required=True)
    parser.add_argument("-w", "--maxweight", type=int, required=True)
    # GA
    parser.add_argument("--ngen", type=int, default=50)
    parser.add_argument("--mu", type=int, default=50)
    parser.add_argument("--lmd", type=int, default=100)
    parser.add_argument("--cxprob", type=float, default=0.5)
    parser.add_argument("--mutprob", type=float, default=0.2)
    # others
    parser.add_argument("-f", "--seed", type=int, default=1)

    args = parser.parse_args()
    n_total = args.ntotal
    random.seed(args.seed)
    np.random.seed(args.seed)

    # preparation
    items = prep_input(args.input)
    creator.create("Fitness", base.Fitness, weights=(-1.0, 1.0))
    ## weights depend on the number of objective functions and its direction
    creator.create("Individual", set, fitness=creator.Fitness)
    toolbox = base.Toolbox()

    # Structure initializers
    toolbox.register(
        "individual", prepInd, creator.Individual, n_total=n_total, n_set=args.nset
        )
    toolbox.register("population", tools.initRepeat, list, toolbox.individual)

    # Prepare algorithm
    toolbox.register("evaluate", eval, items=items, max_weight=args.maxweight)
    toolbox.register("mate", cxSet)
    toolbox.register("mutate", mutSet, n_total=n_total)
    toolbox.register("select", tools.selNSGA2)

    # run
    pop = toolbox.population(n=args.mu)
    hof = tools.ParetoFront()
    stats = tools.Statistics(lambda ind: ind.fitness.values)
    stats.register("avg", np.mean, axis=0)
    stats.register("std", np.std, axis=0)
    stats.register("min", np.min, axis=0)
    stats.register("max", np.max, axis=0)

    algorithms.eaMuPlusLambda(
        pop, toolbox, args.mu, args.lmd, args.cxprob, args.mutprob, args.ngen, stats,halloffame=hof
        )

    # export
    res = []
    for ind in hof:
        res.append([set(ind), ind.fitness.values[0], ind.fitness.values[1]])
        ## individual.fitness.values depends on the number of objectives
    df = pd.DataFrame(res, columns=['set', 'score_obj0', 'score_obj1'])
    df.to_csv(args.output + '/result.txt', sep='\t')

if __name__ == "__main__":
    main()

Discussion