core API

src にあるパッケージが linear-program の公開インターフェースのすべてです。インターフェースファイルは src/pkg.generated.mbti です。パッケージパスは linear-program です。モジュール名にはまだ Luna-Flow/ 名前空間が付いていないため、他のパッケージからは "linear-program" としてインポートします。

係数型 V はジェネリックです。各関数は必要な trait を明示しており、Double はそのすべてを満たします。ソルバーがテストされている型は Double だけです。行列は linear-algebra の mutable パッケージの @mutable.Matrix[V] です。計画行列またはシンプレックス表は、第 0 行に目的関数を、以降の各行に等式制約を 1 つずつ格納し、最終列に右辺を格納します。

このページの例は、次の moon.pkg を持つパッケージ内のテストです。

import {
  "linear-program" @lp,
  "Luna-Flow/linear-algebra/mutable" @la,
  "moonbitlang/core/double",
}

また、次の宣言を共有します。

using @lp {type Variable, type Poly, type Obj_func, type Constraint, type Lp}

fn nonneg(names : Array[String]) -> Array[Variable] {
  names.map(n => Variable::new(n, 0.0, @double.max_value))
}

変数

`Variable`

Variable は下限と上限を持つ名前付きの決定変数です。

type Variable
pub impl @luna-generic.One for Variable
pub impl Compare for Variable
pub impl Eq for Variable
pub impl Show for Variable

Variable は抽象型です。値は Variable::new で作成し、フィールドを直接読むことはできません。2 つの変数は境界に関係なく名前が等しければ等しく、名前順に並びます。境界は保存・表示されますが、ソルバーは読みません。計画のすべての変数は非負として扱われます。

`Variable::new`

Variable::new は名前、下界、上界から変数を作成します。

pub fn Variable::new(String, Double, Double) -> Self

何も制限しない上界には @double.max_value を使います。Lp::to_standard は追加する変数に y1、y2、… と名前を付けるため、自分の変数にはこれらの名前を使わないでください。

`Variable::show_all`

show_all は変数名を境界とともに表示します。

pub fn Variable::show_all(Self) -> String

test "variables" {
  let x = Variable::new("x1", 0.0, 100.0)
  inspect(x, content="x1")
  inspect(x.show_all(), content="x1 : [0, 100]")
  assert_true(x == Variable::new("x1", -5.0, 5.0))
  assert_true(x < Variable::new("x2", 0.0, 1.0))
}

`Variable::to_string`, `Variable::equal`, `Variable::compare`

これらは Show、Eq、Compare のメソッドで、Variable のメソッドとして昇格されたものです。

pub fn Variable::to_string(Self) -> String
pub fn Variable::equal(Self, Self) -> Bool
pub fn Variable::compare(Self, Self) -> Int

to_string は名前を返します。equal と compare は名前だけを見ます。compare は名前の辞書順なので、"x10" は "x2" より前に来ます。新しいコードでは "\{x}"、==、< を使ってください。

`Variable::one`

Variable::one は線形式の定数項を表す無名の変数です。

pub fn Variable::one() -> Self

これは luna-generic の One メソッドをメソッドとして昇格したものです。名前は空で、境界は @double.min_value と @double.max_value です。Poly::from_var_array は対応する変数のない係数をこれに割り当てます。

線形式

`Poly`

Poly[V] は線形式です。変数から非零係数への写像で、変数名の順に保持されます。

type Poly[V] derive(Eq)
pub impl[V : Eq + @luna-generic.Semiring + @luna-generic.Zero + Neg] Neg for Poly[V]
pub impl[V : Eq + Show + @luna-generic.Semiring + @luna-generic.One] Show for Poly[V]

Poly は可変です。add_term_inplace はそれを変更し、共有しているすべての値に変更が見えます。独立した式が必要なときは copy を使います。2 つの式は、同じ変数を等しい係数で持つときに等しくなります。

`Poly::new`, `Poly::from_var_array`

Poly::new は空の式を作成し、Poly::from_var_array は係数から式を作成します。

pub fn[V] Poly::new() -> Self[V]
pub fn[V : Eq + @luna-generic.Semiring] Poly::from_var_array(Array[V], Array[Variable]) -> Self[V]

