fork download
copy 
  1. #!/use/bin/env python
  2. # -*- coding: UTF-8 -*-
  3.  
  4. from numpy import *
  5. from scipy import *
  6. from scipy import linalg
  7. import string
  8. import sys
  9.  
  10. def main():
  11.     args = sys.argv[1:]
  12.     problem = 1
  13.     method = 1
  14.     search = 1
  15.  
  16.     if problem == 1:
  17.         n = 2
  18.         x = array([-36.,114.])
  19.  
  20.         def evalf(x):
  21.             f = x[0]**2 + math.exp(x[0]) + x[1]**4 + x[1]**2 - 2*x[0]*x[1] + 3
  22.             return f
  23.  
  24.         def evalg(x):
  25.             g = array ([2*x[0] +  math.exp(x[0]) - 2*x[1],4*x[1]**3 + 2*x[1] - 2*x[0]])
  26.             return g
  27.  
  28.         def evalh(x):
  29.             H = array([[2 + math.exp(x[0]),-2],[-2,12*x[1]**2 + 2]])
  30.             return H
  31.  
  32.     elif problem == 2:
  33.         n = 2
  34.         x = array([1.,1.])
  35.  
  36.         def evalf(x):
  37.             f = 0.0
  38.             for i in range(0,3):
  39.                 f += evalf_i(i,x)**2
  40.             return f
  41.  
  42.         def evalf_i(i,x):
  43.             y = array([1.5,2.25,2.625])
  44.             f = y[i] - x[0] * (1 - x[1]**(i+1))
  45.             return f
  46.  
  47.         def evalg_i(i,x):
  48.             g = array([x[1]**(i+1) - 1,(i+1) * x[0] * x[1]**i])
  49.             return g
  50.  
  51.         def evalh_i(i,x):
  52.             h = array([[0,(i+1) * x[1]**i],[(i+1) * x[1]**i,(i+1) * i * x[0] * x[1]**(i-1)]])
  53.             return h
  54.  
  55.         def evalg(x):
  56.             g = zeros(2)
  57.             for i in range(0,3):
  58.                 g += 2 * evalf_i(i,x) * evalg_i(i,x)
  59.             return g
  60.  
  61.         def evalh(x):
  62.             H = array([[0.0,0.0],[0.0,0.0]])
  63.             for i in range(0,3):
  64.                 H += 2 * (evalf_i(i,x) * evalh_i(i,x) + dot(evalg_i(i,x).reshape(-1,1),evalg_i(i,x).reshape(1,-1)))
  65.             return H
  66.  
  67.     elif problem == 3:
  68.         n = 20
  69.         Z = random.rand(n,n)
  70.         A = dot(Z.T,Z)
  71.         x = ones(n)
  72.  
  73.         def evalf(x):
  74.             f = inner(dot(A,x),x)
  75.             return f
  76.  
  77.         def evalg(x):
  78.             g = dot(A,x)
  79.             return g
  80.  
  81.         def evalh(x):
  82.             H = A
  83.             return H
  84.  
  85.     elif problem == 4:
  86.         n = 2
  87.         x = array([0.,0.])
  88.         l = array([-5,0])
  89.         u = array([10,15])
  90.  
  91.         def evalf(x):
  92.             f = (x[1] - 5.1*x[0]**2 / (4*math.pi**2) + 5*x[0] / math.pi -6)**2 + 10*(1 - 1/(8*math.pi))*math.cos(x[0]) + 10
  93.             return f
  94.  
  95.         def evalg(x):
  96.             g = array([2*(x[1] - 5.1*x[0]**2 / (4*math.pi**2) + 5*x[0] / math.pi -6)*(-5.1*x[0] / (2*math.pi**2) + 5/math.pi) - 10*(1 - 1/ (8*math.pi))*math.sin(x[0]),2*(x[1] - 5.1*x[0]**2 / (4*math.pi**2) + 5*x[0] / math.pi -6)])
  97.             return g
  98.  
  99.         def evalh(x):
  100.             h = array([[2*(-5.1*x[0] / (2*math.pi**2) + 5/math.pi)**2 -5.1*(x[1] - 5.1*x[0]**2 / (4*math.pi**2) + 5*x[0] / math.pi -6) / math.pi**2,2*(-5.1*x[0] / (2*math.pi**2) + 5/math.pi)],[2*(-5.1*x[0] / (2*math.pi**2) + 5/math.pi),2]])
  101.             return h
  102.  
  103.     elif problem == 5:
  104.         n = 2
  105.         x = array([-1.,-22.])
  106.         l = -32.768 * ones(n)
  107.         u = 32.768 * ones(n)
  108.  
  109.         def evalf1(x):
  110.             f1 = 0
  111.             for i in range(0,n):
  112.                 f1 += x[i]**2
  113.             f1 = -0.2 * sqrt(f1/n)
  114.             return f1
  115.  
  116.         def evalg1(x):
  117.             g1 = zeros(n)
  118.             for i in range(0,n):
  119.                 g1[i] = 0.04 * x[i] / (n * evalf1(x))
  120.             return g1
  121.  
  122.         def evalf2(x):
  123.             f2 = 0
  124.             for i in range(0,n):
  125.                 f2 += math.cos(2 * math.pi * x[i])
  126.             f2 = f2/n
  127.             return f2
  128.  
  129.         def evalg2(x):
  130.             g2 = zeros(n)
  131.             for i in range(0,n):
  132.                 g2[i] = -2 * math.pi * math.sin(2 * math.pi * x[i]) / n
  133.             return g2 
  134.  
  135.         def evalf(x):
  136.             f = -20 * math.exp(evalf1(x)) - math.exp(evalf2(x)) + 20 + e
  137.             return f
  138.  
  139.         def evalg(x):
  140.             g = -20 * math.exp(evalf1(x)) * evalg1(x) - math.exp(evalf2(x)) * evalg2(x)
  141.             return g
  142.  
  143.     else:
  144.         print 'Error!'
  145.         return 0
  146.  
  147.     xi = 1e-4
  148.     rho = 0.5
  149.     epsilon = n * 1e-6
  150.     k = 0
  151.     f_num = 0
  152.  
  153.     while True:
  154.         if k >= 200000:
  155.             print 'k =',k,'\n許容最大反復回数を超えました.'
  156.             return 0
  157.  
  158.  
  159.         if method == 1:
  160.             g = evalg(x)
  161.             if linalg.norm(g) <= epsilon:
  162.                 break
  163.             d = - g
  164.  
  165.         elif method == 2:
  166.             g = evalg(x)
  167.             if linalg.norm(g) <= epsilon:
  168.                 break
  169.             tau = 0
  170.             B = evalh(x)
  171.             I = identity(n)
  172.             while True:
  173.                 try:
  174.                     L = linalg.cho_factor(B + tau * I)
  175.                     d = linalg.cho_solve(L,-g)
  176.                     break
  177.                 except:
  178.                     if tau == 0:
  179.                         tau = 2
  180.                     else:
  181.                         tau = tau * 2
  182.  
  183.         elif method == 3:
  184.             s = 1
  185.             d = fmax(l,fmin(x - s * evalg(x),u)) - x
  186.             if linalg.norm(d) <= epsilon:
  187.                 break
  188.  
  189.         else:
  190.             print 'Error!'
  191.             return 0
  192.  
  193.         t = 1
  194.         f = evalf(x)
  195.         f_num = f_num + 1
  196.         dg = inner(d,evalg(x))
  197.         while True:
  198.             f_step = evalf(x + t * d)
  199.             f_num = f_num + 1
  200.             if f_step <= f + xi * t * dg:
  201.                 break
  202.             if search == 1:
  203.                 t = rho * t
  204.             elif search == 2:
  205.                 t_prev = t
  206.                 t = - dg * t**2 / (2 * (f_step - f - dg * t))
  207.                 if t < 0.1 * t_prev or t > 0.9 * t_prev:
  208.                     t = t_prev / 2.0
  209.             else:
  210.                 print 'Error!'
  211.                 return 0
  212.  
  213.         if method == 3:
  214.             if t * max(abs(d)) <= 1e-16 * max(1,max(abs(x))):
  215.                 break
  216.  
  217.         x = x + t * d
  218.         k = k + 1
  219.         print 'g =',linalg.norm(evalg(x)), '\nx =',x
  220.     print 'f =',evalf(x),'\nx =',x,'\nk =',k,'\n評価回数 =',f_num
  221. if __name__ == '__main__': main()
