⚙️ Utilisation d'un meilleur algorithme d'inversion

Pas pour autant lourd : une fonction de cent lignes

source/engine/traverse.js

@@ -36,6 +36,8 @@ import {
 	applyOrEmpty
 } from './traverse-common-functions'
 
+import uniroot from './uniroot'
+
 let nearley = () => new Parser(Grammar.ParserRules, Grammar.ParserStart)
 
 /*
@@ -403,34 +405,13 @@ export let computeRuleValue = (formuleValue, isApplicable) =>
 		? formuleValue
 		: isApplicable === false ? 0 : formuleValue == 0 ? 0 : null
 
-let computeInversion = (objective, computeGivenInput, currentValue) => {
-	let v = currentValue || objective, // notre première approximation est l'objectif lui-même (on suppose donc qu'ils sont du même ordre de grandeur, ce qui est vrai pour les salaires mais pas forcément pour d'autres variables évidemment)
-		here = computeGivenInput(v)
-
-	console.log('coucou', v, here)
-
-	if (Math.abs(here - objective) < 20 ) {
-		return v
-	}
-
-	let
-		ascend = computeGivenInput(v + 10),
-		descend = computeGivenInput(v - 10)
-
-	if (Math.abs(ascend - objective) < Math.abs(descend - objective))
-		return computeInversion(objective, computeGivenInput, v + 10)
-	else
-		return computeInversion(objective, computeGivenInput, v - 10)
-
-}
-
 
 export let treatRuleRoot = (rules, rule) => {
 	let evaluate = (situationGate, parsedRules, r) => {
 		let inversions = r['inversions possibles']
 		if (inversions) {
 			/*
-			Quel inversion possible est renseignée dans la situation courante ?
+			Quelle inversion possible est renseignée dans la situation courante ?
 			Ex. s'il nous est demandé de calculer le salaire de base, est-ce qu'un candidat à l'inversion, comme
 			le salaire net, a été renseigné ?
 			*/
