経路の計算 - プログラム倉庫

Function ponu(s As Single, mx As Single, nx As Single, u) As Single
Dim m1 As Single
Dim m2 As Single
Dim m3 As Single
Dim n1 As Single
Dim n2 As Single
Dim n3 As Single
Dim pp As Single
Dim dm As Single
Dim dn As Single
Dim nu As Single
pp = 0
m1 = mx
n1 = nx
If m1 < -10 Then pp = 1
If m1 < -10 Then m1 = -10
If m1 > 9 Then pp = 1
If m1 > 9 Then m1 = 9
If n1 < -10 Then pp = 1
If n1 < -10 Then n1 = -10
If n1 > 9 Then pp = 1
If n1 > 9 Then n1 = 9
m2 = Int(m1)
m3 = m2 + 1
n2 = Int(n1)
n3 = n2 + 1
dm = (m1 - m2) * (u(s, m3, n2) - u(s, m2, n2))
dn = (n1 - n2) * (u(s, m2, n3) - u(s, m2, n2))
nu = u(s, m2, n2) + dm + dn
If pp = 1 Then nu = -999
ponu = nu
End Function
Function fastv(n As Single, h As Single, ws, fastu, v) As Single
Dim bxs As Single
Dim uxs As Single
Dim wxs As Single
Dim mxs As Single
Dim nxs As Single
Dim t1 As Single
Dim bxp As Single
Dim t2 As Single
uxs = fastu(n)
wxs = seekw(uxs)
mxs = (wxs - ws(99)) / h
nxs = -0.1 * n * h
fastv = uxs + nearv(99, mxs, nxs, v)
End Function

Function nearfast(nx As Single, fastu) As Single
Dim n1 As Single
Dim n2 As Single
Dim n3 As Single
Dim dn As Single
n1 = nx
pp = 0
If n1 > 9 Then pp = -999
If n1 > 9 Then n1 = 9
If n1 < -9 Then pp = -999
If n1 < -9 Then n1 = -9
n2 = Int(n1)
n3 = n2 + 1
dn = (n1 - n2) * (fastu(n3) - fastu(n2))
nearfast = fastu(n2) + dn + pp
End Function

