ルーレット選択

３分プログラミング

お題

GAとかで使うルーレット選択を書く．いわゆる重み付きランダム選択．重みWがx:y:z=1:2:3の場合，x=1/6, y=2/6, z=3/6の確率で生成する．色んな言語で累積和，２分探索，乱数を扱う良い練習問題なのかな．

python

やり方としては重みの累積和の配列を作っておいて，乱数を１つ生成して，それがどこに入るかを探索する． pythonのnumpyを使うと累積和はcumsumで用意されているのでそれを使うと簡単．また，二部探索もbisectが用意されている．こうゆのはどんどん使う．

import random
import numpy as np
import bisect

def roulette_choice1(w):
    # generate accumulate weight list.
    # It can be also made by using numpy.cumsum()
    tot = []
    acc = 0
    for e in w:
        acc += e
        tot.append(acc)

    r = random.random() * acc
    # generate a random integer and search it in linear manner.
    # It can be also made by using bisect.bisect_right(tot, r)
    for i, e in enumerate(tot):
        if r < e:
            return i

if __name__ == "__main__":
    w = [1, 2, 3, 4, 5]
    test = [0]*len(w)
    for x in range(100000):
        i = roulette_choice1(w)
        test[i] += 1
    print test

C++

C++でpythonの例と同じことをやる場合には，cumsumは無いので適当に作るしかないけど，bisect相当のものとしてupper_bound(またはlower_bound)が使える．upper_boundはソート済みの配列に対して，指定したval以上となる最小の要素のイテレータを返す．探索は２分探索でやってくれる．upper_boundとlower_boundの違いは，”以上”か”より大きい”かの違い．

#include <vector>
#include <random>
#include <algorithm>
#include <iostream>

std::random_device rnd;
std::mt19937 mt(rnd());

template<typename T>
std::vector<T> cumsum(const std::vector<T> v)
{
  std::vector<T> ret;
  double sum = 0;
  for (auto d : v) {
    ret.push_back(sum+d);
    sum += d;
  }
  return ret;
}

template<typename T> 
size_t rouletteSelection(const std::vector<T> weight)
{
  std::vector<T> cumsum_weight = cumsum(weight);
  std::uniform_real_distribution<> rng(0, cumsum_weight.back());

  double r = rng(mt);
  auto it = std::upper_bound(cumsum_weight.begin(), cumsum_weight.end(), r);
  return (it - cumsum_weight.begin());
}


int main()
{
  std::vector<int> v(5);
  std::iota(v.begin(), v.end(), 1);

  std::vector<int> count = {0,0,0,0,0};
  for (int i=0; i<10000; ++i) {
    size_t j = rouletteSelection(v);
    count[j] += 1;
  }
  dumpVector(count); // dump a vector 
  return 0;
}

クロージャ with Python

目的

クロージャの使いドコロがいまいち分からないけれど，動作が（設定パラメタなどで）異なる関数を生成する，という目的で使うのかな？他にはクラス変数っぽい役割で使っているような例があるけどどんな意味があるのか良くわからなかった．

補足説明

pythonの関数はクロージャ機能を持つ．つまり，グローバルスコープ以外で定義された関数は，定義された時の変数情報を記録する．定義された，というのが重要で宣言時じゃない．（また，pythonは普通に関数を引数にしたり，関数を返り値にしたりできる．）

サンプル

genc_funcが関数を生成するもので，funcはグローバルスコープじゃないので，定義時の関数情報を内包する．例では，funcの外側のgen_funcの引数aが記録される．よって，add1はa=1として関数を定義しているので，引数xに1を加える関数となり，add10はa=10は引数xに10を加える関数となる．

def gen_func(a):
    # 'a' is closed.
    def func(x):
        return x+a
    return func

add1 = gen_func(1)
add10 = gen_func(10)

print add1(0)  #=>1
print add10(0) #=>10

pythonでのスレッド

基本事項

  • pythonのthreadは複数コアに跨がらないので，CPU空き時間がないような並列化では意味が無い．そのような場合にはmultiprocessingモジュールを使う．
  • threadingモジュールを使って、threading.Threadのサブクラスを作るか、インスタンスを作る，の２つの方法がある．サブクラスでやる場合に，runに引数は渡せない(runはselfだけを取る)．コンストラクタ(init)で渡しとく必要があることに注意する．
  • 処理はrunメソッドをオーバライドして、startメソッドで実行する。
  • 同期はjoinメソッド。
  • メインスレッド終了時にサブスレッドを終了するにはsetDaemon(True)にする。
  • 外からスレッドを止める(kill, stop)方法は無い(古いVersionには_stop()が隠しメソッドとしてあったが今は無い)．そもそも外からいきなりスレッドをkillことは良くない． 　方法としては，「ストップフラグをスレッドに持たせておいて，それを定期的に確認させる．ストップフラグは外部からいじれるようにしておく」．そして，外部からいじるためにthreading.Eventを使う．Eventはスレッドセーフな（複数のスレッドがアクセスしてもちゃんと同期をとってくれる）変数なのでこれをフラグとして使うんだ．
