Project1

标题: 拿出破轮子:线性规划的单纯形法 [打印本页]
作者: 浮云半仙    时间: 2016-12-15 20:42
标题: 拿出破轮子:线性规划的单纯形法
线性规划能做很多事啊。
听说叶子姐姐会感兴趣然后就拿出了这个破轮子...
具体用法写在注释里面啦....
Inf这次我可用Float::INFINITY表示啦
代码写得丑，大概凑合看吧....
自觉回沟

simplex.zip (1.49 KB, 下载次数: 83)

RUBY 代码复制
  1. =begin
  2.     线性规划的单纯形法
  3.     较为详细的资料见 [url]https://github.com/sxysxy/math/tree/master/simplex[/url]
  4.     
  5.     用法:
  6.     lp = LinearProgramming.new(n, [a1, a2, .. an])
  7.         创建一个线性规划对象，它具有n个变量，目标函数为z = a1*x1 + a2*x2 + ... + an*xn
  8.     lp.addAssert([a1, a2, ..an], b)
  9.         添加一条约束，为 a1*x1 + a2*x2 + ... + an*xn <= b
  10.     默认具有约束 x1 >= 0, x2 >= 0, ..., xn >= 0 (非负约束)
  11.         对于不满足上述形式的约束请使用适当的方法转变为上述形式
  12.     
  13.     lp.solve 求解
  14.         如果无解返回          :no_solution
  15.         如果目标函数发散返回  :infinity
  16.         否则返回              [z, [x1, x2, x3, ..., xn]]，目标函数的最大值和此时变量的取值
  17.     例子:
  18.     lp = LinearProgramming.new(3, [3, 1, 2]) #最大化3x1 + x2 + 2x3
  19.     lp.addAssert [1, 1, 3], 30  #x1 + x2 + 3x3 <= 30
  20.     lp.addAssert [2, 2, 5], 24  #2x1 + 2x2 + 5x3 <= 24
  21.     lp.addAssert [4, 1, 2], 36  #4x1 + x2 + 2x3 <= 36
  22.     print lp.solve  #[28.0, [8.0, 4.0, 0.0]] 目标函数最大值为28，此时x1 = 8, x2 = 4, x3 = 0
  23. =end
  24. class LinearProgramming
  25.     NO_SOLUTION = :no_solution
  26.     INFINITY = :infinity
  27.  
  28.     attr_reader :xs #变量数
  29.     attr_reader :ys #约束条件数
  30.  
  31.     def initialize(xs, f)  #f为目标函数
  32.         @xs = xs
  33.         @ys = 0
  34.         @a = []
  35.         @a << [0.0]+f.map {|e| Float(e)}
  36.         @base = []
  37.         @unbase = []
  38.     end
  39.  
  40.     def addAssert(x, b)
  41.         @a << ([b]+x).map{|e| Float(e)}
  42.         @ys += 1
  43.     end
  44.  
  45.     def solve
  46.         r = simplex
  47.         return NO_SOLUTION if r == 0
  48.         return INFINITY if r == -1
  49.         @unbase.fill 0
  50.         (1..@ys).each {|i| @unbase[@base[i]] = i if @base[i] <= @xs}
  51.         return [@a[0][0], (1..@xs).map{|i| @unbase[i]!=0? @a[@unbase[i]][0]: 0.0}] 
  52.     end
  53.     private
  54.     EPS = 1e-9
  55.     def dcmp(x)
  56.         return -1 if x<-EPS
  57.         return 1 if x>EPS
  58.         return 0
  59.     end
  60.     def swapVar(x, y)
  61.         t = @base[x]
  62.         @base[x] = @unbase[y]
  63.         @unbase[y] = t
  64.     end
  65.     def pivot(x, y)
  66.         swapVar(x, y)
  67.         k = @a[x][y]
  68.         @a[x][y] = 1.0
  69.         @a[x].map! {|e| e/k}
  70.         (0..@ys).each do |i|
  71.             next if x == i || dcmp(@a[i][y]) == 0
  72.             k = @a[i][y]
  73.             @a[i][0] += (i==0? 1:-1)*k*@a[x][0]
  74.             @a[i][y] = 0
  75.             (1..@xs).each do |j|
  76.                 @a[i][j] -= k*@a[x][j]
  77.             end
  78.         end
  79.     end
  80.     def init
  81.         (1..@xs).each{|e| @unbase[e] = e}
  82.         (1..@ys).each{|e| @base[e] = e+@xs}
  83.         x = y = 0
  84.         while true
  85.             (1..@ys).each do |i|
  86.                 x = i if dcmp(@a[i][0]) < 0
  87.             end
  88.             return true if x == 0
  89.             (1..@xs).each do |i|
  90.                 y = i if dcmp(@a[x][i]) < 0
  91.             end
  92.             return false if y == 0
  93.             pivot(x, y)
  94.             x = y = 0
  95.         end
  96.     end
  97.     def simplex
  98.         return 0 unless init
  99.         x = y = 0
  100.         while true
  101.             (1..@xs).each do |i|
  102.                 if dcmp(@a[0][i]) > 0
  103.                     y = i
  104.                     break
  105.                 end
  106.             end
  107.             return 1 if y == 0
  108.             inf = Float::INFINITY
  109.             (1..@ys).each do |i|
  110.                 t = @a[i][0]/@a[i][y]
  111.                 if dcmp(@a[i][y])>0 && (x==0 || t<inf)
  112.                     inf = t
  113.                     x = i
  114.                 end
  115.             end
  116.             return -1 if x == 0
  117.             pivot(x, y)
  118.             x = y = 0
  119.         end
  120.     end
  121. end
  122.  
  123. lp = LinearProgramming.new(3, [3, 1, 2]) #最大化3x1 + x2 + 2x3
  124. lp.addAssert [1, 1, 3], 30  #x1 + x2 + 3x3 <= 30
  125. lp.addAssert [2, 2, 5], 24  #2x1 + 2x2 + 5x3 <= 24
  126. lp.addAssert [4, 1, 2], 36  #4x1 + x2 + 2x3 <= 36
  127. print lp.solve

作者: shitake    时间: 2016-12-15 21:25
嘛 代码写好是慢慢练习的事。ACM的话，因为通常是写“一次性”代码【要么过，要么不过换下一种】。所以对很多地方要求并不高【也没那没多时间去在意】。
但是实际开发恰恰相反，强调代码的可维护性。
ruby的话，可以用用RubyCritic或者rubocop做代码检查。坚持做一段时间的话，基本就没问题啦
还可以参考参考这个：https://github.com/bbatsov/ruby-style-guide
作者: 晴兰    时间: 2016-12-16 08:38
提示: 作者被禁止或删除 内容自动屏蔽



欢迎光临 Project1 (https://rpg.blue/) Powered by Discuz! X3.1