Poly::from_var_array(coeffs, vars) は coeffs[i] と vars[i] を組にし、ゼロ係数を飛ばします。最後の変数より後ろの位置にある係数は定数項 Variable::one() に割り当てられ、そのような係数が複数あるときは最後のものが採用されます。

`Poly::add_term_inplace`

add_term_inplace は式に変数の倍数をその場で加えます。

pub fn[V : Eq + @luna-generic.Semiring + @luna-generic.Zero] Poly::add_term_inplace(Self[V], Variable, V) -> Unit

p.add_term_inplace(x, c) は x の係数に c を加えます。ゼロを加えても何も起こらず、係数がゼロになった項は削除されるため、ゼロ係数が保存されることはありません。

`Poly::copy`

copy は、同じ項を持ち元の式と記憶領域を共有しない式を返します。

pub fn[V : Eq + @luna-generic.Semiring + @luna-generic.Zero] Poly::copy(Self[V]) -> Self[V]

`Poly::to_array`, `Poly::to_vector`

to_array は保存された係数を列挙し、to_vector は与えられた変数ごとに係数を 1 つ列挙します。

pub fn[V] Poly::to_array(Self[V]) -> Array[V]
pub fn[V : @luna-generic.Zero] Poly::to_vector(Self[V], Array[Variable]) -> Array[V]

to_array は非ゼロ係数を変数名の順に、変数なしで返します。to_vector(vars) は vars の順に vars.length() 個の係数を返し、現れない変数にはゼロを入れます。これが計画行列で使われる密な行です。

`Poly::neg`, `Poly::equal`, `Poly::to_string`

これらは Neg、Eq、Show のメソッドで、Poly のメソッドとして昇格されたものです。

pub fn[V : Eq + @luna-generic.Semiring + @luna-generic.Zero + Neg] Poly::neg(Self[V]) -> Self[V]
pub fn[V : Eq] Poly::equal(Self[V], Self[V]) -> Bool
pub fn[V : Eq + Show + @luna-generic.Semiring + @luna-generic.One] Poly::to_string(Self[V]) -> String

-p はすべての係数の符号を反転して新しい式を返します。Show は項を名前順に " + " でつないで表示し、1 に等しい係数は省略します。例: 3.5x1 + x2 + -2x3。

test "linear expressions" {
  let vars = nonneg(["x1", "x2", "x3"])
  let p : Poly[Double] = Poly::from_var_array([3.5, 0.0, -2.0], vars)
  inspect(p, content="3.5x1 + -2x3")
  p.add_term_inplace(vars[1], 1.0)
  p.add_term_inplace(vars[2], 2.0)
  inspect(p, content="3.5x1 + x2")
  inspect(-p, content="-3.5x1 + -1x2")
  debug_inspect(p.to_array(), content="[3.5, 1]")
  debug_inspect(p.to_vector(vars), content="[3.5, 1, 0]")
  assert_true(p.copy() == p)
}

目的関数

`Obj_func`

Obj_func[V] は目的関数です。線形式と、最小化か最大化かの向きからなります。

type Obj_func[V]
pub impl[V : Eq + Show + @luna-generic.Semiring] Show for Obj_func[V]

Show は向きと式を表示します。例: max 3x1 + 4x2。Obj_func::to_string は昇格された Show メソッドです。

pub fn[V : Eq + Show + @luna-generic.Semiring] Obj_func::to_string(Self[V]) -> String

`Obj_func::new`, `Obj_func::from_array`

Obj_func::new は空の式を持つ目的関数を作成し、Obj_func::from_array は式も設定します。

pub fn[V] Obj_func::new(String) -> Self[V]
pub fn[V : Eq + @luna-generic.Semiring] Obj_func::from_array(String, Array[V], Array[Variable]) -> Self[V]

向きは、最小化なら "min"、"Min"、"MIN"、"minimize"、"Minimize" のいずれか、最大化なら "max"、"Max"、"MAX"、"maximize"、"Maximize" のいずれかです。それ以外の文字列は abort します。Obj_func::from_array(sense, coeffs, vars) は Poly::from_var_array で式を作ります。

`Obj_func::set_poly`, `Obj_func::add_term_inplace`

set_poly は式を置き換え、add_term_inplace はその場で項を追加します。