# Threadを継承したクラスを作成する場合
import threading
import time

class MyThread(threading.Thread):
    def __init__(self, i):
        super(MyThread, self).__init__()
        self.i = i
        self.stop_event = threading.Event()
        self.setDaemon(True)
        
    def stop(self):
        self.stop_event.set()

    def run(self):
        print 'START {}'.format(self.i)
        i = 0
        while not self.stop_event.is_set():
            print 'I am no. {}'.format(self.i)
            time.sleep(self.i)
            i += 1
        print 'END {}'.format(self.i)

if __name__ == '__main__':
    th1 = MyThread(1)
    th2 = MyThread(2)
    th1.start()
    th2.start()
    time.sleep(3)
    th1.stop()
    th2.stop()
  • 特定のメソッドをスレッド化したい場合には下記のようにインスタンスを作成する．引数としては，targetでスレッドを指定，nameでスレッド名，argsでメソッドの引数を渡す．
# メソッドをスレッドとして実行
import threading

def say_hello(situation, who):
   print "Good {} {} !".format(situation, who)

if __name__ == '__main__':
   th1 = threading.Thread(target=say_hello, name='morning_thread', args=('morning', 'Nobunaga'))
   th2 = threading.Thread(target=say_hello, name='eveningthread', args=('evening', 'Shingen'))
   th1.start()
   th2.start()

DEAPのメモ

概要

  • 拡張性を最も意識している。データ構造もアルゴリズムもほとんどがカスタマイズ可能。
  • 問題や遺伝子の型を扱う、creator, 遺伝子操作を扱うtoolboxが中心モジュール。
  • インストール
    • 描画用にvpythonも入れておく(pipでコケる場合は面倒なのでaptで入れてしまう．)
  $ pip install deap 
  $ pip install vpython

使い方のメモ

  1. 大枠
    まず，問題表現（個体表現）と適応度を定義する．次に，個体と集団の初期化方法を登録．突然変異，交叉，適応度評価方法，淘汰方法を登録．ここまでできたらGAフレームワークを回すだけ．

  2. 個体表現，適応度の定義

    • Deapでは個体表現，適応度をそれぞれ型(クラス)として定義する．自分で決められるところがDeapのいいところだけど，全部自分でやるのは大変だ．そこで使うのがCreator．Creatorは型定義を簡単にしてくれるクラス製造機みたいなものだ．
    • createの使い方はcreate('type-name', base-type, optional)。optionalは作る型に応じて付ける。つまり，最低限型の名前とそのベースクラスが必要ということだ．オプショナル引数はたぶんインスタンス変数が追加されているんだと思う。
    • 問題の型：creator.create("FitnessMin", base.Fitness, weights=(-1.0,))．weightは単目的であってもタプルで渡すこと．最小化の場合は-1.0, 最大化の場合は1.0にする．多目的の場合はそれぞれの重みを与える．
    • 遺伝子の型：creator.create("Individual", list, fitness=creator.FitnessMin)
    • これをやっておいて、initRepeatコンストラクタなどを活用して別途個体や集団の初期化方法を定義する。initRepeatはinitRepeat(container, constructor, n=num)とすると、containerにnum個のconstructorで生成するオブジェクトを格納する。つまり、コンストラクタを繰り返すもので、遺伝子の初期化などに使えるんだ。
    • initIterate(container, generator)．generatorが生成するオブジェクトをcontainerに代入して返す
    • 木のコンストラクタは3種ある
      • genFull:左右の葉の高さが等しい
      • genGrow:左右の葉の高さが等しくはないことも許す
      • genHalfAndHalf:genFullとgenGrowが半々で発生する。

  3. ToolboxにGAで使用するお決まりの関数の登録

    • toolbox.register("mutate", tools.mutGaussian, mu=0, sigma=1, indpb=0.1)．などとして、mutate, crossoverとかの中身を割り当てる。こうすることで、アルゴリズムのインターフェイスを統一的にしようって狙いだと思われる．統一した名前で突然変異,交叉を定義しておけば，それを回す上位ルーチンはフレームワーク化できる．
    • 突然変異や交叉だけでなく，個体の生成方法，集団の生成方法も定義する．この際に，initRepeatが活躍する．

  4. fitnessのEvaluationを書く

    • これだけは絶対に自分で書く部分

  5. アルゴリズムを書く

    • mutate, crossoverとか、決められた名前で関数をtooboxに定義しておけば algorithm.eaSimpleとかで定形ループは回してくれる。
    • 定形でないことをやりたいなら、ループを自分で書く。
    • 突然変異はtoolの中にいろいろある。注意点は、コピーしないと元のものが変化する。 mutant = toolbox.clone(ind1) ind2, = tools.mutGaussian(mutant, mu=0.0, sigma=0.2, indpb=0.2) del mutant.fitness.values child1, child2 = [toolbox.clone(ind) for ind in (ind1, ind2)] tools.cxBlend(child1, child2, 0.5) del child1.fitness.values del child2.fitness.values

