From e47303662fee7d61b016903a2acb9aef5f341f61 Mon Sep 17 00:00:00 2001
From: David Dorchies <david.dorchies@irstea.fr>
Date: Tue, 9 Jul 2019 14:55:39 +0200
Subject: [PATCH] =?UTF-8?q?#33=20Traitement=20du=20cas=20Z1=20<=20ZDV.=20?=
=?UTF-8?q?=20Fix=20#39=20Remplacement=20de=20la=20dichotomie=20par=20la?=
=?UTF-8?q?=20m=C3=A9thode=20de=20Brent.?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
---
spec/pab/pab.spec.ts | 6 +-
src/dichotomie.ts | 159 ++++++++++++++++++++++++++++++++++++++++---
src/pab/pab.ts | 20 ++++--
src/util/message.ts | 9 ++-
4 files changed, 174 insertions(+), 20 deletions(-)
diff --git a/spec/pab/pab.spec.ts b/spec/pab/pab.spec.ts
index d5dbc1c2..a927c4a8 100644
--- a/spec/pab/pab.spec.ts
+++ b/spec/pab/pab.spec.ts
@@ -205,11 +205,11 @@ describe("Class Pab: ", () => {
describe("when Z1 is too low compared to ZDV", () => {
it("logs should appear if dichotomy does not converge", () => {
pab.calculatedParam = pab.prms.Q;
- pab.prms.Z1.singleValue = 78.26; // @WTF 78.27 in the example above works, 78.26 fails !
+ pab.prms.Z1.singleValue = 76.60;
pab.CalcSerie();
expect(pab.result.hasLog).toBe(true);
expect(pab.result.resultElement.log.messages[0].code)
- .toBe(MessageCode.ERROR_DICHO_CONVERGE); // @TODO or other error message ?
+ .toBe(MessageCode.ERROR_PAB_Z1_LOWER_THAN_UPSTREAM_WALL); // @TODO or other error message ?
});
});
@@ -220,7 +220,7 @@ describe("Class Pab: ", () => {
pab.CalcSerie();
expect(pab.result.nbResultElements).toBe(3);
expect(pab.result.hasLog).toBe(true);
- expect(pab.result.resultElements[0].log.messages[0].code).toBe(MessageCode.ERROR_Z1_LOWER_THAN_Z2);
+ expect(pab.result.resultElements[0].log.messages[0].code).toBe(MessageCode.ERROR_PAB_Z1_LOWER_THAN_Z2);
});
});
diff --git a/src/dichotomie.ts b/src/dichotomie.ts
index 27a1059b..177159c2 100644
--- a/src/dichotomie.ts
+++ b/src/dichotomie.ts
@@ -131,6 +131,11 @@ export class Dichotomie extends Debug {
*/
private analyticalSymbol: string;
+ /**
+ * Targetted value calculated by the function
+ */
+ private _target: number;
+
/**
* Construction de la classe.
* @param nub Noeud de calcul contenant la méthode de calcul Equation
@@ -185,6 +190,7 @@ export class Dichotomie extends Debug {
* @param rInit valeur initiale approximative de x
*/
public Dichotomie(rTarget: number, rTol: number, rInit: number): Result {
+ this._target = rTarget;
// console.log("-----");
// for (let x = 0; x <= 1; x += 0.1)
// console.log(this.CalculX(x).vCalc);
@@ -201,15 +207,34 @@ export class Dichotomie extends Debug {
if (r.ok) {
const inter: SearchInterval = r.intSearch;
this._currentIterations = 0;
- // Dichotomie
- return this.dichoSearch(
- rTarget,
- rTol,
+ // Recherche du zero par la méthode de Brent
+ const s: number = this.uniroot(
+ this.CalcZero,
+ this,
inter.min,
inter.max,
- inter.targets.getVal(inter.minIndex),
- inter.targets.getVal(inter.maxIndex)
+ rTol
);
+ if (s === undefined) {
+ this.debug(
+ "Dichotomie : Non convergence"
+ );
+ return new Result(new Message(MessageCode.ERROR_DICHO_CONVERGE));
+ } else {
+ this.debug(
+ `Dichotomie : ${this.analyticalSymbol}(${this.sVarCalc}=${s}) = ${rTarget}`
+ );
+ return new Result(s);
+ }
+ // Dichotomie
+ // this.dichoSearch(
+ // rTarget,
+ // rTol,
+ // inter.min,
+ // inter.max,
+ // inter.targets.getVal(inter.minIndex),
+ // inter.targets.getVal(inter.maxIndex)
+ // );
} else {
return new Result(r.res);
}
@@ -223,9 +248,9 @@ export class Dichotomie extends Debug {
*/
private Calcul(): Result {
const r: Result = this.nub.Equation(this.analyticalSymbol);
- this.debug(
- "dicho : Calcul(vX=" + this.sVarCalc + "=" + this.vX + ") -> " +
- this.analyticalSymbol + "=" + r.vCalc);
+ // this.debug(
+ // "dicho : Calcul(vX=" + this.sVarCalc + "=" + this.vX + ") -> " +
+ // this.analyticalSymbol + "=" + r.vCalc);
return r;
}
@@ -251,7 +276,8 @@ export class Dichotomie extends Debug {
* @param intMax intervalle maximum (domaine de définition de la variable)
*/
private searchTarget(rTarget: number, intSearch: SearchInterval, intMax: Interval) {
- this.debug("searchTarget : Target " + rTarget + " dans " + intSearch.toString());
+ // tslint:disable-next-line:max-line-length
+ this.debug(`searchTarget debut: Target ${this.analyticalSymbol}=${rTarget} dans ${this.sVarCalc}=${intSearch.toString()}`);
let n = 0;
let ok: boolean = false;
@@ -267,6 +293,8 @@ export class Dichotomie extends Debug {
intSearch.growStep(2);
intSearch.next();
} while (n++ < this._startIntervalMaxSteps);
+ // tslint:disable-next-line:max-line-length
+ this.debug(`searchTarget fin: Target ${this.analyticalSymbol}=${rTarget} dans ${this.sVarCalc}=${intSearch.toString()}`);
return { ok, intSearch };
}
@@ -409,7 +437,8 @@ export class Dichotomie extends Debug {
return r;
}
- this.debug("dichoSearch : yTarget " + Ytarg + " X " + Xmin + "->" + Xmax + " tol " + rTol);
+ // tslint:disable-next-line:max-line-length
+ this.debug(`dichoSearch : yTarget ${this.analyticalSymbol}= ${Ytarg} ${this.sVarCalc}=[${Xmin}; ${Xmax}] tol ${rTol}`);
// tslint:disable-next-line:variable-name
const Xmid = (Xmin + Xmax) / 2;
@@ -437,4 +466,112 @@ export class Dichotomie extends Debug {
}
}
+ private CalcZero(x: number): number {
+ this.vX = x;
+ const r: Result = this.nub.Equation(this.analyticalSymbol);
+ if (r.ok) {
+ // this.debug(`CalcZero ${this.sVarCalc} = ${x} => ${this.analyticalSymbol} = ${r.vCalc}`);
+ return r.vCalc - this._target;
+ } else {
+ return undefined;
+ }
+ }
+
+ /**
+ * 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} func function for which the root is sought.
+ * @param {number} lowerlimit the lower point of the interval to be searched.
+ * @param {number} upperlimit the upper point of the interval to be searched.
+ * @param {number} errorTol the desired accuracy (convergence tolerance).
+ * @param {number} maxIter the maximum number of iterations.
+ * @returns an estimate for the root within accuracy.
+ *
+ */
+ private uniroot<T>(func: (param: number) => number, thisArg: T, lowerLimit: number, upperLimit: number,
+ errorTol: number = 0, maxIter: number = 1000
+ ): number {
+ let a = lowerLimit;
+ let b = upperLimit;
+ let c = a;
+ let fa = func.call(thisArg, a);
+ let fb = func.call(thisArg, b);
+ let fc = fa;
+ let tolAct; // Actual tolerance
+ let newStep; // Step at this iteration
+ let prevStep; // Distance from the last but one to the last approximation
+ let p; // Interpolation step is calculated in the form p/q; division is delayed until the last moment
+ let q;
+
+ while (maxIter-- > 0) {
+
+ prevStep = 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;
+ }
+
+ tolAct = 1e-15 * Math.abs(b) + errorTol / 2;
+ newStep = (c - b) / 2;
+
+ if (Math.abs(newStep) <= tolAct || fb === 0) {
+ return b; // Acceptable approx. is found
+ }
+
+ // Decide if the interpolation can be tried
+ if (Math.abs(prevStep) >= tolAct && Math.abs(fa) > Math.abs(fb)) {
+ // If prev_step was large enough and was in true direction, Interpolatiom may be tried
+ let t1;
+ let cb;
+ let 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(tolAct * q) / 2) &&
+ p < Math.abs(prevStep * q / 2)) {
+ // If (b + p / q) falls in [b,c] and isn't too large it is accepted
+ newStep = p / q;
+ }
+
+ // If p/q is too large then the bissection procedure can reduce [b,c] range to more extent
+ }
+
+ if (Math.abs(newStep) < tolAct) { // Adjust the step to be not less than tolerance
+ newStep = (newStep > 0) ? tolAct : -tolAct;
+ }
+
+ a = b, fa = fb; // Save the previous approx.
+ b += newStep, fb = func.call(thisArg, 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
+ }
+ }
+ return undefined;
+
+ }
+
}
diff --git a/src/pab/pab.ts b/src/pab/pab.ts
index b479e182..b34fef88 100644
--- a/src/pab/pab.ts
+++ b/src/pab/pab.ts
@@ -12,6 +12,7 @@ import { Result } from "../util/result";
import { CloisonAval, ParallelStructureParams } from "./cloison_aval";
import { Cloisons } from "./cloisons";
import { PabParams } from "./pab_params";
+import { ResultElement } from "../util/resultelement";
export { PabParams };
@@ -108,12 +109,22 @@ export class Pab extends Nub {
* @param sVarCalc Calcul
*/
public Calc(sVarCalc: string, rInit?: number): Result {
+ const r: Result = new Result(new ResultElement());
if (!this.downWall.checkVanneLevante()) {
- this.currentResult = new Result(new Message(MessageCode.ERROR_CLOISON_AVAL_UN_OUVRAGE_REGULE));
- return this.result;
+ r.resultElement.addMessage(new Message(MessageCode.ERROR_CLOISON_AVAL_UN_OUVRAGE_REGULE));
+ }
+ if (sVarCalc === "Q") {
+ if (this.prms.Z1.v < this.prms.Z2.v) {
+ r.resultElement.addMessage(new Message(MessageCode.ERROR_PAB_Z1_LOWER_THAN_Z2));
+ }
+ this.children[0].prms.Z1.v = this.prms.Z1.v;
+ this.children[0].prms.Z2.v = -Infinity;
+ if (this.children[0].CalcQ().vCalc < 1E-7) {
+ r.resultElement.addMessage(new Message(MessageCode.ERROR_PAB_Z1_LOWER_THAN_UPSTREAM_WALL));
+ }
}
- if (sVarCalc === "Q" && this.prms.Z1.v < this.prms.Z2.v) {
- this.currentResult = new Result(new Message(MessageCode.ERROR_Z1_LOWER_THAN_Z2));
+ if (r.hasErrorMessages) {
+ this.currentResult = r;
return this.result;
}
return super.Calc(sVarCalc, rInit);
@@ -333,6 +344,7 @@ export class Pab extends Nub {
}
private dbgWall(cl: ParallelStructure) {
+ if (!this.DBG) { return; }
let s: string = "";
for (const p of cl.prms) {
s += p.symbol + " = " + p.v + "; ";
diff --git a/src/util/message.ts b/src/util/message.ts
index e67ca8a4..30682bf0 100644
--- a/src/util/message.ts
+++ b/src/util/message.ts
@@ -119,9 +119,14 @@ export enum MessageCode {
ERROR_CLOISON_AVAL_UN_OUVRAGE_REGULE,
/**
- * La cote amont est plus basse que la cote aval
+ * PAB : La cote amont est plus basse que la cote aval
*/
- ERROR_Z1_LOWER_THAN_Z2,
+ ERROR_PAB_Z1_LOWER_THAN_Z2,
+
+ /**
+ * PAB : La cote amont est plus basse que la cloison amont de la passe
+ */
+ ERROR_PAB_Z1_LOWER_THAN_UPSTREAM_WALL,
/**
* internationalisation : langue non prise en charge
--
GitLab