pub fn[V] Obj_func::set_poly(Self[V], Poly[V]) -> Self[V]
pub fn[V : Eq + @luna-generic.Semiring + @luna-generic.Zero] Obj_func::add_term_inplace(Self[V], Variable, V) -> Unit

set_poly(p) は同じ向きで式が p の新しい目的関数を返します。p はコピーしません。

`Obj_func::judge_max`, `Obj_func::to_min`

judge_max は目的関数が最大化かどうかを判定し、to_min は等価な最小化の目的関数を返します。

pub fn[V] Obj_func::judge_max(Self[V]) -> Bool
pub fn[V : Eq + @luna-generic.Semiring + @luna-generic.Zero + Neg] Obj_func::to_min(Self[V]) -> Self[V]

to_min は最小化の目的関数をそのまま返し、最大化の目的関数では式の符号を反転します。 $\max c^\top x = -\min\,(-c)^\top x$ だからです。

`Obj_func::to_vector`

to_vector は計画行列の目的関数行を返します。

pub fn[V : @luna-generic.Zero] Obj_func::to_vector(Self[V], Int) -> Array[V]

to_vector(n) は to_array の係数を取り、シンプレックス表の 1 行の幅である長さ n + 1 までゼロで埋めます。to_array は非ゼロ係数だけを名前順に列挙するため、結果が変数と揃うのは、すべての変数の係数が非ゼロで、変数が名前順に宣言されている場合だけです。下の Lp::two_stage の警告を参照してください。

test "objectives" {
  let vars = nonneg(["x1", "x2"])
  let obj : Obj_func[Double] = Obj_func::from_array("max", [3.0, 4.0], vars)
  inspect(obj, content="max  3x1 + 4x2")
  assert_true(obj.judge_max())
  inspect(obj.to_min(), content="min  -3x1 + -4x2")
  debug_inspect(obj.to_vector(2), content="[3, 4, 0]")
}

制約

`Constraint`

Constraint[V] は問題の制約の集合です。等式 poly = b と不等式 poly <= b または poly >= b からなります。

type Constraint[V]
pub impl[V : Show + Eq + @luna-generic.Semiring] Show for Constraint[V]
pub fn[V : Show + Eq + @luna-generic.Semiring] Constraint::to_string(Self[V]) -> String

Show は s.t. の後に制約集合を 1 行に 1 つずつ、等式を先にして表示します。計画は Lp::cons_from_Array で制約を受け取ります。単独の Constraint 値は制約集合の構築と表示に便利です。

`Constraint::new`, `Constraint::from_array`

Constraint::new は空の集合を作成し、Constraint::from_array は係数の行から集合を構築します。

pub fn[V] Constraint::new() -> Self[V]
pub fn[V : Eq + @luna-generic.Semiring] Constraint::from_array(eq_array? : Array[(Array[V], V)], ineq_array? : Array[(Array[V], String, V)], Array[Variable]) -> Self[V]

等式の行は (coeffs, b)、不等式の行は (coeffs, relation, b) で、relation は正確に ">=" または "<=" です。それ以外の関係は abort します。どちらの配列も既定値は空です。

`Constraint::add_eqpoly`, `Constraint::add_ineqpoly`

add_eqpoly は等式を、add_ineqpoly は不等式を追加します。

pub fn[V] Constraint::add_eqpoly(Self[V], Poly[V], V) -> Unit
pub fn[V] Constraint::add_ineqpoly(Self[V], Poly[V], String, V) -> Unit

add_ineqpoly は関係の文字列を検査しません。">=" か "<=" を渡してください。

`Constraint::change_eqpoly`, `Constraint::change_ineqpoly`

これらは 1 から数えた位置にある制約を置き換えます。

pub fn[V] Constraint::change_eqpoly(Self[V], Int, Poly[V], V) -> Unit
pub fn[V] Constraint::change_ineqpoly(Self[V], Int, Poly[V], String, V) -> Unit

change_eqpoly(i, poly, b) は 1 から数えて i 番目の等式を置き換えます。i が正でないとき、または等式配列の容量（配列の実装上の詳細）未満でないとき、どちらの関数も abort します。

`Constraint::to_matrix`

to_matrix は等式ごとに密な行 [coefficients..., b] を 1 つ返します。

pub fn[V : @luna-generic.Zero] Constraint::to_matrix(Self[V], Array[Variable]) -> Array[Array[V]]

