From fe928c5eeb4e04f0fab0576f81b4e78457636293 Mon Sep 17 00:00:00 2001 From: Tobias Glaubach Date: Tue, 10 Jul 2018 21:46:12 +0200 Subject: [PATCH 01/18] Implemented a prototype for a polynomial class only based on real values (doubles) --- src/Numerics/Polynomial.cs | 471 +++++++++++++++++++++++++++++++++++++ 1 file changed, 471 insertions(+) create mode 100644 src/Numerics/Polynomial.cs diff --git a/src/Numerics/Polynomial.cs b/src/Numerics/Polynomial.cs new file mode 100644 index 00000000..360cfb27 --- /dev/null +++ b/src/Numerics/Polynomial.cs @@ -0,0 +1,471 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using System.Numerics; +using MathNet.Numerics; +using MathNet.Numerics.LinearAlgebra; +using MathNet.Numerics.LinearAlgebra.Double; +using MathNet.Numerics.Statistics; +using MathNet.Numerics.IntegralTransforms; + +using MathNet.Numerics.LinearAlgebra.Factorization; + + +namespace MathNet.Numerics +{ + /// + /// a class handlin REAL VALUED Polynomials, complex coefficients can not be handled (yet) + /// + public class Polynomial + { + + public double[] Coeffs { get; set; } + + /// + /// indicator if Polynomial was flipped + /// + public bool IsFlipped { get; } + + /// + /// Only needed for the ToString method + /// + public string VarName = "x^"; + + /// + /// Length of Polynomial (max element + 1) e.G x^5 highest element, will give Length = 6 + /// + public int Length + { + get + { + return (Coeffs.Length); + } + } + + /// + /// constructor setting a Polynomial of size n containing only zeros + /// + /// size of Polynomial + public Polynomial(int n) + { + Coeffs = new double[n]; + } + + /// + /// constructor setting Polynomial coefficiens and flipping them if necessary. + /// + /// e.G: + /// var x = new double[] {5, 4, 3, 0, 2}; + /// var xP1 = new Polynomial(x, isFlip:true); + /// var xP2 = new Polynomial(x, isFlip:false); + /// + /// xP1: 5 x^3 + 4 x^2 + 3 x^2 + 0 x^1 + 2 + /// xP2: 2 x^3 + 0 x^2 + 3 x^2 + 4 x^1 + 5 + /// + /// WARNING cut all trailing zeros before, since they would result in zeros at the end + /// + /// Polynomial coefficiens as array + /// use true for flipping + public Polynomial(double[] coeffs, bool isFlip = false) + { + this.Coeffs = new double[Coeffs.Length]; + Array.Copy(coeffs, Coeffs, coeffs.Length); + + if (isFlip) + { + Coeffs = Coeffs.Reverse().ToArray(); + IsFlipped = true; + } + } + + /// + /// constructor setting Polynomial coefficiens + /// + /// just the x^0 part + public Polynomial(double coeff) + { + IsFlipped = false; + this.Coeffs = new double[1]; + Coeffs[0] = coeff; + } + + /// + /// constructor setting Polynomial coefficiens + /// + /// Polynomial coefficiens as array + public Polynomial(double[] Coeffs) + { + this.Coeffs = new double[Coeffs.Length]; + Array.Copy(Coeffs, this.Coeffs, Coeffs.Length); + } + + /// + /// remove all trailing zeros, e.G before: "0.00 x^2 + 1.0 x^1 + 1.00" after: "1.0 x^1 + 1.00" + /// + public void CutTrailZeros() + { + int count = 0; + for (int ii = Length - 1; ii >= 0; ii--) + { + if (Coeffs[ii] == 0.0) + { + count++; + } + else + { + double[] CoeffsHold = new double[Length]; + Coeffs.CopyTo(CoeffsHold, 0); + Array.Resize(ref CoeffsHold, Length - count); + Coeffs = new double[Length - count]; + CoeffsHold.CopyTo(Coeffs, 0); + return; + } + } + } + + #region Operators + + + /// + /// multiplies a Polynomial by a Polynomial using convolution [ASINCO.libs.subfun.conv(a.Coeffs, b.Coeffs)] + /// + /// left Polynomial + /// right Polynomial + /// resulting Polynomial + public static Polynomial operator *( Polynomial a, Polynomial b) + { + // do not cut trailing zeros, since it may corrupt the outcom, if the array is of form 1 + x^-1 + x^-2 + x^-3 + //a.CutTrailZeros(); + //b.CutTrailZeros(); + + double[] ret = conv(a.Coeffs, b.Coeffs); + Polynomial ret_p = new Polynomial(ret); + + //ret_p.CutTrailZeros(); + + return (ret_p); + + } + + /// + /// multiplies a Polynomial by a scalar + /// + /// left Polynomial + /// scalar value + /// resulting Polynomial + public static Polynomial operator *( Polynomial a, double k) + { + for (int ii = 0; ii < a.Length; ii++) + a.Coeffs[ii] *= k; + + return a; + } + + /// + /// adds a scalar to a Polynomial (to the x^0 element) + /// + /// left Polynomial + /// scalar value + /// resulting Polynomial + public static Polynomial operator +( Polynomial a, double k) + { + a.Coeffs[0] += k; + return a; + } + + /// + /// substracs a scalar from a Polynomial (from the x^0 element) + /// + /// left Polynomial + /// scalar value + /// resulting Polynomial + public static Polynomial operator -( Polynomial a, double k) + { + + a.Coeffs[0] -= k; + return a; + } + + /// + /// divide Polynomial by scalar value + /// + /// left Polynomial + /// scalar value + /// resulting Polynomial + public static Polynomial operator /( Polynomial a, double k) + { + for (int ii = 0; ii < a.Length; ii++) + a.Coeffs[ii] /= k; + + return a; + } + + /// + /// Addition of two Polynomials (piecewise) + /// + /// left Polynomial + /// right Polynomial + /// resulting Polynomial + public static Polynomial operator +( Polynomial a, Polynomial b) + { + return add(a, b); + } + + /// + /// substraction of two Polynomials (piecewise) + /// + /// left Polynomial + /// right Polynomial + /// resulting Polynomial + public static Polynomial operator -( Polynomial a, Polynomial b) + { + return substract(a, b); + } + + /// + /// Calculates the complex roots of the Polynomial in the same way as matlab does + /// + /// a vector of complex numbers with the roots + public Complex[] GetRoots() + { + DenseMatrix A = this.GetEigValMatrix(); + Complex[] c_vec; + + if (A == null) + { + if (Coeffs.Length < 2) + { + var val = Coeffs.Length == 1 ? Coeffs[0] : Double.NaN; + c_vec = new Complex[1] { val }; + } + else + c_vec = new Complex[1] { new Complex(-Coeffs[0] / Coeffs[1], 0) }; + } + else + { + Evd eigen = A.Evd(Symmetricity.Asymmetric); + c_vec = eigen.EigenValues.ToArray(); + } + return c_vec; + } + + /// + /// get the eigenvalue matrix A of this Polynomial such that eig(A) = roots of this Polynomial + /// + /// Eigenvalue matrix A + public DenseMatrix GetEigValMatrix() + { + Polynomial pLoc = new Polynomial(this.Coeffs); + pLoc.CutTrailZeros(); + + int n = pLoc.Length - 1; + if (n < 2) + return null; + + double[] p = new double[n]; + + double a0 = pLoc.Coeffs[p.Length]; + + for (int ii = n - 1; ii >= 0; ii--) + p[ii] = -pLoc.Coeffs[ii] / a0; + + DenseMatrix A0 = DenseMatrix.CreateDiagonal(n - 1, n - 1, 1.0); + DenseMatrix A = new DenseMatrix(n); + + A.SetSubMatrix(1, 0, A0); + + A.SetRow(0, p.Reverse().ToArray()); + return A; + } + + /// + /// pointwise division of two Polynomials + /// + /// left Polynomial + /// right Polynomial + /// resulting Polynomial + public static Polynomial DividePointwise( Polynomial a, Polynomial b) + { + if (a.Length != b.Length) + mkSameLength(ref a, ref b); + + int n = a.Length; + double[] res = new double[a.Length]; + + + for (int ii = 0; ii < n; ii++) + { + res[ii] = a.Coeffs[ii] / b.Coeffs[ii]; + } + Polynomial res_poly = new Polynomial(res); + return (res_poly); + } + + /// + /// pointwise multiplication of two Polynomials + /// + /// left Polynomial + /// right Polynomial + /// resulting Polynomial + public static Polynomial MultiplyPointwise( Polynomial a, Polynomial b) + { + if (a.Length != b.Length) + mkSameLength(ref a, ref b); + + int n = a.Length; + double[] res = new double[a.Length]; + + + for (int ii = 0; ii < n; ii++) + { + res[ii] = a.Coeffs[ii] * b.Coeffs[ii]; + } + Polynomial res_poly = new Polynomial(res); + return (res_poly); + } + + /// + /// Addition of two Polynomials (piecewise) + /// + /// left Polynomial + /// right Polynomial + /// resulting Polynomial + public static Polynomial add( Polynomial a, Polynomial b) + { + + if (a.Length != b.Length) + mkSameLength(ref a, ref b); + + int n = a.Length; + double[] res = new double[a.Length]; + + + for (int ii = 0; ii < n; ii++) + { + res[ii] = a.Coeffs[ii] + b.Coeffs[ii]; + } + Polynomial res_poly = new Polynomial(res); + return (res_poly); + } + + /// + /// substraction of two Polynomials (piecewise) + /// + /// left Polynomial + /// right Polynomial + /// resulting Polynomial + public static Polynomial substract( Polynomial a, Polynomial b) + { + + if (a.Length != b.Length) + mkSameLength(ref a, ref b); + + int n = a.Length; + double[] res = new double[a.Length]; + + + for (int ii = 0; ii < n; ii++) + { + res[ii] = a.Coeffs[ii] - b.Coeffs[ii]; + } + Polynomial res_poly = new Polynomial(res); + return (res_poly); + } + + #endregion + + #region Displaying + /// + /// "0.00 x^3 + 0.00 x^2 + 0.00 x^1 + 0.00" like display of this Polynomial + /// + /// string in displayed format + public override string ToString() + { + string strLoc = ""; + + for (int ii = Length - 1; ii >= 0; ii--) + { + + if (ii == 0) + strLoc = String.Concat(strLoc, this.Coeffs[ii].ToString()); + else + strLoc = String.Concat(strLoc, this.Coeffs[ii].ToString(), VarName, ii.ToString(), " + "); + } + return strLoc; + } + + #endregion + + #region Interfacing + + /// + /// This method returns the coefficcients of the Polynomial as an array the "IsFlipped" property, + /// which is set during construction is taken into account automatically. + /// + /// the coefficcients of the Polynomial as an array + public double[] ToArray() + { + if (IsFlipped == true) + return (Coeffs.Reverse().ToArray()); + else + return (Coeffs); + } + + #endregion + + #region Helpers + + private static void mkSameLength(ref Polynomial a, ref Polynomial b) + { + double[] aHold = new double[a.Length]; + double[] bHold = new double[b.Length]; + Array.Copy(a.Coeffs, aHold, a.Length); + Array.Copy(b.Coeffs, bHold, b.Length); + + if (a.Length < b.Length) + { + a.Coeffs = new double[b.Length]; + b.Coeffs = new double[b.Length]; + Array.Copy(aHold, a.Coeffs, aHold.Length); + Array.Copy(bHold, b.Coeffs, bHold.Length); + } + else + { + a.Coeffs = new double[a.Length]; + b.Coeffs = new double[a.Length]; + Array.Copy(aHold, a.Coeffs, aHold.Length); + Array.Copy(bHold, b.Coeffs, bHold.Length); + } + + } + + /// + /// (full) convolution of two arrays + /// + /// left vector + /// right vector + /// convolution of a and b as vector + private static double[] conv(double[] a, double[] b) + { + double[] ret = new double[a.Length + b.Length]; + + for (int i = 0; i < a.Length; i++) + { + for (int j = 0; j < b.Length; j++) + { + ret[i + j] += a[i] * b[j]; + } + } + return ret; + } + #endregion + + } + +} From 5f7f7f0d30eb091c7529cc8c4a82b8142588085b Mon Sep 17 00:00:00 2001 From: Tobias Glaubach Date: Tue, 10 Jul 2018 21:46:36 +0200 Subject: [PATCH 02/18] implemented a prototype for the TransferFunctionDiscrete class --- .../LtiSystems/TransferFunctionDiscrete.cs | 968 ++++++++++++++++++ 1 file changed, 968 insertions(+) create mode 100644 src/Numerics/LtiSystems/TransferFunctionDiscrete.cs diff --git a/src/Numerics/LtiSystems/TransferFunctionDiscrete.cs b/src/Numerics/LtiSystems/TransferFunctionDiscrete.cs new file mode 100644 index 00000000..e2659a6e --- /dev/null +++ b/src/Numerics/LtiSystems/TransferFunctionDiscrete.cs @@ -0,0 +1,968 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Numerics; +using System.Text; +using MathNet.Numerics; + +namespace MathNet.Numerics.LtiSystems +{ + /// Class for LTI discrete transfer functions + public class TransferFunctionDiscrete + { + + private double[] _num; + + private double[] _den; + + /// + /// numberator (input dependent) Polynomial coefficients as array + /// in order + /// => index high ... index low + /// => [n], [n-1], +..., [0] + /// => 1 + q^-1 + ... + q^-n + /// + public double[] num + { + get + { + return _num; + } + set + { + _num = cutTrailingZeros(value); + shiftNumDenIfPossible(); + } + } + + /// den (state dependent) Polynomial coefficients as array + /// in order + /// => index high ... index low + /// => [n], [n-1], +..., [0] + /// => 1 + q^-1 + ... + q^-n + /// + public double[] den + { + get + { + return _den; + } + set + { + _den = cutTrailingZeros(value); + shiftNumDenIfPossible(); + } + } + + + + /// b (input dependent) Polynomial coefficients as array ( + public double[] b + { + get + { + return _num; + } + set + { + _num = cutTrailingZeros(value); + shiftNumDenIfPossible(); + } + } + + /// a (state dependent) Polynomial coefficients as array + public double[] a + { + get + { + return _den; + } + set + { + _den = cutTrailingZeros(value); + shiftNumDenIfPossible(); + } + } + + /// Internal FIR States -> updated in every response calculation + public double[] z_FIR { get; set; } + + /// Internal IIR States -> updated in every response calculation + public double[] z_IIR { get; set; } + + /// any name you want to give this transfer function + public string Name { get; set; } + + /// sampling time of discrete transfer function (default = 1) + public double Ts { get; set; } + + /// variable for transfer function so far all tf's are in the q^-1 (or equivalently z^-1) form. Changing this will have NO nfluence besides displaying te TF + public string variable = "q^-1"; + + + + /// + /// Check if this Transfer Function is stable + /// + /// the tolerance for euclidian distance at which a pole/zero pair is considered to be canceling each other + /// false if system is unsable true if system is stable + public bool IsStable(double numTolerance = 1e-8) + { + + var p = GetPoles(); + + var z = GetZeros(); + var z_isCompensated = new bool[z.Length]; + for (int j = 0; j < p.Length; j++) + { + // check if pole would lead to unstable behaviour + if (p[j].Magnitude > 1.0) + { + + // init some values + double minDistance = Double.PositiveInfinity; + int idxMinDistanceZero = -1; + + // analyze the distance between each zero and the pole now + for (int i = 0; i < z.Length; i++) + { + // check if pole has already been used for compensation + if (z_isCompensated[i]) + continue; + + // calculate geometrical distance between each zero and the pole now and store the closest neighbour + var dist = (p[j] - z[i]).Magnitude; + if (dist < minDistance) + { + minDistance = dist; + idxMinDistanceZero = i; + } + } + + // if closest neighbour is too far away to compensate the unstable pole + if (minDistance >= numTolerance) + return false; // the system is unsable + else + z_isCompensated[idxMinDistanceZero] = true; // if not: mark the zero as already used for compensation + + } + + } + + return true; + + } + + + + /// constructor setting no properties at all + public TransferFunctionDiscrete() + { + this.Ts = 1.0; + } + + /// + /// constructor setting a and b vectors as well as initializing the z_FIR and z_IIR states + /// + public TransferFunctionDiscrete(double b_in, double[] a_in, double Ts_in = 1.0d) + { + if (a_in == null) + throw new ArgumentNullException("a_in"); + + this.a = (double[])a_in.Clone(); + this.b = new double[1]; + this.b[0] = b_in; + this.z_IIR = new double[a_in.Length]; + this.z_FIR = new double[1]; + this.Ts = Ts_in; + } + + /// + /// constructor setting a and b vectors as well as initializing the z_FIR and z_IIR states + /// + public TransferFunctionDiscrete(double[] b_in, double a_in, double Ts_in = 1.0d) + { + if (b_in == null) + throw new ArgumentNullException("b_in"); + + this.a = new double[1]; + this.a[0] = a_in; + this.b = (double[])b_in.Clone(); + this.z_IIR = new double[1]; + this.z_FIR = new double[b_in.Length]; + this.Ts = Ts_in; + } + + /// + /// constructor setting a and b vectors as well as initializing the z_FIR and z_IIR states + /// + public TransferFunctionDiscrete(double b_in, double a_in, double Ts_in = 1.0d) + { + this.a = new double[1]; + this.a[0] = a_in; + this.b = new double[1]; + this.b[0] = b_in; + this.z_IIR = new double[1]; + this.z_FIR = new double[1]; + this.Ts = Ts_in; + } + + /// + /// constructor setting a and b vectors as well as initializing the z_FIR and z_IIR states + /// + public TransferFunctionDiscrete(double[] b_in, double[] a_in, double Ts_in = 1.0d) + { + if (b_in == null) + throw new ArgumentNullException("b_in"); + + if (a_in == null) + throw new ArgumentNullException("a_in"); + + this.a = (double[])a_in.Clone(); + this.b = (double[])b_in.Clone(); + this.z_IIR = new double[a_in.Length]; + this.z_FIR = new double[b_in.Length]; + this.Ts = Ts_in; + } + + + /// + /// Adds delay to the numerator array (shifting the values by d steps) + /// + /// integer value of daly to add to this TF + public void AddDelay(int d) + { + double[] b_new = new double[b.Length + d]; + b.CopyTo(b_new, d); + b = b_new; + } + + #region Helpers + /// + /// if num and den both start at a later step (e.G highest power num = q^-3 and highest power den = q^-4), the whole tf can be shifted by n (3) steps + /// + private void shiftNumDenIfPossible() + { + var offset = 0; + if (_num == null || _den == null) + return; + + var n = Math.Min(_num.Length, _den.Length); + for (int i = 0; i < n; i++) + { + if (num[i] == 0.0d && _den[i] == 0.0d) + offset = i + 1; + else + break; + } + + if (offset > 0) + { + double[] tmp1 = new double[_num.Length - offset]; + Array.Copy(_num, offset, tmp1, 0, tmp1.Length); + double[] tmp2 = new double[_den.Length - offset]; + Array.Copy(_den, offset, tmp2, 0, tmp2.Length); + _num = tmp1; + _den = tmp2; + } + } + + private double[] cutTrailingZeros(double[] vIn) + { + int lengthNew = vIn.Length; + + for (int i = vIn.Length - 1; i >= 0; i--) + { + if (vIn[i] != 0) + { + lengthNew = i + 1; + break; + } + } + + var v = new double[lengthNew]; + Array.Copy(vIn, v, lengthNew); + + return v; + } + + /// + /// checks and adjusts internal states to a and b arrays + /// + private void checkStateSizes() + { + if (this.z_IIR == null) + this.z_IIR = new double[this.a.Length]; + + if (this.z_FIR == null) + this.z_FIR = new double[this.b.Length]; + + if (this.a.Length != this.z_IIR.Length) + this.z_IIR = new double[this.a.Length]; + + if (this.b.Length != this.z_FIR.Length) + this.z_FIR = new double[this.b.Length]; + + if (this.a.Length != this.z_IIR.Length) + this.z_IIR = new double[this.a.Length]; + } + #endregion Helpers + + #region Operators + + /// + /// LTI System theory division of a transfer function object by a scalar + /// + /// transfer function + /// scalar for divison + /// new transfer function object divided by k + public static TransferFunctionDiscrete operator /(TransferFunctionDiscrete G1, double k) + { + Polynomial A1 = new Polynomial(G1.a); + + TransferFunctionDiscrete Gres = new TransferFunctionDiscrete(G1.b, (A1 * k).ToArray(), G1.Ts) + { + Name = G1.Name + }; + return Gres; + } + + /// + /// LTI System theory division of a scalar by a transfer function object + /// + /// scalar value + /// transfer function for division + /// new transfer function object + public static TransferFunctionDiscrete operator /(double k, TransferFunctionDiscrete G1) + { + Polynomial A1 = new Polynomial(G1.a); + + TransferFunctionDiscrete Gres = new TransferFunctionDiscrete((A1 * k).ToArray(), G1.b, G1.Ts) + { + Name = G1.Name + }; + return Gres; + } + + /// + /// LTI System theory multiplication of a transfer function object by a scalar + /// + /// transfer function + /// scalar for multiplication + /// new transfer function object + public static TransferFunctionDiscrete operator *(TransferFunctionDiscrete G1, double k) + { + Polynomial B1 = new Polynomial(G1.b); + + TransferFunctionDiscrete Gres = new TransferFunctionDiscrete((B1 * k).ToArray(), G1.a, G1.Ts) + { + Name = G1.Name + }; + return Gres; + + } + + /// + /// LTI System theory multiplication of a transfer function object by a scalar + /// + /// scalar for multiplication + /// transfer function + /// new transfer function object + public static TransferFunctionDiscrete operator *(double k, TransferFunctionDiscrete G1) + { + Polynomial B1 = new Polynomial(G1.b); + + TransferFunctionDiscrete Gres = new TransferFunctionDiscrete((B1 * k).ToArray(), G1.a, G1.Ts) + { + Name = G1.Name + }; + return Gres; + } + + + /// + /// LTI System theory substraction of a transfer function with a scalar + /// + /// transfer function + /// scalar + /// new transfer function object + public static TransferFunctionDiscrete operator -(TransferFunctionDiscrete G1, double k) + { + Polynomial A1 = new Polynomial(G1.a); + Polynomial B1 = new Polynomial(G1.b); + + Polynomial A_res = (A1); + Polynomial B_res = B1 - (A1 * k); + + TransferFunctionDiscrete Gres = new TransferFunctionDiscrete(B_res.ToArray(), A_res.ToArray(), G1.Ts) + { + Name = G1.Name + }; + return Gres; + } + + /// + /// LTI System theory substraction of a scalar by a transfer function + /// + /// scalar + /// transfer function + /// new transfer function object + public static TransferFunctionDiscrete operator -(double k, TransferFunctionDiscrete G1) + { + Polynomial A1 = new Polynomial(G1.a); + Polynomial B1 = new Polynomial(G1.b); + + Polynomial A_res = (A1); + Polynomial B_res = (A1 * k) - B1; + + TransferFunctionDiscrete Gres = new TransferFunctionDiscrete(B_res.ToArray(), A_res.ToArray(), G1.Ts) + { + Name = G1.Name + }; + return Gres; + } + + /// + /// LTI System theory addition of a transfer function with a scalar + /// + /// transfer function + /// scalar + /// new transfer function object + public static TransferFunctionDiscrete operator +(TransferFunctionDiscrete G1, double k) + { + Polynomial A1 = new Polynomial(G1.a); + Polynomial B1 = new Polynomial(G1.b); + + Polynomial A_res = (A1); + Polynomial B_res = (A1 * k) + B1; + + TransferFunctionDiscrete Gres = new TransferFunctionDiscrete(B_res.ToArray(), A_res.ToArray(), G1.Ts) + { + Name = G1.Name + }; + return Gres; + } + + /// + /// LTI System theory addition of a transfer function with a scalar + /// + /// transfer function + /// scalar + /// new transfer function object + public static TransferFunctionDiscrete operator +(double k, TransferFunctionDiscrete G1) + { + Polynomial A1 = new Polynomial(G1.a); + Polynomial B1 = new Polynomial(G1.b); + + Polynomial A_res = (A1); + Polynomial B_res = B1 + (A1 * k); + + TransferFunctionDiscrete Gres = new TransferFunctionDiscrete(B_res.ToArray(), A_res.ToArray(), G1.Ts) + { + Name = G1.Name + }; + return Gres; + } + + /// + /// LTI System theory addition of two transfer functions + /// + /// transfer function left + /// transfer function right + /// new transfer function object + public static TransferFunctionDiscrete operator +(TransferFunctionDiscrete G1, TransferFunctionDiscrete G2) + { + if (Math.Abs(G1.Ts - G2.Ts) > 1e-12) + throw new ArgumentException(String.Format("The two supplied transfer functions do not have equal sampling times. G1.Ts = {0} G2.Ts = {1}", G1.Ts, G2.Ts)); + + Polynomial A1 = new Polynomial(G1.a); + Polynomial B1 = new Polynomial(G1.b); + + Polynomial A2 = new Polynomial(G2.a); + Polynomial B2 = new Polynomial(G2.b); + + Polynomial A_res = (A1 * A2); + Polynomial B_res = (B1 * A2) + (B2 * A1); + + return new TransferFunctionDiscrete(B_res.ToArray(), A_res.ToArray(), G1.Ts); + } + + /// + /// LTI System theory substraction of two transfer functions + /// + /// transfer function left + /// transfer function right + /// new transfer function object + public static TransferFunctionDiscrete operator -(TransferFunctionDiscrete G1, TransferFunctionDiscrete G2) + { + if (Math.Abs(G1.Ts - G2.Ts) > 1e-12) + throw new ArgumentException(String.Format("The two supplied transfer functions do not have equal sampling times. G1.Ts = {0} G2.Ts = {1}", G1.Ts, G2.Ts)); + + Polynomial A1 = new Polynomial(G1.a); + Polynomial B1 = new Polynomial(G1.b); + + Polynomial A2 = new Polynomial(G2.a); + Polynomial B2 = new Polynomial(G2.b); + + Polynomial A_res = (A1 * A2); + Polynomial B_res = (B1 * A2) - (B2 * A1); + + return new TransferFunctionDiscrete(B_res.ToArray(), A_res.ToArray(), G1.Ts); + } + + /// + /// LTI System theory addition of two transfer functions + /// + /// transfer function left + /// transfer function right + /// new transfer function object + public static TransferFunctionDiscrete operator *(TransferFunctionDiscrete G1, TransferFunctionDiscrete G2) + { + if (Math.Abs(G1.Ts - G2.Ts) > 1e-12) + throw new ArgumentException(String.Format("The two supplied transfer functions do not have equal sampling times. G1.Ts = {0} G2.Ts = {1}", G1.Ts, G2.Ts)); + + Polynomial A1 = new Polynomial(G1.a); + Polynomial B1 = new Polynomial(G1.b); + + Polynomial A2 = new Polynomial(G2.a); + Polynomial B2 = new Polynomial(G2.b); + + Polynomial A_res = A1 * A2; + Polynomial B_res = B1 * B2; + + return new TransferFunctionDiscrete(B_res.ToArray(), A_res.ToArray(), G1.Ts); + } + + /// + /// LTI System theory division of two transfer functions + /// + /// transfer function left + /// transfer function right + /// new transfer function object + public static TransferFunctionDiscrete operator /(TransferFunctionDiscrete G1, TransferFunctionDiscrete G2) + { + if (Math.Abs(G1.Ts - G2.Ts) > 1e-12) + throw new ArgumentException(String.Format("The two supplied transfer functions do not have equal sampling times. G1.Ts = {0} G2.Ts = {1}", G1.Ts, G2.Ts)); + + Polynomial A1 = new Polynomial(G1.a); + Polynomial B1 = new Polynomial(G1.b); + + Polynomial A2 = new Polynomial(G2.a); + Polynomial B2 = new Polynomial(G2.b); + + Polynomial A_res = A1 * B2; + Polynomial B_res = B1 * A2; + + return new TransferFunctionDiscrete(B_res.ToArray(), A_res.ToArray(), G1.Ts); + } + #endregion + + /// calculates y_k = G(q^-1) * x_k for a given x_k array + public IEnumerable CalcResponse(IEnumerable x) + { + return (this.CalcResponse(this.b, this.a, x.ToArray())); + } + + /// calculates y_k = G(q^-1) * x_k for a given x_k array + public double[] CalcResponse(double[] x) + { + // this is basically a two step convoltion and could be replaced by a + // conv implementation. + // however... this code works and replacing it would be more work + + double y_now = 0.0d; + int idx_a = 0; + int idx_b = 0; + double[] y = new double[x.Length]; + + this.checkStateSizes(); + + // Loop all inputs + for (int ii_x = 0; ii_x < x.Length; ii_x++) + { + y_now = 0.0d; + idx_b = 0; + + // loop through b-matrix until end of momentary tempx-array + for (int ii_b = 0; ii_b <= ii_x && idx_b < b.Length; ii_b++) + { + + z_FIR[idx_b] = x[ii_x - ii_b]; + y_now += b[idx_b] * z_FIR[idx_b]; + idx_b++; + } + + + // start at second position, since it's the a-matrix + idx_a = 1; + // loop for a-matrix + for (int ii_a = 0; ii_a <= (ii_x - 1) && idx_a < a.Length; ii_a++) + { + z_IIR[idx_a] = y[(ii_x - 1) - ii_a]; + y_now -= a[idx_a] * z_IIR[idx_a]; + idx_a++; + } + // write result + y[ii_x] = (y_now / a[0]); + z_IIR[0] = y[ii_x]; + } + return (y); + } + + // Todo: Implement FiltFilt + /* + /// + /// A wrapper for the StaticFilters.FiltFilt method using the internal a and b arrays + /// + /// The data to filter + /// initial state coefficients null for aotomatic generation via steady state solution + /// the number of datapoints to pad at each side use less than 0 for Math.Max(a.Length, b.Length) * 3 + /// The filterd data + /// + /// In order to prevent transients at the end or start of the sequence we have to pad it + /// The padding is done by rotating the sequence by 180° at the ends and append it to the data + /// + public double[] FiltFilt(double[] data, double[] zi = null, int padlen = 0) + { + if (this.a == null || this.a.Length == 0) + throw new Exception("This transfer function has no a array with data"); + if (this.b == null || this.b.Length == 0) + throw new Exception("This transfer function has no a array with data"); + + return StaticFilters.FiltFilt(data, this.a, this.b, zi, padlen); + } + */ + + #region Dynamics + + /// + /// returns the impulse response with nSteps for the tf model + /// + /// number of steps for impulse response + /// + public double[] Impulse(int nSteps) + { + + var Inp = new double[nSteps]; + Inp[0] = 1.0; + + var ImpulseResponse = this.CalcResponse(Inp); + return (ImpulseResponse); + } + + /// + /// returns the impulse response with nSettling * 1.3 steps for the tf model + /// + public double[] Impulse() + { + + var nSteps = Convert.ToInt32((double)CalcSettlingSteps() * 1.3); + if (nSteps <= 0) + return null; + var Inp = new double[nSteps]; + Inp[0] = 1.0; + + var ImpulseResponse = this.CalcResponse(Inp); + return (ImpulseResponse); + } + + + public Complex[] Bode(int nPoints = 100) + { + // substituting z = exp(j * omega * Ts) + var omega_vec = Generate.LinearSpaced(nPoints, 0, 2 * Math.PI * 1 / Ts); + + return Bode(omega_vec); + } + + public Complex[] Bode(int nPoints, out double[] omega_vec) + { + // substituting z = exp(j * omega * Ts) + omega_vec = Generate.LinearSpaced(nPoints, 0, 2 * Math.PI * 1 / Ts); + + + return Bode(omega_vec); + } + + public Complex[] Bode(double[] omega_vec) + { + + var nPoints = omega_vec.Length; + + double omega; + double expVal; + Complex zVal; + Complex denVal; + Complex numVal; + + var bodeVal = new Complex[nPoints]; + + for (int idx = 0; idx < nPoints; idx++) + { + + + omega = omega_vec[idx]; + + zVal = new Complex(0.0, 0.0); + + denVal = new Complex(0.0, 0.0); + for (int ii = 0; ii < a.Length; ii++) + { + expVal = ii * omega * Ts; + zVal = new Complex(0.0, expVal); + + denVal += a[ii] * zVal.Exp(); + } + + numVal = new Complex(0.0, 0.0); + for (int ii = 0; ii < b.Length; ii++) + { + expVal = ii * omega * Ts; + zVal = new Complex(0.0, expVal); + + numVal += b[ii] * zVal.Exp(); + } + bodeVal[idx] = numVal / denVal; + } + + return bodeVal; + } + + /// The poles resulting from the denominator Polynomial root + public Complex[] GetPoles() + { + Polynomial a_poly = new Polynomial(a, isFlip:true); + Complex[] r = a_poly.GetRoots(); + return r; + } + + /// The zeros resulting from the nominator Polynomial root + public Complex[] GetZeros() + { + Polynomial b_poly = new Polynomial(b, isFlip:true); + Complex[] r = b_poly.GetRoots(); + return r; + } + + + /// + /// calculate the number of steps the system will need until it can be assumed to be settled + /// + /// tolerance in decimal percent at which to assume that the system is settled (default = 0.3) + /// maximum number of steps to simulate (default = 500000) + /// number of steps at which the system is assumed to be settled, or 0 if unstable + public int CalcSettlingSteps(double tol = 0.03, int n_max = 500000) + { + + // init settling time as zero for never settled + int n_sttl = 0; + + // if the system is unstable return zero since the system will never be settled + if (this.IsStable() == false) + return 0; + + int n_sim = 0; + + double[] dampVals = GetDampings(out double[] EigenFrequencys); + + double dampWorst = dampVals.Min(); + + //for (int ii = 1; ii < dampVals.Length; ii++) + // dampWorst = dampWorst * dampVals[ii]; + + double t_simFull; + + // appromate a settling time based on damping + var t_stlDamp = -Math.Log(tol) / dampWorst; + + // approximate a settling time from time constants + var tau = new double[dampVals.Length]; + for (int ii = 0; ii < tau.Length; ii++) + tau[ii] = 1.0 / (dampVals[ii] * EigenFrequencys[ii]); + + // approx after 5 * biggest time constant + var t_stlTimeConst = tau.Max() * 5; + + // choose bigger approximation + t_simFull = Math.Max(t_stlTimeConst, t_stlDamp); + + // recalculate to number of steps + int nStepsBase = (int)Math.Ceiling(t_simFull / Ts); + + + // simulate impulse responses with n*10*nStepsBase time steps + // incrementing n if necessary until steady state is reached + n_sim = nStepsBase <= 0 ? 5 : nStepsBase; + int count = 0; + while (count < 10 && n_sttl == 0) + { + if (n_sim > n_max) + return n_sttl; + + n_sim = 10 * n_sim; + + double[] dirac_sim = new double[n_sim]; + dirac_sim[0] = 1.0; + var tmp_outp = this.CalcResponse(dirac_sim); + + int idxPos = n_sim - 1; + + // find first step beeing bigger than tolerance + while (idxPos > 0 && n_sttl == 0) + { + if (tmp_outp[idxPos] > tol) + n_sttl = idxPos; + + idxPos--; + } + count++; + } + return n_sttl; + + } + + #endregion Dynamics + + + #region Dampings + /// + /// gets the damping coefficients from this transfer function, + /// since all transfer functions so far are discrete time, + /// these values do not directly translate to lambda. + /// the theoretical recalculation is: + /// Z = -cos(angle(log(lambda))) + /// + /// Array of damping values for this transfer function + public double[] GetDampings() + { + return GetDampings(out double[] f); + } + + /// + /// gets the damping coefficients from this transfer function, + /// since all transfer functions so far are discrete time, + /// these values do not directly translate to lambda. + /// the theoretical recalculation is: + /// Z = -cos(angle(log(lambda))) + /// + /// Array of damping values for this transfer function + public double[] GetDampings(out double[] wn) + { + + var r = GetPoles().Clone() as Complex[]; + var s = new Complex[r.Length]; + var f = new double[r.Length]; + var z = new double[r.Length]; + + for (int idx = 0; idx < r.Length; idx++) + { + s[idx] = Complex.Log(r[idx]) / Ts; + f[idx] = s[idx].Magnitude; + z[idx] = -s[idx].Real / f[idx]; + } + + wn = (double[])f.Clone(); + + return z; + + } + + + #endregion Dampings + + + #region displaying + /// + /// + /// + /// + public string DispTF() + { + return (DispTF(this.b, this.a, this.Name, this.variable.Substring(0, variable.Length - 1))); + } + + public string NumString() + { + var varStr = this.variable.Substring(0, variable.Length - 1); + var num = b.Clone() as double[]; + return getFractString(num, varStr); + } + + public string DenString() + { + var varStr = this.variable.Substring(0, variable.Length - 1); + var den = a.Clone() as double[]; + return getFractString(den, varStr); + } + + + private static string getFractString(double[] num, string varStr) + { + string str1; + string str2; + string strNum = ""; + for (int item = 0; item < num.Length; item++) + { + if (num[item] == 0) + continue; + + //str2 = Math.Abs(num[item]).ToString(); + str2 = Math.Abs(num[item]).ToString("0.######"); + if (item == 0) + { + if (num[item] < 0) + str1 = "-"; + else + str1 = ""; + + strNum = String.Concat(strNum, str1, str2); + } + else + { + if (num[item] > 0) + str1 = " + "; + else + str1 = " - "; + strNum = String.Concat(strNum, str1, str2, varStr, item.ToString()); + } + } + + if (strNum.StartsWith("+") || strNum.StartsWith(" ")) + strNum = strNum.Substring(1); + + return strNum; + } + + + + /// + /// + /// + /// + /// + /// + /// + public static string DispTF(double[] num, double[] den, string name, string varStr = " q^-") + { + + string strNum = getFractString(num, varStr); + string strDen = getFractString(den, varStr); + string strHead = ""; + + if (String.IsNullOrEmpty(name)) + strHead = "TF = "; + else + strHead = name; + + int nbar = Math.Max(strDen.Length, strNum.Length); + + string strBar = new String('-', nbar); + string strOut = String.Concat(strHead, "\n\n", strNum, '\n', strBar, '\n', strDen); + + return (strOut); + } + + #endregion displaying + + } +} From 7754ba52d2c5c25cec5340fbcb5db1836e47b59a Mon Sep 17 00:00:00 2001 From: Tobias Glaubach Date: Tue, 10 Jul 2018 22:47:29 +0200 Subject: [PATCH 03/18] added tests for the Polynomial class --- src/Numerics.Tests/PolynomialTests.cs | 152 ++++++++++++++++++++++++++ 1 file changed, 152 insertions(+) create mode 100644 src/Numerics.Tests/PolynomialTests.cs diff --git a/src/Numerics.Tests/PolynomialTests.cs b/src/Numerics.Tests/PolynomialTests.cs new file mode 100644 index 00000000..0adf3ef6 --- /dev/null +++ b/src/Numerics.Tests/PolynomialTests.cs @@ -0,0 +1,152 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// +// Copyright (c) 2009-2018 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using System; +using System.Linq; +using MathNet.Numerics; +using MathNet.Numerics.LinearRegression; +using MathNet.Numerics.Statistics; +using NUnit.Framework; + +namespace MathNet.Numerics.UnitTests +{ + [TestFixture, Category("Calculus")] + public class PolynomialTests + { + [Test] + public void AddTest() + { + for (int i = 0; i < 5; i++) + { + for (int j = 0; j < 5; j++) + { + var msg = String.Format("At i={0}, j={1}", i, j); + var n = Math.Max(i, j) + 1; + var tgt = new double[n]; + tgt[i] += 1; + tgt[j] += 1; + + var c1 = new double[i + 1]; + var c2 = new double[j + 1]; + + c1[i] = 1.0; + c2[j] = 1.0; + + var p1 = new Polynomial(c1); + var p2 = new Polynomial(c2); + + var p_res = Polynomial.add(p1, p2); + var p_tar = new Polynomial(tgt); + + p_res.CutTrailZeros(); + p_tar.CutTrailZeros(); + + Assert.AreEqual(p_tar.Length, p_res.Length); + for (int k = 0; k < p_res.Length; k++) + { + Assert.AreEqual(p_tar.Coeffs[k], p_res.Coeffs[k], msg); + } + } + } + } + + [Test] + public void SubstractTest() + { + for (int i = 0; i < 5; i++) + { + for (int j = 0; j < 5; j++) + { + var msg = String.Format("At i={0}, j={1}", i, j); + var n = Math.Max(i, j) + 1; + var tgt = new double[n]; + tgt[i] += 1; + tgt[j] -= 1; + + var c1 = new double[i + 1]; + var c2 = new double[j + 1]; + + c1[i] = 1.0; + c2[j] = 1.0; + + var p1 = new Polynomial(c1); + var p2 = new Polynomial(c2); + + var p_res = Polynomial.substract(p1, p2); + var p_tar = new Polynomial(tgt); + + p_res.CutTrailZeros(); + p_tar.CutTrailZeros(); + + Assert.AreEqual(p_tar.Length, p_res.Length); + for (int k = 0; k < p_res.Length; k++) + { + Assert.AreEqual(p_tar.Coeffs[k], p_res.Coeffs[k], msg); + } + } + } + } + + [Test] + public void MultiplyTest() + { + for (int i = 0; i < 5; i++) + { + for (int j = 0; j < 5; j++) + { + var msg = String.Format("At i={0}, j={1}", i, j); + var n = i + j + 1; + var tgt = new double[n]; + tgt[i + j] += 1; + + var c1 = new double[i + 1]; + var c2 = new double[j + 1]; + + c1[i] = 1.0; + c2[j] = 1.0; + + var p1 = new Polynomial(c1); + var p2 = new Polynomial(c2); + + var p_res = p1 * p2; + var p_tar = new Polynomial(tgt); + + p_res.CutTrailZeros(); + p_tar.CutTrailZeros(); + + Assert.AreEqual(p_tar.Length, p_res.Length); + for (int k = 0; k < p_res.Length; k++) + { + Assert.AreEqual(p_tar.Coeffs[k], p_res.Coeffs[k], msg); + } + } + } + } + } +} From 8e4e4ca988b9445860f12d6a28e1638be1e0a1f3 Mon Sep 17 00:00:00 2001 From: Tobias Glaubach Date: Tue, 10 Jul 2018 22:47:47 +0200 Subject: [PATCH 04/18] typo fix and small bugfix --- src/Numerics/LtiSystems/TransferFunctionDiscrete.cs | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/Numerics/LtiSystems/TransferFunctionDiscrete.cs b/src/Numerics/LtiSystems/TransferFunctionDiscrete.cs index e2659a6e..c050eca8 100644 --- a/src/Numerics/LtiSystems/TransferFunctionDiscrete.cs +++ b/src/Numerics/LtiSystems/TransferFunctionDiscrete.cs @@ -560,15 +560,15 @@ namespace MathNet.Numerics.LtiSystems /// calculates y_k = G(q^-1) * x_k for a given x_k array public IEnumerable CalcResponse(IEnumerable x) { - return (this.CalcResponse(this.b, this.a, x.ToArray())); + return this.CalcResponse(x.ToArray()); } /// calculates y_k = G(q^-1) * x_k for a given x_k array public double[] CalcResponse(double[] x) { - // this is basically a two step convoltion and could be replaced by a + // this is basically a two step convolution and could be replaced by a // conv implementation. - // however... this code works and replacing it would be more work + // however... this code works fine and replacing it would be more work double y_now = 0.0d; int idx_a = 0; From 070076eda37e38915409efb500e0c947f2d7900a Mon Sep 17 00:00:00 2001 From: Tobias Glaubach Date: Wed, 11 Jul 2018 10:44:20 +0200 Subject: [PATCH 05/18] gooming, bugfixes and error handling --- src/Numerics/Polynomial.cs | 158 ++++++++++++++++++++++++++++--------- 1 file changed, 122 insertions(+), 36 deletions(-) diff --git a/src/Numerics/Polynomial.cs b/src/Numerics/Polynomial.cs index 360cfb27..841280c2 100644 --- a/src/Numerics/Polynomial.cs +++ b/src/Numerics/Polynomial.cs @@ -2,10 +2,6 @@ using System.Collections.Generic; using System.Linq; using System.Text; -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; using System.Threading.Tasks; using System.Numerics; using MathNet.Numerics; @@ -44,7 +40,7 @@ namespace MathNet.Numerics { get { - return (Coeffs.Length); + return (Coeffs == null ? 0 : Coeffs.Length); } } @@ -54,9 +50,40 @@ namespace MathNet.Numerics /// size of Polynomial public Polynomial(int n) { + if (n < 0) + { + throw new ArgumentOutOfRangeException("n must be postive"); + } Coeffs = new double[n]; } + + /// + /// make Polynomial: e.G 3.0 = 3.0 + 0 x^1 + 0 x^2 + /// + /// just the "x^0" part + public Polynomial(double coeff) + { + IsFlipped = false; + this.Coeffs = new double[1]; + Coeffs[0] = coeff; + } + + /// + /// make Polynomial: e.G new double[] {5, 0, 2} = "5 + 0 x^1 + 2 x^2" + /// + /// Polynomial coefficiens as array + public Polynomial(double[] coeffs) + { + if (coeffs == null) + { + throw new ArgumentNullException("coeffs"); + } + this.Coeffs = new double[coeffs.Length]; + Array.Copy(coeffs, this.Coeffs, coeffs.Length); + } + + /// /// constructor setting Polynomial coefficiens and flipping them if necessary. /// @@ -68,41 +95,23 @@ namespace MathNet.Numerics /// xP1: 5 x^3 + 4 x^2 + 3 x^2 + 0 x^1 + 2 /// xP2: 2 x^3 + 0 x^2 + 3 x^2 + 4 x^1 + 5 /// - /// WARNING cut all trailing zeros before, since they would result in zeros at the end /// /// Polynomial coefficiens as array /// use true for flipping - public Polynomial(double[] coeffs, bool isFlip = false) + public Polynomial(double[] coeffs, bool isFlip) { - this.Coeffs = new double[Coeffs.Length]; + if (coeffs == null) + { + throw new ArgumentNullException("coeffs"); + } + this.Coeffs = new double[coeffs.Length]; Array.Copy(coeffs, Coeffs, coeffs.Length); if (isFlip) { Coeffs = Coeffs.Reverse().ToArray(); IsFlipped = true; - } - } - - /// - /// constructor setting Polynomial coefficiens - /// - /// just the x^0 part - public Polynomial(double coeff) - { - IsFlipped = false; - this.Coeffs = new double[1]; - Coeffs[0] = coeff; - } - - /// - /// constructor setting Polynomial coefficiens - /// - /// Polynomial coefficiens as array - public Polynomial(double[] Coeffs) - { - this.Coeffs = new double[Coeffs.Length]; - Array.Copy(Coeffs, this.Coeffs, Coeffs.Length); + } } /// @@ -129,9 +138,45 @@ namespace MathNet.Numerics } } + #region diff/int + public Polynomial Differentiate() + { + + if (Coeffs.Length == 0) + { + return null; + } + + var t = this.Clone() as Polynomial; + t.CutTrailZeros(); + var cNew = new double[t.Length - 1]; + for (int i = 1; i < t.Coeffs.Length; i++) + { + cNew[i-1] = t.Coeffs[i] * i; + } + var p = new Polynomial(cNew, isFlip: IsFlipped); + p.CutTrailZeros(); + return p; + } + + public Polynomial Integrate() + { + var t = this.Clone() as Polynomial; + t.CutTrailZeros(); + var cNew = new double[t.Length + 1]; + for (int i = 1; i < cNew.Length; i++) + { + cNew[i] = t.Coeffs[i-1] / i; + } + var p = new Polynomial(cNew, isFlip: IsFlipped); + p.CutTrailZeros(); + return p; + } + + #endregion #region Operators - + /// /// multiplies a Polynomial by a Polynomial using convolution [ASINCO.libs.subfun.conv(a.Coeffs, b.Coeffs)] /// @@ -386,20 +431,54 @@ namespace MathNet.Numerics /// /// string in displayed format public override string ToString() + { + return ToString(highestFirst:false); + } + + /// + /// "0.00 x^3 + 0.00 x^2 + 0.00 x^1 + 0.00" like display of this Polynomial + /// + /// string in displayed format + public string ToString(bool highestFirst) { string strLoc = ""; + if (this.Coeffs == null) + { + return "null"; + } + if (this.Coeffs.Length == 0) + { + return ""; + } - for (int ii = Length - 1; ii >= 0; ii--) + if (!highestFirst) { + for (int ii = 0; ii < Length; ii++) + { - if (ii == 0) - strLoc = String.Concat(strLoc, this.Coeffs[ii].ToString()); - else - strLoc = String.Concat(strLoc, this.Coeffs[ii].ToString(), VarName, ii.ToString(), " + "); + if (ii == 0) + strLoc += String.Format("{0} + ", this.Coeffs[ii], VarName, ii); + else if (ii == Length - 1) + strLoc += String.Format("{0}{1}{2}", this.Coeffs[ii], VarName, ii); + else + strLoc += String.Format("{0}{1}{2} + ", this.Coeffs[ii], VarName, ii); + } } + else + { + for (int ii = Length - 1; ii >= 0; ii--) + { + if (ii == 0) + strLoc += this.Coeffs[ii].ToString(); + else + strLoc += String.Format("{0}{1}{2} + ", this.Coeffs[ii], VarName, ii); + } + } + return strLoc; } + #endregion #region Interfacing @@ -464,6 +543,13 @@ namespace MathNet.Numerics } return ret; } + + public object Clone() + { + var p = new double[this.Length]; + Array.Copy(Coeffs, p, Length); + return new Polynomial(p, isFlip: IsFlipped); + } #endregion } From 45c14fb557b7416a63e40e7dac307ec1d29d6460 Mon Sep 17 00:00:00 2001 From: Tobias Glaubach Date: Wed, 11 Jul 2018 10:44:42 +0200 Subject: [PATCH 06/18] adjusted the tets for polynomial --- src/Numerics.Tests/PolynomialTests.cs | 63 +++++++++++++++++++++++++-- 1 file changed, 60 insertions(+), 3 deletions(-) diff --git a/src/Numerics.Tests/PolynomialTests.cs b/src/Numerics.Tests/PolynomialTests.cs index 0adf3ef6..0d1dfee1 100644 --- a/src/Numerics.Tests/PolynomialTests.cs +++ b/src/Numerics.Tests/PolynomialTests.cs @@ -39,6 +39,63 @@ namespace MathNet.Numerics.UnitTests [TestFixture, Category("Calculus")] public class PolynomialTests { + + [TestCase(new double[] { 5, 4, 3, 0, 2 }, "5 + 4x^1 + 3x^2 + 0x^3 + 2x^4")] + [TestCase(new double[0], "")] + [TestCase(new double[] { 0, 4, 3, 0, 0 }, "0 + 4x^1 + 3x^2 + 0x^3 + 0x^4")] + public void ToStringTest(double[] x, string expected) + { + var p = new Polynomial(x); + Assert.AreEqual(expected, p.ToString()); + } + + [TestCase(new double[] { 5, 4, 3, 0, 2 }, new double[] { 4*1, 3*2, 0*3, 2*4 })] + [TestCase(new double[0], null)] + [TestCase(new double[] { 0, 4, 3, 0, 0 }, new double[] { 4*1, 3*2 })] + public void DifferentiateTest(double[] x, double[] expected) + { + var p = new Polynomial(x); + var p_res = p.Differentiate(); + + if (expected == null) + { + Assert.IsNull(p_res); + return; + } + else + { + Assert.AreEqual(expected.Length, p_res.Length, "length mismatch"); + for (int k = 0; k < p_res.Length; k++) + { + Assert.AreEqual(expected[k], p_res.Coeffs[k], "idx: " + k + " mismatch"); + } + } + + } + + [TestCase(new double[] { 5, 4, 3, 0, 2 }, new double[] { 0, 5.0/1.0, 4.0/2.0, 3.0/3.0, 0.0/4.0, 2.0/5.0 })] + [TestCase(new double[0], new double[1] { 0 })] + [TestCase(new double[] { 0, 1, 6, 8 }, new double[] {0, 0.0/1.0, 1.0/2.0, 6.0/3.0, 8.0/4.0})] + public void IntegrateTest(double[] x, double[] expected) + { + var p = new Polynomial(x); + var p_res = p.Integrate(); + + if (expected == null) + { + Assert.IsNull(p_res); + return; + } + else + { + Assert.AreEqual(expected.Length, p_res.Length, "length mismatch"); + for (int k = 0; k < p_res.Length; k++) + { + Assert.AreEqual(expected[k], p_res.Coeffs[k], "idx: " + k + " mismatch"); + } + } + } + [Test] public void AddTest() { @@ -67,7 +124,7 @@ namespace MathNet.Numerics.UnitTests p_res.CutTrailZeros(); p_tar.CutTrailZeros(); - Assert.AreEqual(p_tar.Length, p_res.Length); + Assert.AreEqual(p_tar.Length, p_res.Length, "length mismatch"); for (int k = 0; k < p_res.Length; k++) { Assert.AreEqual(p_tar.Coeffs[k], p_res.Coeffs[k], msg); @@ -104,7 +161,7 @@ namespace MathNet.Numerics.UnitTests p_res.CutTrailZeros(); p_tar.CutTrailZeros(); - Assert.AreEqual(p_tar.Length, p_res.Length); + Assert.AreEqual(p_tar.Length, p_res.Length, "length mismatch"); for (int k = 0; k < p_res.Length; k++) { Assert.AreEqual(p_tar.Coeffs[k], p_res.Coeffs[k], msg); @@ -140,7 +197,7 @@ namespace MathNet.Numerics.UnitTests p_res.CutTrailZeros(); p_tar.CutTrailZeros(); - Assert.AreEqual(p_tar.Length, p_res.Length); + Assert.AreEqual(p_tar.Length, p_res.Length, "length mismatch"); for (int k = 0; k < p_res.Length; k++) { Assert.AreEqual(p_tar.Coeffs[k], p_res.Coeffs[k], msg); From 659477794a65e054ac939eae1145f3c70e8ca755 Mon Sep 17 00:00:00 2001 From: Tobias Glaubach Date: Wed, 11 Jul 2018 16:32:35 +0200 Subject: [PATCH 07/18] renamed "Length" to "Degree" and added the methods "Fit" and "Evaluate" --- src/Numerics.Tests/PolynomialTests.cs | 25 +++--- src/Numerics/Polynomial.cs | 124 +++++++++++++++++--------- 2 files changed, 94 insertions(+), 55 deletions(-) diff --git a/src/Numerics.Tests/PolynomialTests.cs b/src/Numerics.Tests/PolynomialTests.cs index 0d1dfee1..71cba7fb 100644 --- a/src/Numerics.Tests/PolynomialTests.cs +++ b/src/Numerics.Tests/PolynomialTests.cs @@ -64,8 +64,8 @@ namespace MathNet.Numerics.UnitTests } else { - Assert.AreEqual(expected.Length, p_res.Length, "length mismatch"); - for (int k = 0; k < p_res.Length; k++) + Assert.AreEqual(expected.Length, p_res.Degree, "length mismatch"); + for (int k = 0; k < p_res.Degree; k++) { Assert.AreEqual(expected[k], p_res.Coeffs[k], "idx: " + k + " mismatch"); } @@ -88,8 +88,8 @@ namespace MathNet.Numerics.UnitTests } else { - Assert.AreEqual(expected.Length, p_res.Length, "length mismatch"); - for (int k = 0; k < p_res.Length; k++) + Assert.AreEqual(expected.Length, p_res.Degree, "length mismatch"); + for (int k = 0; k < p_res.Degree; k++) { Assert.AreEqual(expected[k], p_res.Coeffs[k], "idx: " + k + " mismatch"); } @@ -118,14 +118,14 @@ namespace MathNet.Numerics.UnitTests var p1 = new Polynomial(c1); var p2 = new Polynomial(c2); - var p_res = Polynomial.add(p1, p2); + var p_res = Polynomial.Add(p1, p2); var p_tar = new Polynomial(tgt); p_res.CutTrailZeros(); p_tar.CutTrailZeros(); - Assert.AreEqual(p_tar.Length, p_res.Length, "length mismatch"); - for (int k = 0; k < p_res.Length; k++) + Assert.AreEqual(p_tar.Degree, p_res.Degree, "length mismatch"); + for (int k = 0; k < p_res.Degree; k++) { Assert.AreEqual(p_tar.Coeffs[k], p_res.Coeffs[k], msg); } @@ -155,14 +155,14 @@ namespace MathNet.Numerics.UnitTests var p1 = new Polynomial(c1); var p2 = new Polynomial(c2); - var p_res = Polynomial.substract(p1, p2); + var p_res = Polynomial.Substract(p1, p2); var p_tar = new Polynomial(tgt); p_res.CutTrailZeros(); p_tar.CutTrailZeros(); - Assert.AreEqual(p_tar.Length, p_res.Length, "length mismatch"); - for (int k = 0; k < p_res.Length; k++) + Assert.AreEqual(p_tar.Degree, p_res.Degree, "length mismatch"); + for (int k = 0; k < p_res.Degree; k++) { Assert.AreEqual(p_tar.Coeffs[k], p_res.Coeffs[k], msg); } @@ -197,13 +197,14 @@ namespace MathNet.Numerics.UnitTests p_res.CutTrailZeros(); p_tar.CutTrailZeros(); - Assert.AreEqual(p_tar.Length, p_res.Length, "length mismatch"); - for (int k = 0; k < p_res.Length; k++) + Assert.AreEqual(p_tar.Degree, p_res.Degree, "length mismatch"); + for (int k = 0; k < p_res.Degree; k++) { Assert.AreEqual(p_tar.Coeffs[k], p_res.Coeffs[k], msg); } } } } + } } diff --git a/src/Numerics/Polynomial.cs b/src/Numerics/Polynomial.cs index 841280c2..2aad4e66 100644 --- a/src/Numerics/Polynomial.cs +++ b/src/Numerics/Polynomial.cs @@ -9,10 +9,9 @@ using MathNet.Numerics.LinearAlgebra; using MathNet.Numerics.LinearAlgebra.Double; using MathNet.Numerics.Statistics; using MathNet.Numerics.IntegralTransforms; - +using MathNet.Numerics.LinearRegression; using MathNet.Numerics.LinearAlgebra.Factorization; - namespace MathNet.Numerics { /// @@ -36,7 +35,7 @@ namespace MathNet.Numerics /// /// Length of Polynomial (max element + 1) e.G x^5 highest element, will give Length = 6 /// - public int Length + public int Degree { get { @@ -120,7 +119,7 @@ namespace MathNet.Numerics public void CutTrailZeros() { int count = 0; - for (int ii = Length - 1; ii >= 0; ii--) + for (int ii = Degree - 1; ii >= 0; ii--) { if (Coeffs[ii] == 0.0) { @@ -128,16 +127,54 @@ namespace MathNet.Numerics } else { - double[] CoeffsHold = new double[Length]; + double[] CoeffsHold = new double[Coeffs.Length]; Coeffs.CopyTo(CoeffsHold, 0); - Array.Resize(ref CoeffsHold, Length - count); - Coeffs = new double[Length - count]; + Array.Resize(ref CoeffsHold, Degree - count); + Coeffs = new double[Degree - count]; CoeffsHold.CopyTo(Coeffs, 0); return; } } } + #region Data Interaction + + /// + /// Least-Squares fitting the points (x,y) to a k-order polynomial y : x -> p0 + p1*x + p2*x^2 + ... + pk*x^k, + /// returning its best fitting parameters as [p0, p1, p2, ..., pk] array, compatible with Evaluate.Polynomial. + /// A polynomial with order/degree k has (k+1) coefficients and thus requires at least (k+1) samples. + /// + public static Polynomial Fit(double[] x, double[] y, int order, DirectRegressionMethod method = DirectRegressionMethod.QR) + { + var pArr = MathNet.Numerics.Fit.Polynomial(x, y, order, method); + return new Polynomial(pArr); + } + + /// + /// Evaluate a polynomial at point x. + /// + /// The location where to evaluate the polynomial at. + public double Evaluate(double z) + { + return MathNet.Numerics.Evaluate.Polynomial(z, Coeffs); + } + + /// + /// Evaluate a polynomial at points z. + /// + /// The locations where to evaluate the polynomial at. + public IEnumerable Evaluate(IEnumerable z) + { + var Lst = new List(); + foreach (var item in z) + { + Lst.Add(Evaluate(item)); + } + return Lst; + } + + #endregion + #region diff/int public Polynomial Differentiate() { @@ -149,7 +186,7 @@ namespace MathNet.Numerics var t = this.Clone() as Polynomial; t.CutTrailZeros(); - var cNew = new double[t.Length - 1]; + var cNew = new double[t.Coeffs.Length - 1]; for (int i = 1; i < t.Coeffs.Length; i++) { cNew[i-1] = t.Coeffs[i] * i; @@ -163,7 +200,7 @@ namespace MathNet.Numerics { var t = this.Clone() as Polynomial; t.CutTrailZeros(); - var cNew = new double[t.Length + 1]; + var cNew = new double[t.Coeffs.Length + 1]; for (int i = 1; i < cNew.Length; i++) { cNew[i] = t.Coeffs[i-1] / i; @@ -174,6 +211,7 @@ namespace MathNet.Numerics } #endregion + #region Operators @@ -206,7 +244,7 @@ namespace MathNet.Numerics /// resulting Polynomial public static Polynomial operator *( Polynomial a, double k) { - for (int ii = 0; ii < a.Length; ii++) + for (int ii = 0; ii < a.Coeffs.Length; ii++) a.Coeffs[ii] *= k; return a; @@ -245,7 +283,7 @@ namespace MathNet.Numerics /// resulting Polynomial public static Polynomial operator /( Polynomial a, double k) { - for (int ii = 0; ii < a.Length; ii++) + for (int ii = 0; ii < a.Coeffs.Length; ii++) a.Coeffs[ii] /= k; return a; @@ -259,7 +297,7 @@ namespace MathNet.Numerics /// resulting Polynomial public static Polynomial operator +( Polynomial a, Polynomial b) { - return add(a, b); + return Add(a, b); } /// @@ -270,7 +308,7 @@ namespace MathNet.Numerics /// resulting Polynomial public static Polynomial operator -( Polynomial a, Polynomial b) { - return substract(a, b); + return Substract(a, b); } /// @@ -309,13 +347,13 @@ namespace MathNet.Numerics Polynomial pLoc = new Polynomial(this.Coeffs); pLoc.CutTrailZeros(); - int n = pLoc.Length - 1; + int n = pLoc.Coeffs.Length - 1; if (n < 2) return null; double[] p = new double[n]; - double a0 = pLoc.Coeffs[p.Length]; + double a0 = pLoc.Coeffs[n]; for (int ii = n - 1; ii >= 0; ii--) p[ii] = -pLoc.Coeffs[ii] / a0; @@ -337,11 +375,11 @@ namespace MathNet.Numerics /// resulting Polynomial public static Polynomial DividePointwise( Polynomial a, Polynomial b) { - if (a.Length != b.Length) + if (a.Coeffs.Length != b.Coeffs.Length) mkSameLength(ref a, ref b); - int n = a.Length; - double[] res = new double[a.Length]; + int n = a.Coeffs.Length; + double[] res = new double[a.Coeffs.Length]; for (int ii = 0; ii < n; ii++) @@ -360,11 +398,11 @@ namespace MathNet.Numerics /// resulting Polynomial public static Polynomial MultiplyPointwise( Polynomial a, Polynomial b) { - if (a.Length != b.Length) + if (a.Coeffs.Length != b.Coeffs.Length) mkSameLength(ref a, ref b); - int n = a.Length; - double[] res = new double[a.Length]; + int n = a.Coeffs.Length; + double[] res = new double[a.Coeffs.Length]; for (int ii = 0; ii < n; ii++) @@ -381,14 +419,14 @@ namespace MathNet.Numerics /// left Polynomial /// right Polynomial /// resulting Polynomial - public static Polynomial add( Polynomial a, Polynomial b) + public static Polynomial Add( Polynomial a, Polynomial b) { - if (a.Length != b.Length) + if (a.Degree != b.Degree) mkSameLength(ref a, ref b); - int n = a.Length; - double[] res = new double[a.Length]; + int n = a.Degree; + double[] res = new double[n]; for (int ii = 0; ii < n; ii++) @@ -405,14 +443,14 @@ namespace MathNet.Numerics /// left Polynomial /// right Polynomial /// resulting Polynomial - public static Polynomial substract( Polynomial a, Polynomial b) + public static Polynomial Substract( Polynomial a, Polynomial b) { - if (a.Length != b.Length) + if (a.Degree != b.Degree) mkSameLength(ref a, ref b); - int n = a.Length; - double[] res = new double[a.Length]; + int n = a.Degree; + double[] res = new double[n]; for (int ii = 0; ii < n; ii++) @@ -453,12 +491,12 @@ namespace MathNet.Numerics if (!highestFirst) { - for (int ii = 0; ii < Length; ii++) + for (int ii = 0; ii < Coeffs.Length; ii++) { if (ii == 0) strLoc += String.Format("{0} + ", this.Coeffs[ii], VarName, ii); - else if (ii == Length - 1) + else if (ii == Coeffs.Length - 1) strLoc += String.Format("{0}{1}{2}", this.Coeffs[ii], VarName, ii); else strLoc += String.Format("{0}{1}{2} + ", this.Coeffs[ii], VarName, ii); @@ -466,7 +504,7 @@ namespace MathNet.Numerics } else { - for (int ii = Length - 1; ii >= 0; ii--) + for (int ii = Coeffs.Length - 1; ii >= 0; ii--) { if (ii == 0) strLoc += this.Coeffs[ii].ToString(); @@ -502,22 +540,22 @@ namespace MathNet.Numerics private static void mkSameLength(ref Polynomial a, ref Polynomial b) { - double[] aHold = new double[a.Length]; - double[] bHold = new double[b.Length]; - Array.Copy(a.Coeffs, aHold, a.Length); - Array.Copy(b.Coeffs, bHold, b.Length); + double[] aHold = new double[a.Coeffs.Length]; + double[] bHold = new double[b.Coeffs.Length]; + Array.Copy(a.Coeffs, aHold, a.Coeffs.Length); + Array.Copy(b.Coeffs, bHold, b.Coeffs.Length); - if (a.Length < b.Length) + if (a.Coeffs.Length < b.Coeffs.Length) { - a.Coeffs = new double[b.Length]; - b.Coeffs = new double[b.Length]; + a.Coeffs = new double[b.Coeffs.Length]; + b.Coeffs = new double[b.Coeffs.Length]; Array.Copy(aHold, a.Coeffs, aHold.Length); Array.Copy(bHold, b.Coeffs, bHold.Length); } else { - a.Coeffs = new double[a.Length]; - b.Coeffs = new double[a.Length]; + a.Coeffs = new double[a.Coeffs.Length]; + b.Coeffs = new double[a.Coeffs.Length]; Array.Copy(aHold, a.Coeffs, aHold.Length); Array.Copy(bHold, b.Coeffs, bHold.Length); } @@ -546,8 +584,8 @@ namespace MathNet.Numerics public object Clone() { - var p = new double[this.Length]; - Array.Copy(Coeffs, p, Length); + var p = new double[this.Coeffs.Length]; + Array.Copy(Coeffs, p, Coeffs.Length); return new Polynomial(p, isFlip: IsFlipped); } #endregion From 2ba34ef919880c59e4601b71f3c0dac55d00e338 Mon Sep 17 00:00:00 2001 From: Tobias Glaubach Date: Wed, 11 Jul 2018 18:08:52 +0200 Subject: [PATCH 08/18] adjusted the xml comments --- src/Numerics/Polynomial.cs | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/Numerics/Polynomial.cs b/src/Numerics/Polynomial.cs index 2aad4e66..1da5cf91 100644 --- a/src/Numerics/Polynomial.cs +++ b/src/Numerics/Polynomial.cs @@ -140,9 +140,7 @@ namespace MathNet.Numerics #region Data Interaction /// - /// Least-Squares fitting the points (x,y) to a k-order polynomial y : x -> p0 + p1*x + p2*x^2 + ... + pk*x^k, - /// returning its best fitting parameters as [p0, p1, p2, ..., pk] array, compatible with Evaluate.Polynomial. - /// A polynomial with order/degree k has (k+1) coefficients and thus requires at least (k+1) samples. + /// Least-Squares fitting the points (x,y) to a k-order polynomial y : x -> p0 + p1*x + p2*x^2 + ... + pk*x^k /// public static Polynomial Fit(double[] x, double[] y, int order, DirectRegressionMethod method = DirectRegressionMethod.QR) { From 8e86f10c3267379f95d3af4cb0695fce84c86e3a Mon Sep 17 00:00:00 2001 From: Tobias Glaubach Date: Wed, 11 Jul 2018 18:09:44 +0200 Subject: [PATCH 09/18] Wrote tests for "GetRoots" method --- src/Numerics.Tests/PolynomialTests.cs | 61 ++++++++++++++++++++++++++- 1 file changed, 60 insertions(+), 1 deletion(-) diff --git a/src/Numerics.Tests/PolynomialTests.cs b/src/Numerics.Tests/PolynomialTests.cs index 71cba7fb..7267d403 100644 --- a/src/Numerics.Tests/PolynomialTests.cs +++ b/src/Numerics.Tests/PolynomialTests.cs @@ -28,7 +28,9 @@ // using System; +using System.Collections.Generic; using System.Linq; +using System.Numerics; using MathNet.Numerics; using MathNet.Numerics.LinearRegression; using MathNet.Numerics.Statistics; @@ -205,6 +207,63 @@ namespace MathNet.Numerics.UnitTests } } } - + + [Test] + public void GetRootsTest() + { + var tol = 1e-14; + var p1 = new Polynomial(1.0); + var r = p1.GetRoots(); + + Assert.AreEqual(1, r.Length, "length mismatch"); + Assert.AreEqual(1.0, r.FirstOrDefault().Real); + + var p2 = new Polynomial(new double[] { 1, 2 }); + + var r2 = p2.GetRoots(); + Assert.AreEqual(1, r2.Length, "length mismatch"); + Assert.AreEqual(-0.5, r2.FirstOrDefault().Real, tol); + + // T.G: the following expected values were generated using + // numpys np.roots(x) method + // which is equivalent to np.polynomial.polynomial.polyroots + + var x_2 = new double[] { -1.0, 1.0 }; + var expected_2 = new List(); + expected_2.Add(new Complex(1.0, 0.0)); + testEqual(x_2, expected_2); + + var x_3 = new double[] { -1.0, 0.0, 1.0 }; + var expected_3 = new List(); + expected_3.Add(new Complex(1.0, 0.0)); + expected_3.Add(new Complex(-1.0, 0.0)); + testEqual(x_3, expected_3); + + var x_4 = new double[] { -1.0, -0.33333333333333337, 0.33333333333333326, 1.0 }; + var expected_4 = new List(); + expected_4.Add(new Complex(0.9999999999999996, 0.0)); + expected_4.Add(new Complex(-0.6666666666666666, 0.7453559924999296)); + expected_4.Add(new Complex(-0.6666666666666666, -0.7453559924999296)); + testEqual(x_4, expected_4); + } + + private void testEqual(double[] x, List eIn) + { + var tol = 1e-10; + var r0 = new Polynomial(x).GetRoots().ToList(); + + var e = eIn.OrderBy(v => v.Real).ToArray(); + var r = r0.OrderBy(v => v.Real).ToArray(); + + Assert.IsNotNull(r); + Assert.AreEqual(e.Length, r.Length, "length mismatch"); + for (int k = 0; k < r.Length; k++) + { + var msg = String.Format("At k={0}", k); + Assert.AreEqual(e[k].Real, r[k].Real, tol, msg); + Assert.AreEqual(e[k].Imaginary, r[k].Imaginary, tol, msg); + } + + } } } From bbd792974e297c5c31f970b58e633f9eafbd8f88 Mon Sep 17 00:00:00 2001 From: Tobias Glaubach Date: Tue, 17 Jul 2018 09:08:10 +0200 Subject: [PATCH 10/18] implemented a first shot of polynomial long division(still seems buggy atm) aswell as tests --- src/Numerics.Tests/PolynomialTests.cs | 115 +++++++++++++++++++-- src/Numerics/Polynomial.cs | 139 +++++++++++++++++++++++--- 2 files changed, 234 insertions(+), 20 deletions(-) diff --git a/src/Numerics.Tests/PolynomialTests.cs b/src/Numerics.Tests/PolynomialTests.cs index 7267d403..41721553 100644 --- a/src/Numerics.Tests/PolynomialTests.cs +++ b/src/Numerics.Tests/PolynomialTests.cs @@ -29,6 +29,7 @@ using System; using System.Collections.Generic; +using System.Diagnostics; using System.Linq; using System.Numerics; using MathNet.Numerics; @@ -123,8 +124,8 @@ namespace MathNet.Numerics.UnitTests var p_res = Polynomial.Add(p1, p2); var p_tar = new Polynomial(tgt); - p_res.CutTrailZeros(); - p_tar.CutTrailZeros(); + p_res.TrimTrailingZeros(); + p_tar.TrimTrailingZeros(); Assert.AreEqual(p_tar.Degree, p_res.Degree, "length mismatch"); for (int k = 0; k < p_res.Degree; k++) @@ -160,8 +161,8 @@ namespace MathNet.Numerics.UnitTests var p_res = Polynomial.Substract(p1, p2); var p_tar = new Polynomial(tgt); - p_res.CutTrailZeros(); - p_tar.CutTrailZeros(); + p_res.TrimTrailingZeros(); + p_tar.TrimTrailingZeros(); Assert.AreEqual(p_tar.Degree, p_res.Degree, "length mismatch"); for (int k = 0; k < p_res.Degree; k++) @@ -196,8 +197,8 @@ namespace MathNet.Numerics.UnitTests var p_res = p1 * p2; var p_tar = new Polynomial(tgt); - p_res.CutTrailZeros(); - p_tar.CutTrailZeros(); + p_res.TrimTrailingZeros(); + p_tar.TrimTrailingZeros(); Assert.AreEqual(p_tar.Degree, p_res.Degree, "length mismatch"); for (int k = 0; k < p_res.Degree; k++) @@ -208,6 +209,82 @@ namespace MathNet.Numerics.UnitTests } } + + public void DivideLongTestScalar(Tuple inVals, Tuple expectedVals) + { + var p1 = new Polynomial(1.0d); + var p2 = new Polynomial(new double[0]); + var tpl = Polynomial.DivideLong(p1, p2); + + + } + + [Test] + public void DivideLongTest() + { + Assert.Throws(typeof(ArgumentOutOfRangeException), () => + { + var p1 = new Polynomial(1.0d); + var p2 = new Polynomial(new double[0]); + var tpl = Polynomial.DivideLong(p1, p2); + }); + Assert.Throws(typeof(ArgumentOutOfRangeException), () => + { + var p1 = new Polynomial(1.0d); + var p2 = new Polynomial(new double[0]); + var tpl = Polynomial.DivideLong(p2, p1); + }); + Assert.Throws(typeof(ArgumentOutOfRangeException), () => + { + var p1 = new Polynomial(new double[0]); + var p2 = new Polynomial(new double[0]); + var tpl = Polynomial.DivideLong(p2, p1); + }); + Assert.Throws(typeof(DivideByZeroException), () => + { + var p1 = new Polynomial(1.0d); + var p2 = new Polynomial(0.0d); + var tpl = Polynomial.DivideLong(p1, p2); + }); + + var p11 = new Polynomial(2.0d); + var p21 = new Polynomial(2.0d); + var tpl1 = Polynomial.DivideLong(p11, p21); + testEqual(new double[] { 1.0 }, tpl1.Item1); + testEqual(new double[] { 0.0 }, tpl1.Item2); + + var p12 = new Polynomial(new double[] { 2.0d, 2.0d }); + var p22 = new Polynomial(2.0d); + var tpl2 = Polynomial.DivideLong(p12, p22); + testEqual(new double[] { 1.0, 1.0 }, tpl2.Item1); + testEqual(new double[] { 0.0 }, tpl2.Item2); + + for (int i = 0; i < 5; i++) + { + for (int j = 0; j < 5; j++) + { + var msg = String.Format("At i={0}, j={1}", i, j); + var ci = new double[Math.Max(2, i+1)]; + var cj = new double[Math.Max(2, j+1)]; + ci[ci.Length - 1] = 2; + ci[ci.Length - 2] = 1; + cj[cj.Length - 1] = 2; + cj[cj.Length - 2] = 1; + + var pi = new Polynomial(ci); + var pj = new Polynomial(cj); + var tgt = Polynomial.Add(pi, pj); + var tpl3 = Polynomial.DivideLong(tgt, pi); + var pquo = tpl3.Item1; + var prem = tpl3.Item2; + var pres = (pquo * pi) + prem; + testEqual(pres, tgt, msg); + } + } + } + + + [Test] public void GetRootsTest() { @@ -263,7 +340,33 @@ namespace MathNet.Numerics.UnitTests Assert.AreEqual(e[k].Real, r[k].Real, tol, msg); Assert.AreEqual(e[k].Imaginary, r[k].Imaginary, tol, msg); } + } + private void testEqual(double[] p_tar, double[] p_res, string msg = null) + { + Assert.AreEqual(p_tar.Length, p_res.Length, "length mismatch"); + for (int k = 0; k < p_res.Length; k++) + { + Assert.AreEqual(p_tar[k], p_res[k], msg); + } + } + + private void testEqual(double[] p_tar, Polynomial p_res, string msg = null) + { + Assert.AreEqual(p_tar.Length, p_res.Degree, "length mismatch"); + for (int k = 0; k < p_res.Degree; k++) + { + Assert.AreEqual(p_tar[k], p_res.Coeffs[k], msg); + } + + } + private void testEqual(Polynomial p_tar, Polynomial p_res, string msg = null) + { + Assert.AreEqual(p_tar.Degree, p_res.Degree, "length mismatch"); + for (int k = 0; k < p_res.Degree; k++) + { + Assert.AreEqual(p_tar.Coeffs[k], p_res.Coeffs[k], msg); + } } } } diff --git a/src/Numerics/Polynomial.cs b/src/Numerics/Polynomial.cs index 1da5cf91..1abd363f 100644 --- a/src/Numerics/Polynomial.cs +++ b/src/Numerics/Polynomial.cs @@ -15,11 +15,13 @@ using MathNet.Numerics.LinearAlgebra.Factorization; namespace MathNet.Numerics { /// - /// a class handlin REAL VALUED Polynomials, complex coefficients can not be handled (yet) + /// a class handling REAL VALUED Polynomials, complex coefficients can not be handled (yet) /// public class Polynomial { - + /// + /// The coefficients of the polynomial in a + /// public double[] Coeffs { get; set; } /// @@ -116,7 +118,7 @@ namespace MathNet.Numerics /// /// remove all trailing zeros, e.G before: "0.00 x^2 + 1.0 x^1 + 1.00" after: "1.0 x^1 + 1.00" /// - public void CutTrailZeros() + public void TrimTrailingZeros() { int count = 0; for (int ii = Degree - 1; ii >= 0; ii--) @@ -183,28 +185,28 @@ namespace MathNet.Numerics } var t = this.Clone() as Polynomial; - t.CutTrailZeros(); + t.TrimTrailingZeros(); var cNew = new double[t.Coeffs.Length - 1]; for (int i = 1; i < t.Coeffs.Length; i++) { cNew[i-1] = t.Coeffs[i] * i; } var p = new Polynomial(cNew, isFlip: IsFlipped); - p.CutTrailZeros(); + p.TrimTrailingZeros(); return p; } public Polynomial Integrate() { var t = this.Clone() as Polynomial; - t.CutTrailZeros(); + t.TrimTrailingZeros(); var cNew = new double[t.Coeffs.Length + 1]; for (int i = 1; i < cNew.Length; i++) { cNew[i] = t.Coeffs[i-1] / i; } var p = new Polynomial(cNew, isFlip: IsFlipped); - p.CutTrailZeros(); + p.TrimTrailingZeros(); return p; } @@ -222,13 +224,13 @@ namespace MathNet.Numerics public static Polynomial operator *( Polynomial a, Polynomial b) { // do not cut trailing zeros, since it may corrupt the outcom, if the array is of form 1 + x^-1 + x^-2 + x^-3 - //a.CutTrailZeros(); - //b.CutTrailZeros(); + //a.TrimTrailingZeros(); + //b.TrimTrailingZeros(); double[] ret = conv(a.Coeffs, b.Coeffs); Polynomial ret_p = new Polynomial(ret); - //ret_p.CutTrailZeros(); + //ret_p.TrimTrailingZeros(); return (ret_p); @@ -310,7 +312,7 @@ namespace MathNet.Numerics } /// - /// Calculates the complex roots of the Polynomial in the same way as matlab does + /// Calculates the complex roots of the Polynomial by eigenvalue decomposition /// /// a vector of complex numbers with the roots public Complex[] GetRoots() @@ -337,13 +339,14 @@ namespace MathNet.Numerics } /// - /// get the eigenvalue matrix A of this Polynomial such that eig(A) = roots of this Polynomial + /// get the eigenvalue matrix A of this Polynomial such that eig(A) = roots of this Polynomial. /// /// Eigenvalue matrix A + /// this matrix is similar to the companion matrix of this polynomial, in such a way, that it's transpose is the columnflip of the companion matrix public DenseMatrix GetEigValMatrix() { Polynomial pLoc = new Polynomial(this.Coeffs); - pLoc.CutTrailZeros(); + pLoc.TrimTrailingZeros(); int n = pLoc.Coeffs.Length - 1; if (n < 2) @@ -459,6 +462,112 @@ namespace MathNet.Numerics return (res_poly); } + /// + /// Division of two polynomials returning the quotient-with-remainder of the two polynomials given + /// + /// left polynomial + /// right polynomial + /// a tuple holding quotient in first and remainder in second + public static Tuple DivideLong(Polynomial a, Polynomial b) + { + if (a == null) + throw new ArgumentNullException("a"); + if (b == null) + throw new ArgumentNullException("b"); + + if (a.Degree <= 0) + throw new ArgumentOutOfRangeException("a Degree must be greater than zero"); + if (b.Degree <= 0) + throw new ArgumentOutOfRangeException("b Degree must be greater than zero"); + + if (b.Coeffs.Last() == 0) + throw new DivideByZeroException("b polynomial ends with zero"); + + var c1 = a.Coeffs.ToArray(); + var c2 = b.Coeffs.ToArray(); + + var n1 = c1.Length; + var n2 = c2.Length; + + double[] quo = null; + double[] rem = null; + + if (n2 == 1) // division by scalar + { + var fact = c2[0]; + quo = new double[n1]; + for (int i = 0; i < n1; i++) + quo[i] = c1[i] / fact; + rem = new double[] { 0 }; + } + else if(n1 < n2) // denominator degree higher than nominator degree + { + // quotient always be 0 and return c1 as remainder + quo = new double[] { 0 }; + rem = c1.ToArray(); + } + else + { + var dn = n1 - n2; + var scl = c2[n2 - 1]; + var c22 = new double[n2 - 1]; + for (int ii = 0; ii < c22.Length; ii++) + c22[ii] = c2[ii] / scl; + + int i = dn; + int j = n1 - 1; + while (i >= 0) + { + var vals = new double[j - i]; + for (int k = 0; k < vals.Length; k++) + vals[k] = c22[k] * c1[j]; + + for (int idx = i; idx < j; idx++) + c1[idx] -= vals[idx - i]; + + i--; + j--; + } + + rem = new double[j + 1]; + quo = new double[n1 - j + 1]; + + Array.Copy(c1, rem, j + 1); + + for (int idx2 = j+1; idx2 < n1; idx2++) + quo[idx2 - j + 1] = c1[idx2] / scl; + + } + + if (rem == null) + throw new NullReferenceException("resulting remainder was null"); + + if (quo == null) + throw new NullReferenceException("resulting quotient was null"); + + + // output mapping + var pQuo = new Polynomial(quo); + var pRem = new Polynomial(rem); + pQuo.TrimTrailingZeros(); + pQuo.TrimTrailingZeros(); + + return new Tuple(pQuo, pRem); + } + + + /// + /// Division of two polynomials returning the quotient-with-remainder of the two polynomials given + /// + /// left polynomial + /// right polynomial + /// a tuple holding quotient in first and remainder in second + public Tuple DivideLong(Polynomial b) + { + return DivideLong(this, b); + } + + #endregion #region Displaying @@ -492,7 +601,9 @@ namespace MathNet.Numerics for (int ii = 0; ii < Coeffs.Length; ii++) { - if (ii == 0) + if (ii == 0 && Coeffs.Length == 1) + strLoc += String.Format("{0}", this.Coeffs[ii], VarName, ii); + else if(ii == 0) strLoc += String.Format("{0} + ", this.Coeffs[ii], VarName, ii); else if (ii == Coeffs.Length - 1) strLoc += String.Format("{0}{1}{2}", this.Coeffs[ii], VarName, ii); From 099d36b7df234e007932d57bb12fc270076628da Mon Sep 17 00:00:00 2001 From: Tobias Glaubach Date: Tue, 17 Jul 2018 09:08:25 +0200 Subject: [PATCH 11/18] added root finding for polynomials of all degrees --- src/Numerics/FindRoots.cs | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/Numerics/FindRoots.cs b/src/Numerics/FindRoots.cs index 3e8af9a7..d8dc8963 100644 --- a/src/Numerics/FindRoots.cs +++ b/src/Numerics/FindRoots.cs @@ -110,6 +110,15 @@ namespace MathNet.Numerics { return RootFinding.Cubic.Roots(d, c, b, a); } + /// + /// Find all roots of a polynomial by calculating the characteristic polynomial of the companion matrix + /// + /// the values for the polynomial in ascending order e.G new double[] {5, 0, 2} = "5 + 0 x^1 + 2 x^2" + /// the roots of the polynomial + public static Complex[] Polynomial(double[] poly) + { + return new Polynomial(poly).GetRoots(); + } /// /// Find all roots of the Chebychev polynomial of the first kind. From 5a2f969ec2c42bc40cd145a9fe450fd67635b27d Mon Sep 17 00:00:00 2001 From: Tobias Glaubach Date: Tue, 17 Jul 2018 20:23:20 +0200 Subject: [PATCH 12/18] renamed TrimTrailingZeros to Trim --- src/Numerics.Tests/PolynomialTests.cs | 12 ++++++------ src/Numerics/Polynomial.cs | 23 ++++++++++++----------- 2 files changed, 18 insertions(+), 17 deletions(-) diff --git a/src/Numerics.Tests/PolynomialTests.cs b/src/Numerics.Tests/PolynomialTests.cs index 41721553..b4b78eba 100644 --- a/src/Numerics.Tests/PolynomialTests.cs +++ b/src/Numerics.Tests/PolynomialTests.cs @@ -124,8 +124,8 @@ namespace MathNet.Numerics.UnitTests var p_res = Polynomial.Add(p1, p2); var p_tar = new Polynomial(tgt); - p_res.TrimTrailingZeros(); - p_tar.TrimTrailingZeros(); + p_res.Trim(); + p_tar.Trim(); Assert.AreEqual(p_tar.Degree, p_res.Degree, "length mismatch"); for (int k = 0; k < p_res.Degree; k++) @@ -161,8 +161,8 @@ namespace MathNet.Numerics.UnitTests var p_res = Polynomial.Substract(p1, p2); var p_tar = new Polynomial(tgt); - p_res.TrimTrailingZeros(); - p_tar.TrimTrailingZeros(); + p_res.Trim(); + p_tar.Trim(); Assert.AreEqual(p_tar.Degree, p_res.Degree, "length mismatch"); for (int k = 0; k < p_res.Degree; k++) @@ -197,8 +197,8 @@ namespace MathNet.Numerics.UnitTests var p_res = p1 * p2; var p_tar = new Polynomial(tgt); - p_res.TrimTrailingZeros(); - p_tar.TrimTrailingZeros(); + p_res.Trim(); + p_tar.Trim(); Assert.AreEqual(p_tar.Degree, p_res.Degree, "length mismatch"); for (int k = 0; k < p_res.Degree; k++) diff --git a/src/Numerics/Polynomial.cs b/src/Numerics/Polynomial.cs index 1abd363f..3213a5b3 100644 --- a/src/Numerics/Polynomial.cs +++ b/src/Numerics/Polynomial.cs @@ -118,7 +118,7 @@ namespace MathNet.Numerics /// /// remove all trailing zeros, e.G before: "0.00 x^2 + 1.0 x^1 + 1.00" after: "1.0 x^1 + 1.00" /// - public void TrimTrailingZeros() + public void Trim() { int count = 0; for (int ii = Degree - 1; ii >= 0; ii--) @@ -138,6 +138,7 @@ namespace MathNet.Numerics } } } + #region Data Interaction @@ -185,28 +186,28 @@ namespace MathNet.Numerics } var t = this.Clone() as Polynomial; - t.TrimTrailingZeros(); + t.Trim(); var cNew = new double[t.Coeffs.Length - 1]; for (int i = 1; i < t.Coeffs.Length; i++) { cNew[i-1] = t.Coeffs[i] * i; } var p = new Polynomial(cNew, isFlip: IsFlipped); - p.TrimTrailingZeros(); + p.Trim(); return p; } public Polynomial Integrate() { var t = this.Clone() as Polynomial; - t.TrimTrailingZeros(); + t.Trim(); var cNew = new double[t.Coeffs.Length + 1]; for (int i = 1; i < cNew.Length; i++) { cNew[i] = t.Coeffs[i-1] / i; } var p = new Polynomial(cNew, isFlip: IsFlipped); - p.TrimTrailingZeros(); + p.Trim(); return p; } @@ -224,13 +225,13 @@ namespace MathNet.Numerics public static Polynomial operator *( Polynomial a, Polynomial b) { // do not cut trailing zeros, since it may corrupt the outcom, if the array is of form 1 + x^-1 + x^-2 + x^-3 - //a.TrimTrailingZeros(); - //b.TrimTrailingZeros(); + //a.Trim(); + //b.Trim(); double[] ret = conv(a.Coeffs, b.Coeffs); Polynomial ret_p = new Polynomial(ret); - //ret_p.TrimTrailingZeros(); + //ret_p.Trim(); return (ret_p); @@ -346,7 +347,7 @@ namespace MathNet.Numerics public DenseMatrix GetEigValMatrix() { Polynomial pLoc = new Polynomial(this.Coeffs); - pLoc.TrimTrailingZeros(); + pLoc.Trim(); int n = pLoc.Coeffs.Length - 1; if (n < 2) @@ -549,8 +550,8 @@ namespace MathNet.Numerics // output mapping var pQuo = new Polynomial(quo); var pRem = new Polynomial(rem); - pQuo.TrimTrailingZeros(); - pQuo.TrimTrailingZeros(); + pQuo.Trim(); + pQuo.Trim(); return new Tuple(pQuo, pRem); } From e0f2378c221b1899e4f1c6668c3a21e6d8a73f8d Mon Sep 17 00:00:00 2001 From: Tobias Glaubach Date: Tue, 17 Jul 2018 20:50:09 +0200 Subject: [PATCH 13/18] changed somne comments. improved the Trim implementation added a UnitTest for Trim --- src/Numerics.Tests/PolynomialTests.cs | 13 ++++++++++++ src/Numerics/Polynomial.cs | 29 +++++++++++---------------- 2 files changed, 25 insertions(+), 17 deletions(-) diff --git a/src/Numerics.Tests/PolynomialTests.cs b/src/Numerics.Tests/PolynomialTests.cs index b4b78eba..342cf67a 100644 --- a/src/Numerics.Tests/PolynomialTests.cs +++ b/src/Numerics.Tests/PolynomialTests.cs @@ -210,6 +210,19 @@ namespace MathNet.Numerics.UnitTests } + [TestCase(new double[] { 5, 4, 0 }, "5 + 4x^1")] + [TestCase(new double[] { 0, 0, 0 }, "0")] + [TestCase(new double[] { 5, 4, 3, 0, 2 }, "5 + 4x^1 + 3x^2 + 0x^3 + 2x^4")] + [TestCase(new double[] { 0, 0, 8, 0, 0 }, "0 + 0x^1 + 8x^2")] + [TestCase(new double[] { 0, 4, 3, 0, 0 }, "0 + 4x^1 + 3x^2")] + public void TrimTest(double[] x, string expected) + { + var p = new Polynomial(x); + p.Trim(); + Assert.AreEqual(expected, p.ToString()); + + } + public void DivideLongTestScalar(Tuple inVals, Tuple expectedVals) { var p1 = new Polynomial(1.0d); diff --git a/src/Numerics/Polynomial.cs b/src/Numerics/Polynomial.cs index 3213a5b3..4992f1d0 100644 --- a/src/Numerics/Polynomial.cs +++ b/src/Numerics/Polynomial.cs @@ -93,8 +93,8 @@ namespace MathNet.Numerics /// var xP1 = new Polynomial(x, isFlip:true); /// var xP2 = new Polynomial(x, isFlip:false); /// - /// xP1: 5 x^3 + 4 x^2 + 3 x^2 + 0 x^1 + 2 - /// xP2: 2 x^3 + 0 x^2 + 3 x^2 + 4 x^1 + 5 + /// xP1: 2 + 0x^1 + 3x^2 + 4x^3 + 5x^4 + /// xP2: 5 + 4x^1 + 3x^2 + 0x^3 + 2x^4 /// /// /// Polynomial coefficiens as array @@ -120,22 +120,17 @@ namespace MathNet.Numerics /// public void Trim() { - int count = 0; - for (int ii = Degree - 1; ii >= 0; ii--) + int i = Degree - 1; + while (i >= 0 && Coeffs[i] == 0.0) + i--; + + if (i < 0) + Coeffs = new double[0]; + else { - if (Coeffs[ii] == 0.0) - { - count++; - } - else - { - double[] CoeffsHold = new double[Coeffs.Length]; - Coeffs.CopyTo(CoeffsHold, 0); - Array.Resize(ref CoeffsHold, Degree - count); - Coeffs = new double[Degree - count]; - CoeffsHold.CopyTo(Coeffs, 0); - return; - } + var hold = new double[i+1]; + Array.Copy(Coeffs, hold, i+1); + Coeffs = hold; } } From 8a34e61275c32c8fff8639a367cd8a91ebf0281c Mon Sep 17 00:00:00 2001 From: Tobias Glaubach Date: Tue, 17 Jul 2018 23:39:14 +0200 Subject: [PATCH 14/18] removed the isflipped Property, improved the clone method, fixed the divideLong --- src/Numerics/Polynomial.cs | 163 ++++++++++++++++++------------------- 1 file changed, 79 insertions(+), 84 deletions(-) diff --git a/src/Numerics/Polynomial.cs b/src/Numerics/Polynomial.cs index 4992f1d0..efc1499b 100644 --- a/src/Numerics/Polynomial.cs +++ b/src/Numerics/Polynomial.cs @@ -24,11 +24,6 @@ namespace MathNet.Numerics /// public double[] Coeffs { get; set; } - /// - /// indicator if Polynomial was flipped - /// - public bool IsFlipped { get; } - /// /// Only needed for the ToString method /// @@ -65,7 +60,6 @@ namespace MathNet.Numerics /// just the "x^0" part public Polynomial(double coeff) { - IsFlipped = false; this.Coeffs = new double[1]; Coeffs[0] = coeff; } @@ -73,59 +67,44 @@ namespace MathNet.Numerics /// /// make Polynomial: e.G new double[] {5, 0, 2} = "5 + 0 x^1 + 2 x^2" /// - /// Polynomial coefficiens as array - public Polynomial(double[] coeffs) + /// Polynomial coefficiens as enumerable + public Polynomial(IEnumerable coeffs) { if (coeffs == null) { throw new ArgumentNullException("coeffs"); } - this.Coeffs = new double[coeffs.Length]; - Array.Copy(coeffs, this.Coeffs, coeffs.Length); + this.Coeffs = coeffs.ToArray(); } - /// - /// constructor setting Polynomial coefficiens and flipping them if necessary. - /// - /// e.G: - /// var x = new double[] {5, 4, 3, 0, 2}; - /// var xP1 = new Polynomial(x, isFlip:true); - /// var xP2 = new Polynomial(x, isFlip:false); - /// - /// xP1: 2 + 0x^1 + 3x^2 + 4x^3 + 5x^4 - /// xP2: 5 + 4x^1 + 3x^2 + 0x^3 + 2x^4 - /// + /// make Polynomial: e.G new double[] {5, 0, 2} = "5 + 0 x^1 + 2 x^2" /// /// Polynomial coefficiens as array - /// use true for flipping - public Polynomial(double[] coeffs, bool isFlip) + public Polynomial(double[] coeffs) { if (coeffs == null) { throw new ArgumentNullException("coeffs"); } this.Coeffs = new double[coeffs.Length]; - Array.Copy(coeffs, Coeffs, coeffs.Length); - - if (isFlip) - { - Coeffs = Coeffs.Reverse().ToArray(); - IsFlipped = true; - } + Array.Copy(coeffs, this.Coeffs, coeffs.Length); } /// - /// remove all trailing zeros, e.G before: "0.00 x^2 + 1.0 x^1 + 1.00" after: "1.0 x^1 + 1.00" + /// remove all trailing zeros, e.G before: "3 + 2 x^1 + 0 x^2" after: "3 + 2 x^1" /// public void Trim() { + if (Degree == 1) + return; + int i = Degree - 1; while (i >= 0 && Coeffs[i] == 0.0) i--; - if (i < 0) - Coeffs = new double[0]; + if (i <= 0) + Coeffs = new double[1] { 0.0 }; else { var hold = new double[i+1]; @@ -187,7 +166,7 @@ namespace MathNet.Numerics { cNew[i-1] = t.Coeffs[i] * i; } - var p = new Polynomial(cNew, isFlip: IsFlipped); + var p = new Polynomial(cNew); p.Trim(); return p; } @@ -201,7 +180,7 @@ namespace MathNet.Numerics { cNew[i] = t.Coeffs[i-1] / i; } - var p = new Polynomial(cNew, isFlip: IsFlipped); + var p = new Polynomial(cNew); p.Trim(); return p; } @@ -219,11 +198,13 @@ namespace MathNet.Numerics /// resulting Polynomial public static Polynomial operator *( Polynomial a, Polynomial b) { + var aa = a.Clone() as Polynomial; + var bb = b.Clone() as Polynomial; // do not cut trailing zeros, since it may corrupt the outcom, if the array is of form 1 + x^-1 + x^-2 + x^-3 //a.Trim(); //b.Trim(); - double[] ret = conv(a.Coeffs, b.Coeffs); + double[] ret = conv(aa.Coeffs, bb.Coeffs); Polynomial ret_p = new Polynomial(ret); //ret_p.Trim(); @@ -240,10 +221,13 @@ namespace MathNet.Numerics /// resulting Polynomial public static Polynomial operator *( Polynomial a, double k) { - for (int ii = 0; ii < a.Coeffs.Length; ii++) - a.Coeffs[ii] *= k; + var aa = a.Clone() as Polynomial; + - return a; + for (int ii = 0; ii < aa.Coeffs.Length; ii++) + aa.Coeffs[ii] *= k; + + return aa; } /// @@ -254,8 +238,10 @@ namespace MathNet.Numerics /// resulting Polynomial public static Polynomial operator +( Polynomial a, double k) { - a.Coeffs[0] += k; - return a; + var aa = a.Clone() as Polynomial; + + aa.Coeffs[0] += k; + return aa; } /// @@ -266,9 +252,10 @@ namespace MathNet.Numerics /// resulting Polynomial public static Polynomial operator -( Polynomial a, double k) { + var aa = a.Clone() as Polynomial; a.Coeffs[0] -= k; - return a; + return aa; } /// @@ -279,10 +266,12 @@ namespace MathNet.Numerics /// resulting Polynomial public static Polynomial operator /( Polynomial a, double k) { - for (int ii = 0; ii < a.Coeffs.Length; ii++) - a.Coeffs[ii] /= k; + var aa = a.Clone() as Polynomial; + + for (int ii = 0; ii < aa.Coeffs.Length; ii++) + aa.Coeffs[ii] /= k; - return a; + return aa; } /// @@ -372,16 +361,19 @@ namespace MathNet.Numerics /// resulting Polynomial public static Polynomial DividePointwise( Polynomial a, Polynomial b) { - if (a.Coeffs.Length != b.Coeffs.Length) - mkSameLength(ref a, ref b); + var aa = a.Clone() as Polynomial; + var bb = b.Clone() as Polynomial; + + if (aa.Coeffs.Length != bb.Coeffs.Length) + mkSameLength(ref aa, ref bb); - int n = a.Coeffs.Length; - double[] res = new double[a.Coeffs.Length]; + int n = aa.Coeffs.Length; + double[] res = new double[aa.Coeffs.Length]; for (int ii = 0; ii < n; ii++) { - res[ii] = a.Coeffs[ii] / b.Coeffs[ii]; + res[ii] = aa.Coeffs[ii] / bb.Coeffs[ii]; } Polynomial res_poly = new Polynomial(res); return (res_poly); @@ -395,16 +387,19 @@ namespace MathNet.Numerics /// resulting Polynomial public static Polynomial MultiplyPointwise( Polynomial a, Polynomial b) { - if (a.Coeffs.Length != b.Coeffs.Length) - mkSameLength(ref a, ref b); + var aa = a.Clone() as Polynomial; + var bb = b.Clone() as Polynomial; + + if (aa.Coeffs.Length != bb.Coeffs.Length) + mkSameLength(ref aa, ref bb); - int n = a.Coeffs.Length; - double[] res = new double[a.Coeffs.Length]; + int n = aa.Coeffs.Length; + double[] res = new double[aa.Coeffs.Length]; for (int ii = 0; ii < n; ii++) { - res[ii] = a.Coeffs[ii] * b.Coeffs[ii]; + res[ii] = aa.Coeffs[ii] * bb.Coeffs[ii]; } Polynomial res_poly = new Polynomial(res); return (res_poly); @@ -418,17 +413,19 @@ namespace MathNet.Numerics /// resulting Polynomial public static Polynomial Add( Polynomial a, Polynomial b) { + var aa = a.Clone() as Polynomial; + var bb = b.Clone() as Polynomial; - if (a.Degree != b.Degree) - mkSameLength(ref a, ref b); + if (aa.Degree != bb.Degree) + mkSameLength(ref aa, ref bb); - int n = a.Degree; + int n = aa.Degree; double[] res = new double[n]; for (int ii = 0; ii < n; ii++) { - res[ii] = a.Coeffs[ii] + b.Coeffs[ii]; + res[ii] = aa.Coeffs[ii] + bb.Coeffs[ii]; } Polynomial res_poly = new Polynomial(res); return (res_poly); @@ -442,17 +439,19 @@ namespace MathNet.Numerics /// resulting Polynomial public static Polynomial Substract( Polynomial a, Polynomial b) { + var aa = a.Clone() as Polynomial; + var bb = b.Clone() as Polynomial; - if (a.Degree != b.Degree) - mkSameLength(ref a, ref b); + if (aa.Degree != bb.Degree) + mkSameLength(ref aa, ref bb); - int n = a.Degree; + int n = aa.Degree; double[] res = new double[n]; for (int ii = 0; ii < n; ii++) { - res[ii] = a.Coeffs[ii] - b.Coeffs[ii]; + res[ii] = aa.Coeffs[ii] - bb.Coeffs[ii]; } Polynomial res_poly = new Polynomial(res); return (res_poly); @@ -476,7 +475,7 @@ namespace MathNet.Numerics if (b.Degree <= 0) throw new ArgumentOutOfRangeException("b Degree must be greater than zero"); - if (b.Coeffs.Last() == 0) + if (b.Coeffs[b.Degree-1] == 0) throw new DivideByZeroException("b polynomial ends with zero"); var c1 = a.Coeffs.ToArray(); @@ -514,24 +513,25 @@ namespace MathNet.Numerics int j = n1 - 1; while (i >= 0) { + var v = c1[j]; var vals = new double[j - i]; - for (int k = 0; k < vals.Length; k++) - vals[k] = c22[k] * c1[j]; - - for (int idx = i; idx < j; idx++) - c1[idx] -= vals[idx - i]; - + for (int k = i; k < j; k++) + c1[k] -= c22[k-i] * v; i--; j--; } - rem = new double[j + 1]; - quo = new double[n1 - j + 1]; + var j1 = j + 1; + var l1 = n1 - j1; - Array.Copy(c1, rem, j + 1); + rem = new double[j1]; + quo = new double[l1]; - for (int idx2 = j+1; idx2 < n1; idx2++) - quo[idx2 - j + 1] = c1[idx2] / scl; + for (int k = 0; k < l1; k++) + quo[k] = c1[k + j1] / scl; + + for (int k = 0; k < j1; k++) + rem[k] = c1[k]; } @@ -545,9 +545,9 @@ namespace MathNet.Numerics // output mapping var pQuo = new Polynomial(quo); var pRem = new Polynomial(rem); + + pRem.Trim(); pQuo.Trim(); - pQuo.Trim(); - return new Tuple(pQuo, pRem); } @@ -633,10 +633,7 @@ namespace MathNet.Numerics /// the coefficcients of the Polynomial as an array public double[] ToArray() { - if (IsFlipped == true) - return (Coeffs.Reverse().ToArray()); - else - return (Coeffs); + return (Coeffs.ToArray()); } #endregion @@ -689,9 +686,7 @@ namespace MathNet.Numerics public object Clone() { - var p = new double[this.Coeffs.Length]; - Array.Copy(Coeffs, p, Coeffs.Length); - return new Polynomial(p, isFlip: IsFlipped); + return new Polynomial(p); } #endregion From f229b266ebd7d8575fda33be8edb214a79d40fb4 Mon Sep 17 00:00:00 2001 From: Tobias Glaubach Date: Tue, 17 Jul 2018 23:53:29 +0200 Subject: [PATCH 15/18] grooming and fixed a test --- src/Numerics.Tests/PolynomialTests.cs | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/src/Numerics.Tests/PolynomialTests.cs b/src/Numerics.Tests/PolynomialTests.cs index 342cf67a..9b96f3be 100644 --- a/src/Numerics.Tests/PolynomialTests.cs +++ b/src/Numerics.Tests/PolynomialTests.cs @@ -39,10 +39,14 @@ using NUnit.Framework; namespace MathNet.Numerics.UnitTests { + /// + /// some of these tests were inspired by numpys tests in python for the Polynomial functions. + /// Thanks to the numpy contributers! + /// [TestFixture, Category("Calculus")] public class PolynomialTests { - + [TestCase(new double[] { 5, 4, 3, 0, 2 }, "5 + 4x^1 + 3x^2 + 0x^3 + 2x^4")] [TestCase(new double[0], "")] [TestCase(new double[] { 0, 4, 3, 0, 0 }, "0 + 4x^1 + 3x^2 + 0x^3 + 0x^4")] @@ -233,7 +237,7 @@ namespace MathNet.Numerics.UnitTests } [Test] - public void DivideLongTest() + public void DivideLongTestWrongInputs() { Assert.Throws(typeof(ArgumentOutOfRangeException), () => { @@ -259,6 +263,11 @@ namespace MathNet.Numerics.UnitTests var p2 = new Polynomial(0.0d); var tpl = Polynomial.DivideLong(p1, p2); }); + } + + [Test] + public void DivideLongTest() + { var p11 = new Polynomial(2.0d); var p21 = new Polynomial(2.0d); @@ -277,8 +286,8 @@ namespace MathNet.Numerics.UnitTests for (int j = 0; j < 5; j++) { var msg = String.Format("At i={0}, j={1}", i, j); - var ci = new double[Math.Max(2, i+1)]; - var cj = new double[Math.Max(2, j+1)]; + var ci = new double[i + 2]; + var cj = new double[j + 2]; ci[ci.Length - 1] = 2; ci[ci.Length - 2] = 1; cj[cj.Length - 1] = 2; @@ -291,13 +300,13 @@ namespace MathNet.Numerics.UnitTests var pquo = tpl3.Item1; var prem = tpl3.Item2; var pres = (pquo * pi) + prem; + pres.Trim(); + testEqual(pres, tgt, msg); } } } - - [Test] public void GetRootsTest() { From 37e405ce2fd62100749f944c28e69887cf79dd7e Mon Sep 17 00:00:00 2001 From: Tobias Glaubach Date: Tue, 17 Jul 2018 23:53:52 +0200 Subject: [PATCH 16/18] small fixes --- src/Numerics/Polynomial.cs | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/Numerics/Polynomial.cs b/src/Numerics/Polynomial.cs index efc1499b..ce2e6cd9 100644 --- a/src/Numerics/Polynomial.cs +++ b/src/Numerics/Polynomial.cs @@ -103,8 +103,10 @@ namespace MathNet.Numerics while (i >= 0 && Coeffs[i] == 0.0) i--; - if (i <= 0) + if (i < 0) Coeffs = new double[1] { 0.0 }; + else if (i == 0) + Coeffs = new double[1] { Coeffs[0] }; else { var hold = new double[i+1]; @@ -686,7 +688,7 @@ namespace MathNet.Numerics public object Clone() { - return new Polynomial(p); + return new Polynomial(this.Coeffs); } #endregion From 76eeb0d445a9e2ca9ad0d90ecb157008a9eeb3c5 Mon Sep 17 00:00:00 2001 From: Tobias Glaubach Date: Tue, 17 Jul 2018 23:57:12 +0200 Subject: [PATCH 17/18] deleted a unwanted file in this branch... will be added in another pull request --- .../LtiSystems/TransferFunctionDiscrete.cs | 968 ------------------ 1 file changed, 968 deletions(-) delete mode 100644 src/Numerics/LtiSystems/TransferFunctionDiscrete.cs diff --git a/src/Numerics/LtiSystems/TransferFunctionDiscrete.cs b/src/Numerics/LtiSystems/TransferFunctionDiscrete.cs deleted file mode 100644 index c050eca8..00000000 --- a/src/Numerics/LtiSystems/TransferFunctionDiscrete.cs +++ /dev/null @@ -1,968 +0,0 @@ -using System; -using System.Collections.Generic; -using System.Linq; -using System.Numerics; -using System.Text; -using MathNet.Numerics; - -namespace MathNet.Numerics.LtiSystems -{ - /// Class for LTI discrete transfer functions - public class TransferFunctionDiscrete - { - - private double[] _num; - - private double[] _den; - - /// - /// numberator (input dependent) Polynomial coefficients as array - /// in order - /// => index high ... index low - /// => [n], [n-1], +..., [0] - /// => 1 + q^-1 + ... + q^-n - /// - public double[] num - { - get - { - return _num; - } - set - { - _num = cutTrailingZeros(value); - shiftNumDenIfPossible(); - } - } - - /// den (state dependent) Polynomial coefficients as array - /// in order - /// => index high ... index low - /// => [n], [n-1], +..., [0] - /// => 1 + q^-1 + ... + q^-n - /// - public double[] den - { - get - { - return _den; - } - set - { - _den = cutTrailingZeros(value); - shiftNumDenIfPossible(); - } - } - - - - /// b (input dependent) Polynomial coefficients as array ( - public double[] b - { - get - { - return _num; - } - set - { - _num = cutTrailingZeros(value); - shiftNumDenIfPossible(); - } - } - - /// a (state dependent) Polynomial coefficients as array - public double[] a - { - get - { - return _den; - } - set - { - _den = cutTrailingZeros(value); - shiftNumDenIfPossible(); - } - } - - /// Internal FIR States -> updated in every response calculation - public double[] z_FIR { get; set; } - - /// Internal IIR States -> updated in every response calculation - public double[] z_IIR { get; set; } - - /// any name you want to give this transfer function - public string Name { get; set; } - - /// sampling time of discrete transfer function (default = 1) - public double Ts { get; set; } - - /// variable for transfer function so far all tf's are in the q^-1 (or equivalently z^-1) form. Changing this will have NO nfluence besides displaying te TF - public string variable = "q^-1"; - - - - /// - /// Check if this Transfer Function is stable - /// - /// the tolerance for euclidian distance at which a pole/zero pair is considered to be canceling each other - /// false if system is unsable true if system is stable - public bool IsStable(double numTolerance = 1e-8) - { - - var p = GetPoles(); - - var z = GetZeros(); - var z_isCompensated = new bool[z.Length]; - for (int j = 0; j < p.Length; j++) - { - // check if pole would lead to unstable behaviour - if (p[j].Magnitude > 1.0) - { - - // init some values - double minDistance = Double.PositiveInfinity; - int idxMinDistanceZero = -1; - - // analyze the distance between each zero and the pole now - for (int i = 0; i < z.Length; i++) - { - // check if pole has already been used for compensation - if (z_isCompensated[i]) - continue; - - // calculate geometrical distance between each zero and the pole now and store the closest neighbour - var dist = (p[j] - z[i]).Magnitude; - if (dist < minDistance) - { - minDistance = dist; - idxMinDistanceZero = i; - } - } - - // if closest neighbour is too far away to compensate the unstable pole - if (minDistance >= numTolerance) - return false; // the system is unsable - else - z_isCompensated[idxMinDistanceZero] = true; // if not: mark the zero as already used for compensation - - } - - } - - return true; - - } - - - - /// constructor setting no properties at all - public TransferFunctionDiscrete() - { - this.Ts = 1.0; - } - - /// - /// constructor setting a and b vectors as well as initializing the z_FIR and z_IIR states - /// - public TransferFunctionDiscrete(double b_in, double[] a_in, double Ts_in = 1.0d) - { - if (a_in == null) - throw new ArgumentNullException("a_in"); - - this.a = (double[])a_in.Clone(); - this.b = new double[1]; - this.b[0] = b_in; - this.z_IIR = new double[a_in.Length]; - this.z_FIR = new double[1]; - this.Ts = Ts_in; - } - - /// - /// constructor setting a and b vectors as well as initializing the z_FIR and z_IIR states - /// - public TransferFunctionDiscrete(double[] b_in, double a_in, double Ts_in = 1.0d) - { - if (b_in == null) - throw new ArgumentNullException("b_in"); - - this.a = new double[1]; - this.a[0] = a_in; - this.b = (double[])b_in.Clone(); - this.z_IIR = new double[1]; - this.z_FIR = new double[b_in.Length]; - this.Ts = Ts_in; - } - - /// - /// constructor setting a and b vectors as well as initializing the z_FIR and z_IIR states - /// - public TransferFunctionDiscrete(double b_in, double a_in, double Ts_in = 1.0d) - { - this.a = new double[1]; - this.a[0] = a_in; - this.b = new double[1]; - this.b[0] = b_in; - this.z_IIR = new double[1]; - this.z_FIR = new double[1]; - this.Ts = Ts_in; - } - - /// - /// constructor setting a and b vectors as well as initializing the z_FIR and z_IIR states - /// - public TransferFunctionDiscrete(double[] b_in, double[] a_in, double Ts_in = 1.0d) - { - if (b_in == null) - throw new ArgumentNullException("b_in"); - - if (a_in == null) - throw new ArgumentNullException("a_in"); - - this.a = (double[])a_in.Clone(); - this.b = (double[])b_in.Clone(); - this.z_IIR = new double[a_in.Length]; - this.z_FIR = new double[b_in.Length]; - this.Ts = Ts_in; - } - - - /// - /// Adds delay to the numerator array (shifting the values by d steps) - /// - /// integer value of daly to add to this TF - public void AddDelay(int d) - { - double[] b_new = new double[b.Length + d]; - b.CopyTo(b_new, d); - b = b_new; - } - - #region Helpers - /// - /// if num and den both start at a later step (e.G highest power num = q^-3 and highest power den = q^-4), the whole tf can be shifted by n (3) steps - /// - private void shiftNumDenIfPossible() - { - var offset = 0; - if (_num == null || _den == null) - return; - - var n = Math.Min(_num.Length, _den.Length); - for (int i = 0; i < n; i++) - { - if (num[i] == 0.0d && _den[i] == 0.0d) - offset = i + 1; - else - break; - } - - if (offset > 0) - { - double[] tmp1 = new double[_num.Length - offset]; - Array.Copy(_num, offset, tmp1, 0, tmp1.Length); - double[] tmp2 = new double[_den.Length - offset]; - Array.Copy(_den, offset, tmp2, 0, tmp2.Length); - _num = tmp1; - _den = tmp2; - } - } - - private double[] cutTrailingZeros(double[] vIn) - { - int lengthNew = vIn.Length; - - for (int i = vIn.Length - 1; i >= 0; i--) - { - if (vIn[i] != 0) - { - lengthNew = i + 1; - break; - } - } - - var v = new double[lengthNew]; - Array.Copy(vIn, v, lengthNew); - - return v; - } - - /// - /// checks and adjusts internal states to a and b arrays - /// - private void checkStateSizes() - { - if (this.z_IIR == null) - this.z_IIR = new double[this.a.Length]; - - if (this.z_FIR == null) - this.z_FIR = new double[this.b.Length]; - - if (this.a.Length != this.z_IIR.Length) - this.z_IIR = new double[this.a.Length]; - - if (this.b.Length != this.z_FIR.Length) - this.z_FIR = new double[this.b.Length]; - - if (this.a.Length != this.z_IIR.Length) - this.z_IIR = new double[this.a.Length]; - } - #endregion Helpers - - #region Operators - - /// - /// LTI System theory division of a transfer function object by a scalar - /// - /// transfer function - /// scalar for divison - /// new transfer function object divided by k - public static TransferFunctionDiscrete operator /(TransferFunctionDiscrete G1, double k) - { - Polynomial A1 = new Polynomial(G1.a); - - TransferFunctionDiscrete Gres = new TransferFunctionDiscrete(G1.b, (A1 * k).ToArray(), G1.Ts) - { - Name = G1.Name - }; - return Gres; - } - - /// - /// LTI System theory division of a scalar by a transfer function object - /// - /// scalar value - /// transfer function for division - /// new transfer function object - public static TransferFunctionDiscrete operator /(double k, TransferFunctionDiscrete G1) - { - Polynomial A1 = new Polynomial(G1.a); - - TransferFunctionDiscrete Gres = new TransferFunctionDiscrete((A1 * k).ToArray(), G1.b, G1.Ts) - { - Name = G1.Name - }; - return Gres; - } - - /// - /// LTI System theory multiplication of a transfer function object by a scalar - /// - /// transfer function - /// scalar for multiplication - /// new transfer function object - public static TransferFunctionDiscrete operator *(TransferFunctionDiscrete G1, double k) - { - Polynomial B1 = new Polynomial(G1.b); - - TransferFunctionDiscrete Gres = new TransferFunctionDiscrete((B1 * k).ToArray(), G1.a, G1.Ts) - { - Name = G1.Name - }; - return Gres; - - } - - /// - /// LTI System theory multiplication of a transfer function object by a scalar - /// - /// scalar for multiplication - /// transfer function - /// new transfer function object - public static TransferFunctionDiscrete operator *(double k, TransferFunctionDiscrete G1) - { - Polynomial B1 = new Polynomial(G1.b); - - TransferFunctionDiscrete Gres = new TransferFunctionDiscrete((B1 * k).ToArray(), G1.a, G1.Ts) - { - Name = G1.Name - }; - return Gres; - } - - - /// - /// LTI System theory substraction of a transfer function with a scalar - /// - /// transfer function - /// scalar - /// new transfer function object - public static TransferFunctionDiscrete operator -(TransferFunctionDiscrete G1, double k) - { - Polynomial A1 = new Polynomial(G1.a); - Polynomial B1 = new Polynomial(G1.b); - - Polynomial A_res = (A1); - Polynomial B_res = B1 - (A1 * k); - - TransferFunctionDiscrete Gres = new TransferFunctionDiscrete(B_res.ToArray(), A_res.ToArray(), G1.Ts) - { - Name = G1.Name - }; - return Gres; - } - - /// - /// LTI System theory substraction of a scalar by a transfer function - /// - /// scalar - /// transfer function - /// new transfer function object - public static TransferFunctionDiscrete operator -(double k, TransferFunctionDiscrete G1) - { - Polynomial A1 = new Polynomial(G1.a); - Polynomial B1 = new Polynomial(G1.b); - - Polynomial A_res = (A1); - Polynomial B_res = (A1 * k) - B1; - - TransferFunctionDiscrete Gres = new TransferFunctionDiscrete(B_res.ToArray(), A_res.ToArray(), G1.Ts) - { - Name = G1.Name - }; - return Gres; - } - - /// - /// LTI System theory addition of a transfer function with a scalar - /// - /// transfer function - /// scalar - /// new transfer function object - public static TransferFunctionDiscrete operator +(TransferFunctionDiscrete G1, double k) - { - Polynomial A1 = new Polynomial(G1.a); - Polynomial B1 = new Polynomial(G1.b); - - Polynomial A_res = (A1); - Polynomial B_res = (A1 * k) + B1; - - TransferFunctionDiscrete Gres = new TransferFunctionDiscrete(B_res.ToArray(), A_res.ToArray(), G1.Ts) - { - Name = G1.Name - }; - return Gres; - } - - /// - /// LTI System theory addition of a transfer function with a scalar - /// - /// transfer function - /// scalar - /// new transfer function object - public static TransferFunctionDiscrete operator +(double k, TransferFunctionDiscrete G1) - { - Polynomial A1 = new Polynomial(G1.a); - Polynomial B1 = new Polynomial(G1.b); - - Polynomial A_res = (A1); - Polynomial B_res = B1 + (A1 * k); - - TransferFunctionDiscrete Gres = new TransferFunctionDiscrete(B_res.ToArray(), A_res.ToArray(), G1.Ts) - { - Name = G1.Name - }; - return Gres; - } - - /// - /// LTI System theory addition of two transfer functions - /// - /// transfer function left - /// transfer function right - /// new transfer function object - public static TransferFunctionDiscrete operator +(TransferFunctionDiscrete G1, TransferFunctionDiscrete G2) - { - if (Math.Abs(G1.Ts - G2.Ts) > 1e-12) - throw new ArgumentException(String.Format("The two supplied transfer functions do not have equal sampling times. G1.Ts = {0} G2.Ts = {1}", G1.Ts, G2.Ts)); - - Polynomial A1 = new Polynomial(G1.a); - Polynomial B1 = new Polynomial(G1.b); - - Polynomial A2 = new Polynomial(G2.a); - Polynomial B2 = new Polynomial(G2.b); - - Polynomial A_res = (A1 * A2); - Polynomial B_res = (B1 * A2) + (B2 * A1); - - return new TransferFunctionDiscrete(B_res.ToArray(), A_res.ToArray(), G1.Ts); - } - - /// - /// LTI System theory substraction of two transfer functions - /// - /// transfer function left - /// transfer function right - /// new transfer function object - public static TransferFunctionDiscrete operator -(TransferFunctionDiscrete G1, TransferFunctionDiscrete G2) - { - if (Math.Abs(G1.Ts - G2.Ts) > 1e-12) - throw new ArgumentException(String.Format("The two supplied transfer functions do not have equal sampling times. G1.Ts = {0} G2.Ts = {1}", G1.Ts, G2.Ts)); - - Polynomial A1 = new Polynomial(G1.a); - Polynomial B1 = new Polynomial(G1.b); - - Polynomial A2 = new Polynomial(G2.a); - Polynomial B2 = new Polynomial(G2.b); - - Polynomial A_res = (A1 * A2); - Polynomial B_res = (B1 * A2) - (B2 * A1); - - return new TransferFunctionDiscrete(B_res.ToArray(), A_res.ToArray(), G1.Ts); - } - - /// - /// LTI System theory addition of two transfer functions - /// - /// transfer function left - /// transfer function right - /// new transfer function object - public static TransferFunctionDiscrete operator *(TransferFunctionDiscrete G1, TransferFunctionDiscrete G2) - { - if (Math.Abs(G1.Ts - G2.Ts) > 1e-12) - throw new ArgumentException(String.Format("The two supplied transfer functions do not have equal sampling times. G1.Ts = {0} G2.Ts = {1}", G1.Ts, G2.Ts)); - - Polynomial A1 = new Polynomial(G1.a); - Polynomial B1 = new Polynomial(G1.b); - - Polynomial A2 = new Polynomial(G2.a); - Polynomial B2 = new Polynomial(G2.b); - - Polynomial A_res = A1 * A2; - Polynomial B_res = B1 * B2; - - return new TransferFunctionDiscrete(B_res.ToArray(), A_res.ToArray(), G1.Ts); - } - - /// - /// LTI System theory division of two transfer functions - /// - /// transfer function left - /// transfer function right - /// new transfer function object - public static TransferFunctionDiscrete operator /(TransferFunctionDiscrete G1, TransferFunctionDiscrete G2) - { - if (Math.Abs(G1.Ts - G2.Ts) > 1e-12) - throw new ArgumentException(String.Format("The two supplied transfer functions do not have equal sampling times. G1.Ts = {0} G2.Ts = {1}", G1.Ts, G2.Ts)); - - Polynomial A1 = new Polynomial(G1.a); - Polynomial B1 = new Polynomial(G1.b); - - Polynomial A2 = new Polynomial(G2.a); - Polynomial B2 = new Polynomial(G2.b); - - Polynomial A_res = A1 * B2; - Polynomial B_res = B1 * A2; - - return new TransferFunctionDiscrete(B_res.ToArray(), A_res.ToArray(), G1.Ts); - } - #endregion - - /// calculates y_k = G(q^-1) * x_k for a given x_k array - public IEnumerable CalcResponse(IEnumerable x) - { - return this.CalcResponse(x.ToArray()); - } - - /// calculates y_k = G(q^-1) * x_k for a given x_k array - public double[] CalcResponse(double[] x) - { - // this is basically a two step convolution and could be replaced by a - // conv implementation. - // however... this code works fine and replacing it would be more work - - double y_now = 0.0d; - int idx_a = 0; - int idx_b = 0; - double[] y = new double[x.Length]; - - this.checkStateSizes(); - - // Loop all inputs - for (int ii_x = 0; ii_x < x.Length; ii_x++) - { - y_now = 0.0d; - idx_b = 0; - - // loop through b-matrix until end of momentary tempx-array - for (int ii_b = 0; ii_b <= ii_x && idx_b < b.Length; ii_b++) - { - - z_FIR[idx_b] = x[ii_x - ii_b]; - y_now += b[idx_b] * z_FIR[idx_b]; - idx_b++; - } - - - // start at second position, since it's the a-matrix - idx_a = 1; - // loop for a-matrix - for (int ii_a = 0; ii_a <= (ii_x - 1) && idx_a < a.Length; ii_a++) - { - z_IIR[idx_a] = y[(ii_x - 1) - ii_a]; - y_now -= a[idx_a] * z_IIR[idx_a]; - idx_a++; - } - // write result - y[ii_x] = (y_now / a[0]); - z_IIR[0] = y[ii_x]; - } - return (y); - } - - // Todo: Implement FiltFilt - /* - /// - /// A wrapper for the StaticFilters.FiltFilt method using the internal a and b arrays - /// - /// The data to filter - /// initial state coefficients null for aotomatic generation via steady state solution - /// the number of datapoints to pad at each side use less than 0 for Math.Max(a.Length, b.Length) * 3 - /// The filterd data - /// - /// In order to prevent transients at the end or start of the sequence we have to pad it - /// The padding is done by rotating the sequence by 180° at the ends and append it to the data - /// - public double[] FiltFilt(double[] data, double[] zi = null, int padlen = 0) - { - if (this.a == null || this.a.Length == 0) - throw new Exception("This transfer function has no a array with data"); - if (this.b == null || this.b.Length == 0) - throw new Exception("This transfer function has no a array with data"); - - return StaticFilters.FiltFilt(data, this.a, this.b, zi, padlen); - } - */ - - #region Dynamics - - /// - /// returns the impulse response with nSteps for the tf model - /// - /// number of steps for impulse response - /// - public double[] Impulse(int nSteps) - { - - var Inp = new double[nSteps]; - Inp[0] = 1.0; - - var ImpulseResponse = this.CalcResponse(Inp); - return (ImpulseResponse); - } - - /// - /// returns the impulse response with nSettling * 1.3 steps for the tf model - /// - public double[] Impulse() - { - - var nSteps = Convert.ToInt32((double)CalcSettlingSteps() * 1.3); - if (nSteps <= 0) - return null; - var Inp = new double[nSteps]; - Inp[0] = 1.0; - - var ImpulseResponse = this.CalcResponse(Inp); - return (ImpulseResponse); - } - - - public Complex[] Bode(int nPoints = 100) - { - // substituting z = exp(j * omega * Ts) - var omega_vec = Generate.LinearSpaced(nPoints, 0, 2 * Math.PI * 1 / Ts); - - return Bode(omega_vec); - } - - public Complex[] Bode(int nPoints, out double[] omega_vec) - { - // substituting z = exp(j * omega * Ts) - omega_vec = Generate.LinearSpaced(nPoints, 0, 2 * Math.PI * 1 / Ts); - - - return Bode(omega_vec); - } - - public Complex[] Bode(double[] omega_vec) - { - - var nPoints = omega_vec.Length; - - double omega; - double expVal; - Complex zVal; - Complex denVal; - Complex numVal; - - var bodeVal = new Complex[nPoints]; - - for (int idx = 0; idx < nPoints; idx++) - { - - - omega = omega_vec[idx]; - - zVal = new Complex(0.0, 0.0); - - denVal = new Complex(0.0, 0.0); - for (int ii = 0; ii < a.Length; ii++) - { - expVal = ii * omega * Ts; - zVal = new Complex(0.0, expVal); - - denVal += a[ii] * zVal.Exp(); - } - - numVal = new Complex(0.0, 0.0); - for (int ii = 0; ii < b.Length; ii++) - { - expVal = ii * omega * Ts; - zVal = new Complex(0.0, expVal); - - numVal += b[ii] * zVal.Exp(); - } - bodeVal[idx] = numVal / denVal; - } - - return bodeVal; - } - - /// The poles resulting from the denominator Polynomial root - public Complex[] GetPoles() - { - Polynomial a_poly = new Polynomial(a, isFlip:true); - Complex[] r = a_poly.GetRoots(); - return r; - } - - /// The zeros resulting from the nominator Polynomial root - public Complex[] GetZeros() - { - Polynomial b_poly = new Polynomial(b, isFlip:true); - Complex[] r = b_poly.GetRoots(); - return r; - } - - - /// - /// calculate the number of steps the system will need until it can be assumed to be settled - /// - /// tolerance in decimal percent at which to assume that the system is settled (default = 0.3) - /// maximum number of steps to simulate (default = 500000) - /// number of steps at which the system is assumed to be settled, or 0 if unstable - public int CalcSettlingSteps(double tol = 0.03, int n_max = 500000) - { - - // init settling time as zero for never settled - int n_sttl = 0; - - // if the system is unstable return zero since the system will never be settled - if (this.IsStable() == false) - return 0; - - int n_sim = 0; - - double[] dampVals = GetDampings(out double[] EigenFrequencys); - - double dampWorst = dampVals.Min(); - - //for (int ii = 1; ii < dampVals.Length; ii++) - // dampWorst = dampWorst * dampVals[ii]; - - double t_simFull; - - // appromate a settling time based on damping - var t_stlDamp = -Math.Log(tol) / dampWorst; - - // approximate a settling time from time constants - var tau = new double[dampVals.Length]; - for (int ii = 0; ii < tau.Length; ii++) - tau[ii] = 1.0 / (dampVals[ii] * EigenFrequencys[ii]); - - // approx after 5 * biggest time constant - var t_stlTimeConst = tau.Max() * 5; - - // choose bigger approximation - t_simFull = Math.Max(t_stlTimeConst, t_stlDamp); - - // recalculate to number of steps - int nStepsBase = (int)Math.Ceiling(t_simFull / Ts); - - - // simulate impulse responses with n*10*nStepsBase time steps - // incrementing n if necessary until steady state is reached - n_sim = nStepsBase <= 0 ? 5 : nStepsBase; - int count = 0; - while (count < 10 && n_sttl == 0) - { - if (n_sim > n_max) - return n_sttl; - - n_sim = 10 * n_sim; - - double[] dirac_sim = new double[n_sim]; - dirac_sim[0] = 1.0; - var tmp_outp = this.CalcResponse(dirac_sim); - - int idxPos = n_sim - 1; - - // find first step beeing bigger than tolerance - while (idxPos > 0 && n_sttl == 0) - { - if (tmp_outp[idxPos] > tol) - n_sttl = idxPos; - - idxPos--; - } - count++; - } - return n_sttl; - - } - - #endregion Dynamics - - - #region Dampings - /// - /// gets the damping coefficients from this transfer function, - /// since all transfer functions so far are discrete time, - /// these values do not directly translate to lambda. - /// the theoretical recalculation is: - /// Z = -cos(angle(log(lambda))) - /// - /// Array of damping values for this transfer function - public double[] GetDampings() - { - return GetDampings(out double[] f); - } - - /// - /// gets the damping coefficients from this transfer function, - /// since all transfer functions so far are discrete time, - /// these values do not directly translate to lambda. - /// the theoretical recalculation is: - /// Z = -cos(angle(log(lambda))) - /// - /// Array of damping values for this transfer function - public double[] GetDampings(out double[] wn) - { - - var r = GetPoles().Clone() as Complex[]; - var s = new Complex[r.Length]; - var f = new double[r.Length]; - var z = new double[r.Length]; - - for (int idx = 0; idx < r.Length; idx++) - { - s[idx] = Complex.Log(r[idx]) / Ts; - f[idx] = s[idx].Magnitude; - z[idx] = -s[idx].Real / f[idx]; - } - - wn = (double[])f.Clone(); - - return z; - - } - - - #endregion Dampings - - - #region displaying - /// - /// - /// - /// - public string DispTF() - { - return (DispTF(this.b, this.a, this.Name, this.variable.Substring(0, variable.Length - 1))); - } - - public string NumString() - { - var varStr = this.variable.Substring(0, variable.Length - 1); - var num = b.Clone() as double[]; - return getFractString(num, varStr); - } - - public string DenString() - { - var varStr = this.variable.Substring(0, variable.Length - 1); - var den = a.Clone() as double[]; - return getFractString(den, varStr); - } - - - private static string getFractString(double[] num, string varStr) - { - string str1; - string str2; - string strNum = ""; - for (int item = 0; item < num.Length; item++) - { - if (num[item] == 0) - continue; - - //str2 = Math.Abs(num[item]).ToString(); - str2 = Math.Abs(num[item]).ToString("0.######"); - if (item == 0) - { - if (num[item] < 0) - str1 = "-"; - else - str1 = ""; - - strNum = String.Concat(strNum, str1, str2); - } - else - { - if (num[item] > 0) - str1 = " + "; - else - str1 = " - "; - strNum = String.Concat(strNum, str1, str2, varStr, item.ToString()); - } - } - - if (strNum.StartsWith("+") || strNum.StartsWith(" ")) - strNum = strNum.Substring(1); - - return strNum; - } - - - - /// - /// - /// - /// - /// - /// - /// - public static string DispTF(double[] num, double[] den, string name, string varStr = " q^-") - { - - string strNum = getFractString(num, varStr); - string strDen = getFractString(den, varStr); - string strHead = ""; - - if (String.IsNullOrEmpty(name)) - strHead = "TF = "; - else - strHead = name; - - int nbar = Math.Max(strDen.Length, strNum.Length); - - string strBar = new String('-', nbar); - string strOut = String.Concat(strHead, "\n\n", strNum, '\n', strBar, '\n', strDen); - - return (strOut); - } - - #endregion displaying - - } -} From 48aa6fe6c4b9fd8ef3237198bca1f3c524c059ba Mon Sep 17 00:00:00 2001 From: Tobias Glaubach Date: Tue, 17 Jul 2018 23:57:23 +0200 Subject: [PATCH 18/18] grooming --- src/Numerics/Polynomial.cs | 1 - 1 file changed, 1 deletion(-) diff --git a/src/Numerics/Polynomial.cs b/src/Numerics/Polynomial.cs index ce2e6cd9..64a57887 100644 --- a/src/Numerics/Polynomial.cs +++ b/src/Numerics/Polynomial.cs @@ -557,7 +557,6 @@ namespace MathNet.Numerics /// /// Division of two polynomials returning the quotient-with-remainder of the two polynomials given /// - /// left polynomial /// right polynomial /// a tuple holding quotient in first and remainder in second public Tuple DivideLong(Polynomial b)