From 1c119d1250e816f549c916c3acac65d1694bb872 Mon Sep 17 00:00:00 2001 From: Antonio De Lucreziis Date: Mon, 20 Jan 2025 17:32:10 +0100 Subject: [PATCH] wow this actually works --- src/MiniMark.tsx | 14 + src/Primale.tsx | 646 ++++++++++++++++++++------------ src/lib-v2/canvas/index.ts | 39 +- src/lib-v2/ro/primal-simplex.ts | 241 ++++++++++++ src/main.tsx | 5 - src/style.css | 31 +- 6 files changed, 721 insertions(+), 255 deletions(-) create mode 100644 src/MiniMark.tsx create mode 100644 src/lib-v2/ro/primal-simplex.ts diff --git a/src/MiniMark.tsx b/src/MiniMark.tsx new file mode 100644 index 0000000..b1bb660 --- /dev/null +++ b/src/MiniMark.tsx @@ -0,0 +1,14 @@ +export const MiniMark = ({ content }: { content: string }) => { + return ( +

/g, '>') + .replace(/\\n/g, '
') + .replace(/\*(.*?)\*/g, '$1') + .replace(/_(.*?)_/g, '$1'), + }} + /> + ) +} diff --git a/src/Primale.tsx b/src/Primale.tsx index 896b044..71d4828 100644 --- a/src/Primale.tsx +++ b/src/Primale.tsx @@ -1,42 +1,128 @@ import { Katex } from './Katex' -import { drawSemiplane } from './lib-v2/canvas' -import { indexSetToLatex, matrixToLatex, rowVectorToLatex, vectorToLatex } from './lib-v2/latex' +import { fillDot, drawSemiplane, drawSimpleArrow } from './lib-v2/canvas' import { range } from './lib-v2/math' import { Matrix } from './lib-v2/matrix' -import { Rational, RationalField } from './lib-v2/rationals' +import { Rational } from './lib-v2/rationals' +import { ProblemComment, computePrimalSimplexSteps } from './lib-v2/ro/primal-simplex' import { Vector } from './lib-v2/vector' +import { MiniMark } from './MiniMark' import { ProblemInput } from './parser-problem' -type Step = { +// type Step = { +// B: number[] +// } + +// 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] : [])) +// } + +const PrimalStep = ({ + iter, + + A, + b, + c, + B, + + x, + xi, + + comments, +}: { + iter: number + + A: Matrix + b: Vector + c: Vector + B: number[] -} + x?: Vector + xi?: Vector -const activeIndices = (input: ProblemInput, x: Vector): number[] => { - const { A, b } = input - const A_x = A.apply(x) + comments: ProblemComment[] +}) => { + return ( +

+
+
+

+ Iterazione {iter + 1} dell'algoritmo +

+
- return A_x.getData().flatMap((a, i) => (a.eq(b.at(i)) ? [i] : [])) + {comments.map(comment => + comment.type === 'formula' ? ( +
+ +
+ ) : ( +
+ +
+ ) + )} +
+
+ +
+
+ ) } export const Primale = ({ input }: { input: ProblemInput }) => { - const steps: Step[] = [{ B: input.B }] + // const steps: Step[] = [{ B: input.B }] + + const problemOutput = computePrimalSimplexSteps({ + A: input.A, + b: input.b, + c: input.c, + B: input.B, + maxIterations: 10, + }) return (
- {steps.map(step => ( - + {problemOutput.steps.map((step, iter) => ( + ))}
) } -const PrimaleCanvas = ({ +const PrimalCanvas = ({ A, b, c, B, x, + xi, }: { A: Matrix b: Vector @@ -44,6 +130,7 @@ const PrimaleCanvas = ({ B: number[] x?: Vector + xi?: Vector }) => { const render = ($canvas: HTMLCanvasElement | null) => { if (!$canvas) { @@ -67,23 +154,61 @@ const PrimaleCanvas = ({ g.lineCap = 'round' g.lineJoin = 'round' - // draw y axis arrow - g.beginPath() - g.moveTo(width / 2, height / 2) - g.lineTo(width / 2, 5) - g.lineTo(width / 2 - 10, 15) - g.moveTo(width / 2, 5) - g.lineTo(width / 2 + 10, 15) - g.stroke() - - // draw x axis arrow - g.beginPath() - g.moveTo(width / 2, height / 2) - g.lineTo(width - 5, height / 2) - g.lineTo(width - 15, height / 2 - 10) - g.moveTo(width - 5, height / 2) - g.lineTo(width - 15, height / 2 + 10) - g.stroke() + g.fillStyle = '#333' + g.font = 'bold 16px sans-serif' + g.textAlign = 'center' + g.textBaseline = 'middle' + + // // draw y axis arrow + // g.beginPath() + // g.moveTo(width / 2, height / 2) + // g.lineTo(width / 2, 5) + // g.lineTo(width / 2 - 10, 15) + // g.moveTo(width / 2, 5) + // g.lineTo(width / 2 + 10, 15) + // g.stroke() + + // // draw x axis arrow + // g.beginPath() + // g.moveTo(width / 2, height / 2) + // g.lineTo(width - 5, height / 2) + // g.lineTo(width - 15, height / 2 - 10) + // g.moveTo(width - 5, height / 2) + // g.lineTo(width - 15, height / 2 + 10) + // g.stroke() + + // draw c vector + const [c1, c2] = c.getData() + const cLen = Math.sqrt(c1.toNumber() ** 2 + c2.toNumber() ** 2) + + // g.save() + g.strokeStyle = 'darkgreen' + g.lineWidth = 2 + drawSimpleArrow( + g, + 50, + height - 50, + 50 + (c1.toNumber() / cLen) * 40, + height - 50 - (c2.toNumber() / cLen) * 40, + 5 + ) + + g.fillStyle = 'darkgreen' + fillDot(g, 50, height - 50, 4) + + g.fillText(`c`, 50 - (c1.toNumber() / cLen) * 20, height - 50 + (c2.toNumber() / cLen) * 20) + + // g.beginPath() + // g.translate(50, height - 50) + // g.rotate(Math.atan2(c2.toNumber(), c1.toNumber())) + // g.moveTo(0, 0) + // g.lineTo(30, 0) + // g.moveTo(30, 0) + // g.lineTo(25, -5) + // g.moveTo(30, 0) + // g.lineTo(25, 5) + // g.stroke() + // g.restore() // g.fillStyle = '#333' // g.font = '16px sans-serif' @@ -108,232 +233,263 @@ const PrimaleCanvas = ({ }) // draw semiplanes in B + B.forEach(i => { const [a1, a2] = A.rowAt(i).getData() const b_i = b.at(i) - drawSemiplane(g, a1.toNumber(), a2.toNumber(), b_i.toNumber(), { lineColor: '#040', lineWidth: 3 }) + drawSemiplane(g, a1.toNumber(), a2.toNumber(), b_i.toNumber(), { + lineColor: '#040', + lineWidth: 3, + }) }) - // draw x + // draw current solution if (x) { const [x1, x2] = x.getData() - g.fillStyle = '#f00' - g.beginPath() - g.ellipse(x1.toNumber(), x2.toNumber(), 0.1, 0.1, 0, 0, 2 * Math.PI) - g.fill() - } - } - 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_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( -
- . -

-

- 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 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 ( -
-
{rows}
-
- -
-
- ) + 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_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( +//
+// . +//

+//

+// 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 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 ( +//
+//
{rows}
+//
+// +//
+//
+// ) +// } 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; + } } } }