wanimaru47's diary

プログラミング等々

線形計画から始める数理計画入門

この記事は信州大学のものづくりサークル「kstm」のアドベントカレンダー11日目の投稿です。

qiita.com

普段は組合せ最適化の研究をしている学生です。組合せ最適化って難しそうだよねって思う人もいると思います。僕もそう思います。しかし、世の中にはそんな僕のような人たちを助けてくれるツールがたくさん存在します。その一つに「数理計画ソルバー」が存在します。

そもそも数理計画ってなに?ってなりますよね。ざっくり説明すると数理計画は目的制約を表す式で記述され、制約の中で目的を最適化(最大化もしくは最小化)します。目的を表す式を目的関数、制約を表す式を制約条件って呼んだりします。

数理計画には様々な分類があります。ここでの分類は、目的関数と制約条件をどのような式で表すかによって分けられています。最も代表的なものは線形計画でしょう。線形計画は、目的関数も制約条件がすべて1次式であるようなものの総称です。線形計画は以下のようなものです。

f:id:wanimaru47:20181210140519p:plain
式1 : 線形計画問題を一般的に表した式
ベクトルx線形計画問題の解を表している。もう少し具体的に見てみましょう。
f:id:wanimaru47:20181210222128p:plain
式2 : 具体的な例
上の例は例1のA, c, b, xの値が次のような場合です。
f:id:wanimaru47:20181210222229p:plain
式3 : 各変数についての具体値
このような式の作り方(定式化)は解きたい問題によって変わってくるのここでは説明しません。

さて、このような問題を線形計画問題へ定式化していかに解くかということが次の問題になります。そこで活躍するのがソルバーです。ソルバーは定式化された問題を広く知られた効率の良いアルゴリズムを用いて解いてくれます。自分でアルゴリズムを設計・実装をしなくていいことが最大の利点でしょう。

ここでは、PuLPと呼ばれるソルバーを用いて実際に問題を解いてみましょう。PuLPPythonで線形計画などの数理計画の問題を解くことのできるライブラリ(ソルバーのラッパー)です。一番手軽に実行できるので選んでみました。上記の式2Pythonプログラムで書くと以下のようになります。

from pulp import *

prob = LpProblem("Sample", LpMaximize)

# 変数の宣言
x1 = LpVariable("x1", 0, None, LpContinuous)
x2 = LpVariable("x2", 0, None, LpContinuous)
x3 = LpVariable("x3", 0, None, LpContinuous)

c = [2, 5, 3]

A = [[1, 1, 0],[0, 1, 1],[1, 0, 0]]

B = [5, 10, 4]

# はじめに追加した式が「目的関数」
prob += lpDot(c, [x1, x2, x3]) # Maximize : 2 * x_1 + 5 * x_2 + 3 * x_3

# 二つ目以降の式は「制約条件」
for i, b in enumerate(B):
    prob += lpDot(A[i], [x1, x2, x3]) <= B[i]
    
prob.solve()

print(prob)

print("x1 =", x1.value())
print("x2 =", x2.value())
print("x3 =", x3.value())

わざわざ解説することもないぐらいに直感的なライブライになっていると思います。内積lpDot関数を用いて表すことが出来ます。

これを実際に実行してみると以下のような結果がえられます。x1, x2, x3の値が解になっているわけです。

Sample:
MAXIMIZE
2*x1 + 5*x2 + 3*x3 + 0
SUBJECT TO
_C1: x1 + x2 <= 5

_C2: x2 + x3 <= 10

_C3: x1 <= 4

VARIABLES
x1 Continuous
x2 Continuous
x3 Continuous

x1 = 0.0
x2 = 5.0
x3 = 5.0

変数の宣言は配列ですることも出来ます。配列で宣言する時の例を以下に書いてみました。加えて、lpSum関数(配列に格納されたPuLP変数の総和の式を返す関数)も使っていました。

from pulp import *

prob = LpProblem("Sample", LpMaximize)

