災害支援
プログラミングとは関係ないけど。
国連世界食糧計画
http://www.wfp.or.jp/
国外の災害支援では国連の組織はもっと活用されてもいいと思うんだ。
既に持続的な活動をして動ける状態のはずだし、クレジットで即時募金できるし。
もちろん特に地震に対してなら日本独自の支援の方法もあると思うけど。
ラムゼーの定理
ラムゼーの定理
http://ja.wikipedia.org/wiki/%E3%83%A9%E3%83%A0%E3%82%BC%E3%83%BC%E3%81%AE%E5%AE%9A%E7%90%86
これで以前のエントリが理解できるわけですね。
グラフ理論すごいな。プログラマは勉強して損はない分野だ。
峠のTaisaiがナイズルを(ry
某MMOとは関係ナイです。
粘菌っぽい動きをする何かを迷路に置いて挙動を観察してみるプログラムです。
プログラムの理論的な根拠はありません。
大雑把に連立拡散方程式で粘菌の伸び縮みと餌信号伝達を時間発展してみただけ。
densityが大きいと粘菌が濃く見えるような感じです。
迷路条件
- 粘菌は初期時刻でスタート位置にいます。
- 餌はスタート位置とゴール位置に等量置かれます。
粘菌の挙動
- 時間が経つごとに伸びていく
- あんまり伸びるとマズイのでそれなりに形を保とうとする(自己組織化っぽく)
- 餌があると信号を体内に伝達していく(信号は時間で減衰する)
class Nenkin class Point attr_accessor :x, :y def initialize(x, y) @x = x; @y = y end def ==(other) other != nil && @x == other.x && @y == other.y end end # 進行方向 @@directions = [Point.new(-1,0), Point.new(1,0), Point.new(0,-1), Point.new(0,1)] # イプシロン @@eps = 1.0e-10 # 進行時間単位 [時間] @@dt = 0.1 # 拡散速度 [時間^-1] @@velocity = 1.0 # 餌減衰係数 < 1.0 @@decay_food = 0.75 # 餌シグナル伝達相対速度 @@trans_food_signal_weight = 100 # 粘菌減衰係数 < 1.0 @@decay_density = 0.999 def initialize # マップ構造 @map = [] # 粘菌密度[無次元]:xy積分で1 @density = [] # 粘菌拡散速度 [時間^-1] @d = [] # 粘菌密度時間発展差分[距離^-2]:Δx=Δy=1.0として陽には現れない @density_delta = [] # 餌シグナル強度 @food = [] # 餌シグナル強度時間発展差分[距離^-2]:Δx=Δy=1.0として陽には現れない @food_delta = [] end def run load_map total_count = 1000000 output_count = 10 (total_count/output_count).times {|i| output_count.times{|j| decay_and_rethrow_food_signal @@trans_food_signal_weight.times{|k| trans_food_signal } move_density decay_and_reborn } printf "[Time=%5.5f]\n", i*output_count*@@dt print_map } end def load_map DATA.read.each_line {|line| @map.push line.chomp.split(//) } @map.each_with_index{|line,i| line.each_with_index{|c,j| case c when 'S' @start_point = Point.new(j,i) @map[i][j] = ' ' when 'G' @goal_point = Point.new(j,i) @map[i][j] = ' ' end } } @map.each{|m| @density.push Array.new(m.size, 0.0) @d.push Array.new(m.size, 0.0) @density_delta.push Array.new(m.size, 0.0) @food.push Array.new(m.size, 0.0) @food_delta.push Array.new(m.size, 0.0) } @d.each_with_index{|line,i| line.each_with_index{|c,j| @d[i][j] = is_passage?(Point.new(j,i)) ? @@velocity : 0.0 } } # 粘菌はスタート位置に配置される @density[@start_point.y][@start_point.x] = 1.0 # 餌はスタート位置とゴール位置に配置される @food[@start_point.y][@start_point.x] += 1.0 @food[@goal_point.y][@goal_point.x] += 1.0 end # 餌シグナルを全体に伝える # 拡散の重みとして粘菌密度比が掛けられる def trans_food_signal for y in 1..@density.size-2 for x in 1..@density[0].size-2 diffusion = 0.0 diffusion += (@food[y][x-1]*@density[y][x]/(@density[y][x]+@density[y][x-1]+@@eps)-@food[y][x]*@density[y][x-1]/(@density[y][x]+@density[y][x-1]+@@eps)) *@d[y][x-1] *@d[y][x] diffusion += (@food[y][x+1]*@density[y][x]/(@density[y][x]+@density[y][x+1]+@@eps)-@food[y][x]*@density[y][x+1]/(@density[y][x]+@density[y][x+1]+@@eps)) *@d[y][x+1] *@d[y][x] diffusion += (@food[y-1][x]*@density[y][x]/(@density[y][x]+@density[y-1][x]+@@eps)-@food[y][x]*@density[y-1][x]/(@density[y][x]+@density[y-1][x]+@@eps)) *@d[y-1][x] *@d[y][x] diffusion += (@food[y+1][x]*@density[y][x]/(@density[y][x]+@density[y+1][x]+@@eps)-@food[y][x]*@density[y+1][x]/(@density[y][x]+@density[y+1][x]+@@eps)) *@d[y+1][x] *@d[y][x] @food_delta[y][x] = @@dt * diffusion * @@velocity end end for y in 1..@density.size-2 for x in 1..@density[0].size-2 @food[y][x] += @food_delta[y][x] end end end # 一定確率で寿命となりスタート地点に再生成される def decay_and_reborn for y in 0..@food.size-1 for x in 0..@food[0].size-1 @density[y][x] *= @@decay_density end end @density[@start_point.y][@start_point.x] += (1.0-@@decay_density) * 1.0 end def decay_and_rethrow_food_signal for y in 0..@food.size-1 for x in 0..@food[0].size-1 @food[y][x] *= @@decay_food end end @food[@start_point.y][@start_point.x] += (1.0-@@decay_food) * 1.0 @food[@goal_point.y][@goal_point.x] += (1.0-@@decay_food) * 1.0 end def normalize_density s = 0.0 for y in 0..@density.size-1 for x in 0..@density[0].size-1 @density[y][x] *= 0.0 if @density[y][x] < 0.0 s += @density[y][x] end end for y in 0..@density.size-1 for x in 0..@density[0].size-1 @density[y][x] /= s end end end # 詳細釣り合い条件を満たす拡散方程式に従い粘菌密度が動く # 拡散の重みとして餌シグナル強度比が掛けられる def move_density r_walk = rand(10000)/100000.0 for y in 1..@density.size-2 for x in 1..@density[0].size-2 # 単純拡散 diffusion = 0.0 diffusion += (@density[y][x-1]-@density[y][x]) *@d[y][x-1] *@d[y][x] diffusion += (@density[y][x+1]-@density[y][x]) *@d[y][x+1] *@d[y][x] diffusion += (@density[y-1][x]-@density[y][x]) *@d[y-1][x] *@d[y][x] diffusion += (@density[y+1][x]-@density[y][x]) *@d[y+1][x] *@d[y][x] # 餌シグナルによる収縮(重み付き拡散方程式) food_shrinkage = 0.0 food_shrinkage += (@density[y][x-1]*@food[y][x]/(@food[y][x]+@food[y][x-1]+@@eps)-@density[y][x]*@food[y][x-1]/(@food[y][x]+@food[y][x-1]+@@eps)) *@d[y][x-1] *@d[y][x] food_shrinkage += (@density[y][x+1]*@food[y][x]/(@food[y][x]+@food[y][x+1]+@@eps)-@density[y][x]*@food[y][x+1]/(@food[y][x]+@food[y][x+1]+@@eps)) *@d[y][x+1] *@d[y][x] food_shrinkage += (@density[y-1][x]*@food[y][x]/(@food[y][x]+@food[y-1][x]+@@eps)-@density[y][x]*@food[y-1][x]/(@food[y][x]+@food[y-1][x]+@@eps)) *@d[y-1][x] *@d[y][x] food_shrinkage += (@density[y+1][x]*@food[y][x]/(@food[y][x]+@food[y+1][x]+@@eps)-@density[y][x]*@food[y+1][x]/(@food[y][x]+@food[y+1][x]+@@eps)) *@d[y+1][x] *@d[y][x] # 自己組織化による収縮(物理的根拠無しな適当方程式) self_shrinkage = 0.0 self_shrinkage -= (@density[y][x-1]*@density[y][x-1]/(@density[y][x]+@density[y][x-1]+@@eps)-@density[y][x]*@density[y][x]/(@density[y][x]+@density[y][x-1]+@@eps)) *@d[y][x-1] *@d[y][x] self_shrinkage -= (@density[y][x+1]*@density[y][x+1]/(@density[y][x]+@density[y][x+1]+@@eps)-@density[y][x]*@density[y][x]/(@density[y][x]+@density[y][x+1]+@@eps)) *@d[y][x+1] *@d[y][x] self_shrinkage -= (@density[y-1][x]*@density[y-1][x]/(@density[y][x]+@density[y-1][x]+@@eps)-@density[y][x]*@density[y][x]/(@density[y][x]+@density[y-1][x]+@@eps)) *@d[y-1][x] *@d[y][x] self_shrinkage -= (@density[y+1][x]*@density[y+1][x]/(@density[y][x]+@density[y+1][x]+@@eps)-@density[y][x]*@density[y][x]/(@density[y][x]+@density[y+1][x]+@@eps)) *@d[y+1][x] *@d[y][x] @density_delta[y][x] = @@dt * (diffusion + food_shrinkage *2.73 + self_shrinkage/2.73) * @@velocity end end for y in 1..@density.size-2 for x in 1..@density[0].size-2 @density[y][x] += @density_delta[y][x] end end normalize_density end def print_map total_density = 0.0 @map.each_with_index{|line,y| line.each_with_index{|c,x| if @start_point.x == x && @start_point.y == y print 'S' elsif @goal_point.x == x && @goal_point.y == y print 'G' elsif c == ' ' if @density[y][x] > 1.0e-1 print '1' elsif @density[y][x] > 1.0e-2 print '2' elsif @density[y][x] > 1.0e-3 print '3' elsif @density[y][x] > 1.0e-4 print '4' elsif @density[y][x] > 1.0e-5 print '5' elsif @density[y][x] > 1.0e-6 print '6' elsif @density[y][x] > 1.0e-7 print '7' elsif @density[y][x] > 1.0e-8 print '8' elsif @density[y][x] > 1.0e-9 print '9' else print ' ' end else print c end total_density += @density[y][x] } print "\n" } printf "-- total density = %5.5f \n", total_density print "==density==\n" @density.each{|line| line.each{|c| printf "%1.5f\t", c } print "\n" } print "==food==\n" @food.each{|line| line.each{|c| printf "%1.5f\t", c } print "\n" } end def is_passage?(point) @map[point.y][point.x] == ' ' end end Nenkin.new.run __END__ ************************** *S* * * * * * * ************* * * * * ************ * * * * ************** *********** * * ** *********************** * * G * * * *********** * * * * ******* * * * * * **************************
実行するとこんな感じの出力が続きます。
[Time=2200.00000] ************************** *S*3*33344444455555566666* *2*3*33*44*************66* *2*233*4444************66* *2222*4444445555566666666* **************5*********** *777777666666666666667777* **7*********************** *776666*65555554444444G44* *77*666666***********4*44* *7766*66666655*******4*44* *7766666*6665555554444444* ************************** -- total density = 1.00000 ==density== 0.00000 0.00000 0.00000 ... 0.00000 0.10715 0.00000 ... ... ==food== 0.00000 0.00000 0.00000 ... ...
迷路の数字は粘菌の密度1(濃い)〜9(薄い)を表わしています。
densityは各点における粘菌密度の実データ(全点足すと1)
foodは各点における餌シグナル強度(全点足すと2のはず)
これが必要なところだけに収縮してくれれば良いのだが。
いろいろいじってみたがちゃんと勉強してから書き直した方がいいみたいだ。
現状は迷路いっぱいに広がるか、スタートゴールの2点に分裂するかになってしまうので、
ループ解消の何らかの機構が必要かも。
一連の流れで
正直この一連のエントリの問題が楽勝なんじゃなかろうか、などと思っていましたが、自分の力不足を感じる結果でした。
若い頃から華々しい成果を出すような才能を持たない凡人なので、日々虚心に努力を続けることが大事だと思ったことが収穫です。
ところで粘菌の動きをプログラミングすることを考えてみたい。
頭の体操
またまた「人生を書き換える者すらいた。」のエントリから
http://okajima.air-nifty.com/b/2010/01/post-c7b3.html
6人でミーティングをする。 どの2人を取っても、初対面か会ったことがあるかのいずれかである。 このとき、「互いに初対面の3人組」か「互いに会ったことのある3人組」の どちらかは存在することを証明せよ
ある人を基準にして知り合いか初対面かを振り分け、それぞれでまた振り分けていく。
つまりニ分木を構成していく。
どこかに左左 or 左右左が出れば「互いに会ったことのある3人組」が抽出される。
どこかに右右 or 右左右が出れば「互いに初対面の3人組」が抽出される。
それで、5人ならばルートノード+左右+右左までは埋められる。これは上記2パターンを満たさない。
しかし、6人目は左左、(左)右右、左右左、右右、(右)左左、右左右のいずれかしか取りえない。
したがって(1+2+2*(2-1)+1で)6人いると「互いに会ったことのある3人組」か「互いに初対面の3人組」のいずれかを抽出できることとなる。
これだと「どの2人を取っても、初対面か会ったことがあるかのいずれかである。」というのがヒントだと理解できる。
「会ったかどうか(酔っ払ってて)覚えてない」を入れると子ノードの数が3になる…「互いに会ったかどうか覚えてない3人組」の存在も考慮すると(1+3+3*(3-1)+3*(3-1)*(3-2)+1で)17人ミーティングで成り立つかな?
区分がnあると必要な人数が 3(n=1), 6(n=2), 17(n=3), 66(n=4), 327(n=5), ... と増えていく。
N(n)=(N(n-1)-1)*n+2; N(1)=3 →同じ区分の人が3人。
分かる人は10秒で説明できるというのはそうなのかも知れない。
ってこれがすぐに分からないとか自分がアレ過ぎる……
お給料もらっててごめんなさい。
【再考察:不完全】
よく見たら間違ってたことに気づく。
必要十分→十分
でした。ある一人は複数のグループを持つことができるので。ただ、上記は十分条件を満たすので証明に使うことはできるはず。
思い込みは危険でした、反省。
ホントに正しいのか自信が無くなってきた。ダメっぽいので再考。画像の色は意味ないです。
三角形を作ってれば「互いに会ったことのある3人組」が出来る。
三角形を作ってなければ「互いに会ったことのある3人組」は出来ない。
このとき「互いに初対面の3人組」について、
任意の3人間で互いに初対面の2人組は出来る。
(a)その2人のどちらかが残りの3人に会ったことがある場合、最初の3人のうちの残りの1人は後者の3人の中で初対面の2人組とは初対面であり、「互いに初対面の3人組」となる。
(b)その2人のどちらも残りの3人の中で会ったことがない人がいればそれが「互いに初対面の3人組」である。
(a),(b)は三角形を作れない制約から成立して十分条件である、はず?
10秒どころか48時間掛かってます。
【1/18:以下間違い】
「どの2人を取っても、初対面か会ったことがあるかのいずれかである。」というのは関連が真偽のいずれかで表現され曖昧さが無い、ということと考えると、
例えば全員が初対面の場合ABCDEF
と表現される。また、
AABCDE
と表現される場合、AAを抽出すると「会ったことがある」ということとする。
つまり、6人は何グループに分類されるか、という問題に帰着される。
「互いに初対面の3人組」
ABCCCCなどのように3グループ以上に分類されることが必要十分条件である。
「互いに会ったことのある3人組」
AAABBCなどのように同一グループの人数が3人以上であることが必要十分条件である。
また、2グループ以下に分類されることはこの条件を満たす。
つまり2グループ以下であることが十分条件である。
ということで、2グループ以下、3グループ以上はいずれかの条件を満たし、この2通りの論理和が全体集合となるため、命題が真であることを示すことができる。
Ruby勉強 -- 追記:統計的なアプローチ(のつもり)
最適解?何それ美味しいの?なコード
解を見つけたり見つけなかったりします。
探索を繰り返して良さそうな経路の評価を高くし、悪そうな経路の評価を低くしていきます。
探索は前回までの評価を元に行いますが、乱数によるフラつきをいれています。
評価、乱数パラメータのバランスで挙動が変わります。また局所解にトラップされると最短経路を見つけることができません。
いろんな代償を払いますがより一般的な問題へ適用可能なはず。(うっかりマヌケな解を出すのが人間ぽくていい)
# http://okajima.air-nifty.com/b/2010/01/post-abc6.html class Maze @@directions = [{'x'=>-1,'y'=>0},{'x'=>1,'y'=>0},{'x'=>0,'y'=>-1},{'x'=>0,'y'=>1}] @@trial_count = 1000 @@swinging = 1.0 @@positive = 1.05 @@negative = 0.99 @@base = 100.0 def run load_map @@trial_count.times {|i| route = Array.new solve_map(@s['y'],@s['x'], 0, route) } if select_route(@g['y'],@g['x']) then @map.each{|line| print line.to_s + "\n" } else print "fail!\n" end end def load_map data = DATA.read @map = [] data.each_line {|line| @map.push line.chomp.split(//) } @cmap = [] @map.each{|m| @cmap.push Array.new(m.size, @@base)} @map.each_with_index{|line,i| line.each_with_index{|c,j| @s = {'x'=>j, 'y'=>i} if c == 'S' @g = {'x'=>j, 'y'=>i} if c == 'G' }} end def solve_map(y,x,c,route) node = {'x'=>x,'y'=>y} # check goal if @map[y][x] == 'G' route.each_with_index{|node,i| @map[node['y']][node['x']] = ' ' @cmap[node['y']][node['x']] *= @@positive*(1.0+1.0/(c*c)) } return end # new node if @map[y][x] == ' ' @map[y][x] = '-' route.push node end # check route walkable = false @@directions.each{|d| walkable |= @map[y+d['y']][x+d['x']] == ' '} if walkable then max_c = -10000.0 ny=-1 nx=-1 @@directions.each{|d| rc = (rand(10000)/5000.0-1.0)*@@swinging if @cmap[y+d['y']][x+d['x']]+rc > max_c && (@map[y+d['y']][x+d['x']] == ' ' || @map[y+d['y']][x+d['x']] == 'G') max_c = @cmap[y+d['y']][x+d['x']]+rc ny = y+d['y'] nx = x+d['x'] end } return solve_map(ny,nx,c+1,route) if nx != -1 && ny != -1 else # abort route.each_with_index{|node,i| @map[node['y']][node['x']] = ' ' @cmap[node['y']][node['x']] *= @@negative } return end end def select_route(y,x) @map[y][x] = '$' if @map[y][x] == ' ' @@directions.each{|d| return true if @map[y+d['y']][x+d['x']] == 'S'} max_c = -10000.0; ny=-1 nx=-1 @@directions.each{|d| if (@cmap[y+d['y']][x+d['x']] > max_c && @map[y+d['y']][x+d['x']] == ' ') max_c = @cmap[y+d['y']][x+d['x']] ny = y+d['y'] nx = x+d['x'] end } return select_route(ny,nx) if nx != -1 && ny != -1 return false end end maze = Maze.new maze.run __END__ ************************** *S* * * * * * * ************* * * * * ************ * * * * ************** *********** * * ** *********************** * * G * * * *********** * * * * ******* * * * * * **************************
Ruby勉強 -- 追記:C言語の場合
C言語だとダメというわけでもないような。
A*が書きづらいというだけかな。しかし配列確保してポインタ付け替えてくだけでいいからむしろC言語の方が簡単なんじゃないかな。
言語の違いは決定的なものではないように思う。Javaでは書きたくないけど。
/* http://okajima.air-nifty.com/b/2010/01/post-abc6.html */ #include <stdio.h> #define MAX_SIZE 100 #define _(X,Y,S) (((Y)*S)+X) void solve_map(int x, int y, int d, char map[], int cmap[], int xs) { if (cmap[_(x,y,xs)]>d && map[_(x,y,xs)]==' ') cmap[_(x,y,xs)]=d; if (cmap[_(x-1,y,xs)]>d+1 && map[_(x-1,y,xs)]==' ') solve_map(x-1,y,d+1,map,cmap,xs); if (cmap[_(x+1,y,xs)]>d+1 && map[_(x+1,y,xs)]==' ') solve_map(x+1,y,d+1,map,cmap,xs); if (cmap[_(x,y-1,xs)]>d+1 && map[_(x,y-1,xs)]==' ') solve_map(x,y-1,d+1,map,cmap,xs); if (cmap[_(x,y+1,xs)]>d+1 && map[_(x,y+1,xs)]==' ') solve_map(x,y+1,d+1,map,cmap,xs); } int select_route(int x, int y, char map[], int cmap[], int xs) { int min_d; if (map[_(x,y,xs)]==' ') map[_(x,y,xs)]='$'; if (map[_(x-1,y,xs)]=='S') return 1; if (map[_(x+1,y,xs)]=='S') return 1; if (map[_(x,y-1,xs)]=='S') return 1; if (map[_(x,y+1,xs)]=='S') return 1; min_d=cmap[_(x,y,xs)]; if (cmap[_(x-1,y,xs)]<min_d && map[_(x-1,y,xs)]==' ') min_d=cmap[_(x-1,y,xs)]; if (cmap[_(x+1,y,xs)]<min_d && map[_(x+1,y,xs)]==' ') min_d=cmap[_(x+1,y,xs)]; if (cmap[_(x,y-1,xs)]<min_d && map[_(x,y-1,xs)]==' ') min_d=cmap[_(x,y-1,xs)]; if (cmap[_(x,y+1,xs)]<min_d && map[_(x,y+1,xs)]==' ') min_d=cmap[_(x,y+1,xs)]; if (cmap[_(x-1,y,xs)]==min_d && map[_(x-1,y,xs)]==' ') return select_route(x-1,y,map,cmap,xs); if (cmap[_(x+1,y,xs)]==min_d && map[_(x+1,y,xs)]==' ') return select_route(x+1,y,map,cmap,xs); if (cmap[_(x,y-1,xs)]==min_d && map[_(x,y-1,xs)]==' ') return select_route(x,y-1,map,cmap,xs); if (cmap[_(x,y+1,xs)]==min_d && map[_(x,y+1,xs)]==' ') return select_route(x,y+1,map,cmap,xs); return 0; } int main(int argc, char *argv[]) { char map[MAX_SIZE*MAX_SIZE]; int cmap[MAX_SIZE*MAX_SIZE]; int xsize, ysize; int sx, sy; /* start */ int gx, gy; /* goal */ FILE *f; char c; int i, x, y, pos; /* load data */ if ((f=fopen(argv[1],"r"))= NULL) exit(1); xsize = 0; for (i=0,pos=0; (c=fgetc(f))!=EOF; i++) { if (c!='\n') { map[pos++]=c; switch (c) { case 'S': sx=(pos-1)%xsize; sy=(pos-1)/xsize; break; case 'G': gx=(pos-1)%xsize; gy=(pos-1)/xsize; break; } } else if (xsize==0) { xsize=i; } } ysize = i/xsize; fclose(f); for (i=0; i<xsize*ysize; i++) cmap[i] = 8192; /* solve */ solve_map(sx, sy, 0, map, cmap, xsize); if (select_route(gx, gy, map, cmap, xsize) == 1) { /* print data */ for (y=0; y<ysize; y++) { for (x=0; x<xsize; x++) putchar(map[_(x,y,xsize)]); putchar('\n'); } } else { printf("fail!\n"); exit(1); } return 0; }