9-2. 遺伝的アルゴリズム（vcopt）でスケジューリング問題を解く（12種のレシピ編）

import numpy as np
import numpy.random as nr
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from vcopt import vcopt

#スケジューリング問題の作成
def create_recipe(recipe_num):
    
    #レシピを入れる配列
    recipe = np.zeros((recipe_num, 4, 25), dtype=int)
    
    #レシピの生成
    for i in range(recipe_num):
        #まず０分目は人間の作業（＝マシンNo.0）
        recipe[i, 0, 0] = i + 1 #レシピIDは１から始める
        #その後、ランダムに遷移しては作業（２～５分）を４回繰り返す
        current_time = 1
        for j in range(4):
            next_machine = nr.randint(4)
            next_time = nr.randint(2, 6)
            #レシピに追加
            recipe[i, next_machine, current_time:current_time+next_time] = i + 1
            current_time += next_time
    
    return recipe

#スケジューリング問題の作成
nr.seed(1)
recipe_num = 3
recipe = create_recipe(recipe_num)
print(recipe)

#評価関数
def recipe_score(para):
    
    #全工程を敷き詰める配列
    process = np.zeros((4, len(para)*25), dtype=int)
    
    #各レシピを、paraの順番に０分から敷き詰める
    for i in range(len(para)):
        #このレシピを
        current_recipe = recipe[para[i]]
        #０分から順に
        for t in range(len(para)*25):
            #すでに敷き詰めたレシピと競合しないかチェックして敷き詰める
            mask = process[:, t:t+25] * current_recipe[:, :]            
            if np.sum(mask) == 0:
                process[:, t:t+25] += current_recipe[:, :]
                break
    
    #最後のプロセスの時刻を返す
    process_marge = np.sum(process, axis=0)
    return np.where(process_marge > 0)[0][-1]

#paraの可視化
def recipe_show_para(para, **info):
    
    #2-opt処理中の諸情報はinfoという辞書に格納されて渡されます
    #これらを受け取って使用することができます
    step_num = info['step_num']
    score = info['score']
    
    #全工程を敷き詰める配列
    process = np.zeros((4, len(para)*25), dtype=int)
    
    #各レシピを、paraの順番に０分から敷き詰める（評価関数と一緒）
    for i in range(len(para)):
        current_recipe = recipe[para[i]]
        for t in range(len(para)*25):
            mask = process[:, t:t+25] * current_recipe[:, :]            
            if np.sum(mask) == 0:
                process[:, t:t+25] += current_recipe[:, :]
                break
    
    #プロット
    plt.figure(figsize=(12, 2))
    for i in range(len(para)):
        #レシピの抽出
        a = np.where(process==i+1)
        #プロット
        plt.plot(a[1], a[0], 'o', color=cm.gist_rainbow(i/len(para)), label=i+1)
        plt.text(a[1][0] - 0.5, a[0][0] - 0.2, str(i+1))
    
    plt.xlabel('time'); plt.ylabel('machine')
    plt.xlim([0, 80]); plt.ylim([3.5, -0.9])
    plt.title('step={}, score={}'.format(step_num, score))
    #plt.savefig('save/{}.png'.format(step_num))
    plt.show()
    print()

#poolの可視化
def recipe_show_pool(pool, **info):
    
    #GA中の諸情報はinfoという辞書に格納されて渡されます
    #これらを受け取って使用することができます
    gen = info['gen']
    best_index = info['best_index']
    best_score = info['best_score']
    mean_score = info['mean_score']
    mean_gap = info['mean_gap']
    time = info['time']
    
    pool_len = len(pool[:100])
    para_len = len(pool[0, :])
    
    #pool（100体まで）の全工程を詰め込む配列
    pool_process = np.zeros((pool_len, 4, para_len*25), dtype=int)
    #詰め込む（伝われ…）
    for n in range(pool_len):
        for i in range(para_len):
            current_recipe = recipe[pool[n, i]]
            for t in range(para_len*25):
                mask = pool_process[n, :, t:t+25] * current_recipe[:, :]            
                if np.sum(mask) == 0:
                    pool_process[n, :, t:t+25] += current_recipe[:, :]
                    break
    
    #pool（100体まで）の最後のプロセスの時刻を縦線でプロット
    plt.figure(figsize=(12, 2))
    for n in range(pool_len):
        process_marge = np.sum(pool_process[n], axis=0)
        last = np.where(process_marge > 0)[0][-1]
        plt.plot([last, last], [-1, 4], '-k', alpha=(2.0/pool_len))
    
    
    #エリートの全工程を敷き詰める
    process = np.zeros((4, para_len*25), dtype=int)
    for i in range(para_len):
        current_recipe = recipe[pool[best_index, i]]
        for t in range(para_len*25):
            mask = process[:, t:t+25] * current_recipe[:, :]            
            if np.sum(mask) == 0:
                process[:, t:t+25] += current_recipe[:, :]
                break
    
    #エリートの全工程をプロット
    for i in range(para_len):
        #レシピの抽出
        a = np.where(process==i+1)
        #プロット
        plt.plot(a[1], a[0], 'o', color=cm.gist_rainbow(i/para_len), label=i+1)
        plt.text(a[1][0] - 0.5, a[0][0] - 0.2, str(i+1))
    
    plt.xlabel('time'); plt.ylabel('machine')
    plt.xlim([0, 80]); plt.ylim([3.5, -0.9])
    plt.title('gen={}, best={} mean={} time={}'.format(gen, best_score, mean_score, time))
    #plt.savefig('save/{}.png'.format(gen_num))
    plt.show()
    print()

#スケジューリング問題の作成
nr.seed(1)
recipe_num = 12
recipe = create_recipe(recipe_num)

#適当に順番を作成
para = np.arange(recipe_num)

#2-opt法で最適化
para, score = vcopt().opt2(para,                             #para
                           recipe_score,                     #score_func
                           0.0,                              #aim
                           show_para_func=recipe_show_para,  #show_para_func=None
                           seed=None)                        #seed=None

#結果の表示
print(para + 1)
print(score)

#スケジューリング問題の作成
nr.seed(1)
recipe_num = 12
recipe = create_recipe(recipe_num)

#パラメータ範囲
para_range = np.arange(recipe_num)

#GAで最適化
para, score = vcopt().tspGA(para_range,                       #para_range
                            recipe_score,                     #score_func
                            0.0,                              #aim
                            show_pool_func=recipe_show_pool,  #show_pool_func=None
                            seed=1)                           #seed=None

#結果の表示
print(para + 1)
print(score)