主要モジュール

  • Creator
    • fitnessや遺伝子などはすべてcreatorを使って、’ある型’として生成する。
    • APIはcreator.create('name-of-type', base-class, optional-args)
    • Fitnessを定義する時は、weightを渡す。特に、weightは多目的な場合のためにタプルで渡す。 最大化の時は正、最小化の時は負の値。 ie) creator.create('FitnessMin', baseFitness, weight=(-1,))
    • 遺伝子(Indivisual)を定義する時は上で定義したfitnessをfitness属性として渡す。 ie) creator.create("Individual", list, fitness=creator.FitnessMin)
  • Toolbox
    • registerで関数と、そのエイリアスを登録する。３つ目以降の引数は登録する関数のパラメタ として渡される。
    • APIは　register('alias-name', function, function-args)
  • initRepeat
    • 3 argな関数。arg1はデータを入れるコンテナ,arg2は個々の生成関数,arg3は個数

GPについて

  • DEAPはGAだけでなくGPも対応している．
  • 終端記号と原始記号を定義する。終端記号は特に、引数(変数)と定数に分類する。
  • トップ関数名や引数の数はgp.PrimitiveSet('top-name', num-of-args)で設定。
  • 引数はデフォルトはARG0, ARG1,..になっている．必要ならrenameArgumentsで名前を付ける。
  • 定数はすべてを指定する必要はなく、例えばこの範囲の整数、とかはEphemeralConstantとして 定義する。下の例は,-1,0,1のどれかの定数が入る。
pset = gp.PrimitiveSet("MAIN", 1)  
pset.renameArguments(ARG0='x')
pset.addPrimitive(operator.add, 2)
pset.addEphemeralConstant("rand101", lambda: random.randint(-1,1))
  • 木のパースのときにスタックが溢れるので、90以上の木の高さは扱えない
  • 木の描画
    • networkxとgraphvizの例があるけど、networkxも結局graphvizで描画するっぽい。 ならばgraphvizを入れようということで、
      • sudo pip install pygraphviz だとimportErrorでコケる。
      • pip install pygraphviz --install-option="--include-path=/usr/include/graphviz" --install-option="--library-path=/usr/lib/graphviz/"
    • sudo apt-get install libgraphviz-dev を先に入れておく。

サンプル

DEAPはマルチプロセッシング含めて，たくさんのサンプルが提供されている．それを見ながら理解するのが一番分かりやすいと思う．そこで，OneMaxはドキュメントにも説明があるので，ナップザック問題のサンプルを例題にコメントを追加する．元のサンプルには最大アイテム数の制約があるがそれは無くしている．また，ナップザックの大きさを明に与えずに多目的関数として定義している．

import random
import numpy

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

IND_INIT_SIZE = 5 
MAX_WEIGHT = 50
NBR_ITEMS = 20

# 決まった乱数を出すために種は指定しておく
random.seed(64)

# 問題の入力（アイテムの重さと価値）定義
# ここでは乱数でアイテムを生成する．
items = {}
for i in range(NBR_ITEMS):
    items[i] = (random.randint(1, 10), random.uniform(0, 100))

# 問題定義（適応度と個体）
creator.create("Fitness", base.Fitness, weights=(-1.0, 1.0))
creator.create("Individual", set, fitness=creator.Fitness)

toolbox = base.Toolbox()

# 遺伝子生成ルーチン登録（サンプルでは遺伝子をアトリビュートと呼んでいる）
# 一つの個体は複数の遺伝子から構成される．つまり個体の各ビットのこと．
toolbox.register("attr_item", random.randrange, NBR_ITEMS)

# 個体の生成ルーチン登録
# 個体はナップザックに含まれるアイテムの集合として定義している．
toolbox.register("individual", tools.initRepeat, creator.Individual, 
    toolbox.attr_item, IND_INIT_SIZE)

# 集団の生成ルーチン登録
# 何個生成するのかの引数は明に与えない．
toolbox.register("population", tools.initRepeat, list, toolbox.individual)

# 適応度の計算ルーチン定義
def evalKnapsack(individual):
    weight = 0.0
    value = 0.0
    for item in individual:
        weight += items[item][0]
        value += items[item][1]
    return weight, value

# 交叉ルーチン定義
def cxSet(ind1, ind2):
    .... # 省略
    return ind1, ind2
    