# 配列で宣言をした時
X = LpVariable.dicts("x", range(3),  0, None, LpContinuous)

c = [2, 5, 3]

A = [[1, 1, 0],[0, 1, 1],[1, 0, 0]]

B = [5, 10, 4]

# はじめに追加した式が「目的関数」
prob += lpSum([c[i] * X[i] for i in range(3)]) # Maximize : 2 * x_1 + 5 * x_2 + 3 * x_3

# 二つ目以降の式は「制約条件」
for i, b in enumerate(B):
    prob += lpSum([A[i][j] * X[j] for j in range(3)]) <= B[i]
    
prob.solve()

print(prob)

print("x1 =", x1.value())
print("x2 =", x2.value())
print("x3 =", x3.value())

こんな感じ線形計画問題PuLPを使って簡単に解くことが出来ます。PuLPは他にも混合整数計画を解くことが出来ます。整数計画は変数のとる値が整数値であり、目的関数と制約条件が1次式で記述された問題です。混合は有利数値と整数値が混ざった問題を表しています。実に便利です。

落とし所がわからなかったのでこの辺で終了とします。皆さんもソルバーライフを楽しんでみてください。以下はソルバーの紹介と文献紹介になっています。

ソルバーの例

developers.google.com

www.ibm.com

www.gurobi.com

www.wolframalpha.com

ACM-ICPC2017 アジア地区予選

qiita.com
この記事は、20日目の記事です。


ICPCのアジア地区予選に参加してきました。
結果は、3完32位でした。

1日目

JAXA見学

一日目は、JAXA見学をしてきました。f:id:wanimaru47:20171216125603j:plain

普段は見ることのできない宇宙開発をみれて楽しかったです。

練習フェーズ

普段はHHKBでコーディングをしているので、普通の英字キーボードだとCtrlキーの位置が違って全然コーディングができませんでした。
宿に戻ってから色々調べてCtrlとCapsを入れ替える方法を調べました。
(ただ、戻し方を調べずに本番の時に苦労しました。)

2日目

2日目は、本番でした。自分の担当は、問題読解と解法を考える役でした。
コーディングは、 @goryudyuma
@o0p1q が担当しました。

A問題とB問題は、初めからコーダーの2人に任せて、D問題から後を読みました。
あとで、C問題が難解だと言われて読みました。

A問題とB問題は、コーダーの2人がACしてくれたのでとても助かりました。
C問題は、2ACしたところで自分が読みました。問題文が並列性に関して誤読していて、理解に時間がかかってしまいました。
読解したあとは、簡単に解けそうだったのでコーディングは任せました。

その後、順位表からI問題、F問題が簡単だとわかったので、解法を考え始めました。F問題を自分がやって、I問題をもう2人にやってもらいました。
F問題はしばらくして、解法が思いつきましたが、"SOSO"と"SAD"の判定が詰められなくて死んでいました。
I問題は、最大の方がわかっていなかったらしく自分も参加して3分ほどで解法は思いつきましたが、その時点で残り10分になってしまっていてACできませんでした。

まだまだ、鍛錬が足りていなかったので、取り組む問題を間違えたり、解法が詰められなかったりして、もっと頑張りたかったなと思いました。

久しぶりの長時間コンテストだったの楽しめてよかったです。
あと、C問題以降はほぼ全て自分が問題文を読んだので、それだけでめっちゃ疲れました。

気合いを入れる!!

新年度が始まり時間が経ってしまったが今年度の抱負を書きたいと思う。

去年はICPCで国内予選を突破でき、とてもいい結果が残せたと思う。いい結果が残せた一方で今一歩競技プログラミングの力を伸ばせなかった。国内予選が終わった後の喪失感みたいな燃え尽き症候群に襲われ、時間を無駄にとかってしまったと思う。


でも、そんな中でも今までほとんど挑戦してこなかったインターンに参加して、開発とか企画とか色んなことを学べて有意義な時間を過ごせた。