@@ -439,14 +420,18 @@ export let treatRuleRoot = (rules, rule) => {
 			if (fixedObjective != null) {
 				let
 					fixedObjectiveRule = findRuleByName(parsedRules, fixedObjective),
-					nodeValue = computeInversion(
-						situationGate(fixedObjective),
-						i =>
-							evaluateNode(
-								n => (r.name === n || n === 'sys.filter') ? i : situationGate(n),
-								parsedRules,
-								fixedObjectiveRule
-							).nodeValue
+					fx = x => evaluateNode(
+						n => (r.name === n || n === 'sys.filter') ? x : situationGate(n), //TODO pourquoi doit-on nous préoccuper de sys.filter ?
+						parsedRules,
+						fixedObjectiveRule
+					).nodeValue,
+					tolerancePercentage = 0.00001,
+					nodeValue = uniroot(
+						x => fx(x) - situationGate(fixedObjective),
+						0,
+						1000000000,
+						tolerancePercentage * situationGate(fixedObjective),
+						100
 					)
 
 				return {nodeValue}

source/engine/uniroot.js

@@ -0,0 +1,116 @@
+/**
+ * Searches the interval from <tt>lowerLimit</tt> to <tt>upperLimit</tt>
+ * for a root (i.e., zero) of the function <tt>func</tt> with respect to
+ * its first argument using Brent's method root-finding algorithm.
+ *
+ * Translated from zeroin.c in http://www.netlib.org/c/brent.shar.
+ *
+ * Copyright (c) 2012 Borgar Thorsteinsson <borgar@borgar.net>
+ * MIT License, http://www.opensource.org/licenses/mit-license.php
+ *
+ * @param {function} function for which the root is sought.
+ * @param {number} the lower point of the interval to be searched.
+ * @param {number} the upper point of the interval to be searched.
+ * @param {number} the desired accuracy (convergence tolerance).
+ * @param {number} the maximum number of iterations.
+ * @returns an estimate for the root within accuracy.
+ *
+ */
+export default function uniroot(func, lowerLimit, upperLimit, errorTol, maxIter) {
+	var a = lowerLimit,
+		b = upperLimit,
+		c = a,
+		fa = func(a),
+		fb = func(b),
+		fc = fa,
+		tol_act, // Actual tolerance
+		new_step, // Step at this iteration
+		prev_step, // Distance from the last but one to the last approximation
+		p, // Interpolation step is calculated in the form p/q; division is delayed until the last moment
+		q
+
+	errorTol = errorTol || 0
+	maxIter = maxIter || 1000
+
+	while (maxIter-- > 0) {
+		prev_step = b - a
+
+		if (Math.abs(fc) < Math.abs(fb)) {
+			// Swap data for b to be the best approximation
+			(a = b), (b = c), (c = a)
+			;(fa = fb), (fb = fc), (fc = fa)
+		}
+
+		tol_act = 1e-15 * Math.abs(b) + errorTol / 2
+		new_step = (c - b) / 2
+
+		if (Math.abs(new_step) <= tol_act || fb === 0) {
+			return b // Acceptable approx. is found
+		}
+
+		// Decide if the interpolation can be tried
+		if (Math.abs(prev_step) >= tol_act && Math.abs(fa) > Math.abs(fb)) {
+			// If prev_step was large enough and was in true direction, Interpolatiom may be tried
+			var t1, cb, t2
+			cb = c - b
+			if (a === c) {
+				// If we have only two distinct points linear interpolation can only be applied
+				t1 = fb / fa
+				p = cb * t1
+				q = 1.0 - t1
+			} else {
+				// Quadric inverse interpolation
+				(q = fa / fc), (t1 = fb / fc), (t2 = fb / fa)
+				p = t2 * (cb * q * (q - t1) - (b - a) * (t1 - 1))
+				q = (q - 1) * (t1 - 1) * (t2 - 1)
+			}
+
+			if (p > 0) {
+				q = -q // p was calculated with the opposite sign; make p positive
+			} else {
+				p = -p // and assign possible minus to q
+			}
+
+			if (
+				p < 0.75 * cb * q - Math.abs(tol_act * q) / 2 &&
+				p < Math.abs(prev_step * q / 2)
+			) {
+				// If (b + p / q) falls in [b,c] and isn't too large it is accepted
+				new_step = p / q
+			}
+
+			// If p/q is too large then the bissection procedure can reduce [b,c] range to more extent
+		}
+
+		if (Math.abs(new_step) < tol_act) {
+			// Adjust the step to be not less than tolerance
+			new_step = new_step > 0 ? tol_act : -tol_act
+		}
+
+		(a = b), (fa = fb) // Save the previous approx.
+		;(b += new_step), (fb = func(b)) // Do step to a new approxim.
+
+		if ((fb > 0 && fc > 0) || (fb < 0 && fc < 0)) {
+			(c = a), (fc = fa) // Adjust c for it to have a sign opposite to that of b
+		}
+	}
+}
+
+/*
+var test_counter;
+function f1 (x) { test_counter++; return (Math.pow(x,2)-1)*x - 5; }
+function f2 (x) { test_counter++; return Math.cos(x)-x; }
+function f3 (x) { test_counter++; return Math.sin(x)-x; }
+function f4 (x) { test_counter++; return (x + 3) * Math.pow(x - 1, 2); }
+[
+  [f1, 2, 3],
+  [f2, 2, 3],
+  [f2, -1, 3],
+  [f3, -1, 3],
+  [f4, -4, 4/3]
+].forEach(function (args) {
+  test_counter = 0;
+  var root = uniroot.apply( pv, args );
+  ;;;console.log( 'uniroot:', args.slice(1), root, test_counter );
+})
+*/

test/inversion.test.js

@@ -42,12 +42,12 @@ describe("inversions", () => {
 
 				- nom: brut
 				  format: euro
-				  inversions possibles: 
+				  inversions possibles:
 				    - net
 			`,
 			rules = yaml.safeLoad(rawRules).map(enrichRule),
 			analysis = analyseSituation(rules, "brut")(stateSelector)
 
-		expect(analysis.nodeValue).to.be.closeTo(2570, 0.001)
+		expect(analysis.nodeValue).to.be.closeTo(2000/(77/100), 0.0001*2000)
 	})
 })