Function gotow(s As Single, m As Single, n As Single, h As Single, ws, u, v) As Single
Dim bxs As Single
Dim bx1 As Single
Dim bx2 As Single
Dim uxs As Single
Dim wxs As Single
Dim mxs As Single
Dim nxs As Single
Dim ux1 As Single
Dim wx1 As Single
Dim mx1 As Single
Dim nx1 As Single
Dim ux2 As Single
Dim wx2 As Single
Dim mx2 As Single
Dim nx2 As Single
Dim j As Single
Dim t1 As Single
Dim bxp As Single
Dim t2 As Single
j = 3
bxs = 0
uxs = ponu(s, m, bxs, u)
wxs = seekw(uxs)
mxs = (wxs - ws(s - 1)) / h
nxs = n - bxs
vxs = uxs + nearv(s - 1, mxs, nxs, v)
bxp = 888
t2 = 0
Do Until t2 > 10
t1 = 0
Do Until t1 > 100
bx1 = bxs + j
ux1 = ponu(s, m, bx1, u)
wx1 = seekw(ux1)
mx1 = (wx1 - ws(s - 1)) / h
nx1 = n - bx1
vx1 = ux1 + nearv(s - 1, mx1, nx1, v)
bx2 = bxs - j
ux2 = ponu(s, m, bx2, u)
wx2 = seekw(ux2)
mx2 = (wx2 - ws(s - 1)) / h
nx2 = n - bx2
vx2 = ux2 + nearv(s - 1, mx2, nx2, v)
If vx1 > vxs Then bxs = bx1
If vx1 > vxs Then vxs = vx1
If vx2 > vxs Then bxs = bx2
If vx2 > vxs Then vxs = vx2
If (bxp - bxs) ^ 2 < 10 ^ (-9) Then t1 = 1000
bxp = bxs
t1 = t1 + 1
Loop
t2 = t2 + 1
j = j / 2
Loop
uxs = ponu(s, m, bxs, u)
wxs = seekw(uxs)
mxs = (wxs - ws(s - 1)) / h
gotow = mxs
End Function
Function gotob(s As Single, m As Single, n As Single, h As Single, ws, u, v) As Single
Dim bxs As Single
Dim bx1 As Single
Dim bx2 As Single
Dim uxs As Single
Dim wxs As Single
Dim mxs As Single
Dim nxs As Single
Dim ux1 As Single
Dim wx1 As Single
Dim mx1 As Single
Dim nx1 As Single
Dim ux2 As Single
Dim wx2 As Single
Dim mx2 As Single
Dim nx2 As Single
Dim j As Single
Dim t1 As Single
Dim bxp As Single
Dim t2 As Single
j = 3
bxs = 0
uxs = ponu(s, m, bxs, u)
wxs = seekw(uxs)
mxs = (wxs - ws(s - 1)) / h
nxs = n - bxs
vxs = uxs + nearv(s - 1, mxs, nxs, v)
bxp = 888
t2 = 0
Do Until t2 > 10
t1 = 0
Do Until t1 > 100
bx1 = bxs + j
ux1 = ponu(s, m, bx1, u)
wx1 = seekw(ux1)
mx1 = (wx1 - ws(s - 1)) / h
nx1 = n - bx1
vx1 = ux1 + nearv(s - 1, mx1, nx1, v)
bx2 = bxs - j
ux2 = ponu(s, m, bx2, u)
wx2 = seekw(ux2)
mx2 = (wx2 - ws(s - 1)) / h
nx2 = n - bx2
vx2 = ux2 + nearv(s - 1, mx2, nx2, v)
If vx1 > vxs Then bxs = bx1
If vx1 > vxs Then vxs = vx1
If vx2 > vxs Then bxs = bx2
If vx2 > vxs Then vxs = vx2
If (bxp - bxs) ^ 2 < 10 ^ (-9) Then t1 = 1000
bxp = bxs
t1 = t1 + 1
Loop
t2 = t2 + 1
j = j / 2
Loop
gotob = n - bxs
End Function
Function seekfastu(th, bp As Single) As Single
Dim ls As Single
Dim cs As Single
Dim us As Single
Dim ys As Single
Dim ws As Single
Dim l1 As Single
Dim c1 As Single
Dim u1 As Single
Dim l2 As Single
Dim y2 As Single
Dim y1 As Single
Dim c2 As Single
Dim u2 As Single
Dim h As Single
Dim lp As Single
Dim t1 As Single
Dim t2 As Single
Dim e As Single
e = 10 ^ (-5)
h = 0.1
ls = 0.5
ys = th(100) * ls
cs = ys - bp
us = Log(cs) + Log(1 - ls)
t2 = 0
Do Until t2 > 10
t1 = 0
Do Until t1 > 100
l1 = ls + h
c1 = th(100) * l1 - bp
y1 = th(100) * l1
u1 = Log(c1) + Log(1 - l1)
l2 = ls - h
c2 = th(100) * l2 - bp
z = 0
u2 = Log(c2) + Log(1 - l2)
If u1 > us Then ls = l1
If u1 > us Then us = u1
If u2 > us Then ls = l2
If u2 > us Then us = u2
If (lp - ls) ^ 2 < e Then t1 = 1000
lp = ls
t1 = t1 + 1
Loop
h = h / 2
t2 = t2 + 1
Loop
seekfastu = us
End Function
Function seekv(s As Single, m As Single, n As Single, h As Single, ws, u, v) As Single
Dim bxs As Single
Dim bx1 As Single
Dim bx2 As Single
Dim uxs As Single
Dim wxs As Single
Dim mxs As Single
Dim nxs As Single
Dim ux1 As Single
Dim wx1 As Single
Dim mx1 As Single
Dim nx1 As Single
Dim ux2 As Single
Dim wx2 As Single
Dim mx2 As Single
Dim nx2 As Single
Dim j As Single
Dim t1 As Single
Dim bxp As Single
Dim t2 As Single
j = 3
bxs = 0
uxs = nearu(s, m, bxs, u)
wxs = seekw(uxs)
mxs = (wxs - ws(s - 1)) / h
nxs = n - bxs
vxs = uxs + nearv(s - 1, mxs, nxs, v)
bxp = 888
t2 = 0
Do Until t2 > 10
t1 = 0
Do Until t1 > 100
bx1 = bxs + j
ux1 = nearu(s, m, bx1, u)
wx1 = seekw(ux1)
mx1 = (wx1 - ws(s - 1)) / h
nx1 = n - bx1
vx1 = ux1 + nearv(s - 1, mx1, nx1, v)
bx2 = bxs - j
ux2 = nearu(s, m, bx2, u)
wx2 = seekw(ux2)
mx2 = (wx2 - ws(s - 1)) / h
nx2 = n - bx2
vx2 = ux2 + nearv(s - 1, mx2, nx2, v)
If vx1 > vxs Then bxs = bx1
If vx1 > vxs Then vxs = vx1
If vx2 > vxs Then bxs = bx2
If vx2 > vxs Then vxs = vx2
If (bxp - bxs) ^ 2 < 10 ^ (-9) Then t1 = 1000
bxp = bxs
t1 = t1 + 1
Loop
t2 = t2 + 1
j = j / 2
Loop
seekv = vxs
End Function