インターンのときに食べた、とてもうま辛い麻婆豆腐


今年はやっと研究室にも配属され本格的にアルゴリズムを学べるような気がしている。新年度が始まって一ヶ月が経ちようやく「研究」のイメージがぼやぼや見え始めた。そしてあと一ヶ月もすれば今年もICPCがやってくる。今年も国内予選を突破したい。でも去年ほど真面目に競技プログラミングに打ち込めていない自分がいて残念な感じがある。でも、まだ一ヶ月はあるのでどんどん問題解いて国内予選を当たり前のように突破したい。
(去年、ギリギリだった奴が何を言っているんじゃいって感じだ…)

今年もアジア地区いくぞ!!

そして、修士課程も見越して今年はたくさん論文とかアルゴリズムの本に触れて、今後の3年間を有意義なものにしたい。研究はハフマン
木生成問題について拡張を行うことをテーマにした。論文を読み始めて知らないアルゴリズムとか定理とか出てきてもっと成長できそう。今年はハフマン木のプロになるぞw

ACM-ICPC2015 つくば大会

(かなり時間が経ってしまったけど参加記書きます。)

ICPC参加三年目にしてようやく地区大会に駒を進めることができました。
結果は予選順位とほぼ変わらず40位という結果に終わってしまいました。

1日目はJava Challengeとレセプションが行われました。
今年のJava ChallengeはC++のサンプルコードがあったのでC++で参加しました。
Java ChallengeなのにJavaじゃなかったw)
レセプション時は、他大学の人たちはチーム紹介の際に、
かなり笑いを取りに行っていて面白かったです。
(来年は自分のチームもちょっと遊んなだ感じにしてもいいかも)

2日目は本番です。
A問題見て解けそうだったので、書いてみました。WAでした。
かなり、バグらせて結局チームメイトに交代…AC
(あとで、解説聞いてもっと考察するべきだったと反省)
B問題の問題概要を聞いてまだ解法が定めっていなかったので、
解法を考えてコーディングしてもらうことに。
C問題はチームメイト任せてD問題を読むことに、
解法を考えつつ、B問題のコーディングが詰まっていたのでサポート。
(自分が良かれと思って入れた)条件式が邪魔をしていたらしく
if文消して無条件に更新対象としてコーディングしなおしたらAC
結局、この日は2完でした。
(自分がかなり足を引っ張ってペナルティーを出してしまったので申し訳なかった。)

3日目
高エネルギー加速器研究所とサイバーダイン社の展示場の見学に行ってきました。
サイバーダインでサイボーグスーツ的なものの操縦体験ができたので満足でした。


反省:
予選に比べてAOJ−ICPC埋めが行えずに、
完全にチームメイトの足を引っ張ってしまいました。
最近、全然問題解いてなくてかなり競プロ力が低下してきているので
今一度気合を入れ直して頑張りたいと思います。

ACM-ICPC2015 国内予選参加記(真面目)

ACM−ICPC2015国内予選が今月26日に行われました。
今年は去年と同じメンバ同じチーム名で参加しました。
結果は全体で41位、予選通過順位31位でした。

順位と結果 | ACM-ICPC 2015 Tsukuba


今年自分がやったこと
・AOJ−ICPC埋め
・勉強会を催した
  (自分の知識のアウトプットとモチーベーションの維持向上)

 今年は自分が中心となり勉強会を4月中旬から予選がある週まで週一で実施。初めの方は自分がスライドを作成してアルゴリズムと実装方法の紹介。5月の第3週に大学のOBの方に勉強会にお越しいただきいい刺激になった。(多分他の参加者も刺激になったのかと)その後は過去問で自分が面白かったと思った問題を出題もしくは解けない問題を解く時間にしてペアプロとかやりながらチーム力を向上した。(そう思ってる。そう思いたい。)あと、模擬国内の開催を例年大学で行わないものの大学で会場を設置して本番環境でやった。事前に会場の勝手が分かったので当日慌てることは特になくよかった。

 AOJ-ICPC埋めは春休みあたりから真剣にやり始めた。春休み開始時は十数問しか通してなかった。結果的に88問を通して国内予選当日に。100問ぐらいは埋めたかったものの全然間に合わずかなり不安だったが結果的にいい方向に。数日前に解いた問題を持ちにC問を(自分の中では)速攻で通すことができた。その問題が以下の2問。