Success #stdin #stdout 0.07s 26800KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
g = 1195678.86248 
x = [-35.99084473 -66.86169434]
g = 1003.18899876 
x = [-35.99461314   6.11675176]
g = 82.5821342539 
x = [-35.33662306  -1.69299267]
g = 91.9919743385 
x = [-33.23389616  -3.1891575 ]
g = 91.2432314816 
x = [-29.47830383   1.16426162]
g = 82.8941483651 
x = [-25.64798315  -3.06059971]
g = 74.4062517155 
x = [-22.82456022   1.28334375]
g = 36.7818080337 
x = [-19.81107222  -2.25855179]
g = 146.292838491 
x = [-2.25855179  3.23092697]
g = 8.25356886148 
x = [-1.91872505 -1.32807649]
g = 4.00210506822 
x = [-1.78941216 -0.3045165 ]
g = 4.09367201056 
x = [-1.08872892 -1.01872642]
g = 1.11947536454 
x = [-1.11330881 -0.50760811]
g = 0.830207348466 
x = [-0.89257604 -0.67966511]
g = 0.589188384651 
x = [-0.88852039 -0.4721529 ]
g = 0.347020989267 
x = [-0.7831526  -0.57508037]
g = 0.210909095536 
x = [-0.79335726 -0.48892738]
g = 0.116932966172 
x = [-0.75422324 -0.52426424]
g = 0.0661016468296 
x = [-0.7568377  -0.49514814]
g = 0.036595597756 
x = [-0.74327984 -0.50459662]
g = 0.0206770950296 
x = [-0.74282613 -0.49545898]
g = 0.0118688930714 
x = [-0.73808441 -0.49751749]
g = 0.0174013163254 
x = [-0.73653186 -0.49178973]
g = 0.00896249341193 
x = [-0.73385366 -0.49521794]
g = 0.00465067219438 
x = [-0.73454966 -0.49308815]
g = 0.00242481753719 
x = [-0.73374927 -0.49393146]
g = 0.0012795801447 
x = [-0.73386676 -0.49333675]
g = 0.0006862586574 
x = [-0.73361405 -0.49353289]
g = 0.000376723047772 
x = [-0.73361609 -0.49336134]
g = 0.000212856227995 
x = [-0.73353109 -0.4934019 ]
g = 0.000321318850996 
x = [-0.73350706 -0.49329822]
g = 0.000165256999167 
x = [-0.7334581  -0.49336191]
g = 8.53728931285e-05 
x = [-0.73347135 -0.49332277]
g = 4.44130480264e-05 
x = [-0.73345681 -0.4933384 ]
g = 2.3357795694e-05 
x = [-0.73345911 -0.49332754]
g = 1.24837584106e-05 
x = [-0.73345455 -0.49333118]
g = 6.82282708861e-06 
x = [-0.73345464 -0.49332806]
g = 3.83612276092e-06 
x = [-0.73345311 -0.49332883]
g = 2.22737270034e-06 
x = [-0.73345291 -0.49332789]
g = 3.02709890153e-06 
x = [-0.73345182 -0.49332812]
g = 1.56245154447e-06 
x = [-0.73345207 -0.49332741]
f = 3.59713802496 
x = [-0.73345207 -0.49332741] 
k = 41 
評価回数 = 206
https://ideone.com/AuO23i
language:
Python (cpython 2.7.16)
created:
7 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!