Function nearu(s As Single, mx As Single, nx As Single, u) As Single
Dim m1 As Single
Dim m2 As Single
Dim m3 As Single
Dim n1 As Single
Dim n2 As Single
Dim n3 As Single
Dim pp As Single
Dim dm As Single
Dim dn As Single
Dim nu As Single
pp = 0
m1 = mx
n1 = nx
If n1 < -10 Then pp = 1
If n1 < -10 Then n1 = -10
If n1 > 9 Then pp = 1
If n1 > 9 Then n1 = 9
n2 = Int(n1)
n3 = n2 + 1
dn = (n1 - n2) * (u(s, m2, n3) - u(s, m2, n2))
nu = u(s, m2, n2) + dn
If pp = 1 Then nu = -999
nearu = nu
End Function

Function nearv(s As Single, mx As Single, nx As Single, v) As Single
Dim m1 As Single
Dim m2 As Single
Dim m3 As Single
Dim n1 As Single
Dim n2 As Single
Dim n3 As Single
Dim pp As Single
Dim dm As Single
Dim dn As Single
Dim nv As Single
pp = 0
m1 = mx
n1 = nx
If m1 < -9 Then pp = 1
If m1 < -9 Then m1 = -9
If m1 > 9 Then pp = 1
If m1 > 9 Then m1 = 9
If n1 < -29 Then pp = 1
If n1 < -29 Then n1 = 9
If n1 > 29 Then pp = 1
If n1 > 29 Then n1 = 9
m2 = Int(m1)
m3 = m2 + 1
n2 = Int(n1)
n3 = n2 + 1
dm = (m1 - m2) * (v(s, m3, n2) - v(s, m2, n2))
dn = (n1 - n2) * (v(s, m2, n3) - v(s, m2, n2))
nv = v(s, m2, n2) + dm + dn
If pp = 1 Then nv = -999
nearv = nv
End Function


Function seeku(s As Single, th, wp As Single, bp As Single) As Single
Dim ls As Single
Dim cs As Single
Dim us As Single
Dim ys As Single
Dim ws As Single
Dim l1 As Single
Dim c1 As Single
Dim u1 As Single
Dim l2 As Single
Dim y2 As Single
Dim y1 As Single
Dim c2 As Single
Dim u2 As Single
Dim h As Single
Dim lp As Single
Dim t1 As Single
Dim t2 As Single
Dim e As Single
e = 10 ^ (-5)
h = 0.1
ls = (bp + 0.1) / th(s)
ys = th(s) * ls
cs = ys - bp
us = Log(cs) + Log(1 - ls)
ws = Log(cs) + Log(1 - ys / th(s + 1))
If ws > wp Then us = -999
t2 = 0
Do Until t2 > 10
t1 = 0
Do Until t1 > 100
l1 = ls + h
If l1 > 0.99 Then l1 = ls
c1 = th(s) * l1 - bp
y1 = th(s) * l1
u1 = Log(c1) + Log(1 - l1)
w1 = Log(c1) + Log(1 - y1 / th(s + 1))
If w1 > wp Then u1 = -999
l2 = ls - h
If l2 < 0.01 Then l2 = ls
c2 = th(s) * l2 - bp
If c2 < 0.01 Then l2 = ls
c2 = th(s) * l2 - bp
y2 = th(s) * l2
u2 = Log(c2) + Log(1 - l2)
w2 = Log(c2) + Log(1 - y2 / th(s + 1))
If w2 > wp Then u2 = -999
If u1 > us Then ls = l1
If u1 > us Then us = u1
If u2 > us Then ls = l2
If u2 > us Then us = u2
If (lp - ls) ^ 2 < e Then t1 = 1000
lp = ls
t1 = t1 + 1
Loop
h = h / 2
t2 = t2 + 1
Loop
seeku = us
End Function
Function seekw(wp As Single) As Single
Dim x1 As Single
Dim x2 As Single
Dim x3 As Single
Dim w1 As Single
Dim w2 As Single
Dim t As Single
x1 = 0.3
x2 = 0.7
Do Until t > 100
w1 = 2 * Log(x1)
w2 = 2 * Log(x2)
x3 = x2 + (wp - w2) * (x2 - x1) / (w2 - w1)
x1 = x2
x2 = x3
If (wp - w2) ^ 2 < 10 ^ (-5) Then t = 1000
t = t + 1
Loop
seekw = x2
End Function


