simplex.rb

#--
# MAIN MODULE -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
#++
module Simplex

  #--
  # UTILITIES -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
  #++
  module Util

    #
    # Converts number to Rational
    #
    def rationalize(num)
      case num
      when Rational
        num
      when Integer, String
        num.to_r
      when Array
        num.map { |n| rationalize(n) }
      else
        rationalize(num.to_s) # Float
      end
    end

    #
    # Tells if a number is numeric.
    #
    def numeric?(num)
      num.is_a?(Numeric)
    end

  end


  #--
  # BIG 'M' CLASS -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
  #++
  class BigM
    include Comparable # For <, >, ==
    include Util       # For utilities

    #
    # Public attributes
    #
    attr_reader :m, :c

    #
    # Constructs a new BigM.
    #
    def initialize(num = nil)
      @m = @c = 0r

      case num
      when BigM
        @m, @c = num.to_a
      when Numeric
        @c = rationalize(num)
      when Array
        tm, tc = rationalize(num)

        @m = tm if tm
        @c = tc if tc
      when String
        num = num.strip

        unless num.empty?
          # Remove other than [0-9], +, -, or M
          num = num.upcase.gsub(/[^\d\-\+M]/, '')

          if num.include?('M')
            tm, tc = num.split('M')

            if tm
              @m = rationalize(tm == '-' ? -1 : (tm.empty? ? 1 : tm))
            else
              @m = 1r
            end
          else
            @c = rationalize(num)
          end
        end
      end if num
    end

    #
    # Unary plus operation.
    #
    def +@
      make(m, c)
    end

    #
    # Unary minus operation.
    #
    def -@
      make(-m, -c)
    end

    #
    # Plus operation.
    #
    def +(num)
      tm, tc = values(num)

      make(m + tm, c + tc)
    end

    #
    # Minus operation.
    #
    def -(num)
      self + -num
    end

    #
    # Multiplication operation.
    #
    def *(num)
      raise Error, "Can't multiply by #{num.class}." unless numeric?(num)

      make(m * num, c * num)
    end

    #
    # Division operation.
    #
    def /(num)
      raise Error, "Can't divide by  #{num.class}." unless numeric?(num)
      raise Error, "Can't divide by zero." if num == 0

      make(m / num, c / num)
    end

    #
    # Comparator with other numbers.
    #
    def <=>(num)
      case num
      when BigM
        # Compare entities
        m == num.m ? (c <=> num.c) : (m < num.m ? -1 : 1)
      when Numeric
        # If zero, compare integers
        # Else, negative M or positive M wins
        m == 0 ? (c <=> num) : (m < 0 ? -1 : 1)
      else
        raise Error, "Can't compare with #{num.class}."
      end
    end

    #
    # Returns array representation.
    #
    def to_a
      [m, c]
    end

    #
    # Returns string representation.
    #
    def to_s
      res = ''

      # If m is not 0
      unless m == 0
        am = m.abs
        res += '-' if m < 0
        res += am.to_s unless am == 1
        res += 'M'
      end

      # If c is not 0
      unless c == 0
        res += '+' if c > 0 and !res.empty?
        res += c.to_s
      end

      # If neither one, fallback c
      res += c.to_s if res.empty?

      res
    end

    #
    # Coerce num to BigM.
    #
    def coerce(num)
      [make(0, num), self]
    end

    private

    #
    # Creates a new BigM.
    #
    def make(m, c)
      self.class.new([m, c])
    end

    #
    # Parses values into array.
    #
    def values(num)
      case num
      when BigM
        num.to_a
      when Numeric
        [0, num]
      else
        raise Error, "Invalid values type #{num.class}."
      end
    end


    #
    # BigM useful constants.
    #
    ZERO = new(0)
    ONE  = new([1])


    #
    # A generic BigM Error.
    #
    class Error < StandardError; end

  end


  #--
  # MATRIX UTILS -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
  #++
  class Matrix

    #
    # Singleton Matrix.
    #
    class << self

      #
      # Creates an identity matrix.
      #
      def identity(n)
        Array.new(n) do |i|
          tmp = Array.new(n, 0)
          tmp[i] = 1
          tmp
        end
      end

      #
      # Multiply two matrices.
      #
      def multiply(a, b)
        rca = a.length
        cca = a[0].length
        rcb = b.length
        ccb = b[0].length

        raise Error, 'Matrix dimension mismatch.' if cca != rcb

        Array.new(rca) do |i|
          Array.new(ccb) do |j|
            (0...cca).inject(0) do |val, k|
              val + a[i][k] * b[k][j]
            end
          end
        end
      end

    end


    #
    # A generic Matrix Error.
    #
    class Error < StandardError; end

  end


  #--
  # PROBLEM CLASSES -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
  #++
  class Problem

    #
    # Public attributes
    #
    attr_reader :objective, :constraints, :order

    #
    # Constructs a new Problem.
    #
    def initialize(objective, constraints)
      @objective = objective
      @constraints = constraints

      @order = objective.order if objective.respond_to?(:order)

      @standard = false
      @was_min = false

      validate
    end

    #
    # Tells if Problem is in standard form or not.
    #
    def standard?
      @standard
    end

    #
    # Tells if the Problem was initially MIN.
    #
    def was_min?
      @was_min
    end

    #
    # Converts Problem to standard form.
    #
    def standarize!
      unless standard?
        # Change objective from MIN to MAX
        if objective.type == Objective::Type::MIN
          @was_min = true

          objective.coefficients.map! { |coefficient| -coefficient }
        end

        # For each constraint we have
        constraints.each_with_index do |constraint, idx1|
          # Get slack or artificial coefficient
          slack, artif = case constraint.type
                         when Constraint::Type::GREATER
                           [-1r, 1r]
                         when Constraint::Type::EQUAL
                           [nil, 1r]
                         when Constraint::Type::LESS
                           [1r, nil]
                         end

          # Update constraint itself
          constraint.lhs << slack if slack
          constraint.lhs << artif if artif

          # Update objective function
          objective.coefficients << 0r         if slack
          objective.coefficients << -BigM::ONE if artif

          # For each constraint we have
          constraints.each_with_index do |second, idx2|
            # Except current one
            next if idx1 == idx2

            # Update other constraint
            second.lhs << 0r if slack
            second.lhs << 0r if artif
          end
        end

        # Mark as standarized
        @standard = true
      end
    end

    private

    #
    # Validates constructor arguments.
    #
    def validate
      raise Error, 'Objective must be instance of Objective.' unless objective.is_a?(Objective)
      raise Error, 'Constraints must be array.' unless constraints.is_a?(Array)

      constraints.each do |constraint|
        raise Error, 'Constraint must be instance of Constraint.' unless constraint.is_a?(Constraint)
        raise Error, 'Constraint order does not match.' unless constraint.order == order
      end
    end


    #
    # A generic Problem Error.
    #
    class Error < StandardError; end


    #
    # Class for a Problem Objective.
    #
    class Objective
      include Util # For utilities

      #
      # Public attributes.
      #
      attr_reader :type, :coefficients, :order

      #
      # Constructs a new Problem Objective.
      #
      def initialize(type, coefficients)
        @type = type
        @coefficients = rationalize(coefficients)

        @order = coefficients.length if coefficients.respond_to?(:length)

        validate
      end

      private

      #
      # Validates constructor arguments.
      #
      def validate
        sanitize_type

        raise Error, 'Coefficients must be array.' unless coefficients.is_a?(Array)
        raise Error, 'Invalid objective variables length.' if order < 2

        coefficients.each do |coefficient|
          raise Error, 'Coefficient must be a number.' unless numeric?(coefficient)
        end
      end

      #
      # Converts Objective type to constant.
      #
      def sanitize_type
        case type
        when Fixnum
          raise Error, 'Unrecognized objective type.' unless [Type::MAX, Type::MIN].include?(type)
        when String
          case type.downcase[0..2]
          when 'max'
            @type = Type::MAX
          when 'min'
            @type = Type::MIN
          else
            raise Error, "Unrecognized objective type: #{type}."
          end
        else
          raise Error, 'Invalid objective type.'
        end
      end


      #
      # Problem Objective type constants.
      #
      module Type
        MAX = 0x15
        MIN = 0x30
      end


      #
      # A generic Problem Objective Error.
      #
      class Error < StandardError; end

    end


    #
    # Class for a Problem Constraint.
    #
    class Constraint
      include Util # For utilities

      #
      # Public attributes.
      #
      attr_reader :lhs, :type, :rhs, :order

      #
      # Constructs a new Problem Constraint.
      #
      def initialize(lhs, type, rhs)
        @lhs = rationalize(lhs)
        @type = type
        @rhs = rationalize(rhs)

        @order = lhs.length if lhs.respond_to?(:length)

        validate
      end

      private

      #
      # Validates constructor arguments.
      #
      def validate
        sanitize_type

        raise Error, 'Constraint LHS must be an array.' unless lhs.is_a?(Array)
        raise Error, "Invalid constraint LHS length: #{lhs.length}." if lhs.length < 2

        lhs.each do |coefficient|
          raise Error, 'Coefficient must be a number.' unless numeric?(coefficient)
        end

        raise Error, 'Invalid constraint RHS.' unless numeric?(rhs)
      end

      #
      # Converts Constraint type to constant.
      #
      def sanitize_type
        case type
        when Fixnum
          raise Error, 'Unrecognized constraint type.' unless [Type::GREATER, Type::EQUAL, Type::LESS].include?(type)
        when String
          case type.strip
          when '>='
            @type = Type::GREATER
          when '=', '=='
            @type = Type::EQUAL
          when '<='
            @type = Type::LESS
          else
            raise Error, "Unrecognized constraint type: #{type}."
          end
        else
          raise Error, 'Invalid constraint type.'
        end
      end


      #
      # Problem Constraint type constants.
      #
      module Type
        GREATER = 0x12
        EQUAL   = 0x24
        LESS    = 0x36
      end


      #
      # A generic Problem Constraint Error.
      #
      class Error < StandardError; end

    end

  end


  #--
  # SIMPLEX SOLVER -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
  #++
  class Solver

    #
    # Public attributes.
    #
    attr_reader :problem

    #
    # Constructs a new Solver.
    #
    def initialize(problem)
      @problem = problem

      validate
      prepare
    end

    def solve
      unless is_optimal?
        # Get in and out variables
        vin, vout = get_inout
        voutIdx = @Vb.transpose[0].index(vout)
        vinIdx = @x.transpose[0].index(vin)

        # Update Vb and Cb
        @Vb[voutIdx][0] = vin
        @Cb[0][voutIdx] = @c[0][vinIdx]

        @E = Matrix.identity(@B.length).map { |row| row.map(&:to_r) }

        # Get eta column for exchange
        piv = @Aj[voutIdx][0]
        eta = Array.new(@Aj.length)
        @Aj.each_with_index do |row, idx2|
          tmp = row[0]

          eta[idx2] = [(idx2 == voutIdx ? 1 : -tmp) / piv]
        end

        # Exchange eta into (E)
        @E.each_with_index do |row, r|
          row.each_with_index do |col, c|
            next unless c == voutIdx

          @E[r][c] = eta[r][0]
          end
        end

        # Update (B-1) as (E)*(B)
        @B = Matrix.multiply(@E, @B)

        # Update (xb) as (B)*(b)
        @xb = Matrix.multiply(@B, @b)

        # Iterate again
        return solve
      end

      res = {}

      tmp = @Vb.transpose[0]
      @x.transpose[0].each_with_index do |var, idx|
        fld = tmp.index(var)

        res[var + '*'] = fld ? @xb[fld][0] : 0
      end

      res['Z*'] = Matrix.multiply(@Cb, @xb)[0][0]
      res['Z*'] *= -1 if problem.was_min?

      res
    end

    private

    #
    # Prepares variables for first iteration.
    #
    def prepare
      problem.standarize!

      objective = problem.objective
      constraints = problem.constraints

      # Calculate coefficients (c)
      @c = [objective.coefficients]

      # Calculate table (A)
      @A = []

      # Calculate LHS (b)
      @b = []

      # Fill (A) and (b)
      constraints.each do |constraint|
        @A << constraint.lhs
        @b << [constraint.rhs]
      end

      # Get B-1 (B)
      @B = Matrix.identity(constraints.length).map { |row| row.map(&:to_r) }

      # Calculate basic LHS (xb)
      @xb = @b.clone

      # Fill all variables (x)
      @x = objective.coefficients.length.times.map { |i| ['x' + (i + 1).to_s] }

      # Find basic variables (Vb)
      @Vb = find_vb

      # Find basic coefficients (Cb)
      @Cb = [@Vb.map { |arr| @c[0][@x.transpose[0].index(arr[0])] }]
    end

    #
    # Calculates Zj - Cj.
    #
    def zjcj
      @nonVb = @x.transpose[0] - @Vb.transpose[0]

      @zjcj = @nonVb.map do |var|
        idx = @x.transpose[0].index(var)

        a = [@A.transpose[idx]].transpose
        tmp = Matrix.multiply(@Cb, @B)
        tmp = Matrix.multiply(tmp, a)

        tmp[0][0] - @c[0][idx]
      end
    end

    #
    # Checks if we have an optimal solution.
    #
    def is_optimal?
      zjcj.select { |val| val < 0 }.length == 0
    end

    #
    # Gets in and out variables.
    #
    def get_inout
      vin = @nonVb[@zjcj.index(@zjcj.min)]
      idx = @x.transpose[0].index(vin)

      @Aj = Matrix.multiply(@B, [@A.transpose[idx]].transpose)
      vout = nil
      min = Float::MAX
      @xb.transpose[0].each_with_index do |val, idx2|
        next if @Aj[idx2][0] <= 0

        tmp = val / @Aj[idx2][0]
        if tmp < min
          min = tmp
          vout = @Vb[idx2][0]
        end
      end

      [vin, vout]
    end

    #
    # Finds basic variables in (A).
    #
    def find_vb
      tmp = @A.transpose

      @B.map { |row| [@x[tmp.index(row)][0]] }
    end

    #
    # Validates constructor arguments.
    #
    def validate
      raise Error, 'Solver expects a Problem.' unless problem.is_a?(Problem)
    end


    #
    # A generic Solver Error.
    #
    class Error < StandardError; end

  end


  #--
  # CONSOLE CLIENT -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
  #++
  class CLI

    #
    # Public attributes.
    #
    attr_reader :order

    #
    # Starts asking for problem.
    #
    def start
      @order = ask('Number of variables').to_i
      raise Error, 'Invalid number of variables.' if @order < 2

      objective = ask_objective

      contr = ask('Number of constraints').to_i
      raise Error, 'Invalid number of constraints.' if contr == 0

      constraints = []
      contr.times do |i|
        constraints << ask_constraint(i)
      end

      problem = Problem.new(objective, constraints)
      problem.standarize!

      solver = Solver.new(problem)

      res = solver.solve

      puts 'Result:'
      res.each { |key, value| puts "#{key} = #{value}" }
    end

    private

    #
    # Asks a question and returns answer.
    #
    def ask(msg)
      print "#{msg}: "
      gets.chomp.strip
    end

    #
    # Asks and builds a Problem Objective.
    #
    def ask_objective
      objfn = ask('Objective function').gsub(/[^\w\+\-\=\(\)]/, '')
      match = objfn.match(/^(Max|Min)\(([A-Za-z]+)\)\=(.*)$/)
      raise Error, 'Invalid objective function.' unless match

      type = match[1]
      terms = parse_terms(match[3])
      coefficients = make_coefficients(terms)

      Problem::Objective.new(type, coefficients)
    end

    #
    # Asks and builds a Problem Constraint.
    #
    def ask_constraint(i)
      contr = ask("Constraint #{i + 1}").gsub(/[^x\d\+\-\<\>\=]/, '')
      match = contr.match(/^([^\<\>\=]+)([\<\>\=]?\=)(\d+)$/)
      raise Error, 'Invalid constraint.' unless match

      lhs = make_coefficients(parse_terms(match[1]))
      sign = match[2]
      rhs = match[3].to_r

      Problem::Constraint.new(lhs, sign, rhs)
    end

    #
    # Parses a sequence of algebraic terms.
    #
    def parse_terms(terms)
      match = terms.scan(/(\+|\-)?(\d+)?x(\d+)/)

      match.map do |term|
        sign, coefficient, variable = term
        sign = '+' unless sign
        coefficient = 1 unless coefficient

        coefficient = coefficient.to_r
        coefficient = -coefficient if sign == '-'

        [coefficient, variable.to_i]
      end
    end

    #
    # Turns term arrays into coefficients.
    #
    def make_coefficients(terms)
      res = Array.new(order + 1, 0)
      terms.each do |term|
        res[term[1]] = term[0]
      end
      res.shift

      res
    end


    #
    # A generic CLI Error.
    #
    class Error < StandardError; end

  end

end


#--
# BUILT-IN OVERRIDES -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
#++
class Rational

  #
  # Overrides to_s to remove '/1' portions.
  #
  alias_method :old_to_s, :to_s
  private :old_to_s
  def to_s
    return numerator.to_s if denominator == 1
    old_to_s
  end

  #
  # Overrides comparison to work with BigM.
  #
  alias_method :old_cmp, :<=>
  private :old_cmp
  def <=>(num)
    if num.is_a?(Simplex::BigM)
      return -(num <=> self)
    end
    old_cmp(num)
  end

end


# Start program
cli = Simplex::CLI.new
cli.start