diff --git a/spec/regime_uniforme_puissance.spec.ts b/spec/regime_uniforme_puissance.spec.ts index e8d46167c9c693f91435b9abc3c7750676d3a302..d8b5d98d4e4a4a29e75929ab54045f87a57bd596 100644 --- a/spec/regime_uniforme_puissance.spec.ts +++ b/spec/regime_uniforme_puissance.spec.ts @@ -1,6 +1,6 @@ /// <reference path="../node_modules/@types/jasmine/index.d.ts" /> -import { Result } from "../src/base"; +import { Result, ResultCode } from "../src/base"; import { RegimeUniforme } from "../src/regime_uniforme"; import { ParamsSectionPuiss, cSnPuiss } from "../src/section/section_puissance"; import { precDigits, precDist } from "./nubtest"; @@ -18,9 +18,7 @@ describe('Class RegimeUniforme / section puissance :', () => { 1 // YB= hauteur de berge // YCL=Condition limite en cote à l'amont ou à l'aval ); - let sect = new cSnPuiss(undefined, prms); - let ru = new RegimeUniforme(sect, false); expect(ru.Calc("k").vCalc).toBeCloseTo(0.635, precDigits); @@ -37,9 +35,7 @@ describe('Class RegimeUniforme / section puissance :', () => { 1 // YB= hauteur de berge // YCL=Condition limite en cote à l'amont ou à l'aval ); - let sect = new cSnPuiss(undefined, prms); - let ru = new RegimeUniforme(sect, false); expect(ru.Calc("LargeurBerge").vCalc).toBeCloseTo(3.474, 2); @@ -56,9 +52,7 @@ describe('Class RegimeUniforme / section puissance :', () => { 1 // YB= hauteur de berge // YCL=Condition limite en cote à l'amont ou à l'aval ); - let sect = new cSnPuiss(undefined, prms); - let ru = new RegimeUniforme(sect, false); expect(ru.Calc("Ks", 1e-8).vCalc).toBeCloseTo(33.774, 2); @@ -76,9 +70,7 @@ describe('Class RegimeUniforme / section puissance :', () => { 1 // YB= hauteur de berge // YCL=Condition limite en cote à l'amont ou à l'aval ); - let sect = new cSnPuiss(undefined, prms); - let ru = new RegimeUniforme(sect, false); expect(ru.Calc("If").vCalc).toBeCloseTo(0.002, 2); @@ -95,9 +87,7 @@ describe('Class RegimeUniforme / section puissance :', () => { 1 // YB= hauteur de berge // YCL=Condition limite en cote à l'amont ou à l'aval ); - let sect = new cSnPuiss(undefined, prms); - let ru = new RegimeUniforme(sect, false); expect(ru.Calc("Q").vCalc).toBeCloseTo(1.421, precDigits); @@ -114,9 +104,7 @@ describe('Class RegimeUniforme / section puissance :', () => { 1 // YB= hauteur de berge // YCL=Condition limite en cote à l'amont ou à l'aval ); - let sect = new cSnPuiss(undefined, paramCnl); - let ru = new RegimeUniforme(sect, false); expect(ru.Calc("Y").vCalc).toBeCloseTo(0.742, precDigits); @@ -124,7 +112,7 @@ describe('Class RegimeUniforme / section puissance :', () => { }); describe('débordement :', () => { - it('k should be 0.635', () => { + it('k should be undefined', () => { let prms = new ParamsSectionPuiss(undefined, // coef 2, // tirant d'eau 4, // largeur de berge @@ -135,12 +123,29 @@ describe('Class RegimeUniforme / section puissance :', () => { 1 // YB= hauteur de berge // YCL=Condition limite en cote à l'amont ou à l'aval ); - let sect = new cSnPuiss(undefined, prms); + let ru = new RegimeUniforme(sect, false); + + let res: Result = ru.Calc("k"); + expect(res.vCalc).toBeUndefined(); + expect(res.code).toBe(ResultCode.ERROR_DICHO_INIT_DOMAIN); + }); + it('k should be 0.635', () => { + let prms = new ParamsSectionPuiss(undefined, // coef + 2, // tirant d'eau + 4, // largeur de berge + 40, // Ks=Strickler + 10, // Q=Débit + 0.001, // If=pente du fond + precDist, // précision + 1 // YB= hauteur de berge + // YCL=Condition limite en cote à l'amont ou à l'aval + ); + let sect = new cSnPuiss(undefined, prms); let ru = new RegimeUniforme(sect, false); - expect(ru.Calc("k").vCalc).toBeCloseTo(0.635, precDigits); + expect(ru.Calc("k").vCalc).toBeCloseTo(0.933, precDigits); }); it('LargeurBerge should be 0.721', () => { @@ -160,7 +165,7 @@ describe('Class RegimeUniforme / section puissance :', () => { expect(ru.Calc("LargeurBerge").vCalc).toBeCloseTo(0.721, precDigits); }); - it('Strickler should be 4.366', () => { + it('Strickler should be 4.367', () => { let prms = new ParamsSectionPuiss(0.5, // coef 2, // tirant d'eau 4, // largeur de berge @@ -171,14 +176,51 @@ describe('Class RegimeUniforme / section puissance :', () => { 1 // YB= hauteur de berge // YCL=Condition limite en cote à l'amont ou à l'aval ); - let sect = new cSnPuiss(undefined, prms); - let ru = new RegimeUniforme(sect, false); - expect(ru.Calc("Ks", 1e-8).vCalc).toBeCloseTo(4.366, precDigits); + expect(ru.Calc("Ks", 1e-8).vCalc).toBeCloseTo(4.367, precDigits); + }); + + it('Ks should be undefined', () => { + let prms = new ParamsSectionPuiss(0.5, // coef + 2, // tirant d'eau + 4, // largeur de berge + undefined, // Ks=Strickler + 1.2, // Q=Débit + 0.001, // If=pente du fond + precDist, // précision + 1 // YB= hauteur de berge + // YCL=Condition limite en cote à l'amont ou à l'aval + ); + let sect = new cSnPuiss(undefined, prms); + let ru = new RegimeUniforme(sect); + ru.dichoStartIntervalMaxSteps = 3; + + let res: Result = ru.Calc("Ks", 5000000); + expect(res.vCalc).toBeUndefined(); + expect(res.code).toBe(ResultCode.ERROR_DICHO_INITVALUE_HIGH); }); + it('Ks should be undefined', () => { + let prms = new ParamsSectionPuiss(0.5, // coef + 2, // tirant d'eau + 4, // largeur de berge + undefined, // Ks=Strickler + 1.2, // Q=Débit + 0.001, // If=pente du fond + precDist, // précision + 1 // YB= hauteur de berge + // YCL=Condition limite en cote à l'amont ou à l'aval + ); + let sect = new cSnPuiss(undefined, prms); + let ru = new RegimeUniforme(sect); + ru.dichoStartIntervalMaxSteps = 1; + + let res: Result = ru.Calc("Ks", 1e-8); + expect(res.vCalc).toBeUndefined(); + expect(res.code).toBe(ResultCode.ERROR_DICHO_INITVALUE_LOW); + }); it('If should be 0.001', () => { let prms = new ParamsSectionPuiss(0.5, // coef @@ -191,9 +233,7 @@ describe('Class RegimeUniforme / section puissance :', () => { 1 // YB= hauteur de berge // YCL=Condition limite en cote à l'amont ou à l'aval ); - let sect = new cSnPuiss(undefined, prms); - let ru = new RegimeUniforme(sect, false); expect(ru.Calc("If").vCalc).toBeCloseTo(0.001, precDigits); @@ -210,9 +250,7 @@ describe('Class RegimeUniforme / section puissance :', () => { 1 // YB= hauteur de berge // YCL=Condition limite en cote à l'amont ou à l'aval ); - let sect = new cSnPuiss(undefined, prms); - let ru = new RegimeUniforme(sect, false); expect(ru.Calc("Q").vCalc).toBeCloseTo(10.993, precDigits); @@ -229,9 +267,7 @@ describe('Class RegimeUniforme / section puissance :', () => { 1 // YB= hauteur de berge // YCL=Condition limite en cote à l'amont ou à l'aval ); - let sect = new cSnPuiss(undefined, paramCnl); - let ru = new RegimeUniforme(sect, false); expect(ru.Calc("Y").vCalc).toBeCloseTo(0.742, precDigits); diff --git a/spec/regime_uniforme_rect.spec.ts b/spec/regime_uniforme_rect.spec.ts index d1fe6496666aa1250a3b766eef1682d16f82c0cf..e490f81dcb23beff8f21b5fd1454a3592cc20078 100644 --- a/spec/regime_uniforme_rect.spec.ts +++ b/spec/regime_uniforme_rect.spec.ts @@ -26,7 +26,7 @@ describe('Class RegimeUniforme / section rectangulaire :', () => { expect(ru.Calc("LargeurBerge").vCalc).toBeCloseTo(2.5, precDigits); }); - it('Strickler should be 30.618', () => { + it('Strickler should be 30.619', () => { let prms = new ParamsSectionRectang(0.8, // tirant d'eau 2.5, // largeur de fond undefined, // Ks=Strickler @@ -40,7 +40,7 @@ describe('Class RegimeUniforme / section rectangulaire :', () => { let sect = new cSnRectang(undefined, prms); let ru = new RegimeUniforme(sect); - expect(ru.Calc("Ks", 1e-8).vCalc).toBeCloseTo(30.618, precDigits); + expect(ru.Calc("Ks", 1e-8).vCalc).toBeCloseTo(30.619, precDigits); }); it('If should be 0.001', () => { @@ -135,7 +135,7 @@ describe('Class RegimeUniforme / section rectangulaire :', () => { expect(ru.Calc("LargeurBerge").vCalc).toBeCloseTo(2.5, precDigits); }); - it('Strickler should be 9.04', () => { + it('Strickler should be 9.041', () => { let prms = new ParamsSectionRectang(2, // tirant d'eau 2.5, // largeur de fond undefined, // Ks=Strickler @@ -149,7 +149,7 @@ describe('Class RegimeUniforme / section rectangulaire :', () => { let sect = new cSnRectang(undefined, prms); let ru = new RegimeUniforme(sect); - expect(ru.Calc("Ks", 1e-8).vCalc).toBeCloseTo(9.04, precDigits); + expect(ru.Calc("Ks", 1e-8).vCalc).toBeCloseTo(9.041, precDigits); }); it('If should be 0.001', () => { diff --git a/spec/regime_uniforme_trapeze.spec.ts b/spec/regime_uniforme_trapeze.spec.ts index 3795359d3e9c07d8a023d47e95797b8399a85ea8..4625d803e24ad1979d404ddd6643221454385e74 100644 --- a/spec/regime_uniforme_trapeze.spec.ts +++ b/spec/regime_uniforme_trapeze.spec.ts @@ -153,7 +153,7 @@ describe('Class RegimeUniforme / section trapèze :', () => { expect(ru.Calc("Fruit").vCalc).toBeCloseTo(0.56, precDigits); }); - it('Ks should be 5.744', () => { + it('Ks should be 5.745', () => { let prms = new ParamsSectionTrapez(2.5, // largeur de fond 0.56, // fruit 2, // tirant d'eau @@ -167,7 +167,7 @@ describe('Class RegimeUniforme / section trapèze :', () => { let sect = new cSnTrapez(undefined, prms); let ru = new RegimeUniforme(sect); - expect(ru.Calc("Ks", 1e-8).vCalc).toBeCloseTo(5.744, precDigits); + expect(ru.Calc("Ks", 1e-8).vCalc).toBeCloseTo(5.745, precDigits); }); it('If should be 0.001', () => { diff --git a/src/base.ts b/src/base.ts index ff3902082344a76a75b70d8f9e21a8f4c47db1b8..e20284f8063f7a0ae379aed6f0ccf38d93761b96 100644 --- a/src/base.ts +++ b/src/base.ts @@ -1,3 +1,31 @@ +export enum ResultCode { + OK = 0, + + /** + * la dichotomie n'a pas pu trouver automatiquement d'intervalle de départ + * car la valeur initiale de la variable est trop haute + */ + ERROR_DICHO_INITVALUE_HIGH = -1, + + /** + * la dichotomie n'a pas pu trouver automatiquement d'intervalle de départ + * car la valeur initiale de la variable est trop basse + */ + ERROR_DICHO_INITVALUE_LOW = -2, + + /** + * la dichotomie n'a pas pu trouver automatiquement d'intervalle de départ + * car la valeur cible de la fonction n'existe pas pour des valeurs de la + * variable dans son domaine de définition, cad il n'existe pas de solution + */ + ERROR_DICHO_INIT_DOMAIN = -3, + + /** + * la dichotomie n'a pas pu converger + */ + ERROR_DICHO_CONVERGE = -4, +} + /** * Résultat de calcul comprenant la valeur du résultat et des calculs annexes (flag, calculs intermédiaires...) */ @@ -5,15 +33,23 @@ export class Result { /** Valeur calculée */ private _vCalc: number; + /** + * code d'erreur + */ + private _code: ResultCode; + /** Variables intermédiaires, flags d'erreur */ public extraVar: {}; - constructor(v: number) { + constructor(v: number, c: ResultCode = ResultCode.OK) { this._vCalc = v + this._code = c; this.extraVar = {} }; get vCalc() { return this._vCalc; } + + get code() { return this._code; } } /** @@ -56,9 +92,16 @@ export abstract class Debug { * Méthode simulant l'opérateur booléen xor * @see http://www.howtocreate.co.uk/xor.html */ -export function XOR(a: boolean, b: boolean) { +export function XOR(a: boolean, b: boolean): boolean { return (a || b) && !(a && b); } +/** + * Méthode simulant l'opérateur booléen identité (vrai si a=b) + */ +export function BoolIdentity(a: boolean, b: boolean): boolean { + return (a && b) || (!a && !b); +} + export class UndefinedError extends Error { } diff --git a/src/dichotomie.ts b/src/dichotomie.ts index 2c0512bfe47618c585abe7abd7e15ac7e662603a..2c49318ec294289f6b8a3223ed8947e39801b950 100644 --- a/src/dichotomie.ts +++ b/src/dichotomie.ts @@ -1,15 +1,195 @@ -import { XOR, Debug, Result } from "./base"; +import { XOR, BoolIdentity, Debug, Result, ResultCode, UndefinedError } from "./base"; import { Nub } from "./nub"; +import { ParamDefinition, ParamDomain, ParamDomainValue } from "./param" + +/** + * couple de valeurs non ordonnées + */ +class Pair { + protected _val1: number; + protected _val2: number; + + constructor(v1: number, v2: number) { + this.setValues(v1, v2); + } + + get val1() { + return this._val1; + } + + get val2() { + return this._val2; + } + + setValues(v1: number, v2: number) { + this._val1 = v1; + this._val2 = v2; + } + + setPair(p: Pair) { + this.setValues(p._val1, p._val2); + } + + undefine() { + this._val1 = undefined; + this._val2 = undefined; + } + + get min() { + return Math.min(this._val1, this._val2); + } + + get max() { + return Math.max(this._val1, this._val2); + } + + intervalHasValue(v: number) { + return this.min <= v && v <= this.max; + } +} + +/** + * couple de valeurs ordonnées + */ +class Interval extends Pair { + constructor(bound1, bound2) { + super(bound1, bound2); + } + + checkValue(v: number) { + if (v == undefined) + throw new UndefinedError("Interval.checkValue() : invalid undefined value"); + + if (!this.intervalHasValue(v)) + throw new RangeError("Interval : value " + v + " is outside of " + this.toString()); + } + + get length(): number { + return this.max - this.min; + } + + intersect(i: Interval): Interval { + let min = Math.max(this.min, i.min); + let max = Math.min(this.max, i.max); + return new Interval(min, max); + } + + setBounds(i: Pair) { + this.setPair(i); + } + + toString(): string { + return "[" + this.min + "," + this.max + "]"; + } +} + +class SearchInterval extends Interval { + private _step; + + private _dicho: Dichotomie; + + private _targets: Pair; + + constructor(d: Dichotomie, min, max, s) { + super(min, max); + if (s == 0) + throw new RangeError("SearchInterval : invalid null step"); + this._step = s; + this._dicho = d; + } + + reverse() { + this._step = -this._step; + } + + growStep(k: number) { + if (k == 0) + throw new RangeError("SearchInterval : invalid null step grow"); + this._step *= k; + } + + get step() { + return this._step; + } + + get targets() { + if (this._targets == undefined) + this._targets = new Pair(undefined, undefined); + return this._targets; + } + + private updateTargets() { + let t1 = this.targets.val1; + if (t1 == undefined) + t1 = this._dicho.CalculX(this._val1).vCalc + + let t2 = this.targets.val2; + if (t2 == undefined) + t2 = this._dicho.CalculX(this._val2).vCalc; + this.targets.setValues(t1, t2); + } + + /** + * get target value for lower bound + */ + get targetLower() { + this.updateTargets(); + return this.targets.val1; + } + + /** + * get target value for upper bound + */ + get targetUpper() { + this.updateTargets(); + return this.targets.val2; + } + + next() { + if (this._step > 0) { + this._val1 = this._val2; + this._val2 += this._step; + this.targets.setValues(this.targets.val2, undefined); + } + else { + this._val2 = this._val1; + this._val1 += this._step; + this.targets.setValues(undefined, this.targets.val1); + } + } + + hasTargetValue(targ: number) { + this.updateTargets(); + return this.targets.intervalHasValue(targ); + } + + // non fiable : certaines fonctions (coef section puiss par ex) sont globalement monotones mais peuvent inverser localement leur sens de variation + // isIncreasingFunction(): boolean { + // return this.targetUpper > this.targetLower; + // } +} + /** * calcul par dichotomie */ export class Dichotomie extends Debug { /** Pas de parcours de l'intervalle pour initialisation dichotomie */ - private readonly IDEFINT = 100; + // private readonly IDEFINT = 100; /** Nombre d'itérations maximum de la dichotomie */ - private readonly IDICMAX = 100; + // private readonly IDICMAX = 100; + + /** + * définition de la variable de la fonction + */ + private _paramX: ParamDefinition; + + + /** + * nombre d'étapes de recherche de l'intervalle de départ + */ + private _startIntervalMaxSteps = 100; /** * Construction de la classe. @@ -18,17 +198,33 @@ export class Dichotomie extends Debug { */ constructor(private nub: Nub, private sVarCalc: string, dbg: boolean = false) { super(dbg); + this._paramX = this.nub.getParameter(this.sVarCalc); } /** * Valeur inconnue à rechercher */ get vX(): number { - return this.nub.getParameter(this.sVarCalc).v; + return this.paramX.v; } set vX(vCalc: number) { - this.nub.getParameter(this.sVarCalc).v = vCalc; + this.paramX.v = vCalc; + } + + get paramX() { + return this._paramX; + } + + /** + * étapes de recherche de l'intervalle de départ + */ + get startIntervalMaxSteps() { + return this._startIntervalMaxSteps; + } + + set startIntervalMaxSteps(n: number) { + this._startIntervalMaxSteps = n; } /** @@ -37,128 +233,267 @@ export class Dichotomie extends Debug { * On utilise la première variable du tableau des variables pouvant être calculée analytiquement * Il faudra s'assurer que cette première variable correspond à la méthode de calcul la plus rapide */ - private Calcul() { - //return this.nub.Equation(this.nub.getFirstAnalyticalParameter().symbol); - let targetParam = this.nub.getFirstAnalyticalParameter(); - let r = this.nub.Equation(targetParam.symbol); + private Calcul(): Result { + let targetParam: ParamDefinition = this.nub.getFirstAnalyticalParameter(); + let r: Result = this.nub.Equation(targetParam.symbol); this.debug("dicho : Calcul(vX=" + this.sVarCalc + "=" + this.vX + ") -> " + targetParam.symbol + "=" + r.vCalc); return r; } - private CalculX(x: number) { + CalculX(x: number): Result { this.vX = x return this.Calcul(); } /** - * Calcul à l'ouvrage - * @param rTarget Valeur cible - * @param rtol Précision attendue du résultat - * @param rInit Valeur initiale de l'inconnue à rechercher - */ - Dichotomie(rTarget: number, rTol: number, rInit: number): Result { - let res: Result; + * Détermine si une fonction est croissante ou décroissante. + * On fait ça en sondant plusieurs points sur un intervalle : on ne peut pas + * prendre 2 valeurs de la fonction pour 2 valeurs successives de la variable car certaines + * fonctions (Q=f(coef de la section puissance) par ex) ne sont pas strictement monotones (globalement + * monotones mais peuvent inverser localement leur sens de variation) + * @param d domaine de définition de la variable + */ + private isIncreasingFunction(d: ParamDomain): boolean { + let n = 20; + let variation = 0; + let min, max, step; + + switch (d.domain) { + case ParamDomainValue.INTERVAL: + min = d.minValue; + max = d.maxValue; + break; + + case ParamDomainValue.ANY: + case ParamDomainValue.NOT_NULL: + min = -100; + max = 100; + break; - let XminInit: number = 1E-8; - let XmaxInit: number = Math.max(1, rInit) * 100; - - let DX = (XmaxInit - XminInit) / this.IDEFINT; - let nIterMax = Math.floor(Math.max(XmaxInit - rInit, rInit - XminInit) / DX + 1); - let Xmin = rInit; - let Xmax = rInit; - let X1 = rInit; - let X2 = rInit; - this.debug("dicho : rInit(" + this.sVarCalc + ")=" + rInit + " target=" + rTarget); - let v = this.CalculX(rInit).vCalc; - - if (isNaN(v)) { - rInit += 1e-8; - v = this.CalculX(rInit).vCalc; + case ParamDomainValue.POS: + min = 1e-8; + max = 200; + break; + + case ParamDomainValue.POS_NULL: + min = 0; + max = 200; + break; + + default: + throw "unsupported parameter domain value" } - let v1 = v; - let v2 = v; - - this.debug("dicho : recherche intervalle"); - - //** @todo : Chercher en dehors de l'intervalle en le décalant à droite ou à gauche en fonction de la valeur */ - let nIter: number; - for (nIter = 1; nIter < nIterMax; nIter++) { - //Ouverture de l'intervalle des deux côtés : à droite puis à gauche - Xmax = Xmax + DX; - if (XOR(Xmax > XmaxInit, DX <= 0)) - Xmax = XmaxInit; - this.debug("dicho interv. : iter=" + nIter); - v = this.CalculX(Xmax).vCalc; - if (XOR(v1 < v2, v <= v2)) { - v2 = v; - X2 = Xmax; - } - if (XOR(v1 < v2, v >= v1)) { - v1 = v; - X1 = Xmax; - } + let x = min; + step = (max - min) / n; + let lastY; + for (let i = 0; i < n; i++) { + if (x == 0 && d.domain == ParamDomainValue.NOT_NULL) + x = 1e-8; + let y = this.CalculX(x).vCalc; - Xmin = Xmin - DX; - if (XOR(Xmin < XminInit, DX <= 0)) { - Xmin = XminInit; + if (i > 0) { + if (y > lastY) + variation++; + else + variation--; } + lastY = y; - this.debug("dicho interv. : iter=" + nIter); - v = this.CalculX(Xmin).vCalc; - if (XOR(v1 < v2, v <= v2)) { - v2 = v; - X2 = Xmin; - } - if (XOR(v1 < v2, v >= v1)) { - v1 = v; - X1 = Xmin; - } + x += step; + } + + if (variation == 0) + throw "unable to determinate function direction of variation"; - if (XOR(rTarget > v1, rTarget >= v2)) + return variation > 0; + } + + /** + * Détermine l'intervalle de recherche initial + * @param rTarget valeur cible de la fonction + * @param rInit valeur initiale de la variable + */ + getStartInterval(rTarget: number, rInit: number): any { + let prmDom: ParamDomain = this.paramX.getDomain(); + let min, max, step; + + switch (prmDom.domain) { + case ParamDomainValue.INTERVAL: + min = prmDom.minValue; + max = prmDom.maxValue; + step = (max - min) / 100.0; + break; + + case ParamDomainValue.ANY: + min = -Infinity; + max = Infinity; + step = 1; + break; + + case ParamDomainValue.NOT_NULL: + min = -Infinity; + max = Infinity; + step = 1; + if (rInit == 0) + rInit = 1e-8; + break; + + case ParamDomainValue.POS: + min = 1e-8; + max = Infinity; + step = 1; + break; + + case ParamDomainValue.POS_NULL: + min = 0; + max = Infinity; + step = 1; break; + + default: + throw "unsupported parameter domain value"; } - this.debug("intervalle initial " + X1 + ", " + X2); - - // Gestion de l'absence de solution dans l'intervalle de recherche - if (nIter >= this.IDEFINT) { - this.debug("nIter >= this.IDEFINT"); - this.debug("dicho interv. : iter=" + nIter); - if (v2 < rTarget && v1 < rTarget) { - // Cote de l'eau trop basse pour passer le débit il faut ouvrir un autre ouvrage - res = this.CalculX(XmaxInit); - } - else { - // Cote de l'eau trop grande il faut fermer l'ouvrage - res = this.CalculX(XminInit); + + let intMax: Interval = new Interval(min, max); + intMax.checkValue(rInit); + + let intSearch: SearchInterval = new SearchInterval(this, rInit, rInit + step, step); + + // let inc = intSearch.isIncreasingFunction(); + let inc = this.isIncreasingFunction(prmDom); + + let initTarget = this.CalculX(rInit).vCalc; + // if (initTarget > rTarget && inc || initTarget < rTarget && !inc) + if (BoolIdentity(initTarget > rTarget, inc)) + intSearch.reverse(); + + let n = 0; + let ok: boolean = false; + do { + var inters: Interval = intSearch.intersect(intMax); + if (inters.length == 0) + break; + intSearch.setBounds(inters); + if (intSearch.hasTargetValue(rTarget)) { + ok = true; + break; } - res.extraVar["flag"] = -1; // la valeur cible n'est pas dans l'intervalle de recherche - return res; + intSearch.growStep(2); + intSearch.next(); + } while (n++ < this._startIntervalMaxSteps); + + // intervalle trouvé + + if (ok) + return { ok, intSearch }; + + // la valeur cible de la fonction est elle trouvable ? + + let aps = this.nub.getFirstAnalyticalParameter().symbol; + let m; + + let res: Result; + let errDomain = false; + switch (prmDom.domain) { + case ParamDomainValue.INTERVAL: + let si: SearchInterval = new SearchInterval(this, intMax.min, intMax.max, 1); + errDomain = !si.hasTargetValue(rTarget); + break; + + case ParamDomainValue.POS: + case ParamDomainValue.POS_NULL: + let y = this.CalculX(1e-8).vCalc; + errDomain = inc && (rTarget < y); + break; + + case ParamDomainValue.ANY: + break; + + default: + throw "unsupported parameter domain value"; } - // Dichotomie - this.debug("start dicho"); - let X = rInit; - for (nIter = 1; nIter <= this.IDICMAX; nIter++) { - this.debug("dicho : iter=" + nIter); - v = this.CalculX(X).vCalc; - if (rTarget != 0 && Math.abs(X1 - X2) <= rTol) { break; } - if (XOR(rTarget < v, v1 <= v2)) { - // rTarget < IQ et v(X1) > v(X2) ou pareil en inversant les inégalités - X1 = this.vX; - } - else { - // rTarget < IQ et v(X1) < v(X2) ou pareil en inversant les inégalités - X2 = this.vX; - } - X = (X2 + X1) * 0.5; - //this.debug((X)); + if (errDomain) { + // m = "Dichotomy : target " + aps + "=" + rTarget + " does not exist in interval " + intMax.toString() + " for variable " + this.paramX.symbol; + res = new Result(undefined, ResultCode.ERROR_DICHO_INIT_DOMAIN); + res.extraVar["targetSymbol"] = aps; // symbole de la variable calculée par la fonction + res.extraVar["targetValue"] = rTarget; // valeur cible pour la fonction + res.extraVar["variableInterval"] = intMax.toString(); // valeur de la fonction pour la valeur initiale de la variable + res.extraVar["variableSymbol"] = this.paramX.symbol; // symbole de la variable de la fonction } - if (nIter == this.IDICMAX) { - res = this.Calcul(); - res.extraVar["flag"] = -1; // la valeur cible n'est pas dans l'intervalle de recherche - return res; + else { + // m = "Dichotomy : initial value " + this.paramX.symbol + "=" + rInit; + // if (intSearch.step < 0) + // m += " is to high"; + // else + // m += " is to low"; + // m += " (target is " + aps + "=" + rTarget + ", " + aps + "(" + this.paramX.symbol + "=" + rInit + ")=" + initTarget + ")"; + if (intSearch.step < 0) + res = new Result(undefined, ResultCode.ERROR_DICHO_INITVALUE_HIGH); + else + res = new Result(undefined, ResultCode.ERROR_DICHO_INITVALUE_LOW); + + res.extraVar["variableSymbol"] = this.paramX.symbol; // symbole de la variable de la fonction + res.extraVar["variableInitValue"] = rInit; // valeur initiale de la variable + res.extraVar["targetSymbol"] = aps; // symbole de la variable calculée par la fonction + res.extraVar["targetValue"] = rTarget; // valeur cible pour la fonction + res.extraVar["initTarget"] = initTarget; // valeur de la fonction pour la valeur initiale de la variable } - return new Result(X); + + return { ok, res }; + } + + /** + * @param Xmin valeur min de l'intervalle + * @param Xmax valeur max de l'intervalle + * @param Y valeur cible de la fonction + * @param Ymin valeur de la fonction pour Xmin + * @param Ymax valeur de la fonction pour Xmax + */ + private dichoSearch(Ytarg: number, rTol: number, Xmin: number, Xmax: number, Ymin: number, Ymax: number) { + let Xmid = (Xmin + Xmax) / 2; + + if (Math.abs(Xmin - Xmax) <= rTol) + return new Result(Xmid); + + let Ymid = this.CalculX(Xmid).vCalc; + if (Ymin < Ymax) { + // fonction croissante + + if (Ytarg > Ymid) + return this.dichoSearch(Ytarg, rTol, Xmid, Xmax, Ymid, Ymax); // intervalle supérieur + return this.dichoSearch(Ytarg, rTol, Xmin, Xmid, Ymin, Ymid); // intervalle inférieur + } + else { + // fonction décroissante + + if (Ytarg > Ymid) + return this.dichoSearch(Ytarg, rTol, Xmin, Xmid, Ymin, Ymid); // intervalle inférieur + return this.dichoSearch(Ytarg, rTol, Xmid, Xmax, Ymid, Ymax); // intervalle supérieur + } + } + + Dichotomie(rTarget: number, rTol: number, rInit: number): Result { + let res: Result; + + // console.log("-----"); + // for (let x = 0; x <= 1; x += 0.1) + // console.log(this.CalculX(x).vCalc); + // console.log("-----"); + + // recherche de l'intervalle de départ + + this.debug("intervalle de départ"); + let r = this.getStartInterval(rTarget, rInit); + if (r.ok) + var interv: SearchInterval = r.intSearch; + else + return r.res; + + // Dichotomie + + this.debug("start dicho"); + return this.dichoSearch(rTarget, rTol, interv.min, interv.max, interv.targetLower, interv.targetUpper); } } diff --git a/src/nub.ts b/src/nub.ts index c5f1223e5be0a77609e0f642d83a42c215ff05ca..e77c81c9e61517d24150ab33dc17385baaea3c55 100644 --- a/src/nub.ts +++ b/src/nub.ts @@ -6,6 +6,20 @@ import { ComputeNode, ParamDefinition, IParamsEquation } from "./param" * Classe abstraite de Noeud de calcul : classe de base pour tous les calculs */ export abstract class Nub extends ComputeNode { + private _dichoStartIntervalMaxSteps: number = 100; + + /* + * paramétrage de la dichotomie + */ + + /** + * étapes de recherche de l'intervalle de départ + */ + set dichoStartIntervalMaxSteps(n: number) { + this._dichoStartIntervalMaxSteps = n; + } + + private checkParametersCalculability() { let res = []; @@ -37,7 +51,6 @@ export abstract class Nub extends ComputeNode { return this.Solve(sVarCalc, rInit, rPrec); } - CalcSerie(svarCalc: string, serie: Serie): Result[] { /** @todo faire une boucle pour appeler this.Calc avec chaque valeur de serie.values * @@ -47,7 +60,7 @@ export abstract class Nub extends ComputeNode { return results; } - public getParameter(name: string) { + public getParameter(name: string): ParamDefinition { for (let ps in this._prms) { let p: ParamDefinition = this._prms[ps]; if (p.symbol == name) @@ -74,6 +87,7 @@ export abstract class Nub extends ComputeNode { */ Solve(sVarCalc: string, rInit: number, rPrec: number): Result { let dicho: Dichotomie = new Dichotomie(this, sVarCalc, this.DBG); + dicho.startIntervalMaxSteps = this._dichoStartIntervalMaxSteps; var target = this.getFirstAnalyticalParameter(); return dicho.Dichotomie(target.v, rPrec, rInit); } diff --git a/src/param.ts b/src/param.ts index f164f4a6eacfafe4636aa14b7fa07275b14f58b1..601741b9c2f32aa0531b4ba30fe9bbf4bfe11022 100644 --- a/src/param.ts +++ b/src/param.ts @@ -38,7 +38,7 @@ export class ParamDomain { private _maxValue: number; - constructor(d: ParamDomainValue, min: number = 0, max: number = -1) { + constructor(d: ParamDomainValue, min: number = undefined, max: number = undefined) { this.checkValue(d, min, max); this._domain = d; this._minValue = min; @@ -48,7 +48,7 @@ export class ParamDomain { private checkValue(val, min, max) { switch (val) { case ParamDomainValue.INTERVAL: - if (min > max) + if (min == undefined || max == undefined || min > max) throw "invalid " + min + "/" + max + " min/max boundaries for 'interval' parameter definition domain" break; @@ -133,10 +133,6 @@ export class ParamDefinition { this._value = val; } - get() { - return this._value; - } - /** * getter symbole */ @@ -144,6 +140,13 @@ export class ParamDefinition { return this._symbol; } + /** + * méthodes de domaine + */ + getDomain(): ParamDomain { + return this._domain; + } + /* * méthodes de calculabilité */ diff --git a/src/regime_uniforme.ts b/src/regime_uniforme.ts index a92d7434dff869b86796e9493ebdd4575b3c7356..35a9d0d0650b38351c3353a78029db010cd574e9 100644 --- a/src/regime_uniforme.ts +++ b/src/regime_uniforme.ts @@ -29,8 +29,8 @@ export class RegimeUniforme extends Nub { if (this.Sn.prms.If.v <= 0) return 0; - this.debug("RU: Calc_Qn : "); - this.debug(this.Sn.prms); + // this.debug("RU: Calc_Qn : "); + // this.debug(this.Sn.prms); return this.Sn.prms.Ks.v * Math.pow(this.Sn.Calc("R", this.Sn.prms.Y.v), 2 / 3) * this.Sn.Calc("S", this.Sn.prms.Y.v) * Math.sqrt(this.Sn.prms.If.v); }