diff --git a/.travis.yml b/.travis.yml index d250ed9..a582fbd 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,5 +1,10 @@ language: ruby rvm: - - 1.9.3 - - 2.0.0 - - jruby + - 2.2 + - 2.3 + - 2.4 + - ruby + - ruby-head + - jruby-9.1.9 +install: gem install minitest +script: rake diff --git a/Gemfile b/Gemfile deleted file mode 100644 index 75db47d..0000000 --- a/Gemfile +++ /dev/null @@ -1,6 +0,0 @@ -source "https://rubygems.org" -gem "rake" - -group :test do - gem "test-unit" -end diff --git a/Gemfile.lock b/Gemfile.lock deleted file mode 100644 index 9b744b9..0000000 --- a/Gemfile.lock +++ /dev/null @@ -1,12 +0,0 @@ -GEM - remote: https://rubygems.org/ - specs: - rake (10.1.0) - test-unit (2.0.0.0) - -PLATFORMS - ruby - -DEPENDENCIES - rake - test-unit diff --git a/Rakefile b/Rakefile index cb70dcb..9014442 100644 --- a/Rakefile +++ b/Rakefile @@ -1,7 +1,19 @@ require 'rake/testtask' -task :default => :test - Rake::TestTask.new do |t| - t.test_files = FileList['test/**/*_test.rb'] + t.pattern = ['test/*.rb'] + t.warning = true +end + +Rake::TestTask.new(bench: :loadavg) do |t| + t.pattern = ['test/bench/*.rb'] + t.warning = true + t.description = "Run benchmarks" end + +desc "Show current system load" +task :loadavg do + puts "/proc/loadavg %s" % (File.read("/proc/loadavg") rescue "Unavailable") +end + +task :default => :test diff --git a/lib/simplex.rb b/lib/simplex.rb index 98c893f..573b678 100644 --- a/lib/simplex.rb +++ b/lib/simplex.rb @@ -1,192 +1,164 @@ -require 'matrix' - -class Vector - public :[]= -end - class Simplex DEFAULT_MAX_PIVOTS = 10_000 - class UnboundedProblem < StandardError - end + class Error < RuntimeError; end + class UnboundedProblem < Error; end + class SanityCheck < Error; end + class TooManyPivots < Error; end attr_accessor :max_pivots + # c - coefficients of objective function; size: num_vars + # a - inequality lhs coefficients; 2dim size: num_inequalities, num_vars + # b - inequality rhs constants size: num_inequalities def initialize(c, a, b) - @pivot_count = 0 + num_vars = c.size + num_inequalities = b.size + raise(ArgumentError, "a doesn't match b") unless a.size == num_inequalities + raise(ArgumentError, "a doesn't match c") unless a.first.size == num_vars + @max_pivots = DEFAULT_MAX_PIVOTS - # Problem dimensions - @num_non_slack_vars = a.first.length - @num_constraints = b.length + # Problem dimensions; these never change + @num_non_slack_vars = num_vars + @num_constraints = num_inequalities @num_vars = @num_non_slack_vars + @num_constraints # Set up initial matrix A and vectors b, c - @c = Vector[*c.map {|c1| -1*c1 } + [0]*@num_constraints] - @a = a.map {|a1| Vector[*(a1.clone + [0]*@num_constraints)]} - @b = Vector[*b.clone] - - unless @a.all? {|a| a.size == @c.size } and @b.size == @a.length - raise ArgumentError, "Input arrays have mismatched dimensions" - end - - 0.upto(@num_constraints - 1) {|i| @a[i][@num_non_slack_vars + i] = 1 } + @c = c.map { |flt| -1 * flt } + Array.new(@num_constraints, 0) + @a = a.map.with_index { |ary, i| + if ary.size != @num_non_slack_vars + raise ArgumentError, "a is inconsistent" + end + ary + Array.new(@num_constraints) { |ci| ci == i ? 1 : 0 } + } + @b = b # set initial solution: all non-slack variables = 0 - @x = Vector[*([0]*@num_vars)] @basic_vars = (@num_non_slack_vars...@num_vars).to_a - update_solution - end - - def solution - solve - current_solution - end - - def current_solution - @x.to_a[0...@num_non_slack_vars] + self.update_solution end + # does not modify vector / matrix def update_solution - 0.upto(@num_vars - 1) {|i| @x[i] = 0 } + @x = Array.new(@num_vars, 0) + + @basic_vars.each { |basic_var| + idx = nil + @num_constraints.times { |i| + if @a[i][basic_var] == 1 + idx =i + break + end + } + raise(SanityCheck, "no idx for basic_var #{basic_var} in a") unless idx + @x[basic_var] = @b[idx] + } + end - @basic_vars.each do |basic_var| - row_with_1 = row_indices.detect do |row_ix| - @a[row_ix][basic_var] == 1 - end - @x[basic_var] = @b[row_with_1] - end + def solution + self.solve + self.current_solution end def solve - while can_improve? - @pivot_count += 1 - raise "Too many pivots" if @pivot_count > max_pivots - pivot + count = 0 + while self.can_improve? + count += 1 + raise(TooManyPivots, count.to_s) unless count < @max_pivots + self.pivot end end - def can_improve? - !!entering_variable + def current_solution + @x[0...@num_non_slack_vars] end - def variables - (0...@c.size).to_a + def can_improve? + !self.entering_variable.nil? end + # idx of @c's minimum negative value + # nil when no improvement is possible + # def entering_variable - variables.select { |var| @c[var] < 0 }. - min_by { |var| @c[var] } + (0...@c.size).select { |i| @c[i] < 0 }.min_by { |i| @c[i] } end def pivot - pivot_column = entering_variable - pivot_row = pivot_row(pivot_column) - raise UnboundedProblem unless pivot_row - leaving_var = basic_variable_in_row(pivot_row) - replace_basic_variable(leaving_var => pivot_column) + pivot_column = self.entering_variable or return nil + pivot_row = self.pivot_row(pivot_column) or raise UnboundedProblem + leaving_var = nil + @a[pivot_row].each_with_index { |a, i| + if a == 1 and @basic_vars.include?(i) + leaving_var = i + break + end + } + raise(SanityCheck, "no leaving_var") if leaving_var.nil? + + @basic_vars.delete(leaving_var) + @basic_vars.push(pivot_column) + @basic_vars.sort! pivot_ratio = Rational(1, @a[pivot_row][pivot_column]) # update pivot row - @a[pivot_row] *= pivot_ratio - @b[pivot_row] = pivot_ratio * @b[pivot_row] + @a[pivot_row] = @a[pivot_row].map { |val| val * pivot_ratio } + @b[pivot_row] = @b[pivot_row] * pivot_ratio # update objective - @c -= @c[pivot_column] * @a[pivot_row] + # @c -= @c[pivot_column] * @a[pivot_row] + @c = @c.map.with_index { |val, i| + val - @c[pivot_column] * @a[pivot_row][i] + } # update A and B - (row_indices - [pivot_row]).each do |row_ix| - r = @a[row_ix][pivot_column] - @a[row_ix] -= r * @a[pivot_row] - @b[row_ix] -= r * @b[pivot_row] - end + @num_constraints.times { |i| + next if i == pivot_row + r = @a[i][pivot_column] + @a[i] = @a[i].map.with_index { |val, j| val - r * @a[pivot_row][j] } + @b[i] = @b[i] - r * @b[pivot_row] + } - update_solution - end - - def replace_basic_variable(hash) - from, to = hash.keys.first, hash.values.first - @basic_vars.delete(from) - @basic_vars << to - @basic_vars.sort! + self.update_solution end def pivot_row(column_ix) - row_ix_a_and_b = row_indices.map { |row_ix| - [row_ix, @a[row_ix][column_ix], @b[row_ix]] - }.reject { |_, a, b| - a == 0 - }.reject { |_, a, b| - (b < 0) ^ (a < 0) # negative sign check - } - row_ix, _, _ = *last_min_by(row_ix_a_and_b) { |_, a, b| - Rational(b, a) + min_ratio = nil + idx = nil + @num_constraints.times { |i| + a, b = @a[i][column_ix], @b[i] + next if a == 0 or (b < 0) ^ (a < 0) + ratio = Rational(b, a) + idx, min_ratio = i, ratio if min_ratio.nil? or ratio <= min_ratio } - row_ix - end - - def basic_variable_in_row(pivot_row) - column_indices.detect do |column_ix| - @a[pivot_row][column_ix] == 1 and @basic_vars.include?(column_ix) - end - end - - def row_indices - (0...@a.length).to_a - end - - def column_indices - (0...@a.first.size).to_a + idx end def formatted_tableau - if can_improve? - pivot_column = entering_variable - pivot_row = pivot_row(pivot_column) + if self.can_improve? + pivot_column = self.entering_variable + pivot_row = self.pivot_row(pivot_column) else pivot_row = nil end - num_cols = @c.size + 1 - c = formatted_values(@c.to_a) - b = formatted_values(@b.to_a) - a = @a.to_a.map {|ar| formatted_values(ar.to_a) } + c = @c.to_a.map { |flt| "%2.3f" % flt } + b = @b.to_a.map { |flt| "%2.3f" % flt } + a = @a.to_a.map { |vec| vec.to_a.map { |flt| "%2.3f" % flt } } if pivot_row a[pivot_row][pivot_column] = "*" + a[pivot_row][pivot_column] end max = (c + b + a + ["1234567"]).flatten.map(&:size).max result = [] - result << c.map {|c| c.rjust(max, " ") } + result << c.map { |str| str.rjust(max, " ") } a.zip(b) do |arow, brow| - result << (arow + [brow]).map {|a| a.rjust(max, " ") } + result << (arow + [brow]).map { |val| val.rjust(max, " ") } result.last.insert(arow.length, "|") end - lines = result.map {|b| b.join(" ") } + lines = result.map { |ary| ary.join(" ") } max_line_length = lines.map(&:length).max lines.insert(1, "-"*max_line_length) lines.join("\n") end - - def formatted_values(array) - array.map {|c| "%2.3f" % c } - end - - # like Enumerable#min_by except if multiple values are minimum - # it returns the last - def last_min_by(array) - best_element, best_value = nil, nil - array.each do |element| - value = yield element - if !best_element || value <= best_value - best_element, best_value = element, value - end - end - best_element - end - - def assert(boolean) - raise unless boolean - end - end - diff --git a/lib/simplex/parse.rb b/lib/simplex/parse.rb new file mode 100644 index 0000000..00f360e --- /dev/null +++ b/lib/simplex/parse.rb @@ -0,0 +1,123 @@ +class Simplex + module Parse + class Error < RuntimeError; end + class InvalidExpression < Error; end + class InvalidInequality < Error; end + class InvalidTerm < Error; end + + # coefficient concatenated with a single letter variable, e.g. "-1.23x" + TERM_RGX = %r{ + \A # starts with + (-)? # possible negative sign + (\d+(?:\.\d*)?)? # possible float (optional) + ([a-zA-Z]) # single letter variable + \z # end str + }x + + # a float or integer, possibly negative + CONSTANT_RGX = %r{ + \A # starts with + -? # possible negative sign + \d+ # integer portion + (?:\.\d*)? # possible decimal portion + \z # end str + }x + + def self.inequality(str) + lhs, rhs = str.split('<=') + if lhs.nil? or lhs.empty? or rhs.nil? or rhs.empty? + raise(InvalidInequality, "#{str}") + end + rht = self.tokenize(rhs) + raise(InvalidInequality, "#{str}; bad rhs: #{rhs}") unless rht.size == 1 + c = rht.first + raise(InvalidInequality, "bad rhs: #{rhs}") if !c.match CONSTANT_RGX + return self.expression(lhs), c.to_f + end + + # ignore leading and trailing spaces + # ignore multiple spaces + def self.tokenize(str) + str.strip.split(/\s+/) + end + + # rules: variables are a single letter + # may have a coefficient (default: 1.0) + # only sum and difference operations allowed + # normalize to all sums with possibly negative coefficients + # valid inputs: + # 'x + y' => [1.0, 1.0], [:x, :y] + # '2x - 5y' => [2.0, -5.0], [:x, :y] + # '-2x - 3y + -4z' => [-2.0, -3.0, -4.0], [:x, :y, :z] + def self.expression(str) + terms = self.tokenize(str) + negative = false + coefficients = {} + while !terms.empty? + # consume plus and minus operations + term = terms.shift + if term == '-' + negative = true + term = terms.shift + elsif term == '+' + negative = false + term = terms.shift + end + + coefficient, variable = self.term(term) + raise("double variable: #{str}") if coefficients.key?(variable) + coefficients[variable] = negative ? coefficient * -1 : coefficient + end + coefficients + end + + def self.term(str) + matches = str.match TERM_RGX + raise(InvalidTerm, str) unless matches + flt = (matches[2] || 1).to_f * (matches[1] ? -1 : 1) + sym = matches[3].to_sym # consider matches[3].downcase.to_sym + return flt, sym + end + end + + def self.problem(maximize: nil, constraints: [], **kwargs) + if maximize + obj, maximize = maximize, true + elsif kwargs[:minimize] + obj, maximize = kwargs[:minimize], false + else + raise(ArgumentError, "one of maximize/minimize expected") + end + unless obj.is_a?(String) + raise(ArgumentError, "bad expr: #{expr} (#{expr.class})") + end + obj_cof = Parse.expression(obj) + + c = [] # coefficients of objective expression + a = [] # array (per constraint) of the inequality's lhs coefficients + b = [] # rhs (constant) for the inequalities / constraints + + # this determines the order of coefficients + letter_vars = obj_cof.keys + letter_vars.each { |v| c << obj_cof[v] } + + constraints.each { |str| + unless str.is_a?(String) + raise(ArgumentError, "bad constraint: #{str} (#{str.class})") + end + cofs = [] + ineq_cofs, rhs = Parse.inequality(str) + letter_vars.each { |v| + raise("constraint #{str} is missing var #{v}") unless ineq_cofs.key?(v) + cofs << ineq_cofs[v] + } + a.push cofs + b.push rhs + } + self.new(c, a, b) + end + + def self.maximize(expression, *ineqs) + self.problem(maximize: expression, constraints: ineqs).solution + end +end diff --git a/metrics/bench b/metrics/bench new file mode 100644 index 0000000..32918eb --- /dev/null +++ b/metrics/bench @@ -0,0 +1,5 @@ +/proc/loadavg 0.01 0.01 0.00 1/77 1173 +Warming up -------------------------------------- + Simplex Array 84.000 i/100ms +Calculating ------------------------------------- + Simplex Array 849.092 (± 1.3%) i/s - 2.604k in 3.067337s diff --git a/test/bench/simplex.rb b/test/bench/simplex.rb new file mode 100644 index 0000000..10728e0 --- /dev/null +++ b/test/bench/simplex.rb @@ -0,0 +1,98 @@ +require 'simplex' +require 'benchmark/ips' + +Benchmark.ips do |b| + b.config time: 3, warmup: 0.5 + + b.report("Simplex Array") { + Simplex.new([1, 1], + [[2, 1], + [1, 2]], + [4, 3]).solution + + Simplex.new([3, 4], + [[1, 1], + [2, 1]], + [4, 5]).solution + + Simplex.new([2, -1], + [[1, 2], + [3, 2],], + [6, 12]).solution + + Simplex.new([60, 90, 300], + [[1, 1, 1], + [1, 3, 0], + [2, 0, 1]], + [600, 600, 900]).solution + + Simplex.new([70, 210, 140], + [[1, 1, 1], + [5, 4, 4], + [40, 20, 30]], + [100, 480, 3200]).solution + + Simplex.new([2, -1, 2], + [[2, 1, 0], + [1, 2, -2], + [0, 1, 2]], + [10, 20, 5]).solution + + Simplex.new([11, 16, 15], + [[1, 2, Rational(3, 2)], + [Rational(2, 3), Rational(2, 3), 1], + [Rational(1, 2), Rational(1, 3), Rational(1, 2)]], + [12_000, 4_600, 2_400]).solution + + Simplex.new([5, 4, 3], + [[2, 3, 1], + [4, 1, 2], + [3, 4, 2]], + [5, 11, 8]).solution + + Simplex.new([3, 2, -4], + [[1, 4, 0], + [2, 4,-2], + [1, 1,-2]], + [5, 6, 2]).solution + + Simplex.new([2, -1, 8], + [[2, -4, 6], + [-1, 3, 4], + [0, 0, 2]], + [3, 2, 1]).solution + + Simplex.new([100_000, 40_000, 18_000], + [[20, 6, 3], + [0, 1, 0], + [-1,-1, 1], + [-9, 1, 1]], + [182, 10, 0, 0]).solution + + Simplex.new([1, 2, 1, 2], + [[1, 0, 1, 0], + [0, 1, 0, 1], + [1, 1, 0, 0], + [0, 0, 1, 1]], + [1, 4, 2, 2]).solution + + Simplex.new([10, -57, -9, -24], + [[0.5, -5.5, -2.5, 9], + [0.5, -1.5, -0.5, 1], + [ 1, 0, 0, 0]], + [0, 0, 1]).solution + + Simplex.new([25, 20], + [[20, 12], + [1, 1]], + [1800, 8*15]).solution + } + +#b.report("Simplex Array") { +#} + +#b.report("Simplex Matrix") { +#} + +# b.compare! +end diff --git a/test/parse.rb b/test/parse.rb new file mode 100644 index 0000000..73f816b --- /dev/null +++ b/test/parse.rb @@ -0,0 +1,92 @@ +require 'simplex/parse' +require 'minitest/autorun' + +describe Simplex::Parse do + P = Simplex::Parse + + describe "Parse.term" do + it "must parse valid terms" do + { "-1.2A" => [-1.2, :A], + "99x" => [99.0, :x], + "z" => [ 1.0, :z], + "-b" => [-1.0, :b] }.each { |valid, expected| + P.term(valid).must_equal expected + } + end + + it "must reject invalid terms" do + ["3xy", "24/7x", "x17", "2*x"].each { |invalid| + proc { P.term(invalid) }.must_raise RuntimeError + } + end + end + + describe "Parse.expression" do + it "must parse valid expressions" do + { "x + y" => { x: 1.0, y: 1.0 }, + "2x - 5y" => { x: 2.0, y: -5.0 }, + "-2x - 3y + -42.7z" => { x: -2.0, y: -3.0, z: -42.7 }, + " -5y + -x " => { y: -5.0, x: -1.0 }, + "a - -b" => { a: 1.0, b: 1.0 }, + "a A b" => { a: 1.0, :A => 1.0, b: 1.0 }, + }.each { |valid, expected| + P.expression(valid).must_equal expected + } + end + + it "must reject invalid expressions" do + ["a2 + b2 = c2", + "x + xy", + "x * 2"].each { |invalid| + proc { P.expression(invalid) }.must_raise P::Error + } + end + end + + describe "Parse.tokenize" do + it "ignores leading or trailing whitespace" do + P.tokenize(" 5x + 2.9y ").must_equal ["5x", "+", "2.9y"] + end + + it "ignores multiple spaces" do + P.tokenize("5x + 2.9y").must_equal ["5x", "+", "2.9y"] + end + end + + describe "Parse.inequality" do + it "must parse valid inequalities" do + { "x + y <= 4" => [{ x: 1.0, y: 1.0 }, 4.0], + "0.94a - 22.1b <= -14.67" => [{ a: 0.94, b: -22.1 }, -14.67], + "x <= 0" => [{ x: 1.0 }, 0], + }.each { |valid, expected| + P.inequality(valid).must_equal expected + } + end + + it "must reject invalid inequalities" do + ["x + y >= 4", + "0.94a - 22.1b <= -14.67c", + "x < 0", + ].each { |invalid| + proc { P.inequality(invalid) }.must_raise P::Error + } + end + end +end + +describe "Simplex.maximize" do + it "must problem stuff" do + prob = Simplex.problem(maximize: 'x + y', + constraints: ['2x + y <= 4', + 'x + 2y <= 3']) + sol = prob.solution + sol.must_equal [Rational(5, 3), Rational(2, 3)] + end + + it "must maximize stuff" do + Simplex.maximize('x + y', + '2x + y <= 4', + 'x + 2y <= 3').must_equal [Rational(5, 3), + Rational(2, 3)] + end +end diff --git a/test/simplex_test.rb b/test/simplex.rb similarity index 92% rename from test/simplex_test.rb rename to test/simplex.rb index 6b406e8..125e67c 100644 --- a/test/simplex_test.rb +++ b/test/simplex.rb @@ -1,8 +1,7 @@ -require 'test/unit' -$:.push(File.expand_path("../../lib", __FILE__)) +require 'minitest/autorun' require 'simplex' -class SimplexTest < Test::Unit::TestCase +class SimplexTest < Minitest::Test def test_2x2 result = Simplex.new( [1, 1], @@ -180,14 +179,14 @@ def test_cycle2 ], [0, 0] ) - assert_raise Simplex::UnboundedProblem do + assert_raises Simplex::UnboundedProblem do simplex.solution end end def test_error_mismatched_dimensions - assert_raise ArgumentError do - result = Simplex.new( + assert_raises ArgumentError do + Simplex.new( [10, -57, -9], [ [0.5, -5.5, -2.5, 9], @@ -198,8 +197,8 @@ def test_error_mismatched_dimensions ) end - assert_raise ArgumentError do - result = Simplex.new( + assert_raises ArgumentError do + Simplex.new( [10, -57, -9, 2], [ [0.5, -5.5, 9, 4], @@ -210,8 +209,8 @@ def test_error_mismatched_dimensions ) end - assert_raise ArgumentError do - result = Simplex.new( + assert_raises ArgumentError do + Simplex.new( [10, -57, -9, 2], [ [0.5, -5.5, 9, 4], @@ -263,16 +262,16 @@ def test_cup_factory # [1, -2] # ) # while simplex.can_improve? - # puts + # puts # puts simplex.formatted_tableau # simplex.pivot # end # p :done - # puts + # puts # puts simplex.formatted_tableau #end - + def test_unbounded simplex = Simplex.new( [1, 1, 1], @@ -282,7 +281,7 @@ def test_unbounded ], [5, 7] ) - assert_raise Simplex::UnboundedProblem do + assert_raises Simplex::UnboundedProblem do simplex.solution end end