In [ ]:
from sympy import *
from IPython.display import display
init_printing(use_latex='mathjax')
from data import * #data ディレクトリ(パッケージ)からSD.py(供給/需要マトリックス)をよぶ

X = SD2.X4 # X をインポートしたデータに設定
Y = SD2.Y4

print("供給行列")
display(X)
print("需要行列")
display(Y)

A=Y-X # Xp=Yp を解く

p1, p2, p3, p4, p5 = symbols('p1, p2, p3, p4, p5') #5人5財まで想定
ps = (0,0,0,
      Matrix([[p1],[p2],[p3]]),
      Matrix([[p1],[p2],[p3],[p4]]),
      Matrix([[p1],[p2],[p3],[p4],[p5]])
     )
p=ps[A.shape[1]]#価格ベクトルをAのコラム数にあわせて ps タルプから選択
eqivprice = solve(A*p,p)#Ap=0を,pについて解く。戻り値はpiをキーにもつdictionary.

#switch case 文がないpython
print("各自の収入=支出となるような価格ベクトル =")
if A.shape[1] == 3:
    ep = [eqivprice[p1]/p3,eqivprice[p2]/p3,p3/p3]
    display(ep)
    print("交換比率は、価格の逆数")
    display((1,ep[0]/ep[1],ep[0]/ep[2]))
elif A.shape[1] == 4:
    ep = [eqivprice[p1]/p4,eqivprice[p2]/p4,eqivprice[p3]/p4,p4/p4]
    display(ep)
    print("交換比率は、価格の逆数")
    display((1,ep[0]/ep[1],ep[0]/ep[2],ep[0]/ep[3]))
elif A.shape[1] == 5:
    ep = [eqivprice[p1]/p5,eqivprice[p2]/p5,eqivprice[p3]/p5,eqivprice[p4]/p5,p5/p5]
    display(ep)
    print("交換比率は、価格の逆数")
    display((1,ep[0]/ep[1],ep[0]/ep[2],ep[0]/ep[3],ep[0]/ep[4]))
else:
    print("Matrixのサイズは 3x3,4X4,5X5 のみ")


def exchange(w1,i,j,l,m):# (i,j)が、(l,m)から、m財をw1だけ取得
    man = ['A','B','C','D','E']
    ware = ['リンネル','上衣','茶','砂糖','鉄']
    xrate = ep[m]/ep[j]
    w2 = w1*xrate #(i,j) が交換で相手に与えるj財の量
    X[i,m] += w1 # i が需要する
    X[l,m] += -w1 #lが供給する
    X[i,j] += -w2 #iが供給
    X[l,j] += w2 #lが需要 
    print("{}{}:{}単位が、{}{}:{}単位と、1:{}の比率で交換される".format(man[i],ware[j],w2,man[l],ware[m],w1,xrate))
    display(X)

#Demand 型の交換
# (i,j)が、(l,m)から、m財をw1だけ取得して、j財w2を与える
def exchangeD(w1,i,j,l,m):
    xrate = ep[m]/ep[j]
    w2 = w1*xrate #(i,j) が交換で相手に与えるj財の量
    X[i,m] += w1 # i が需要する
    X[l,m] += -w1 #lが供給する
    X[i,j] += -w2 #iが供給
    X[l,j] += w2 #lが需要 
    man = ['A','B','C','D','E']
    ware = ['リンネル','上衣','茶','砂糖','鉄']
    print("{}{}:{}単位が、{}{}:{}単位と、1:{}の比率で交換される".format(man[i],ware[j],w2,man[l],ware[m],w1,xrate))
    display(X)
    
#Suply 型の交換
# (i,j)が、(l,m)に、j財をw1だけ与えてm財 w2を取得する。
def exchangeS(w1,i,j,l,m):
    xrate = ep[j]/ep[m]
    w2 = w1*xrate #(i,j) が交換で相手から得られるm財の量
    X[i,j] += -w1 # i が供給
    X[l,j] += w1 #lが需要
    X[i,m] += w2 #iが需要
    X[l,m] += -w2 #lが供給
    man = ['A','B','C','D','E']
    ware = ['リンネル','上衣','茶','砂糖','鉄']
    print("{}{}:{}単位が、{}{}:{}単位と、1:{}の比率で交換される".format(man[i],ware[j],w1,man[l],ware[m],w2,xrate))
    display(X)

print("\n<<交換開始>>")

'''
交換のプロセスをプログラミング
    A,B,C,... と順に自分の需要をみたす
        - B->A という交換はなし B->C,D,...のみ
        - 自分の余剰の状態をチェックして、その財kを記憶する。
            --TODO:初期値では対角要素だが、交換のプロセスで別のように変わる可能性があるので注意
            --TODO:変わっても余剰財は1種類だと仮定してコーディングしてある。ここ、バグの可能性?
            --TODO:対策 余剰財のリストをつくっておいて、A[i,k1]で足りなければA[i,k2].... 
            --TODO:とすればよいが面倒
            --Done:supply配列に納めて処理
        - 自分の不足分A[i,m]をA[0,m] から下に向かって探してゆき
        - 余剰をもっている人 l  が見つかったら exchange 関数をよぶ。
        - exchange(A[i,m],i,k, l,m) :不足分A[i,m]を、(i,k)財を渡して、(l,m)財と得る
'''
(r,c) = Y.shape
for i in range(r): # A,B,C,.. の順に需要を満たしてゆく(「最後は予定調和的に満たされる」=>r-1と思ったのは誤り)
    A = Y - X # 過不足をしめす行列をreset
    #供給量のチェック
    supply=[] #供給可能な財番号の配列
    for j in range(c):
        if A[i,j] < 0:# 余剰がある財(初期値では対角要素だが、ズレる可能性があるのでチェック。注意!)
            supply.append(j)
    for m in range(c): #リンネル、上衣、茶...... の順に みていって         
        if A[i,m] > 0:# もしm 財の交換が必要なら
            for l in range(r):# m財をもっている人たちの保有状態を順にあたってゆき
                if A[l,m] < 0 : #マイナスなら余剰をもっているということになるので、
                    #まずsupply[0]で交換、つぎにsupply[1]で...
                    #以下のアルゴリズムは不完全:バグあり注意
                    for s in supply:
                        w = -A[i,s]*ep[s]/ep[m] # sを与えて得られるm財の量
                        if w < A[i,m]:
                            exchangeS(-A[i,s],i,s,l,m)
                            A= Y-X # 交換の結果をAに反映させる
                        else:
                            exchangeD(A[i,m],i,s,l,m) #自分の余剰財(i,k) を、この人の余剰財 (l,m) と交換して、A[i,m]を手に入れる。
                            A = Y-X # 交換の結果をAに反映させる
print("<<完了>>")