拿出破轮子:线性规划的单纯形法

浮云半仙

线性规划能做很多事啊。
听说叶子姐姐会感兴趣然后就拿出了这个破轮子...
具体用法写在注释里面啦....
Inf这次我可用Float::INFINITY表示啦
代码写得丑，大概凑合看吧....
自觉回沟

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

2016-12-15 20:40 上传

点击文件名下载附件

[pre lang="ruby"]=begin
    线性规划的单纯形法
    较为详细的资料见 https://github.com/sxysxy/math/tree/master/simplex
   
    用法:
    lp = LinearProgramming.new(n, [a1, a2, .. an])
        创建一个线性规划对象，它具有n个变量，目标函数为z = a1*x1 + a2*x2 + ... + an*xn
    lp.addAssert([a1, a2, ..an], b)
        添加一条约束，为 a1*x1 + a2*x2 + ... + an*xn <= b
    默认具有约束 x1 >= 0, x2 >= 0, ..., xn >= 0 (非负约束)
        对于不满足上述形式的约束请使用适当的方法转变为上述形式
   
    lp.solve 求解
        如果无解返回          :no_solution
        如果目标函数发散返回  :infinity
        否则返回              [z, [x1, x2, x3, ..., xn]]，目标函数的最大值和此时变量的取值
    例子:
    lp = LinearProgramming.new(3, [3, 1, 2]) #最大化3x1 + x2 + 2x3
    lp.addAssert [1, 1, 3], 30  #x1 + x2 + 3x3 <= 30
    lp.addAssert [2, 2, 5], 24  #2x1 + 2x2 + 5x3 <= 24
    lp.addAssert [4, 1, 2], 36  #4x1 + x2 + 2x3 <= 36
    print lp.solve  #[28.0, [8.0, 4.0, 0.0]] 目标函数最大值为28，此时x1 = 8, x2 = 4, x3 = 0
=end
class LinearProgramming
    NO_SOLUTION = :no_solution
    INFINITY = :infinity
   
    attr_reader :xs #变量数
    attr_reader :ys #约束条件数
   
    def initialize(xs, f)  #f为目标函数
        @xs = xs
        @ys = 0
        @a = []
        @a << [0.0]+f.map {|e| Float(e)}
        @base = []
        @unbase = []
    end
   
    def addAssert(x, b)
        @a << (+x).map{|e| Float(e)}
        @ys += 1
    end
   
    def solve
        r = simplex
        return NO_SOLUTION if r == 0
        return INFINITY if r == -1
        @unbase.fill 0
        (1..@ys).each {|i| @unbase[@base] = i if @base <= @xs}
        return [@a[0][0], (1..@xs).map{|i| @unbase!=0? @a[@unbase][0]: 0.0}]
    end
    private
    EPS = 1e-9
    def dcmp(x)
        return -1 if x<-EPS
        return 1 if x>EPS
        return 0
    end
    def swapVar(x, y)
        t = @base[x]
        @base[x] = @unbase[y]
        @unbase[y] = t
    end
    def pivot(x, y)
        swapVar(x, y)
        k = @a[x][y]
        @a[x][y] = 1.0
        @a[x].map! {|e| e/k}
        (0..@ys).each do |i|
            next if x == i || dcmp(@a[y]) == 0
            k = @a[y]
            @a[0] += (i==0? 1:-1)*k*@a[x][0]
            @a[y] = 0
            (1..@xs).each do |j|
                @a[j] -= k*@a[x][j]
            end
        end
    end
    def init
        (1..@xs).each{|e| @unbase[e] = e}
        (1..@ys).each{|e| @base[e] = e+@xs}
        x = y = 0
        while true
            (1..@ys).each do |i|
                x = i if dcmp(@a[0]) < 0
            end
            return true if x == 0
            (1..@xs).each do |i|
                y = i if dcmp(@a[x]) < 0
            end
            return false if y == 0
            pivot(x, y)
            x = y = 0
        end
    end
    def simplex
        return 0 unless init
        x = y = 0
        while true
            (1..@xs).each do |i|
                if dcmp(@a[0]) > 0
                    y = i
                    break
                end
            end
            return 1 if y == 0
            inf = Float::INFINITY
            (1..@ys).each do |i|
                t = @a[0]/@a[y]
                if dcmp(@a[y])>0 && (x==0 || t<inf)
                    inf = t
                    x = i
                end
            end
            return -1 if x == 0
            pivot(x, y)
            x = y = 0
        end
    end
end

lp = LinearProgramming.new(3, [3, 1, 2]) #最大化3x1 + x2 + 2x3
lp.addAssert [1, 1, 3], 30  #x1 + x2 + 3x3 <= 30
lp.addAssert [2, 2, 5], 24  #2x1 + 2x2 + 5x3 <= 24
lp.addAssert [4, 1, 2], 36  #4x1 + x2 + 2x3 <= 36
print lp.solve
[/pre]

shitake

嘛 代码写好是慢慢练习的事。ACM的话，因为通常是写“一次性”代码【要么过，要么不过换下一种】。所以对很多地方要求并不高【也没那没多时间去在意】。
但是实际开发恰恰相反，强调代码的可维护性。
ruby的话，可以用用RubyCritic或者rubocop做代码检查。坚持做一段时间的话，基本就没问题啦
还可以参考参考这个：https://github.com/bbatsov/ruby-style-guide