列は与えられた変数の順に並びます。不等式は除外されるため、標準形の計画の制約に対して呼び出してください。

test "constraints" {
  let vars = nonneg(["x1", "x2"])
  let c : Constraint[Double] = Constraint::from_array(
    eq_array=[([1.0, 1.0], 4.0)],
    ineq_array=[([2.0, 1.0], "<=", 6.0)],
    vars,
  )
  inspect(
    c,
    content=(
      #|s.t. x1 + x2 = 4
      #|     2x1 + x2 <= 6
      #|    
    ),
  )
  debug_inspect(c.to_matrix(vars), content="[[1, 1, 4]]")
}

線形計画問題

`Lp`

Lp[V] は線形計画です。変数、目的関数、制約、およびそれらから作られた計画行列からなります。

type Lp[V]
pub impl[V : Show + Eq + @luna-generic.Semiring] Show for Lp[V]
pub fn[V : Show + Eq + @luna-generic.Semiring] Lp::to_string(Self[V]) -> String

計画を作る関数は新しい Lp を返し、引数は変更しません。Show は目的関数、制約、および各変数とその境界を表示します。

`Lp::new`, `Lp::cons_from_Array`, `Lp::reset_obj_byarray`

Lp::new は制約のない計画を作成します。残りの 2 つはその制約または目的関数の係数を置き換えます。

pub fn[V : @luna-generic.Zero] Lp::new(Obj_func[V], Array[Variable]) -> Self[V]
pub fn[V : Eq + @luna-generic.Semiring] Lp::cons_from_Array(Self[V], eq_array? : Array[(Array[V], V)], ineq_array? : Array[(Array[V], String, V)]) -> Self[V]
pub fn[V : Compare + @luna-generic.Semiring + @luna-generic.Zero] Lp::reset_obj_byarray(Self[V], Array[V]) -> Self[V]

cons_from_Array は計画の変数を使って Constraint::from_array と同様にすべての制約を置き換えます。reset_obj_byarray(coeffs) は向きを保ったまま新しい目的関数の係数を設定します。

`Lp::from_matrix`

Lp::from_matrix は計画行列から等式制約の計画を作成します。

pub fn[V : Eq + @luna-generic.Semiring] Lp::from_matrix(@mutable.Matrix[V], Array[Variable], String) -> Self[V]

Lp::from_matrix(m, vars, sense) は第 0 行（最後の要素を除く）から目的関数を、以降の各行 $[a_i \mid b_i]$ から等式 $a_i \cdot x = b_i$ を読み取ります。計画は m 自体をコピーせずに自分の行列として保持します。

`Lp::to_standard`

to_standard は標準形の等価な計画を返します。

pub fn[V : Eq + @luna-generic.Semiring + Compare + Neg + @luna-generic.Zero] Lp::to_standard(Self[V]) -> Self[V]

結果は最小化で、等式制約だけを持ち、右辺は非負です。

最大化の目的関数は符号を反転します（Obj_func::to_min）。
$b < 0$ の等式には $-1$ を掛けます。
各不等式には新しい変数 y1、y2、…（不等式の順に数える）が与えられます。必要なら $b \ge 0$ となるよう行に $-1$ を掛けたうえで、<= の行には係数 $+1$ 、>= の行には係数 $-1$ で加えます。

新しい変数は計画の変数の後ろに追加されます。各場合の導出は設計ノートにあります。

`Lp::get_coeff_matrix`, `Lp::get_objfunc_vector`, `Lp::get_b_vector`

これらは計画行列から $A$ 、 $c$ 、 $b$ を読み取ります。

pub fn[V] Lp::get_coeff_matrix(Self[V]) -> @mutable.Matrix[V]
pub fn[V] Lp::get_objfunc_vector(Self[V]) -> Array[V]
pub fn[V] Lp::get_b_vector(Self[V]) -> Array[V]

計画行列には目的関数と等式制約しか含まれないため、不等式を持つ計画については to_standard の後でなければ、これらのアクセサは計画を表しません。