JAG-channel | Aizu Online Judge
Broken Cipher Generator | Aizu Online Judge

以上の問題で結果が出せました。JAGの方々には大変感謝しています。また、今回の大会に向けて沢山の先輩方から練習問題を選出して頂いたりして質の高い練習が行えました。大変感謝しています。

[大会当時]
今年は去年とは違いA,B問題に関してはほとんどNOタッチでチームメイトに任せて、C,D問題の解法を考えるこちに専念しました。A,B問題がACしたところでC問題の説明をしてコーディングして貰うことに、しかしちょっと時間がかかりそうで今年は早めにD問題のコーディングに入りたかったので自分がコーダーになるこちに、バグらせることなくACして暫定31位に。今年はチームメイトもバグを出さずに3完したのでペナルティをかなり少なくすることができました。

[アジア地区に向けて]
今回の国内予選ではD問題が解けるレベルではなかったし、解けない問題を避ける判断力などがたりなかった気がします。まずはアルゴリズムと実装力の強化を行いたいと思います。あとは後期になったら出来る限りチームの連携を強化したいと思います。

code thanks festival 2014 B日程

B日程参加してきました。

AからEまでやるだけでF、G、Hの三問が本命の問題という感じがしましたね。
F問題まではけっこう順調に進んで、焼く肉食べれるかと思いましたが、
G問題とH問題が解けなくて残念な結果になりました。

オンサイトでは13位で、オンライン合わせると70位ぐらいで
今まででは、なかなかいい感じの結果が残せた方なので
楽しめた数少ないコンテストの一つになった気がしますw

F問題だけ載せておきます。

#include <iostream>
#include <sstream>
#include <cstdio>
#include <map>
#include <vector>
#include <string>
#include <algorithm>
#include <functional>
#include <list>
#include <stack>
#include <queue>
#include <cstdlib>
#include <tuple>
#define INF 1e9
#define EPS 1e-9
#define MOD 1000000007
using namespace std;

typedef list<int> L;
typedef pair <int,int> P;
typedef vector<int> V;
typedef queue<int> Q;
typedef stack<int> S;
typedef map<string,int> M;

string s, v, t;

int main()
{
    cin >> s >> v >> t;
    V res(s.size() + 1, 0);

    bool flag = true;
    for (int i = 0; i < v.size(); i++) {
        if (v[i] != s[i]) flag = false;
    }
    if (flag) res[v.size()] = 1;

    flag = true;
    for (int i = 0; i < t.size(); i++) {
        if (t[i] != s[i]) flag = false;
    }
    if (flag) res[t.size()] = 1;


    for (int i = 1; i < s.size(); i++) {
        if (res[i]) {
            int j = 0;
            flag = true;
            while (j < v.size()) {
                if (s[i+j] != v[j]) flag = false;
                j++;
            }
            if (flag) res[i+v.size()] += res[i];
            res[i+v.size()] %= MOD;
            j = 0;
            flag = true;
            while (j < t.size()) {
                if (s[i+j] != t[j]) flag = false;
                j++;
            }
            if (flag) res[i+t.size()] += res[i];
            res[i+t.size()] %= MOD;
        }
    }
    cout << res[s.size()] << endl;

    return 0;
}

ICPC2014 国内予選 反省

ICPC2014国内予選が終わって一週間経ってようやく反省書きます。
問題はA問、B問の二問ACでした。
去年の0問と比べればいい結果でしたが、一問目に時間を割きすぎてしまって
残念な結果になってしまいました。
また、来年に向けて頑張りたいです。