def mutSet(individual):
    ... # 省略
    return individual,

# GA基本オペレータの登録（適応度計算，交叉，変異，選択）
toolbox.register("evaluate", evalKnapsack)
toolbox.register("mate", cxSet)
toolbox.register("mutate", mutSet)
toolbox.register("select", tools.selNSGA2)

def main():
    random.seed(64)
    NGEN = 50
    MU = 50
    LAMBDA = 100
    CXPB = 0.7
    MUTPB = 0.2
    
    pop = toolbox.population(n=MU)
    hof = tools.ParetoFront()　　　# 多目的最適化としているから

    # ログ関数の登録
    stats = tools.Statistics(lambda ind: ind.fitness.values)
    stats.register("avg", numpy.mean, axis=0)
    stats.register("std", numpy.std, axis=0)
    stats.register("min", numpy.min, axis=0)
    stats.register("max", numpy.max, axis=0)
    
    algorithms.eaMuPlusLambda(pop, toolbox, MU, LAMBDA, CXPB, MUTPB, NGEN, stats,
                              halloffame=hof)
    
    return pop, stats, hof
                 
if __name__ == "__main__":
    main()

Gazebo memo

  • worldファイルにモデル、センサなどすべての要素が記載される。フォーマットはSDF。
  • model(ODEの剛体と幾何両方含むもの)もSDFで定義される。worldファイルからインクルードする。
  • gzserver <world_file>でシミュレーションを実行する。
  • GUIツールとしてgzclientが使える．ビュワーとしてもインタラクティブな反映ツールとしても使える。
  • プロセス間通信はgoogle::protobufsとソケットの結合。
  • ~/.gazeboにモデルを置けば呼び出せたり、スクリーンキャプチャもここに置かれる。

モデルの作り方

  • ~/.gazebo/models/にディレクトリを作って、その中にmodel.configとmy_model.sdfを作成する。

vpythonのメモ

インストール

pip版はエラーでうまく使えなかった．aptでインストールする．

$ sudo apt-get install python-visual

　メモ

  • 座標系は右手系（xは→，yは↑，zは手前方向）（ODEも右手系なので悩む必要はない）
  • 右マウスぐりぐりで回転．両方押しでスパン．
  • cylinder, arrow, cone, pyramid のposは一つの端面，対してbox, sphere, ring の posは中心．
  • curveでラインをかけるけど、順次curveオブジェクトに点を追加(append)することで直線を伸ばしていくことができる．
  • rate(100)は前回のrateから「最大」1/100秒待つ。「最大」というのが重要で、sleep(100)みたいに1/100[sec]待つわけじゃなくて、あくまでも前のrate()呼び出しから最大で待つ時間。だから、例えば他の処理ですでに0.005秒使っていたら、残りは0.005しか待たないし、そもそも0.01をすでに使いきっていたら止まらない。
  • テキストを入れたいときはlabel．

　例

サンプルとして公式のチュートリアルの例が下記． ボールが箱のなかで跳ね返るアニメーション．進む向きと大きさ（速度ベクトルの軸と大きさ）と軌跡を描く．注意点としてはsphereにverocity, trailなどのプロパティは無くて，その場で定義している．

import random
from visual import *


def draw_box(w):
    box(pos=(w/2, 0,0),  size=(0.1, w, w), color=color.red)
    box(pos=(-w/2, 0,0), size=(0.1, w, w), color=color.red)
    box(pos=(0,w/2,0), size=(w, 0.1, w), color=color.green)
    box(pos=(0,-w/2,0), size=(w, 0.1, w), color=color.green)
    box(pos=(0,0,-w/2), size=(w, w, 0.1), color=color.blue)

box_size=10
draw_box(box_size)

# ball
ball = sphere(color=color.white, radius=0.2, pos=(0,0,0))
ball.v = vector(random.random(), random.random(), random.random())
dt = 0.3
v_arrow = arrow(pos=ball.pos, axis=ball.v*3, color=color.yellow)
ball.trail = curve(color=ball.color)
l = label(pos=(-box_size/2, -box_size/2, box_size/2), text='pos=', box=0)

while 1:
    rate(100)
    ball.pos = ball.pos + ball.v*dt
    v_arrow.pos, v_arrow.axis= ball.pos, ball.v*3
    ball.trail.append(pos=ball.pos)
    if not (-box_size/2 < ball.x < box_size/2):
        ball.v.x = -ball.v.x
    if not (-box_size/2 < ball.y < box_size/2):
        ball.v.y = -ball.v.y
    if not (-box_size/2 < ball.z < box_size/2):
        ball.v.z = -ball.v.z
    l.text = 'pos=%s' % ball.pos

    if scene.kb.keys:
        key = scene.kb.getkey()
        if key == 'q' or key == 'Q':
            display.exit = 1
            exit()