test "standard form" {
  let vars = nonneg(["x1", "x2"])
  let lp : Lp[Double] = Lp::new(Obj_func::from_array("max", [3.0, 4.0], vars), vars)
    .cons_from_Array(ineq_array=[([2.0, 1.0], "<=", 12.0), ([1.0, 3.0], ">=", 2.0)])
  let s = lp.to_standard()
  inspect(
    s,
    content=(
      #|min  -3x1 + -4x2
      #|s.t. 2x1 + x2 + y1 = 12
      #|     x1 + 3x2 + -1y2 = 2
      #|      x1 : [0, 1.7976931348623157e+308] x2 : [0, 1.7976931348623157e+308] y1 : [0, 1.7976931348623157e+308] y2 : [0, 1.7976931348623157e+308]
      #|
    ),
  )
  debug_inspect(s.get_objfunc_vector(), content="[-3, -4, 0, 0]")
  debug_inspect(s.get_b_vector(), content="[12, 2]")
  inspect(
    s.get_coeff_matrix(),
    content=(
      #||2, 1, 1, 0|
      #||1, 3, 0, -1|
    ),
  )
}

求解

`Lp::two_stage`

two_stage は標準形の問題を二段階単体法で解きます。

pub fn[V : @luna-generic.Zero + @luna-generic.One + Compare + Mul + Div + Sub + ApproximatelyZero + Show + Neg] Lp::two_stage(Self[V]) -> (Array[V], V)

to_standard の結果に対して呼び出してください。標準形のすべての変数（追加された y 変数を含む）の値を計画の変数の順に返し、さらに標準形の最適値を返します。標準形は常に最小化なので、最大化の計画では最大値は返り値の符号反転です。

計画が実行不可能（第 1 段階が非ゼロの人工目的で終わる）なとき、および非有界なとき、two_stage は abort します。各シンプレックスの実行は 1000 反復で停止し、終わっていなければ Maximum iterations reached, may not have converged を表示します。その場合も結果はエラーなしで返されます。

test "solve" {
  let vars = nonneg(["x1", "x2"])
  let lp : Lp[Double] = Lp::new(Obj_func::from_array("max", [3.0, 4.0], vars), vars)
    .cons_from_Array(ineq_array=[([2.0, 1.0], "<=", 12.0), ([1.0, 3.0], "<=", 10.0)])
  let (x, value) = lp.to_standard().two_stage()
  debug_inspect(x, content="[5.199999999999999, 1.6000000000000005, 0, 0]")
  inspect(-value, content="22")
}

`Lp::pivot`

pivot は単体表に対して Gauss–Jordan のピボット操作をその場で一回行います。

pub fn[V : @luna-generic.Zero + Compare + Div + Mul + Sub] Lp::pivot(@mutable.Matrix[V], Int, Int) -> Unit

Lp::pivot(t, r, s) は第 r 行をピボット $t_{rs}$ で割り、第 0 行を含む他のすべての行 i から新しい第 r 行の $t_{is}$ 倍を引いて、第 s 列を単位ベクトル $e_r$ にします。位置が行列の外にあるとき、またはピボットがちょうどゼロのときは abort します。

`Lp::simplex_iteration`

simplex_iteration は単体表に対してその場で単体法を実行します。

pub fn[V : @luna-generic.Zero + Compare + Sub + Mul + Div + Show] Lp::simplex_iteration(@mutable.Matrix[V], max_iterations? : Int, debug? : Bool) -> Unit

シンプレックス表は正準形でなければなりません。すなわち、第 0 行は被約費用を持ち最終列は $-z$ 、各基底列は単位ベクトルで、右辺は非負です。各反復では、第 0 行で最も負の要素を持つ列を基底に入れ（同値なら最も左）、 $\bar a_{is} > 0$ の行のうち比 $\bar b_i / \bar a_{is}$ が最小の行を基底から出します（同値なら最も上）。第 0 行に負の要素がなくなると停止し、入る列に正の要素がなければ Problem is unbounded で abort し、max_iterations 回（既定は 1000）でメッセージを表示して停止します。debug=true のときは各ピボットとその後の表を表示します。

