module Unodos::Solver

Constants

NAMED_BASES
TOLERANCE

Public Class Methods

find_solve(vectors, bvector, &block) click to toggle source
# File lib/unodos/solver.rb, line 72
def self.find_solve(vectors, bvector, &block)
  (1...vectors.size).each do |i|
    least_square_solve [vectors[0], vectors[i]], bvector do |vs, pos|
      block.call vs, pos.map { |c| c == 1 ? i : 0 }
    end
  end
end
least_square_solve(vectors, bvector, &block) click to toggle source
# File lib/unodos/solver.rb, line 93
def self.least_square_solve(vectors, bvector, &block)
  mat = Matrix[*vectors.transpose]
  tmat = mat.transpose
  m = tmat * mat
  b = tmat * Vector[*bvector]
  vs = lup_solve m.lup, b.to_a
  return unless vs
  max_diff = bvector.each_with_index.map do |bv, i|
    (vectors.zip(vs).sum { |a, v| v * a[i] } - bv).abs
  end.max
  block.call vs, (0...vectors.size).to_a if max_diff < TOLERANCE
end
lup_solve(lup, b) click to toggle source
# File lib/unodos/solver.rb, line 13
def self.lup_solve(lup, b)
  size = b.size
  m = b.values_at(*lup.pivots)
  mat_l = lup.l
  mat_u = lup.u
  size.times do |k|
    (k + 1).upto(size - 1) do |i|
      m[i] -= m[k] * mat_l[i, k]
    end
  end
  (size - 1).downto(0) do |k|
    next if m[k] == 0
    return nil if mat_u[k, k].zero?
    m[k] = m[k].quo mat_u[k, k]
    k.times do |i|
      m[i] -= m[k] * mat_u[i, k]
    end
  end
  m
end
match_vector(vector, bvector) click to toggle source
# File lib/unodos/solver.rb, line 65
def self.match_vector(vector, bvector)
  v, b = vector.zip(bvector).max_by { |a,| a.abs }
  a = b.quo v
  err = vector.zip(bvector).map { |v, b| (v * a - b).abs }.max
  a if err < TOLERANCE
end
recursive_solve(vectors, bvector, first_required, &block) click to toggle source
# File lib/unodos/solver.rb, line 106
def self.recursive_solve(vectors, bvector, first_required, &block)
  size = vectors.size
  out_size = bvector.size
  lup = Matrix[*vectors.transpose].lup
  mat_l = lup.l
  mat_u = lup.u.to_a
  bvector = bvector.values_at(*lup.pivots)
  out_size.times do |k|
    (k + 1).upto(out_size - 1) do |i|
      bvector[i] -= bvector[k] * mat_l[i, k]
    end
  end
  solved = lambda do |u, selected|
    b = bvector.dup
    (selected.size - 1).downto 0 do |i|
      j = selected[i]
      next if b[i] == 0
      return if u[i][j] == 0
      b[i] = b[i].quo u[i][j]
      i.times do |k|
        b[k] -= b[i] * u[k][j]
      end
    end
    block.call b, selected
  end
  solve = lambda do |u, selected, index|
    return solved.call u, selected if selected.size == out_size
    j = selected.size
    restore_index = (j .. [index, out_size - 1].min).max_by do |k|
      u[k][index].abs
    end
    u[j], u[restore_index] = u[restore_index], u[j]
    bvector[j], bvector[restore_index] = bvector[restore_index], bvector[j]
    restore = (j + 1 .. [index, out_size - 1].min).map do |k|
      v = u[j][index] == 0 ? 0 : u[k][index].quo(u[j][index])
      (index + 1 ... size).each do |l|
        u[k][l] -= v * u[j][l]
      end
      bvector[k] -= v * bvector[j]
      [k, v]
    end
    selected.push index
    solve.call u, selected, index + 1
    selected.pop
    restore.reverse_each do |k, v|
      (index + 1 ... size).each do |l|
        u[k][l] += v * u[j][l]
      end
      bvector[k] += v * bvector[j]
    end
    u[j], u[restore_index] = u[restore_index], u[j]
    bvector[j], bvector[restore_index] = bvector[restore_index], bvector[j]
    solve.call u, selected, index + 1 if size - index > out_size - selected.size
  end
  if first_required
    solve.call mat_u, [0], 1
  else
    solve.call mat_u, [], 0
  end
end
select_solve(vectors, bvector, max_items, first_required, &block) click to toggle source
# File lib/unodos/solver.rb, line 80
def self.select_solve(vectors, bvector, max_items, first_required, &block)
  if first_required && max_items == 1
    a = match_vector vectors[0], bvector
    block.call [a], [0] if a
  elsif vectors.size < bvector.size
    least_square_solve vectors, bvector, &block
  elsif first_required && max_items == 2
    find_solve vectors, bvector, &block
  else
    recursive_solve vectors, bvector, first_required, &block
  end
end
solve(list, differential_level, min_cost) click to toggle source
# File lib/unodos/solver.rb, line 34
def self.solve(list, differential_level, min_cost)
  base_cost = differential_level
  max_items = min_cost - base_cost - 1
  return nil if max_items <= 0
  result = nil
  vector_max_bases = NAMED_BASES.map do |base|
    vector = (differential_level..list.size-1).map do |i|
      base.proc.call(i)
    end
    [vector, vector.map(&:abs).max, base]
  end
  if differential_level > 0
    (1..differential_level).each do |level|
      vector = list.take(list.size - level).drop(differential_level - level)
      vector_max_bases.unshift [vector, vector.map(&:abs).max, Unodos::DifferentialBase.new(level)]
    end
    list = list.drop differential_level
  end
  select_solve vector_max_bases.map(&:first), list, max_items, differential_level > 0 do |vs, pos|
    rs = vector_max_bases.values_at(*pos).zip(vs)
    rs.each { |r| r[1] = 0 if (r[0][1] * r[1]).abs < TOLERANCE }
    next if differential_level > 0 && rs[0][1] == 0
    cost = rs.sum { |(_, _, base), v| v == 0 ? 0 : 1 } + base_cost
    if cost < min_cost
      min_cost = cost
      result = rs
    end
  end
  [min_cost, result] if result
end