Function lx(s As Single, th, tl As Single, tr As Single) As Single
Dim ls As Single
ls = ((1 - tl) * th(s) - tr) / (2 * (1 - tl) * th(s))
If ls < 0 Then ls = 0
lx = ls
End Function
Function cx(s As Single, th, tl As Single, tr As Single) As Single
Dim ls As Single
ls = ((1 - tl) * th(s) - tr) / (2 * (1 - tl) * th(s))
If ls < 0 Then ls = 0
cx = (1 - tl) * th(s) * ls + tr
End Function

Function tls(th) As Single
Dim maxw As Single
Dim tl As Single
Dim tr As Single
Dim tp As Single
Dim w1 As Single
Dim n As Single
maxw = -999
For n = 1 To 400
tl = 0.001 * n
tr = trs(th, tl)
w1 = wel(th, tl, tr)
If w1 > maxw Then tp = tl
If w1 > maxw Then maxw = w1
Next
tls = tp
End Function

Function trs(th, tl As Single) As Single
Dim tr1 As Single
Dim tr2 As Single
Dim tr3 As Single
Dim b1 As Single
Dim b2 As Single
Dim t As Single
tr1 = 0
tr2 = 0.1
t = 0
Do Until t > 100
b1 = bud(th, tl, tr1)
b2 = bud(th, tl, tr2)
tr3 = tr2 - b2 * (tr2 - tr1) / (b2 - b1)
tr1 = tr2
tr2 = tr3
If (tr1 - tr2) ^ 2 < 10 ^ (-5) Then t = 1000
t = t + 1
Loop
trs = tr2
End Function
Function wel(th, tl As Single, tr As Single) As Single
Dim w1 As Single
Dim s As Single
Dim ls As Single
Dim cs As Single
w1 = 0
For s = 1 To 100
ls = lx(s, th, tl, tr)
If ls < 0 Then ls = 0
cs = cx(s, th, tl, tr)
w1 = w1 + Log(cs) + Log(1 - ls)
Next
wel = w1
End Function
Function bud(th, tl As Single, tr As Single) As Single
Dim w1 As Single
Dim s As Single
Dim ls As Single
Dim cs As Single
Dim ys As Single
ys = 0
cs = 0
For s = 1 To 100
ls = lx(s, th, tl, tr)
If ls < 0 Then ls = 0
ys = ys + th(s) * ls
cs = cs + cx(s, th, tl, tr)
Next
bud = ys - cs
End Function
Private Sub Command1_Click()
Dim th(1 To 100) As Single
Dim s As Single
Dim m As Single
Dim n As Single
Dim tl As Single
Dim tr As Single
Dim bs(1 To 100) As Single
Dim ws(1 To 99) As Single
Dim fastu(-100 To 100) As Single
Dim u(1 To 99, -10 To 10, -10 To 10) As Single
Dim v(1 To 99, -10 To 10, -30 To 30) As Single
Dim w1 As Single
Dim wp As Single
Dim bp As Single
Dim h As Single
For s = 1 To 100
th(s) = 0.5 + 0.01 * s
Next
tl = tls(th)
tr = trs(th, tl)
Debug.Print tl, tr
For s = 1 To 100
bs(s) = th(s) * lx(s, th, tl, tr) - cx(s, th, tl, tr)
Next
For s = 1 To 99
w1 = Log(cx(s, th, tl, tr)) + Log(1 - th(s) * lx(s, th, tl, tr) / th(s + 1))
ws(s) = seekw(w1)
Next
h = 10 ^ (-3)
For s = 1 To 99
For m = -10 To 10
For n = -10 To 10
wp = 2 * Log(ws(s) + h * m)
bp = bs(s) + h * n
u(s, m, n) = seeku(s, th, wp, bp)
Next
Next
Next
For s = 1 To 99
For m = -10 To 10
For n = -30 To 30
v(s, m, n) = -999
Next
Next
Next
s = 1
For m = -10 To 10
For n = -10 To 10
v(s, m, n) = u(s, m, n)
Next
Next
For s = 2 To 99
For m = -10 To 10
For n = -10 To 10
v(s, m, n) = seekv(s, m, n, h, ws, u, v)
Next
Next
Next
For n = -100 To 100
fastu(n) = seekfastu(th, bs(100) + n * 0.1 * h)
Next
Dim opw(1 To 99) As Single
Dim opb(1 To 99) As Single
opw(99) = (seekw(fastu(-100)) - ws(99)) / h
opb(99) = 10
For s1 = 1 To 98
s = 100 - s1
opb(s - 1) = gotob(s, opw(s), opb(s), h, ws, u, v)
opw(s - 1) = gotow(s, opw(s), opb(s), h, ws, u, v)
Debug.Print opb(s - 1), opw(s - 1)
Next
End Sub