test "simplex by hand" {
  // max 3x1 + 4x2, 2x1 + x2 + y1 = 12, x1 + 3x2 + y2 = 10; y1, y2 basic
  let t : @la.Matrix[Double] = @la.Matrix::from_2d_array([
    [-3.0, -4.0, 0.0, 0.0, 0.0],
    [2.0, 1.0, 1.0, 0.0, 12.0],
    [1.0, 3.0, 0.0, 1.0, 10.0],
  ])
  Lp::pivot(t, 2, 1) // x2 enters, y2 leaves
  inspect(t[0][4], content="13.333333333333334")
  Lp::simplex_iteration(t)
  inspect(t[0][4], content="22")
}

`Lp::phase_1`, `Lp::phase_2`

phase_1 と phase_2 は two_stage の前半と後半です。

pub fn[V : @luna-generic.Zero + @luna-generic.One + Compare + Sub + Compare + Div + Mul + Show] Lp::phase_1(@mutable.Matrix[V], Array[Int]) -> @mutable.Matrix[V]
pub fn[V : @luna-generic.Zero + @luna-generic.One + ApproximatelyZero + Div + Sub + Mul + Show + Compare] Lp::phase_2(Array[Variable], Array[V], @mutable.Matrix[V]) -> (Array[V], V)

Lp::phase_1(t, artificial_index) は人工列を追加したシンプレックス表を受け取ります。その第 0 行は各人工列で $1$ 、それ以外で $0$ です。artificial_index[i] は制約 i の人工変数の列で、人工変数がなければ 0 です。人工変数を持つ各制約行を第 0 行から引いて表を正準形にし、simplex_iteration を実行して同じ行列を返します。

Lp::phase_2(vars, c, t) は第 1 段階がゼロに達したことを確認し、先頭 vars.length() 個の変数の列と右辺を新しい表にコピーし、目的関数 c を第 0 行に書き込み、基底列（ApproximatelyZero の範囲で単位ベクトルに等しい列）を第 0 行から消去し、simplex_iteration を実行して、基底解と第 0 行の最後の要素（ $-z$ ）を返します。第 1 段階が非ゼロの目的で終わった場合は W* from Phase1 isn't zero, Lp doesn't have solution で abort します。

計画から第 1 段階の表と artificial_index を作るヘルパーは非公開です。自分で表を作るのでなければ two_stage を使ってください。

数値の許容誤差

`ApproximatelyZero`

ApproximatelyZero は、第 2 段階で値をゼロとみなすかどうかを決めます。

pub trait ApproximatelyZero {
  fn is_zero_eps(Self) -> Bool
}
pub impl ApproximatelyZero for Int
pub impl ApproximatelyZero for Double

Int は 0 に等しいときだけゼロです。Double は $|x| < 10^{-15}$ のときゼロで、これは絶対的なしきい値です。独自の係数型でこの trait を実装すれば、two_stage で使えます。

`double_equal_to_zero`

double_equal_to_zero は is_zero_eps を呼び出します。

pub fn[V : ApproximatelyZero] double_equal_to_zero(V) -> Bool

test "tolerance" {
  assert_true(@lp.double_equal_to_zero(1.5e-17))
  assert_false(@lp.double_equal_to_zero(1.0e-12))
  assert_true(@lp.double_equal_to_zero(0))
}

非推奨

以下の昇格メソッドはソース互換性のために残されています。インターフェースファイルからは隠されており、他のパッケージから呼ぶと警告が出ます。

メソッド	代替
`Variable::not_equal`, `Poly::not_equal`	`a != b`
`Variable::op_lt`, `op_le`, `op_gt`, `op_ge`	`<`, `<=`, `>`, `>=`
`Variable::output`, `Poly::output`, `Obj_func::output`, `Constraint::output`, `Lp::output`	`to_string()` または `"\{x}"`

core API

変数

Variable

Variable::new

Variable::show_all

Variable::to_string, Variable::equal, Variable::compare

Variable::one

線形式

Poly

Poly::new, Poly::from_var_array

Poly::add_term_inplace

Poly::copy

Poly::to_array, Poly::to_vector

Poly::neg, Poly::equal, Poly::to_string

目的関数

Obj_func

Obj_func::new, Obj_func::from_array

Obj_func::set_poly, Obj_func::add_term_inplace

Obj_func::judge_max, Obj_func::to_min

Obj_func::to_vector

制約

Constraint

Constraint::new, Constraint::from_array

Constraint::add_eqpoly, Constraint::add_ineqpoly

Constraint::change_eqpoly, Constraint::change_ineqpoly

Constraint::to_matrix