[0] = { A, b, c, B: step.B }
-
- const A_B = A.slice({ rows: step.B })
- const A_B_inverse = A_B.inverse2x2()
- const b_B = b.slice(step.B)
-
- rows.push(
-
-
-
- )
-
- const x = A_B_inverse.apply(b_B)
- canvasOptions.x = x
-
- rows.push(
-
-
-
- )
-
- const y_B = A_B_inverse.transpose().apply(c)
-
- rows.push(
-
-
-
- )
-
- const y_Zero = Vector.zero(RationalField, A.rows)
-
- const y = y_Zero.with(step.B, y_B)
-
- rows.push(
-
-
-
- )
-
- const I_x = activeIndices(input, x)
-
- rows.push(
-
-
-
- )
-
- const isDegenerate = I_x.length < A_B.rows
- const isDualAdmissible = y_B.getData().every(y => y.geq(RationalField.zero))
- const isDualDegenerate = y_B.getData().some(y => y.eq(RationalField.zero))
-
- rows.push(
-
-
- La soluzione primale è {isDegenerate ? 'degenere' : 'non degenere'}.
-
-
- La soluzione duale è {isDualAdmissible ? 'ammissibile' : 'non ammissibile'} e{' '}
- {isDualDegenerate ? 'degenere' : 'non degenere'}.
-
-
- )
-
- if (!isDualAdmissible) {
- const h = Math.min(...y.getData().flatMap((y, i) => (y.lt(RationalField.zero) ? [i] : [])))
-
- rows.push(
-
-
-
- )
-
- const e_h = Vector.oneHot(RationalField, A.rows, h).slice(step.B)
- console.log(e_h)
-
- // const xi = Vec.neg(Mat.apply(A_B_inverse, e_h))
- const xi = A_B_inverse.apply(e_h).neg()
-
- rows.push(
-
-
-
- )
-
- const N = range(0, A.rows).filter(i => !step.B.includes(i))
- const A_N = A.slice({ rows: N })
-
- const A_N__xi = A_N.apply(xi)
-
- rows.push(
-
-
-
- )
-
- if (!A_N__xi.getData().every(x => x.leq(RationalField.zero))) {
- const positiveIndices = N.filter(i => A_N__xi.at(A_N.forwardRowIndices[i]).gt(RationalField.zero))
-
- const [k, lambda] = positiveIndices
- .map<[number, Rational]>(i => [i, b.at(i).sub(A.rowAt(i).dot(x)).div(A.rowAt(i).dot(xi))])
- .reduce(([i1, lambda1], [i2, lambda2]) => (lambda2.lt(lambda1) ? [i2, lambda2] : [i1, lambda1]))
-
- rows.push(
-
- 0 \\right\\}`,
- `${lambda}`,
- ].join(' = ')}
- />
-
- )
- rows.push(
-
- 0 \\right\\}`,
- `${k + 1}`,
- ].join(' = ')}
- />
-
- )
- rows.push(
-
- i !== h), k].toSorted())}`,
- ].join(' = ')}
- />
-
- )
- } else {
- rows.push(
-
-
- La soluzione duale è illimitata.
-
-
+ g.lineWidth = 30 / g.canvas.width
+ g.strokeStyle = 'darkgreen'
+ drawSimpleArrow(
+ g,
+ x1.toNumber(),
+ x2.toNumber(),
+ x1.toNumber() + c1.toNumber() / cLen,
+ x2.toNumber() + c2.toNumber() / cLen,
+ 0.125
)
+
+ // draw xi
+ if (xi) {
+ const [xi1, xi2] = xi.getData()
+ const xiLen = Math.sqrt(xi1.toNumber() ** 2 + xi2.toNumber() ** 2)
+
+ g.strokeStyle = '#44d'
+ g.lineWidth = 50 / g.canvas.width
+ drawSimpleArrow(
+ g,
+ x1.toNumber(),
+ x2.toNumber(),
+ x1.toNumber() + xi1.toNumber() / xiLen,
+ x2.toNumber() + xi2.toNumber() / xiLen,
+ 0.125
+ )
+ }
+
+ g.fillStyle = '#d44'
+ fillDot(g, x1.toNumber(), x2.toNumber(), 0.2)
}
}
- return (
-
- )
+ return
}
+
+// export const PrimaleStep = ({ input, step }: { input: ProblemInput; step: Step }) => {
+// const { A, b, c } = input
+
+// const rows = []
+// const canvasOptions: Parameters[0] = { A, b, c, B: step.B }
+
+// const A_B = A.slice({ rows: step.B })
+// const A_B_inverse = A_B.inverse2x2()
+// const b_B = b.slice(step.B)
+
+// rows.push(
+//
+//
+//
+// )
+
+// const x = A_B_inverse.apply(b_B)
+// canvasOptions.x = x
+
+// rows.push(
+//
+//
+//
+// )
+
+// const y_B = A_B_inverse.transpose().apply(c)
+
+// rows.push(
+//
+//
+//
+// )
+
+// const y_Zero = Vector.zero(RationalField, A.rows)
+
+// const y = y_Zero.with(step.B, y_B)
+
+// rows.push(
+//
+//
+//
+// )
+
+// const I_x = activeIndices(input, x)
+
+// rows.push(
+//
+//
+//
+// )
+
+// const isDegenerate = I_x.length < A_B.rows
+// const isDualAdmissible = y_B.getData().every(y => y.geq(RationalField.zero))
+// const isDualDegenerate = y_B.getData().some(y => y.eq(RationalField.zero))
+
+// rows.push(
+//
+//
+// La soluzione primale è {isDegenerate ? 'degenere' : 'non degenere'}.
+//
+//
+// La soluzione duale è {isDualAdmissible ? 'ammissibile' : 'non ammissibile'} e{' '}
+// {isDualDegenerate ? 'degenere' : 'non degenere'}.
+//
+//
+// )
+
+// if (!isDualAdmissible) {
+// const h = Math.min(...y.getData().flatMap((y, i) => (y.lt(RationalField.zero) ? [i] : [])))
+
+// rows.push(
+//
+//
+//
+// )
+
+// const e_h = Vector.oneHot(RationalField, A.rows, h).slice(step.B)
+// console.log(e_h)
+
+// // const xi = Vec.neg(Mat.apply(A_B_inverse, e_h))
+// const xi = A_B_inverse.apply(e_h).neg()
+
+// rows.push(
+//
+//
+//
+// )
+
+// const N = range(0, A.rows).filter(i => !step.B.includes(i))
+// const A_N = A.slice({ rows: N })
+
+// const A_N__xi = A_N.apply(xi)
+
+// rows.push(
+//
+//
+//
+// )
+
+// if (!A_N__xi.getData().every(x => x.leq(RationalField.zero))) {
+// const positiveIndices = N.filter(i => A_N__xi.at(A_N.forwardRowIndices[i]).gt(RationalField.zero))
+
+// const [k, lambda] = positiveIndices
+// .map<[number, Rational]>(i => [i, b.at(i).sub(A.rowAt(i).dot(x)).div(A.rowAt(i).dot(xi))])
+// .reduce(([i1, lambda1], [i2, lambda2]) => (lambda2.lt(lambda1) ? [i2, lambda2] : [i1, lambda1]))
+
+// rows.push(
+//
+// 0 \\right\\}`,
+// `${lambda}`,
+// ].join(' = ')}
+// />
+//
+// )
+// rows.push(
+//
+// 0 \\right\\}`,
+// `${k + 1}`,
+// ].join(' = ')}
+// />
+//
+// )
+// rows.push(
+//
+// i !== h), k].toSorted())}`,
+// ].join(' = ')}
+// />
+//
+// )
+// } else {
+// rows.push(
+//
+//
+// La soluzione duale è illimitata.
+//
+//
+// )
+// }
+// }
+
+// return (
+//
+// )
+// }
diff --git a/src/lib-v2/canvas/index.ts b/src/lib-v2/canvas/index.ts
index 02c916d..2c7e1d3 100644
--- a/src/lib-v2/canvas/index.ts
+++ b/src/lib-v2/canvas/index.ts
@@ -30,7 +30,7 @@ export function drawSemiplane(
p2 = b / a2
}
- const normalize = Math.sqrt(a1 ** 2 + a2 ** 2) * 1.5
+ const normalize = Math.sqrt(a1 ** 2 + a2 ** 2) * 1.25
const gradient = g.createLinearGradient(p1, p2, p1 - a1 / normalize, p2 - a2 / normalize)
gradient.addColorStop(0, gradientAccent)
@@ -81,3 +81,40 @@ export function drawSemiplane(
}
g.stroke()
}
+
+export function drawSimpleArrow(
+ g: CanvasRenderingContext2D,
+ x1: number,
+ y1: number,
+ x2: number,
+ y2: number,
+ size: number
+) {
+ const arrowLength = Math.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2)
+
+ g.save()
+ g.beginPath()
+ g.translate(x1, y1)
+ g.rotate(Math.atan2(y2 - y1, x2 - x1))
+ g.moveTo(0, 0)
+ g.lineTo(arrowLength, 0)
+
+ g.moveTo(arrowLength - size, -size)
+ g.lineTo(arrowLength, 0)
+ g.moveTo(arrowLength - size, +size)
+ g.lineTo(arrowLength, 0)
+ g.stroke()
+ g.restore()
+}
+
+export function fillDot(g: CanvasRenderingContext2D, x: number, y: number, radius: number) {
+ g.beginPath()
+ g.arc(x, y, radius, 0, 2 * Math.PI)
+ g.fill()
+}
+
+export function strokeDot(g: CanvasRenderingContext2D, x: number, y: number, radius: number) {
+ g.beginPath()
+ g.arc(x, y, radius, 0, 2 * Math.PI)
+ g.stroke()
+}
diff --git a/src/lib-v2/ro/primal-simplex.ts b/src/lib-v2/ro/primal-simplex.ts
new file mode 100644
index 0000000..364fae6
--- /dev/null
+++ b/src/lib-v2/ro/primal-simplex.ts
@@ -0,0 +1,241 @@
+import { indexSetToLatex, matrixToLatex, rowVectorToLatex, vectorToLatex } from '../latex'
+import { range } from '../math'
+import { Matrix } from '../matrix'
+import { Rational, RationalField } from '../rationals'
+import { Vector } from '../vector'
+
+export type ProblemInput = {
+ A: Matrix
+ b: Vector
+ c: Vector
+
+ B: number[]
+
+ maxIterations?: number
+}
+
+export type ProblemComment = { type: 'formula' | 'text'; content: string }
+
+export type ProblemStep = {
+ B: number[]
+
+ x?: Vector
+ xi?: Vector
+
+ comments: ProblemComment[]
+}
+
+export type ProblemStatus =
+ | { result: 'optimal'; x: Vector }
+ | { result: 'unbounded'; xi: Vector }
+ | { result: 'infeasible' }
+ | { result: 'too-many-iterations' }
+
+export type ProblemOutput = {
+ status: ProblemStatus
+ steps: ProblemStep[]
+}
+
+const activeIndices = (input: ProblemInput, x: Vector): number[] => {
+ const { A, b } = input
+ const A_x = A.apply(x)
+
+ return A_x.getData().flatMap((a, i) => (a.eq(b.at(i)) ? [i] : []))
+}
+
+export function computePrimalSimplexSteps(input: ProblemInput): ProblemOutput {
+ input.maxIterations ??= 10
+
+ const { A, b, c } = input
+
+ let B = input.B
+ let stepResult_B: number[] | undefined = undefined
+ let stepResult_xi: Vector | undefined = undefined
+
+ let status: ProblemStatus | null = null
+
+ const steps: ProblemStep[] = []
+
+ while (status === null && steps.length < input.maxIterations) {
+ const comments: ProblemComment[] = []
+
+ const A_B = A.slice({ rows: B })
+ const A_B_inverse = A_B.inverse2x2()
+ const b_B = b.slice(B)
+
+ comments.push({
+ type: 'formula',
+ content: [
+ `A_B = ${matrixToLatex(A_B)}`,
+ `b_B = ${vectorToLatex(b_B)}`,
+ `c^t = ${rowVectorToLatex(c)}`,
+ ].join(' \\qquad '),
+ })
+
+ const x = A_B_inverse.apply(b_B)
+
+ comments.push({
+ type: 'formula',
+ content: [
+ String.raw`\bar{x}`,
+ `A_B^{-1} b_B`,
+ `${matrixToLatex(A_B)}^{-1} ${vectorToLatex(b_B)}`,
+ `${matrixToLatex(A_B_inverse)} ${vectorToLatex(b_B)}`,
+ `${vectorToLatex(x)}`,
+ ].join(' = '),
+ })
+
+ const y_B = A_B_inverse.transpose().apply(c)
+
+ comments.push({
+ type: 'formula',
+ content: [
+ String.raw`\bar{y}_B^t`,
+ `c_B^t A_B^{-1}`,
+ `${rowVectorToLatex(c)} ${matrixToLatex(A_B_inverse)}`,
+ `${rowVectorToLatex(y_B)}`,
+ ].join(' = '),
+ })
+
+ const y_Zero = Vector.zero(RationalField, A.rows)
+
+ const y = y_Zero.with(B, y_B)
+
+ comments.push({ type: 'formula', content: String.raw`\implies \bar{y}^t = ${rowVectorToLatex(y)}` })
+
+ const I_x = activeIndices(input, x)
+
+ comments.push({
+ type: 'formula',
+ content: [
+ String.raw`I(\bar{x})`,
+ String.raw`\{ i \in \{1, \dots, m\} \mid A_i \bar{x}_i = b_i \}`,
+ indexSetToLatex(I_x),
+ ].join(' = '),
+ })
+
+ const isDegenerate = I_x.length < A_B.rows
+ const isDualAdmissible = y_B.getData().every(y => y.geq(RationalField.zero))
+ const isDualDegenerate = y_B.getData().some(y => y.eq(RationalField.zero))
+
+ comments.push({
+ type: 'text',
+ content: `La soluzione primale è *${isDegenerate ? 'degenere' : 'non degenere'}*.`,
+ })
+
+ comments.push({
+ type: 'text',
+ content: `La soluzione duale è *${isDualAdmissible ? 'ammissibile' : 'non ammissibile'}* e *${
+ isDualDegenerate ? 'degenere' : 'non degenere'
+ }*.`,
+ })
+
+ if (!isDualAdmissible) {
+ const h = Math.min(...y.getData().flatMap((y, i) => (y.lt(RationalField.zero) ? [i] : [])))
+
+ comments.push({
+ type: 'formula',
+ content: [
+ // prettier-ignore
+ `h`,
+ String.raw`\min \{ i \in B \mid \bar{y}_i < 0 \}`,
+ `${h + 1}`,
+ ].join(' = '),
+ })
+
+ const e_h = Vector.oneHot(RationalField, A.rows, h).slice(B)
+ console.log(e_h)
+
+ const xi = A_B_inverse.apply(e_h).neg()
+ stepResult_xi = xi
+
+ comments.push({
+ type: 'formula',
+ content: [
+ // prettier-ignore
+ `\\xi`,
+ `-A_B^{-1} u_{B(h)}`,
+ `${vectorToLatex(xi)}`,
+ ].join(' = '),
+ })
+
+ const N = range(0, A.rows).filter(i => !B.includes(i))
+ const A_N = A.slice({ rows: N })
+
+ const A_N__xi = A_N.apply(xi)
+
+ comments.push({
+ type: 'formula',
+ content: [
+ [`N = \\{1, \\dots, m\\} \\setminus B = ${indexSetToLatex(N)}`],
+ [`A_N \\xi`, `${matrixToLatex(A_N)} ${vectorToLatex(xi)}`, `${vectorToLatex(A_N__xi)}`].join(' = '),
+ ].join(' \\qquad '),
+ })
+
+ if (!A_N__xi.getData().every(x => x.leq(RationalField.zero))) {
+ const [k, lambda] = N.filter(i => A_N__xi.at(A_N.forwardRowIndices[i]).gt(RationalField.zero))
+ .map<[number, Rational]>(i => [i, b.at(i).sub(A.rowAt(i).dot(x)).div(A.rowAt(i).dot(xi))])
+ .reduce(([i1, lambda1], [i2, lambda2]) => (lambda1.lt(lambda2) ? [i1, lambda1] : [i2, lambda2]))
+
+ comments.push({
+ type: 'formula',
+ content: [
+ `\\bar\\lambda`,
+ `\\min_i \\left\\{ \\frac{b_i - A_i \\bar{x}}{A_i \\xi} \\; \\middle| \\; i \\in N, A_i \\xi > 0 \\right\\}`,
+ `${lambda}`,
+ ].join(' = '),
+ })
+
+ comments.push({
+ type: 'formula',
+ content: [
+ `k`,
+ `\\argmin_i \\left\\{ \\bar\\lambda_i \\; \\middle| \\; i \\in N, A_i \\xi > 0 \\right\\}`,
+ `${k + 1}`,
+ ].join(' = '),
+ })
+
+ comments.push({
+ type: 'formula',
+ content: [
+ `B'`,
+ `B \\setminus \\{h + 1\\} \\cup \\{k + 1\\}`,
+ `${indexSetToLatex([...B.filter(i => i !== h), k].toSorted())}`,
+ ].join(' = '),
+ })
+
+ stepResult_B = [...B.filter(i => i !== h), k].toSorted()
+ } else {
+ comments.push({ type: 'text', content: 'La soluzione è *illimitata*' })
+ status = { result: 'unbounded', xi }
+ }
+ } else {
+ comments.push({ type: 'text', content: 'La soluzione è *ottima*' })
+ status = { result: 'optimal', x }
+ }
+
+ steps.push({
+ B,
+ x,
+ xi: stepResult_xi,
+
+ comments,
+ })
+
+ if (stepResult_B) {
+ B = stepResult_B
+
+ stepResult_B = undefined
+ stepResult_xi = undefined
+ }
+ }
+
+ if (status === null) {
+ status = { result: 'too-many-iterations' }
+ }
+
+ return {
+ status,
+ steps,
+ }
+}
diff --git a/src/main.tsx b/src/main.tsx
index 326cceb..08136df 100644
--- a/src/main.tsx
+++ b/src/main.tsx
@@ -54,11 +54,6 @@ const App = () => {
Svolgimento
{'result' in problemValuesResult && }
-
- Debug
-
- {JSON.stringify(problemValuesResult, null, 4)}
-
>
)
}
diff --git a/src/style.css b/src/style.css
index 0a2427a..dd48cc8 100644
--- a/src/style.css
+++ b/src/style.css
@@ -9,9 +9,12 @@
box-sizing: border-box;
}
+ html {
+ height: 100%;
+ }
+
html,
body {
- height: 100%;
min-height: 100%;
}
@@ -45,6 +48,10 @@
color: #333;
}
+ p {
+ line-height: 1.75;
+ }
+
pre,
code {
font-family: 'JetBrains Mono', monospace;
@@ -66,7 +73,7 @@
color: #333;
- padding: 3rem 1rem;
+ padding: 6rem 1rem;
> * {
display: block;
@@ -91,6 +98,7 @@
display: grid;
grid-auto-flow: row;
justify-items: center;
+ grid-template-columns: 1fr 1fr;
border: 1px solid #ddd;
box-shadow: 0 0 1rem rgba(0, 0, 0, 0.1);
@@ -101,10 +109,13 @@
max-width: fit-content;
> .step {
+ grid-column: span 2;
+
display: grid;
- grid-template-columns: 1fr 1fr;
+ grid-template-columns: subgrid;
+
justify-items: center;
- align-items: center;
+ /* align-items: center; */
gap: 0;
@media (width < 768px) {
@@ -112,17 +123,29 @@
}
> .algebraic-step {
+ justify-self: stretch;
+
/* border: 1px solid #ddd; */
/* box-shadow: 0 0 1rem rgba(0, 0, 0, 0.1); */
border-radius: 1rem;
padding: 0 1rem;
+
+ > * + * {
+ margin-top: 0.5rem;
+ }
}
> .geometric-step {
/* padding: 1rem; */
min-width: 30rem;
}
+
+ padding: 2rem 0;
+
+ &:not(:first-child) {
+ border-top: 1px solid #ddd;
+ }
}
}
}