123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692169316941695169616971698169917001701170217031704170517061707170817091710171117121713171417151716171717181719172017211722172317241725172617271728172917301731173217331734173517361737173817391740174117421743174417451746174717481749175017511752175317541755175617571758175917601761176217631764176517661767176817691770177117721773177417751776177717781779178017811782178317841785178617871788178917901791179217931794179517961797179817991800180118021803180418051806180718081809181018111812181318141815181618171818181918201821182218231824182518261827182818291830183118321833183418351836183718381839184018411842184318441845184618471848184918501851185218531854185518561857185818591860186118621863186418651866186718681869187018711872187318741875187618771878187918801881188218831884188518861887188818891890189118921893189418951896189718981899190019011902190319041905190619071908190919101911191219131914191519161917191819191920192119221923192419251926192719281929193019311932193319341935193619371938193919401941194219431944194519461947194819491950195119521953 |
- /* Polynomial.cpp
- *
- * Copyright (C) 1993-2018 David Weenink
- *
- * This code is free software; you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation; either version 2 of the License, or (at
- * your option) any later version.
- *
- * This code is distributed in the hope that it will be useful, but
- * WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this work. If not, see <http://www.gnu.org/licenses/>.
- */
- /*
- djmw 20020813 GPL header
- djmw 20030619 Added SVD_compute before SVD_solve
- djmw 20060510 Polynomial_to_Roots: changed behaviour. All roots found are now saved.
- In previous version a nullptr pointer was returned. New error messages.
- djmw 20071012 Added: o_CAN_WRITE_AS_ENCODING.h
- djmw 20071201 Melder_warning<n>
- djmw 20080122 float -> double
- djmw 20110304 Thing_new
- */
- #include "Polynomial.h"
- #include "SVD.h"
- #include "NUMclapack.h"
- #include "TableOfReal_extensions.h"
- #include "NUMmachar.h"
- #include "oo_DESTROY.h"
- #include "Polynomial_def.h"
- #include "oo_COPY.h"
- #include "Polynomial_def.h"
- #include "oo_EQUAL.h"
- #include "Polynomial_def.h"
- #include "oo_CAN_WRITE_AS_ENCODING.h"
- #include "Polynomial_def.h"
- #include "oo_WRITE_TEXT.h"
- #include "Polynomial_def.h"
- #include "oo_WRITE_BINARY.h"
- #include "Polynomial_def.h"
- #include "oo_READ_TEXT.h"
- #include "Polynomial_def.h"
- #include "oo_READ_BINARY.h"
- #include "Polynomial_def.h"
- #include "oo_DESCRIPTION.h"
- #include "Polynomial_def.h"
- /* Evaluate polynomial and derivative jointly
- c [1..n] -> degree n-1 !!
- */
- void Polynomial_evaluateWithDerivative (Polynomial me, double x, double *out_f, double *out_df) {
- longdouble p = my coefficients [my numberOfCoefficients], dp = 0.0, xc = x;
- for (integer i = my numberOfCoefficients - 1; i > 0; i --) {
- dp = dp * xc + p;
- p = p * xc + my coefficients [i];
- }
- if (out_f) *out_f = (double) p;
- if (out_df) *out_df = (double) dp;
- }
- /* Get value and derivative */
- static void Polynomial_evaluateWithDerivative_z (Polynomial me, dcomplex *in_z, dcomplex *out_p, dcomplex *out_dp) {
- longdouble pr = my coefficients [my numberOfCoefficients], pi = 0.0;
- longdouble dpr = 0.0, dpi = 0.0, x = in_z -> re, y = in_z -> im;
- for (integer i = my numberOfCoefficients - 1; i > 0; i --) {
- longdouble tr = dpr;
- dpr = dpr * x - dpi * y + pr;
- dpi = tr * y + dpi * x + pi;
- tr = pr;
- pr = pr * x - pi * y + my coefficients [i];
- pi = tr * y + pi * x;
- }
- if (out_p) *out_p = { (double) pr, (double) pi };
- if (out_dp) *out_dp = { (double) dpr, (double) dpi };
- }
- void Polynomial_evaluateDerivatives (Polynomial me, double x, double *derivatives, integer numberOfDerivatives) {
- /* Evaluate polynomial c [1]+c [2]*x+...degree*x^degree in derivative [0] and derivatives [1..numberOfDerivatives] */
- integer degree = my numberOfCoefficients - 1;
- numberOfDerivatives = numberOfDerivatives > degree ? degree : numberOfDerivatives;
-
- derivatives [0] = my coefficients [my numberOfCoefficients];
- for (integer j = 1; j <= numberOfDerivatives; j ++)
- derivatives [j] = 0.0;
- for (integer i = degree - 1; i >= 0; i--) {
- integer n =
- numberOfDerivatives < degree - i ? numberOfDerivatives : degree - i;
- for (integer j = n; j >= 1; j --)
- derivatives [j] = derivatives [j] * x + derivatives [j - 1];
- derivatives [0] = derivatives [0] * x + my coefficients [i + 1]; // evaluate polynomial (Horner)
- }
- double fact = 1.0;
- for (integer j = 2; j <= numberOfDerivatives; j ++) {
- fact *= j;
- derivatives [j] *= fact;
- }
- }
- /*
- void polynomial_divide (double *u, integer m, double *v, integer n, double *q, double *r);
- Purpose:
- Find the quotient q(x) and the remainder r(x) polynomials that result from the division of
- the polynomial u(x) = u [1] + u [2]*x^1 + u [3]*x^2 + ... + u [m]*x^(m-1) by the
- polynomial v(x) = v [1] + v [2]*x^1 + v [3]*x^2 + ... + v [n]*x^(n-1), such that
- u(x) = v(x)*q(x) + r(x).
- The arrays u, v, q and r have to be dimensioned as u [1...m], v [1..n], q [1...m] and r [1...m],
- respectively.
- On return, the q [1..m-n] and r [1..n-1] contain the quotient and the remainder
- polynomial coefficients, repectively.
- See Knuth, The Art of Computer Programming, Volume 2: Seminumerical algorithms,
- Third edition, Section 4.6.1 Algorithm D (the algorithm as described has been modified
- to prevent overwriting of the u-polynomial).
- */
- static void polynomial_divide (double *u, integer m, double *v, integer n, double *q, double *r) {
- // Copy u [1..m] into r [1..n] to prevent overwriting of u.
- // Put the q coefficients to zero for cases n > m.
- for (integer k = 1; k <= m; k ++) {
- r [k] = u [k];
- q [k] = 0.0;
- }
- for (integer k = m - n + 1; k > 0; k --) { /* D1 */
- q [k] = r [n + k - 1] / v [n]; /* D2 with u -> r*/
- for (integer j = n + k - 1; j >= k; j --) {
- r [j] -= q [k] * v [j - k + 1];
- }
- }
- }
- static void Polynomial_polish_realroot (Polynomial me, double *x, integer maxit) {
- double xbest = *x, pmin = 1e308;
- if (! NUMfpp)
- NUMmachar ();
- for (integer i = 1; i <= maxit; i ++) {
- double p, dp;
- Polynomial_evaluateWithDerivative (me, *x, & p, & dp);
- double fabsp = fabs (p);
- if (fabsp > pmin || fabs (fabsp - pmin) < NUMfpp -> eps) {
- /*
- We stop, because the approximation is getting worse or we cannot get any closer.
- Return the previous (hitherto best) value for x.
- */
- *x = xbest;
- return;
- }
- pmin = fabsp;
- xbest = *x;
- if (fabs (dp) == 0.0)
- return;
- double dx = p / dp; // Newton-Raphson
- *x -= dx;
- }
- // Melder_throw (U"Maximum number of iterations exceeded.");
- }
- static void Polynomial_polish_complexroot_nr (Polynomial me, dcomplex *z, integer maxit) {
- dcomplex zbest = *z;
- double pmin = 1e308;
- if (! NUMfpp)
- NUMmachar ();
- for (integer i = 1; i <= maxit; i ++) {
- dcomplex p, dp;
- Polynomial_evaluateWithDerivative_z (me, z, &p, &dp);
- double fabsp = dcomplex_abs (p);
- if (fabsp > pmin || fabs (fabsp - pmin) < NUMfpp -> eps) {
- /*
- We stop, because the approximation is getting worse.
- Return the previous (hitherto best) value for z.
- */
- *z = zbest;
- return;
- }
- pmin = fabsp;
- zbest = *z;
- if (dcomplex_abs (dp) == 0.0)
- return;
- dcomplex dz = dcomplex_div (p, dp); // Newton-Raphson
- *z = dcomplex_sub (*z, dz);
- }
- // Melder_throw (U"Maximum number of iterations exceeded.");
- }
- /*
- Symbolic evaluation of polynomial coefficients.
- Recurrence: P [n] = (a [n] * x + b [n]) P [n-1] + c [n] P [n-2],
- where P [n] is any orthogonal polynomial of degree n.
- P [n] is an array of coefficients p [k] representing: p [1] + p [2] x + ... p [n+1] x^n.
- Preconditions:
- degree > 1;
- pnm1 : polynomial of degree n - 1
- pnm2 : polynomial of degree n - 2
- */
- static void NUMpolynomial_recurrence (double *pn, integer degree, double a, double b, double c, double *pnm1, double *pnm2) {
- Melder_assert (degree > 1);
- pn [1] = b * pnm1 [1] + c * pnm2 [1];
- for (integer i = 2; i <= degree - 1; i ++) {
- pn [i] = a * pnm1 [i - 1] + b * pnm1 [i] + c * pnm2 [i];
- }
- pn [degree] = a * pnm1 [degree - 1] + b * pnm1 [degree];
- pn [degree + 1] = a * pnm1 [degree];
- }
- /* frozen [1..ma] */
- static void svdcvm (double **v, integer mfit, integer ma, int *frozen, double *w, double **cvm) {
- autoNUMvector<double> wti (1, mfit);
- for (integer i = 1; i <= mfit; i ++) {
- if (w [i] != 0.0) {
- wti [i] = 1.0 / (w [i] * w [i]);
- } else {
- ; // TODO: write up an explanation for why it is not necessary to do anything if w [i] is zero
- }
- }
- for (integer i = 1; i <= mfit; i ++) {
- for (integer j = 1; j <= i; j ++) {
- longdouble sum = 0.0;
- for (integer k = 1; k <= mfit; k ++) {
- sum += v [i] [k] * v [j] [k] * wti [k];
- }
- cvm [j] [i] = cvm [i] [j] = (double) sum;
- }
- }
- for (integer i = mfit + 1; i <= ma; i ++) {
- for (integer j = 1; j <= i; j ++) {
- cvm [j] [i] = cvm [i] [j] = 0.0;
- }
- }
- integer k = mfit;
- for (integer j = ma; j > 0; j --) {
- // if (! frozen || ! frozen [i]) why i?? TODO
- if (! frozen || ! frozen [j]) {
- for (integer i = 1; i <= ma; i ++) {
- double t = cvm [i] [k];
- cvm [i] [k] = cvm [i] [j];
- cvm [i] [j] = t;
- }
- for (integer i = 1; i <= ma; i ++) {
- double t = cvm [k] [i];
- cvm [k] [i] = cvm [j] [i];
- cvm [j] [i] = t;
- }
- k --;
- }
- }
- }
- /********* FunctionTerms *********************************************/
- Thing_implement (FunctionTerms, Function, 0);
- double structFunctionTerms :: v_evaluate (double /* x */) {
- return undefined;
- }
- dcomplex structFunctionTerms :: v_evaluate_z (dcomplex /* z */) {
- return { undefined, undefined };
- }
- void structFunctionTerms :: v_evaluateTerms (double /* x */, double terms []) {
- for (integer i = 1; i <= numberOfCoefficients; i ++) {
- terms [i] = undefined;
- }
- }
- integer structFunctionTerms :: v_getDegree () {
- return numberOfCoefficients - 1;
- }
- void structFunctionTerms :: v_getExtrema (double x1, double x2, double *p_xmin, double *p_ymin, double *p_xmax, double *p_ymax) { // David, geen aparte naam hier nodig: ???
- integer numberOfPoints = 1000;
- // Melder_warning (L"defaultGetExtrema: extrema calculated by sampling the interval");
- double x = x1, dx = (x2 - x1) / (numberOfPoints - 1);
- double xmn = x, xmx = xmn, ymn = v_evaluate (x), ymx = ymn; // xmin, xmax .. would shadow member
- for (integer i = 2; i <= numberOfPoints; i ++) {
- x += dx;
- double y = v_evaluate (x);
- if (y > ymx) {
- ymx = y; xmx = x;
- } else if (y < ymn) {
- ymn = y; xmn = x;
- }
- }
- if (p_xmin)
- *p_xmin = xmn;
- if (p_xmax)
- *p_xmax = xmx;
- if (p_ymin)
- *p_ymin = ymn;
- if (p_ymax)
- *p_ymax = ymx;
- }
- static inline void FunctionTerms_extendCapacityIf (FunctionTerms me, integer minimum) {
- if (my _capacity < minimum) {
- NUMvector_append <double> (& my coefficients, 1, & minimum);
- my _capacity = minimum;
- }
- }
- void FunctionTerms_init (FunctionTerms me, double xmin, double xmax, integer numberOfCoefficients) {
- my coefficients = NUMvector<double> (1, numberOfCoefficients);
- my numberOfCoefficients = numberOfCoefficients;
- my _capacity = numberOfCoefficients;
- my xmin = xmin; my xmax = xmax;
- }
- autoFunctionTerms FunctionTerms_create (double xmin, double xmax, integer numberOfCoefficients) {
- try {
- autoFunctionTerms me = Thing_new (FunctionTerms);
- FunctionTerms_init (me.get(), xmin, xmax, numberOfCoefficients);
- return me;
- } catch (MelderError) {
- Melder_throw (U"FunctionTerms not created.");
- }
- }
- void FunctionTerms_initFromString (FunctionTerms me, double xmin, double xmax, conststring32 s, bool allowTrailingZeros) {
- integer numberOfCoefficients;
- autoVEC numbers = VEC_createFromString (s);
- if (! allowTrailingZeros) {
- while (numbers [numbers.size] == 0 && numberOfCoefficients > 1)
- numberOfCoefficients --;
- }
- FunctionTerms_init (me, xmin, xmax, numberOfCoefficients);
- NUMvector_copyElements (numbers.at, my coefficients, 1, numberOfCoefficients);
- }
- integer FunctionTerms_getDegree (FunctionTerms me) {
- return my v_getDegree ();
- }
- void FunctionTerms_setDomain (FunctionTerms me, double xmin, double xmax) {
- my xmin = xmin; my xmax = xmax;
- }
- double FunctionTerms_evaluate (FunctionTerms me, double x) {
- return my v_evaluate (x);
- }
- dcomplex FunctionTerms_evaluate_z (FunctionTerms me, dcomplex z) {
- return my v_evaluate_z (z);
- }
- void FunctionTerms_evaluateTerms (FunctionTerms me, double x, double terms []) {
- my v_evaluateTerms (x, terms);
- }
- void FunctionTerms_getExtrema (FunctionTerms me, double x1, double x2, double *xmin, double *ymin, double *xmax, double *ymax) {
- if (x2 <= x1) {
- x1 = my xmin;
- x2 = my xmax;
- }
- my v_getExtrema (x1, x2, xmin, ymin, xmax, ymax);
- }
- double FunctionTerms_getMinimum (FunctionTerms me, double x1, double x2) {
- double xmin, xmax, ymin, ymax;
- FunctionTerms_getExtrema (me, x1, x2, &xmin, &ymin, &xmax, &ymax);
- return ymin;
- }
- double FunctionTerms_getXOfMinimum (FunctionTerms me, double x1, double x2) {
- double xmin, xmax, ymin, ymax;
- FunctionTerms_getExtrema (me, x1, x2, &xmin, &ymin, &xmax, &ymax);
- return xmin;
- }
- double FunctionTerms_getMaximum (FunctionTerms me, double x1, double x2) {
- double xmin, xmax, ymin, ymax;
- FunctionTerms_getExtrema (me, x1, x2, &xmin, &ymin, &xmax, &ymax);
- return ymax;
- }
- double FunctionTerms_getXOfMaximum (FunctionTerms me, double x1, double x2) {
- double xmin, xmax, ymin, ymax;
- FunctionTerms_getExtrema (me, x1, x2, &xmin, &ymin, &xmax, &ymax);
- return xmax;
- }
- static void Graphics_polyline_clipTopBottom (Graphics g, double *x, double *y, integer numberOfPoints, double ymin, double ymax) {
- integer index = 0;
- if (numberOfPoints < 2)
- return;
- double x1 = x [0], y1 = y [0];
- double xb = x1, yb = y1;
- for (integer i = 1; i < numberOfPoints; i ++) {
- double x2 = x [i], y2 = y [i];
- if (! ( (y1 > ymax && y2 > ymax) || (y1 < ymin && y2 < ymin))) {
- double dxy = (x2 - x1) / (y1 - y2);
- double xcros_max = x1 - (ymax - y1) * dxy;
- double xcros_min = x1 - (ymin - y1) * dxy;
- if (y1 > ymax && y2 < ymax) {
- // Line enters from above: start new segment. Save start values.
- xb = x [i - 1];
- yb = y [i - 1];
- index = i - 1;
- y [i - 1] = ymax;
- x [i - 1] = xcros_max;
- }
- if (y1 > ymin && y2 < ymin) {
- // Line leaves at bottom: draw segment. Save end values and restore them
- // Origin of arrays for Graphics_polyline are at element 0 !!!
- double xe = x [i], ye = y [i];
- y [i] = ymin;
- x [i] = xcros_min;
- Graphics_polyline (g, i - index + 1, x + index, y + index);
- x [index] = xb;
- y [index] = yb;
- x [i] = xe;
- y [i] = ye;
- }
- if (y1 < ymin && y2 > ymin) {
- // Line enters from below: start new segment. Save start values
- xb = x [i - 1];
- yb = y [i - 1];
- index = i - 1;
- y [i - 1] = ymin;
- x [i - 1] = xcros_min;
- }
- if (y1 < ymax && y2 > ymax) {
- // Line leaves at top: draw segment. Save and restore
- double xe = x [i], ye = y [i];
- y [i] = ymax;
- x [i] = xcros_max;
- Graphics_polyline (g, i - index + 1, x + index, y + index);
- x [index] = xb;
- y [index] = yb;
- x [i] = xe;
- y [i] = ye;
- }
- } else {
- index = i;
- }
- y1 = y2;
- x1 = x2;
- }
- if (index < numberOfPoints - 1) {
- Graphics_polyline (g, numberOfPoints - index, x + index, y + index);
- x [index] = xb; y [index] = yb;
- }
- }
- void FunctionTerms_draw (FunctionTerms me, Graphics g, double xmin, double xmax, double ymin, double ymax, int extrapolate, int garnish) {
- integer numberOfPoints = 1000;
- autoNUMvector<double> y (1, numberOfPoints + 1);
- autoNUMvector<double> x (1, numberOfPoints + 1);
- if (xmax <= xmin) {
- xmin = my xmin;
- xmax = my xmax;
- }
- double fxmin = xmin, fxmax = xmax;
- if (! extrapolate) {
- if (xmax < my xmin || xmin > my xmax)
- return;
- if (xmin < my xmin)
- fxmin = my xmin;
- if (xmax > my xmax)
- fxmax = my xmax;
- }
- if (ymax <= ymin) {
- double x1, x2;
- FunctionTerms_getExtrema (me, fxmin, fxmax, &x1, &ymin, &x2, &ymax);
- }
- Graphics_setInner (g);
- Graphics_setWindow (g, xmin, xmax, ymin, ymax);
- // Draw only the parts within [fxmin, fxmax] X [ymin, ymax].
- double dx = (fxmax - fxmin) / (numberOfPoints - 1);
- for (integer i = 1; i <= numberOfPoints; i ++) {
- x [i] = fxmin + (i - 1.0) * dx;
- y [i] = FunctionTerms_evaluate (me, x [i]);
- }
- //Graphics_polyline_clipTopBottom (g, x+1, y+1, numberOfPoints, ymin, ymax);
- Graphics_polyline_clipTopBottom (g, &x [1], &y [1], numberOfPoints, ymin, ymax);
- Graphics_unsetInner (g);
- if (garnish) {
- Graphics_drawInnerBox (g);
- Graphics_marksBottom (g, 2, true, true, false);
- Graphics_marksLeft (g, 2, true, true, false);
- }
- }
- void FunctionTerms_drawBasisFunction (FunctionTerms me, Graphics g, integer index, double xmin, double xmax, double ymin, double ymax, int extrapolate, int garnish) {
- if (index < 1 || index > my numberOfCoefficients)
- return;
- autoFunctionTerms thee = Data_copy (me);
- for (integer i = 1; i <= my numberOfCoefficients; i ++)
- thy coefficients [i] = 0.0;
- thy coefficients [index] = 1.0;
- thy numberOfCoefficients = index;
- FunctionTerms_draw (thee.get(), g, xmin, xmax, ymin, ymax, extrapolate, garnish);
- }
- void FunctionTerms_setCoefficient (FunctionTerms me, integer index, double value) {
- Melder_require (index > 0 && index <= my numberOfCoefficients, U"Index out of range [1, ", my numberOfCoefficients, U"].");
- Melder_require (value == 0.0 && index < my numberOfCoefficients, U"You should not remove the highest degree term.");
- my coefficients [index] = value;
- }
- /********** Polynomial ***********************************************/
- Thing_implement (Polynomial, FunctionTerms, 1);
- double structPolynomial :: v_evaluate (double x) {
- longdouble p = coefficients [numberOfCoefficients];
- for (integer i = numberOfCoefficients - 1; i > 0; i --)
- p = p * x + coefficients [i];
- return (double) p;
- }
- dcomplex structPolynomial :: v_evaluate_z (dcomplex z) {
- longdouble x = z.re, y = z.im;
- longdouble pr = coefficients [numberOfCoefficients];
- longdouble pi = 0.0;
- for (integer i = numberOfCoefficients - 1; i > 0; i --) {
- longdouble prtmp = pr;
- pr = pr * x - pi * y + coefficients [i];
- pi = prtmp * y + pi * x;
- }
- return { (double) pr, (double) pi };
- }
- void structPolynomial :: v_evaluateTerms (double x, double terms []) {
- terms [1] = 1;
- for (integer i = 2; i <= numberOfCoefficients; i ++)
- terms [i] = terms [i - 1] * x;
- }
- void structPolynomial :: v_getExtrema (double x1, double x2, double *p_xmin, double *p_ymin, double *p_xmax, double *p_ymax) {
- try {
- integer degree = numberOfCoefficients - 1;
- double xmn = x1, ymn = v_evaluate (x1);
- double xmx = x2, ymx = v_evaluate (x2);
- if (ymn > ymx) {
- std::swap (ymn, ymx);
- std::swap (xmn, xmx);
- }
- if (degree < 2)
- return;
- autoPolynomial d = Polynomial_getDerivative (this);
- autoRoots r = Polynomial_to_Roots (d.get());
- for (integer i = 1; i <= degree - 1; i ++) {
- double x = (r -> v [i]).re;
- if (x > x1 && x < x2) {
- double y = v_evaluate (x);
- if (y > ymx) {
- ymx = y;
- xmx = x;
- } else if (y < ymn) {
- ymn = y;
- xmn = x;
- }
- }
- }
- if (p_xmin)
- *p_xmin = xmn;
- if (p_xmax)
- *p_xmax = xmx;
- if (p_ymin)
- *p_ymin = ymn;
- if (p_ymax)
- *p_ymax = ymx;
- } catch (MelderError) {
- structFunctionTerms :: v_getExtrema (x1, x2, p_xmin, p_ymin, p_xmax, p_ymax);
- Melder_clearError ();
- }
- }
- autoPolynomial Polynomial_create (double xmin, double xmax, integer degree) {
- try {
- autoPolynomial me = Thing_new (Polynomial);
- FunctionTerms_init (me.get(), xmin, xmax, degree + 1);
- return me;
- } catch (MelderError) {
- Melder_throw (U"Polynomial not created.");
- }
- }
- autoPolynomial Polynomial_createFromString (double lxmin, double lxmax, conststring32 s) {
- try {
- autoPolynomial me = Thing_new (Polynomial);
- FunctionTerms_initFromString (me.get(), lxmin, lxmax, s, 0);
- return me;
- } catch (MelderError) {
- Melder_throw (U"Polynomial not created from string.");
- }
- }
- void Polynomial_scaleCoefficients_monic (Polynomial me) {
- double cn = my coefficients [my numberOfCoefficients];
- if (cn == 1 || my numberOfCoefficients <= 1)
- return;
- for (integer i = 1; i < my numberOfCoefficients; i ++)
- my coefficients [i] /= cn;
- my coefficients [my numberOfCoefficients] = 1.0;
- }
- /*
- Transform the polynomial as if the domain were [xmin, xmax].
- Some polynomials (Legendre) are defined on the domain [-1,1]. The domain
- for x may be extended to [xmin, xmax] by a transformation such as
- x' = (2 * x - (xmin + xmax)) / (xmax - xmin) -1 < x' < x.
- This procedure transforms x' back to x.
- */
- autoPolynomial Polynomial_scaleX (Polynomial me, double xmin, double xmax) {
- try {
- Melder_assert (xmin < xmax);
- autoPolynomial thee = Polynomial_create (xmin, xmax, my numberOfCoefficients - 1);
- thy coefficients [1] = my coefficients [1];
- if (my numberOfCoefficients == 1)
- return thee;
- // x = a x + b
- // Constraints:
- // my xmin = a xmin + b; a = (my xmin - my xmax) / (xmin - xmax);
- // my xmax = a xmax + b; b = my xmin - a * xmin
- double a = (my xmin - my xmax) / (xmin - xmax);
- double b = my xmin - a * xmin;
- thy coefficients [2] = my coefficients [2] * a;
- thy coefficients [1] += my coefficients [2] * b;
- if (my numberOfCoefficients == 2)
- return thee;
- autoNUMvector<double> buf (1, 3 * my numberOfCoefficients);
- double *pn = buf.peek();
- double *pnm1 = pn + my numberOfCoefficients;
- double *pnm2 = pnm1 + my numberOfCoefficients;
- // Start the recursion: P [1] = a x + b; P [0] = 1;
- pnm1 [2] = a;
- pnm1 [1] = b;
- pnm2 [1] = 1;
- for (integer n = 2; n <= my numberOfCoefficients - 1; n ++) {
- double *t1 = pnm1, *t2 = pnm2;
- NUMpolynomial_recurrence (pn, n, a, b, 0, pnm1, pnm2);
- if (my coefficients [n + 1] != 0) {
- for (integer j = 1; j <= n + 1; j ++)
- thy coefficients [j] += my coefficients [n + 1] * pn [j];
- }
- pnm1 = pn;
- pnm2 = t1;
- pn = t2;
- }
- return thee;
- } catch (MelderError) {
- Melder_throw (me, U": cannot scale.");
- }
- }
- double Polynomial_evaluate (Polynomial me, double x) {
- return my v_evaluate (x);
- }
- dcomplex Polynomial_evaluate_z (Polynomial me, dcomplex z) {
- return my v_evaluate_z (z);
- }
- static void Polynomial_evaluate_z_cart (Polynomial me, double r, double phi, double *re, double *im) {
- double rn = 1;
- *re = my coefficients [1];
- *im = 0.0;
- if (r == 0.0)
- return;
- for (integer i = 2; i <= my numberOfCoefficients; i ++) {
- rn *= r;
- double arg = (i - 1) * phi;
- *re += my coefficients [i] * rn * cos (arg);
- *im += my coefficients [i] * rn * sin (arg);
- }
- }
- autoPolynomial Polynomial_getDerivative (Polynomial me) {
- try {
- if (my numberOfCoefficients == 1)
- return Polynomial_create (my xmin, my xmax, 0);
- autoPolynomial thee = Polynomial_create (my xmin, my xmax, my numberOfCoefficients - 2);
- for (integer i = 1; i <= thy numberOfCoefficients; i ++)
- thy coefficients [i] = i * my coefficients [i + 1];
- return thee;
- } catch (MelderError) {
- Melder_throw (me, U": no derivative created.");
- }
- }
- autoPolynomial Polynomial_getPrimitive (Polynomial me, double constant) {
- try {
- autoPolynomial thee = Polynomial_create (my xmin, my xmax, my numberOfCoefficients);
- for (integer i = 1; i <= my numberOfCoefficients; i ++)
- thy coefficients [i + 1] = my coefficients [i] / i;
- thy coefficients [1] = constant;
- return thee;
- } catch (MelderError) {
- Melder_throw (me, U": no primitive created.");
- }
- }
- /* P(x)= (x-roots [1])*(x-roots [2])*..*(x-roots [numberOfRoots]) */
- void Polynomial_initFromRealRoots (Polynomial me, constVEC roots) {
- try {
- FunctionTerms_extendCapacityIf (me, roots.size + 1);
- double *c = & my coefficients [1];
- integer n = 1;
- c [0] = - roots [1];
- c [1] = 1.0;
- for (integer i = 2; i <= roots.size; i ++) {
- c [n + 1] = c [n];
- for (integer j = n; j >= 1; j --)
- c [j] = c [j - 1] - c [j] * roots [i];
- c [0] *= -roots [i];
- n ++;
- }
- my numberOfCoefficients = n + 1;
- } catch (MelderError) {
- Melder_throw (me, U": not initalized from real roots.");
- }
- }
- autoPolynomial Polynomial_createFromRealRootsString (double xmin, double xmax, conststring32 s) {
- try {
- autoPolynomial me = Thing_new (Polynomial);
- autoVEC roots = VEC_createFromString (s);
- FunctionTerms_init (me.get(), xmin, xmax, roots.size + 1);
- Polynomial_initFromRealRoots (me.get(), roots.get());
- return me;
- } catch (MelderError) {
- Melder_throw (U"Polynomial not created from roots.");
- }
-
- }
- /* Product (i=1; a.size; (1 + a[i]*x + x^2)
- * Postcondition : my numberOfCoeffcients = 2*a.size+1
- */
- void Polynomial_initFromProductOfSecondOrderTerms (Polynomial me, constVEC a) {
- FunctionTerms_extendCapacityIf (me, 2 * a.size + 1);
- my coefficients [1] = my coefficients [3] = 1.0;
- my coefficients [2] = a [1];
- integer ncoef = 3;
- for (integer i = 2; i <= a.size; i ++) {
- my coefficients [ncoef + 1] = a [i] * my coefficients [ncoef] + my coefficients [ncoef - 1];
- my coefficients [ncoef + 2] = my coefficients [ncoef];
- for (integer j = ncoef; j > 2; j --)
- my coefficients [j] += a [i] * my coefficients [j - 1] + my coefficients [j - 2];
- my coefficients [2] += a [i]; // a [i] * my coefficients [1]
- ncoef += 2;
- }
- my numberOfCoefficients = ncoef;
- }
- autoPolynomial Polynomial_createFromProductOfSecondOrderTermsString (double xmin, double xmax, conststring32 s) {
- try {
- autoPolynomial me = Thing_new (Polynomial);
- autoVEC a = VEC_createFromString (s);
- FunctionTerms_init (me.get(), xmin, xmax, 2 * a.size + 1);
- Polynomial_initFromProductOfSecondOrderTerms (me.get(), a.get());
- return me;
- } catch (MelderError) {
- Melder_throw (U"Polynomial not created from second order terms string.");
- }
- }
- double Polynomial_getArea (Polynomial me, double xmin, double xmax) {
- if (xmax >= xmin) {
- xmin = my xmin;
- xmax = my xmax;
- }
- autoPolynomial p = Polynomial_getPrimitive (me, 0);
- double area = FunctionTerms_evaluate (p.get(), xmax) - FunctionTerms_evaluate (p.get(), xmin);
- return area;
- }
- /* P(x) * (x-a)
- * Postcondition: my numberOfCoefficients = old_numberOfCoefficients + 1
- */
- void Polynomial_multiply_firstOrderFactor (Polynomial me, double factor) {
- integer n = my numberOfCoefficients;
- FunctionTerms_extendCapacityIf (me, n + 1);
-
- my coefficients [n + 1] = my coefficients [n];
- for (integer j = n; j >= 2; j --)
- my coefficients [j] = my coefficients [j - 1] - my coefficients [j] * factor;
- my coefficients [1] *= -factor;
- my numberOfCoefficients += 1;
- }
- /* P(x) * (x^2 - a)
- * Postcondition: my numberOfCoefficients = old_numberOfCoefficients + 2
- */
- void Polynomial_multiply_secondOrderFactor (Polynomial me, double factor) {
- integer n = my numberOfCoefficients;
- FunctionTerms_extendCapacityIf (me, n + 2);
- my coefficients [n + 2] = my coefficients [n];
- my coefficients [n + 1] = my coefficients [n - 1];
- for (integer j = n; j >= 3; j --)
- my coefficients [j] = my coefficients [j - 2] - factor * my coefficients [j];
- my coefficients [2] *= - factor;
- my coefficients [1] *= - factor;
- my numberOfCoefficients += 2;
- }
- autoPolynomial Polynomials_multiply (Polynomial me, Polynomial thee) {
- try {
- integer n1 = my numberOfCoefficients, n2 = thy numberOfCoefficients;
- Melder_require (my xmax > thy xmin && my xmin < thy xmax, U"Domains should overlap.");
-
- double xmin = my xmin > thy xmin ? my xmin : thy xmin;
- double xmax = my xmax < thy xmax ? my xmax : thy xmax;
- autoPolynomial him = Polynomial_create (xmin, xmax, n1 + n2 - 2);
- for (integer i = 1; i <= n1; i ++) {
- for (integer j = 1; j <= n2; j ++) {
- integer k = i + j - 1;
- his coefficients [k] += my coefficients [i] * thy coefficients [j];
- }
- }
- return him;
- } catch (MelderError) {
- Melder_throw (U"Polynomials not multiplied.");
- }
- }
- void Polynomials_divide (Polynomial me, Polynomial thee, autoPolynomial *q, autoPolynomial *r) {
- integer degree, m = my numberOfCoefficients, n = thy numberOfCoefficients;
- if (! q && ! r)
- return;
- autoNUMvector<double> qc (1, m);
- autoNUMvector<double> rc (1, m);
- autoPolynomial aq, ar;
- polynomial_divide (my coefficients, m, thy coefficients, n, qc.peek(), rc.peek());
- if (q) {
- degree = std::max (m - n, (integer) 0);
- aq = Polynomial_create (my xmin, my xmax, degree);
- if (degree >= 0)
- NUMvector_copyElements (qc.peek(), aq -> coefficients, 1, degree + 1);
- *q = aq.move();
- }
- if (r) {
- degree = n - 2;
- if (m >= n)
- degree --;
- if (degree < 0)
- degree = 0;
- while (degree > 1 && rc [degree] == 0.0)
- degree --;
- ar = Polynomial_create (my xmin, my xmax, degree);
- NUMvector_copyElements (rc.peek(), ar -> coefficients, 1, degree + 1);
- *r = ar.move();
- }
- }
- /******** LegendreSeries ********************************************/
- Thing_implement (LegendreSeries, FunctionTerms, 0);
- double structLegendreSeries :: v_evaluate (double x) {
- double p = our coefficients [1];
- // Transform x from domain [xmin, xmax] to domain [-1, 1]
- if (x < our xmin || x > our xmax)
- return undefined;
- double pim1 = x = (2 * x - our xmin - our xmax) / (our xmax - our xmin);
- if (numberOfCoefficients > 1) {
- double pim2 = 1, twox = 2.0 * x, f2 = x, d = 1.0;
- p += our coefficients [2] * pim1;
- for (integer i = 3; i <= our numberOfCoefficients; i ++) {
- double f1 = d ++;
- f2 += twox;
- double pi = (f2 * pim1 - f1 * pim2) / d;
- p += our coefficients [i] * pi;
- pim2 = pim1; pim1 = pi;
- }
- }
- return p;
- }
- void structLegendreSeries :: v_evaluateTerms (double x, double terms []) {
- if (x < our xmin || x > our xmax) {
- for (integer i = 1; i <= numberOfCoefficients; i ++)
- terms [i] = undefined;
- return;
- }
- // Transform x from domain [xmin, xmax] to domain [-1, 1]
- x = (2.0 * x - our xmin - our xmax) / (our xmax - our xmin);
- terms [1] = 1;
- if (our numberOfCoefficients > 1) {
- double twox = 2.0 * x, f2 = x, d = 1.0;
- terms [2] = x;
- for (integer i = 3; i <= numberOfCoefficients; i ++) {
- double f1 = d ++;
- f2 += twox;
- terms [i] = (f2 * terms [i - 1] - f1 * terms [i - 2]) / d;
- }
- }
- }
- void structLegendreSeries :: v_getExtrema (double x1, double x2, double *p_xmin, double *p_ymin, double *p_xmax, double *p_ymax) {
- try {
- autoPolynomial p = LegendreSeries_to_Polynomial (this);
- FunctionTerms_getExtrema (p.get(), x1, x2, p_xmin, p_ymin, p_xmax, p_ymax);
- } catch (MelderError) {
- structFunctionTerms :: v_getExtrema (x1, x2, p_xmin, p_ymin, p_xmax, p_ymax);
- Melder_clearError ();
- }
- }
- autoLegendreSeries LegendreSeries_create (double xmin, double xmax, integer numberOfPolynomials) {
- try {
- autoLegendreSeries me = Thing_new (LegendreSeries);
- FunctionTerms_init (me.get(), xmin, xmax, numberOfPolynomials);
- return me;
- } catch (MelderError) {
- Melder_throw (U"LegendreSeries not created.");
- }
- }
- autoLegendreSeries LegendreSeries_createFromString (double xmin, double xmax, conststring32 s) {
- try {
- autoLegendreSeries me = Thing_new (LegendreSeries);
- FunctionTerms_initFromString (me.get(), xmin, xmax, s, 0);
- return me;
- } catch (MelderError) {
- Melder_throw (U"LegendreSeries not created from string.");
- }
- }
- autoLegendreSeries LegendreSeries_getDerivative (LegendreSeries me) {
- try {
- autoLegendreSeries thee = LegendreSeries_create (my xmin, my xmax, my numberOfCoefficients - 1);
- for (integer n = 1; n <= my numberOfCoefficients - 1; n ++) {
- // P [n]'(x) = Sum (k=0..nonNegative, (2n - 4k - 1) P [n-2k-1](x))
- integer n2 = n - 1;
- for (integer k = 0; n2 >= 0; k ++, n2 -= 2)
- thy coefficients [n2 + 1] += (2 * n - 4 * k - 1) * my coefficients [n + 1];
- }
- return thee;
- } catch (MelderError) {
- Melder_throw (me, U": no derivative created.");
- }
- }
- autoPolynomial LegendreSeries_to_Polynomial (LegendreSeries me) {
- try {
- double xmin = -1, xmax = 1;
- autoPolynomial thee = Polynomial_create (xmin, xmax, my numberOfCoefficients - 1);
- thy coefficients [1] = my coefficients [1]; /* * p [1] */
- if (my numberOfCoefficients == 1)
- return thee;
- thy coefficients [2] = my coefficients [2]; /* * p [2] */
- if (my numberOfCoefficients > 2) {
- autoNUMvector<double> buf (1, 3 * my numberOfCoefficients);
- double *pn = buf.peek();
- double *pnm1 = pn + my numberOfCoefficients;
- double *pnm2 = pnm1 + my numberOfCoefficients;
- // Start the recursion: P [1] = x; P [0] = 1;
- pnm1 [2] = 1;
- pnm2 [1] = 1;
- for (integer n = 2; n <= my numberOfCoefficients - 1; n ++) {
- double a = (2.0 * n - 1.0) / (double) n;
- double c = - (n - 1.0) / (double) n;
- double *t1 = pnm1, *t2 = pnm2;
- NUMpolynomial_recurrence (pn, n, a, 0, c, pnm1, pnm2);
- if (my coefficients [n + 1] != 0.0) {
- for (integer j = 1; j <= n + 1; j ++)
- thy coefficients [j] += my coefficients [n + 1] * pn [j];
- }
- pnm1 = pn; pnm2 = t1; pn = t2;
- }
- }
- if (my xmin != xmin || my xmax != xmax)
- thee = Polynomial_scaleX (thee.get(), my xmin, my xmax);
- return thee;
- } catch (MelderError) {
- Melder_throw (me, U": not converted to Polynomial.");
- }
- }
- /********* Roots ****************************************************/
- Thing_implement (Roots, ComplexVector, 0);
- autoRoots Roots_create (integer numberOfRoots) {
- try {
- autoRoots me = Thing_new (Roots);
- ComplexVector_init (me.get(), 1, numberOfRoots);
- return me;
- } catch (MelderError) {
- Melder_throw (U"Roots not created.");
- }
- }
- void Roots_fixIntoUnitCircle (Roots me) {
- dcomplex z10 { 1.0, 0.0 };
- for (integer i = my min; i <= my max; i ++) {
- if (dcomplex_abs (my v [i]) > 1.0)
- my v [i] = dcomplex_div (z10, dcomplex_conjugate (my v [i]));
- }
- }
- static void NUMdcvector_extrema_re (dcomplex v [], integer lo, integer hi, double *min, double *max) {
- *min = *max = v [lo].re;
- for (integer i = lo + 1; i <= hi; i ++) {
- if (v [i].re < *min) {
- *min = v [i].re;
- } else if (v [i].re > *max) {
- *max = v [i].re;
- }
- }
- }
- static void NUMdcvector_extrema_im (dcomplex v [], integer lo, integer hi, double *min, double *max) {
- *min = *max = v [lo].im;
- for (integer i = lo + 1; i <= hi; i ++) {
- if (v [i]. im < *min) {
- *min = v [i]. im;
- } else if (v [i]. im > *max) {
- *max = v [i]. im;
- }
- }
- }
- void Roots_draw (Roots me, Graphics g, double rmin, double rmax, double imin, double imax,
- conststring32 symbol, int fontSize, bool garnish)
- {
- int oldFontSize = Graphics_inqFontSize (g);
- double eps = 1e-6;
- if (rmax <= rmin) {
- NUMdcvector_extrema_re (my v, 1, my max, & rmin, & rmax);
- }
- double denominator = fabs (rmax) > fabs (rmin) ? fabs (rmax) : fabs (rmin);
- if (denominator == 0.0)
- denominator = 1.0;
- if (fabs ((rmax - rmin) / denominator) < eps) {
- rmin -= 1.0;
- rmax += 1.0;
- }
- if (imax <= imin) {
- NUMdcvector_extrema_im (my v, 1, my max, & imin, & imax);
- }
- denominator = fabs (imax) > fabs (imin) ? fabs (imax) : fabs (imin);
- if (denominator == 0.0)
- denominator = 1.0;
- if (fabs ((imax - imin) / denominator) < eps) {
- imin -= 1;
- imax += 1;
- }
- Graphics_setInner (g);
- Graphics_setWindow (g, rmin, rmax, imin, imax);
- Graphics_setFontSize (g, fontSize);
- Graphics_setTextAlignment (g, Graphics_CENTRE, Graphics_HALF);
- for (integer i = 1; i <= my max; i ++) {
- double re = my v [i].re, im = my v [i].im;
- if (re >= rmin && re <= rmax && im >= imin && im <= imax)
- Graphics_text (g, re, im, symbol);
- }
- Graphics_setFontSize (g, oldFontSize);
- Graphics_unsetInner (g);
- if (garnish) {
- Graphics_drawInnerBox (g);
- if (rmin * rmax < 0.0)
- Graphics_markLeft (g, 0.0, true, true, true, U"0");
- if (imin * imax < 0.0)
- Graphics_markBottom (g, 0.0, true, true, true, U"0");
- Graphics_marksLeft (g, 2, true, true, false);
- Graphics_textLeft (g, true, U"Imaginary part");
- Graphics_marksBottom (g, 2, true, true, false);
- Graphics_textBottom (g, true, U"Real part");
- }
- }
- autoRoots Polynomial_to_Roots (Polynomial me) {
- try {
- integer np1 = my numberOfCoefficients, n = np1 - 1, n2 = n * n;
- Melder_require (n > 0, U"Cannot find roots of a constant function.");
-
- // Allocate storage for Hessenberg matrix (n * n) plus real and imaginary
- // parts of eigenvalues wr [1..n] and wi [1..n].
- autoNUMvector <double> hes (1, n2 + n + n);
- double *wr = & hes [n2];
- double *wi = & hes [n2 + n];
- // Fill the upper Hessenberg matrix (storage is Fortran)
- // C: [i] [j] -> Fortran: (j-1)*n + i
- for (integer i = 1; i <= n; i ++) {
- hes [(i - 1) * n + 1] = - (my coefficients [np1 - i] / my coefficients [np1]);
- if (i < n) {
- hes [(i - 1) * n + 1 + i] = 1;
- }
- }
- // Find out the working storage needed
- char job = 'E', compz = 'N';
- integer ilo = 1, ihi = n, ldh = n, ldz = n, lwork = -1, info;
- double *z = 0, wt [1];
- NUMlapack_dhseqr (&job, &compz, &n, &ilo, &ihi, &hes [1], &ldh, &wr [1], &wi [1], z, &ldz, wt, &lwork, &info);
- if (info != 0) {
- if (info < 0)
- Melder_throw (U"Programming error. Argument ", info, U" in NUMlapack_dhseqr has illegal value.");
- }
- lwork = Melder_ifloor (wt [0]);
- autoNUMvector <double> work (1, lwork);
- // Find eigenvalues.
- NUMlapack_dhseqr (&job, &compz, &n, &ilo, &ihi, &hes [1], &ldh, &wr [1], &wi [1], z, &ldz, &work [1], &lwork, &info);
- integer nrootsfound = n;
- integer ioffset = 0;
- if (info > 0) {
- // if INFO = i, NUMlapack_dhseqr failed to compute all of the eigenvalues. Elements i+1:n of
- // WR and WI contain those eigenvalues which have been successfully computed
- nrootsfound -= info;
- Melder_require (nrootsfound > 0, U"No roots found.");
- Melder_warning (U"Calculated only ", nrootsfound, U" roots.");
- ioffset = info;
- } else if (info < 0) {
- Melder_throw (U"Programming error. Argument ", info, U" in NUMlapack_dhseqr has illegal value.");
- }
- autoRoots thee = Roots_create (nrootsfound);
- for (integer i = 1; i <= nrootsfound; i ++) {
- (thy v [i]).re = wr [ioffset + i];
- (thy v [i]).im = wi [ioffset + i];
- }
- Roots_Polynomial_polish (thee.get(), me);
- return thee;
- } catch (MelderError) {
- Melder_throw (me, U": no roots can be calculated.");
- }
- }
- void Roots_sort (Roots me) {
- (void) me;
- }
- // Precondition: complex roots occur in pairs (a,bi), (a,-bi) with b>0
- void Roots_Polynomial_polish (Roots me, Polynomial thee) {
- integer i = my min, maxit = 80;
- while (i <= my max) {
- double im = my v [i].im, re = my v [i].re;
- if (im != 0.0) {
- Polynomial_polish_complexroot_nr (thee, & my v [i], maxit);
- if (i < my max && im == -my v [i + 1].im && re == my v [i + 1].re) {
- my v [i + 1].re = my v [i].re; my v [i + 1].im = -my v [i].im;
- i ++;
- }
- } else {
- Polynomial_polish_realroot (thee, & (my v [i].re), maxit);
- }
- i ++;
- }
- }
- autoPolynomial Roots_to_Polynomial (Roots me, bool rootsAreReal) {
- try {
- (void) me;
- autoPolynomial thee;
- if (! rootsAreReal)
- throw MelderError();
- return thee;
- } catch (MelderError) {
- Melder_throw (U"Not implemented yet");
- }
- }
- static void dpoly_nr (double x, double *f, double *df, void *closure) {
- Polynomial_evaluateWithDerivative ((Polynomial) closure, x, f, df);
- }
- double Polynomial_findOneSimpleRealRoot_nr (Polynomial me, double xmin, double xmax) {
- double root;
- NUMnrbis (dpoly_nr, xmin, xmax, me, &root);
- return root;
- }
- static double dpoly_r (double x, void *closure) {
- return Polynomial_evaluate ((Polynomial) closure, x);
- }
- double Polynomial_findOneSimpleRealRoot_ridders (Polynomial me, double xmin, double xmax) {
- return NUMridders (dpoly_r, xmin, xmax, me);
- }
- void Polynomial_divide_firstOrderFactor (Polynomial me, double factor, double *p_remainder) { // P(x)/(x-a)
- double remainder = undefined;
- if (my numberOfCoefficients > 1) {
- remainder = my coefficients [my numberOfCoefficients];
- for (integer j = my numberOfCoefficients - 1; j > 0; j --) {
- double tmp = my coefficients [j];
- my coefficients [j] = remainder;
- remainder = tmp + remainder * factor;
- }
- my numberOfCoefficients --;
- } else {
- my coefficients [1] = 0.0;
- }
- if (p_remainder)
- *p_remainder = remainder;
- }
- void Polynomial_divide_secondOrderFactor (Polynomial me, double factor) {
- if (my numberOfCoefficients > 2) {
- integer n = my numberOfCoefficients;
- /* c [1]+c [2]*x...c [n+1]*x^n / (x^2 - a) = r [1]+r [2]*x+...r [n-1]x^(n-2) + possible remainder a [1]+a [2]*x)
- * r [j] = c [j+2]+factor*r [j+2] */
- double cjp2 = my coefficients [n];
- double cjp1 = my coefficients [n - 1];
- my coefficients [n] = my coefficients [n - 1] = 0.0;
- for (integer j = n - 2; j > 0; j --) {
- double cj = my coefficients [j];
- my coefficients [j] = cjp2 + factor * my coefficients [j + 2];
- cjp2 = cjp1;
- cjp1 = cj;
- }
- my numberOfCoefficients -= 2;
- } else {
- my numberOfCoefficients = 1;
- my coefficients [1] = 0.0;
- }
- }
- void Roots_setRoot (Roots me, integer index, double re, double im) {
- Melder_require (index >= my min && index <= my max, U"Index should be in interval [1, ", my max, U"].");
- my v [index].re = re;
- my v [index].im = im;
- }
- dcomplex Roots_evaluate_z (Roots me, dcomplex z) {
- dcomplex result = {1, 0};
- for (integer i = my min; i <= my max; i ++) {
- dcomplex t = dcomplex_sub (z, my v [i]);
- result = dcomplex_mul (result, t);
- }
- return result;
- }
- autoSpectrum Roots_to_Spectrum (Roots me, double nyquistFrequency, integer numberOfFrequencies, double radius) {
- try {
- Melder_require (numberOfFrequencies > 1,
- U"Number of frequencies should be greater than 1.");
- autoSpectrum thee = Spectrum_create (nyquistFrequency, numberOfFrequencies);
- double phi = NUMpi / (numberOfFrequencies - 1);
- dcomplex z;
- for (integer i = 1; i <= numberOfFrequencies; i ++) {
- z.re = radius * cos ( (i - 1) * phi);
- z.im = radius * sin ( (i - 1) * phi);
- dcomplex s = Roots_evaluate_z (me, z);
- thy z [1] [i] = s.re; thy z [2] [i] = s.im;
- }
- return thee;
- } catch (MelderError) {
- Melder_throw (me, U": no Spectrum calculated.");
- }
- }
- integer Roots_getNumberOfRoots (Roots me) {
- return my max;
- }
- dcomplex Roots_getRoot (Roots me, integer index) {
- Melder_require (index > 0 && index <= my max,
- U"Root index out of range.");
- return my v [index];
- }
- /* Can be speeded up by doing a FFT */
- autoSpectrum Polynomial_to_Spectrum (Polynomial me, double nyquistFrequency, integer numberOfFrequencies, double radius) {
- try {
- Melder_require (numberOfFrequencies > 1,
- U"Number of frequencies should be greater than 1.");
- autoSpectrum thee = Spectrum_create (nyquistFrequency, numberOfFrequencies);
- double phi = NUMpi / (numberOfFrequencies - 1);
- for (integer i = 1; i <= numberOfFrequencies; i ++) {
- double re, im;
- Polynomial_evaluate_z_cart (me, radius, (i - 1) * phi, &re, &im);
- thy z [1] [i] = re;
- thy z [2] [i] = im;
- }
- return thee;
- } catch (MelderError) {
- Melder_throw (me, U": no Spectrum calculated.");
- }
- }
- /****** ChebyshevSeries ******************************************/
- Thing_implement (ChebyshevSeries, FunctionTerms, 0);
- /*
- p(x) = sum (k=1..numberOfCoefficients, c [k]*T [k](x')) - c [1] / 2;
- Numerical evaluation via Clenshaw's recurrence equation (NRC: 5.8.11)
- d [m+1] = d [m] = 0;
- d [j] = 2 x' d [j+1] - d [j+2] + c [j];
- p(x) = d [0] = x' d [1] - d [2] + c [0] / 2;
- x' = (2 * x - xmin - xmax) / (xmax - xmin)
- */
- double structChebyshevSeries :: v_evaluate (double x) {
- if (x < our xmin || x > our xmax)
- return undefined;
- double d1 = 0.0, d2 = 0.0;
- if (numberOfCoefficients > 1) {
- // Transform x from domain [xmin, xmax] to domain [-1, 1]
- x = (2.0 * x - our xmin - our xmax) / (our xmax - our xmin);
- double x2 = 2.0 * x;
- for (integer i = our numberOfCoefficients; i > 1; i --) {
- double tmp = d1;
- d1 = x2 * d1 - d2 + our coefficients [i];
- d2 = tmp;
- }
- }
- return x * d1 - d2 + our coefficients [1];
- }
- /*
- T [n](x) = 2*x*T [n-1] - T [n-2](x) n >= 2
- */
- void structChebyshevSeries :: v_evaluateTerms (double x, double *terms) {
- if (x < our xmin || x > our xmax) {
- for (integer i = 1; i <= our numberOfCoefficients; i ++)
- terms [i] = undefined;
- return;
- }
- terms [1] = 1.0;
- if (numberOfCoefficients > 1) {
- // Transform x from domain [xmin, xmax] to domain [-1, 1]
- terms [2] = x = (2.0 * x - xmin - xmax) / (xmax - xmin);
- for (integer i = 3; i <= numberOfCoefficients; i ++)
- terms [i] = 2.0 * x * terms [i - 1] - terms [i - 2];
- }
- }
- void structChebyshevSeries :: v_getExtrema (double x1, double x2, double *p_xmin, double *p_ymin, double *p_xmax, double *p_ymax) {
- try {
- autoPolynomial p = ChebyshevSeries_to_Polynomial (this);
- FunctionTerms_getExtrema (p.get(), x1, x2, p_xmin, p_ymin, p_xmax, p_ymax);
- } catch (MelderError) {
- Melder_throw (this, U"Extrema cannot be calculated");
- }
- }
- autoChebyshevSeries ChebyshevSeries_create (double lxmin, double lxmax, integer numberOfPolynomials) {
- try {
- autoChebyshevSeries me = Thing_new (ChebyshevSeries);
- FunctionTerms_init (me.get(), lxmin, lxmax, numberOfPolynomials);
- return me;
- } catch (MelderError) {
- Melder_throw (U"ChebyshevSeries not created.");
- }
- }
- autoChebyshevSeries ChebyshevSeries_createFromString (double lxmin, double lxmax, conststring32 s) {
- try {
- autoChebyshevSeries me = Thing_new (ChebyshevSeries);
- FunctionTerms_initFromString (me.get(), lxmin, lxmax, s, 0);
- return me;
- } catch (MelderError) {
- Melder_throw (U"ChebyshevSeries not created from string.");
- };
- }
- autoPolynomial ChebyshevSeries_to_Polynomial (ChebyshevSeries me) {
- try {
- double xmin = -1, xmax = 1;
- autoPolynomial thee = Polynomial_create (xmin, xmax, my numberOfCoefficients - 1);
- thy coefficients [1] = my coefficients [1] /* * p [1] */;
- if (my numberOfCoefficients == 1)
- return thee;
- thy coefficients [2] = my coefficients [2];
- if (my numberOfCoefficients > 2) {
- autoNUMvector<double> buf (1, 3 * my numberOfCoefficients);
- double *pn = buf.peek();
- double *pnm1 = pn + my numberOfCoefficients;
- double *pnm2 = pnm1 + my numberOfCoefficients;
- // Start the recursion: T [1] = x; T [0] = 1;
- pnm1 [2] = 1.0;
- pnm2 [1] = 1.0;
- double a = 2.0, b = 0.0, c = -1.0;
- for (integer n = 2; n <= my numberOfCoefficients - 1; n ++) {
- double *t1 = pnm1, *t2 = pnm2;
- NUMpolynomial_recurrence (pn, n, a, b, c, pnm1, pnm2);
- if (my coefficients [n + 1] != 0) {
- for (integer j = 1; j <= n + 1; j ++)
- thy coefficients [j] += my coefficients [n + 1] * pn [j];
- }
- pnm1 = pn;
- pnm2 = t1;
- pn = t2;
- }
- }
- if (my xmin != xmin || my xmax != xmax)
- thee = Polynomial_scaleX (thee.get(), my xmin, my xmax);
- return thee;
- } catch (MelderError) {
- Melder_throw (me, U"; not converted to Polynomial.");
- };
- }
- void FunctionTerms_RealTier_fit (FunctionTerms me, RealTier thee, int freeze [], double tol, int ic, autoCovariance *c) {
- try {
- integer numberOfData = thy points.size;
- integer numberOfParameters = my numberOfCoefficients;
- integer numberOfFreeParameters = numberOfParameters;
- Melder_require (numberOfData > 1, U"The number of data point should be larger than 1.");
- autoFunctionTerms frozen = Data_copy (me);
- autoVEC terms = VECzero (my numberOfCoefficients);
- autoVEC p = VECzero (numberOfParameters);
- autoVEC y_residual = VECzero (numberOfData);
- autoCovariance ac;
- if (ic)
- ac = (Covariance_create (numberOfParameters));
- integer k = 1;
- for (integer j = 1; j <= my numberOfCoefficients; j ++) {
- if (freeze && freeze [j]) {
- numberOfFreeParameters--;
- } else {
- p [k] = my coefficients [j]; k ++;
- frozen -> coefficients [j] = 0;
- }
- }
-
- Melder_require (numberOfFreeParameters > 0, U"No free parameters left.");
- autoSVD svd = SVD_create (numberOfData, numberOfFreeParameters);
- double sigma = RealTier_getStandardDeviation_points (thee, my xmin, my xmax);
-
- Melder_require (isdefined (sigma), U"Not enough data points in fit interval.");
- for (integer i = 1; i <= numberOfData; i ++) {
- // Only 'residual variance' must be explained by the model
- // Evaluate only with the frozen parameters
- RealPoint point = thy points.at [i];
- double x = point -> number;
- double y = point -> value;
- double y_frozen = ( numberOfFreeParameters == numberOfParameters ? 0.0 :
- FunctionTerms_evaluate (frozen.get(), x));
- y_residual [i] = (y - y_frozen) / sigma;
- // Data matrix
- FunctionTerms_evaluateTerms (me, x, terms.at);
- k = 0;
- for (integer j = 1; j <= my numberOfCoefficients; j ++) {
- if (! freeze || ! freeze [j]) {
- k ++;
- svd -> u [i] [k] = terms [j] / sigma;
- }
- }
- }
- // SVD and evaluation of the singular values
- if (tol > 0.0)
- SVD_setTolerance (svd.get(), tol);
- SVD_compute (svd.get());
- autoVEC result = SVD_solve (svd.get(), y_residual.get());
- // Put fitted values at correct position
- k = 1;
- for (integer j = 1; j <= my numberOfCoefficients; j ++) {
- if (! freeze || ! freeze [j])
- my coefficients [j] = result [k ++];
- }
- if (ic)
- svdcvm (svd -> v.at, numberOfFreeParameters, numberOfParameters, freeze, svd -> d.at, ac -> data.at);
- if (c) *c = ac.move();
- } catch (MelderError) {
- Melder_throw (me, U" & ", thee, U": no fit.");
- }
- }
- autoPolynomial RealTier_to_Polynomial (RealTier me, integer degree, double tol, int ic, autoCovariance *cvm) {
- try {
- autoPolynomial thee = Polynomial_create (my xmin, my xmax, degree);
- FunctionTerms_RealTier_fit (thee.get(), me, 0, tol, ic, cvm);
- return thee;
- } catch (MelderError) {
- Melder_throw (me, U": no Polynomial fitted.");
- }
- }
- autoLegendreSeries RealTier_to_LegendreSeries (RealTier me, integer degree, double tol, int ic, autoCovariance *cvm) {
- try {
- autoLegendreSeries thee = LegendreSeries_create (my xmin, my xmax, degree);
- FunctionTerms_RealTier_fit (thee.get(), me, 0, tol, ic, cvm);
- return thee;
- } catch (MelderError) {
- Melder_throw (me, U": no LegendreSeries fitted.");
- }
- }
- autoChebyshevSeries RealTier_to_ChebyshevSeries (RealTier me, integer degree, double tol, int ic, autoCovariance *cvm) {
- try {
- autoChebyshevSeries thee = ChebyshevSeries_create (my xmin, my xmax, degree);
- FunctionTerms_RealTier_fit (thee.get(), me, 0, tol, ic, cvm);
- return thee;
- } catch (MelderError) {
- Melder_throw (me, U":no ChebyshevSeries fitted.");
- };
- }
- /******* Splines *************************************************/
- /*
- Functions for calculating an mspline and an ispline. These functions should replace
- the functions in NUM2.c. Before we can do that we first have to adapt the spline-
- dependencies in MDS.c.
- Formally nKnots == order + numberOfInteriorKnots + order.
- We forget about the multiple knots at start and end.
- Our point-sequece xmin < interiorKont [1] < ... < interiorKnot [n] < xmax.
- nKnots is now numberOfinteriorKnots + 2.
- */
- static double NUMmspline2 (constVEC points, integer order, integer index, double x) {
- integer numberOfSplines = points.size + order - 2;
- double m [Spline_MAXIMUM_DEGREE + 2];
- Melder_assert (points.size > 2 && order > 0 && index > 0);
- if (index > numberOfSplines)
- return undefined;
- /*
- Find the range/interval where x is located.
- M-splines of order k have degree k-1.
- M-splines are zero outside interval [knot [i], knot [i+order]).
- First and last 'order' knots are equal, i.e.,
- knot [1] = ... = knot [order] && knot [nKnots-order+1] = ... knot [nKnots].
- */
- integer index_b = index - order + 1;
- index_b = std::max (index_b, (integer) 1);
- if (x < points [index_b])
- return 0;
- integer index_e = index_b + std::min (index, order);
- index_e = std::min (points.size, index_e);
- if (x > points [index_e])
- return 0;
- // Calculate M [i](x|1,t) according to eq.2.
- for (integer k = 1; k <= order; k ++) {
- integer k1 = index - order + k, k2 = k1 + 1;
- m [k] = 0;
- if (k1 > 0 && k2 <= points.size && x >= points [k1] && x < points [k2])
- m [k] = 1 / (points [k2] - points [k1]);
- }
- // Iterate to get M [i](x|k,t)
- for (integer k = 2; k <= order; k ++) {
- for (integer j = 1; j <= order - k + 1; j ++) {
- integer k1 = index - order + j, k2 = k1 + k;
- if (k2 > 1 && k1 < 1) {
- k1 = 1;
- } else if (k2 > points.size && k1 < points.size) {
- k2 = points.size;
- }
- if (k1 > 0 && k2 <= points.size) {
- double p1 = points [k1], p2 = points [k2];
- m [j] = k * ((x - p1) * m [j] + (p2 - x) * m [j + 1]) /
- ((k - 1) * (p2 - p1));
- }
- }
- }
- return m [1];
- }
- static double NUMispline2 (constVEC points, integer order, integer index, double x) {
- Melder_assert (points.size > 2 && order > 0 && index > 0);
- integer index_b = index - order + 1;
- index_b = std::max (index_b, (integer) 1);
- if (x < points [index_b])
- return 0.0;
- integer index_e = index_b + std::min (index, order);
- index_e = std::min (points.size, index_e);
- if (x > points [index_e])
- return 1.0;
- integer j;
- for (j = index_e - 1; j >= index_b; j--) {
- if (x > points [j]) {
- break;
- }
- }
- // Equation 5 in Ramsay's article contains some errors!!!
- // 1. the interval selection must be 'j-k <= i <= j' instead of
- // 'j-k+1 <= i <= j'
- // 2. the summation index m starts at 'i+1' instead of 'i'
- double y = 0.0;
- for (integer m = index + 1; m <= j + order; m ++) {
- integer km = m - order, kmp = km + order + 1;
- km = std::max (km, (integer) 1);
- kmp = std::min (kmp, points.size);
- y += (points [kmp] - points [km]) * NUMmspline2 (points, order + 1, m, x);
- }
- return y /= (order + 1);
- }
- Thing_implement (Spline, FunctionTerms, 0);
- double structSpline :: v_evaluate (double /* x */) {
- return 0.0;
- }
- integer structSpline :: v_getDegree () {
- return degree;
- }
- integer structSpline :: v_getOrder () {
- return degree + 1;
- }
- /* Precondition: FunctionTerms part inited + degree */
- static void Spline_initKnotsFromString (Spline me, integer degree, conststring32 interiorKnots_string) {
- Melder_require (degree <= Spline_MAXIMUM_DEGREE, U"Degree should be <= ", Spline_MAXIMUM_DEGREE, U".");
-
- autoVEC interiorKnots = VEC_createFromString (interiorKnots_string);
- VECsort_inplace (interiorKnots.get());
- if (interiorKnots [1] <= my xmin || interiorKnots [interiorKnots.size] > my xmax)
- Melder_throw (U"Knots should be inside domain.");
- my degree = degree;
- integer order = Spline_getOrder (me); /* depends on spline type !! */
- integer n = interiorKnots.size + order;
- Melder_require (my numberOfCoefficients == n, U"Number of coefficients should equal ", n, U".");
- my numberOfKnots = interiorKnots.size + 2;
- my knots = VECzero (my numberOfKnots);
- for (integer i = 1; i <= interiorKnots.size; i ++) {
- my knots [i + 1] = interiorKnots [i];
- }
- my knots [1] = my xmin;
- my knots [my numberOfKnots] = my xmax;
- }
- void Spline_init (Spline me, double xmin, double xmax, integer degree, integer numberOfCoefficients, integer numberOfKnots) {
- Melder_require (degree <= Spline_MAXIMUM_DEGREE, U"Degree should be <= ", Spline_MAXIMUM_DEGREE, U".");
-
- FunctionTerms_init (me, xmin, xmax, numberOfCoefficients);
- my knots = VECzero (numberOfKnots);
- my degree = degree;
- my numberOfKnots = numberOfKnots;
- my knots [1] = xmin;
- my knots [numberOfKnots] = xmax;
- }
- void Spline_drawKnots (Spline me, Graphics g, double xmin, double xmax, double ymin, double ymax, int garnish) {
- integer order = Spline_getOrder (me);
- if (xmax <= xmin) {
- xmin = my xmin;
- xmax = my xmax;
- }
- if (xmax < my xmin || xmin > my xmax)
- return;
- if (ymax <= ymin) {
- double x1, x2;
- FunctionTerms_getExtrema (me, xmin, xmax, &x1, &ymin, &x2, &ymax);
- }
- Graphics_setWindow (g, xmin, xmax, ymin, ymax);
- if (my knots [1] >= xmin && my knots [1] <= xmax) {
- Graphics_markTop (g, my knots [1], false, true, true,
- ! garnish ? U"" :
- order == 1 ? U"t__1_" :
- order == 2 ? U"{t__1_, t__2_}" :
- Melder_cat (U"{t__1_..t__", order, U"_}")
- );
- }
- for (integer i = 2; i <= my numberOfKnots - 1; i ++) {
- if (my knots [i] >= xmin && my knots [i] <= xmax) {
- Graphics_markTop (g, my knots [i], false, true, true,
- ! garnish ? U"" :
- Melder_cat (U"t__", i + order - 1, U"_")
- );
- }
- }
- if (my knots [my numberOfKnots] >= xmin && my knots [my numberOfKnots] <= xmax) {
- integer numberOfKnots = ! garnish ? 0 : my numberOfKnots + 2 * (order - 1);
- Graphics_markTop (g, my knots [my numberOfKnots], false, true, true,
- ! garnish ? U"" :
- order == 1 ? Melder_cat (U"t__", numberOfKnots, U"_") :
- order == 2 ? Melder_cat (U"{t__", numberOfKnots - 1, U"_, t__", numberOfKnots, U"_}") :
- Melder_cat (U"{t__", numberOfKnots - order + 1, U"_..t__", numberOfKnots, U"_}")
- );
- }
- }
- integer Spline_getOrder (Spline me) {
- return my v_getOrder ();
- }
- autoSpline Spline_scaleX (Spline me, double xmin, double xmax) {
- try {
- Melder_assert (xmin < xmax);
- autoSpline thee = Data_copy (me);
- thy xmin = xmin;
- thy xmax = xmax;
- // x = a x + b
- // Constraints:
- // my xmin = a xmin + b; a = (my xmin - my xmax) / (xmin - xmax);
- // my xmax = a xmax + b; b = my xmin - a * xmin
- double a = (xmin - xmax) / (my xmin - my xmax);
- double b = xmin - a * my xmin;
- for (integer i = 1; i <= my numberOfKnots; i ++) {
- thy knots [i] = a * my knots [i] + b;
- }
- return thee;
- } catch (MelderError) {
- Melder_throw (U"Scaled Spline not created.");
- }
- }
- /********* MSplines ************************************************/
- double structMSpline :: v_evaluate (double x) {
- if (x < our xmin || x > our xmax)
- return 0.0;
- double result = 0.0;
- for (integer i = 1; i <= numberOfCoefficients; i ++) {
- if (coefficients [i] != 0.0)
- result += coefficients [i] * NUMmspline2 (knots.get(), degree + 1, i, x);
- }
- return result;
- }
- void structMSpline :: v_evaluateTerms (double x, double *terms) {
- if (x < our xmin || x > our xmax)
- return;
- for (integer i = 1; i <= numberOfCoefficients; i ++)
- terms [i] = NUMmspline2 (knots.get(), degree + 1, i, x);
- }
- Thing_implement (MSpline, Spline, 0);
- autoMSpline MSpline_create (double xmin, double xmax, integer degree, integer numberOfInteriorKnots) {
- try {
- autoMSpline me = Thing_new (MSpline);
- integer numberOfCoefficients = numberOfInteriorKnots + degree + 1;
- integer numberOfKnots = numberOfCoefficients + degree + 1;
- Spline_init (me.get(), xmin, xmax, degree, numberOfCoefficients, numberOfKnots);
- return me;
- } catch (MelderError) {
- Melder_throw (U"MSpline not created.");
- }
- }
- autoMSpline MSpline_createFromStrings (double xmin, double xmax, integer degree, conststring32 coef, conststring32 interiorKnots) {
- try {
- Melder_require (degree <= Spline_MAXIMUM_DEGREE,
- U"Degree should be <= ", Spline_MAXIMUM_DEGREE, U".");
- autoMSpline me = Thing_new (MSpline);
- FunctionTerms_initFromString (me.get(), xmin, xmax, coef, 1);
- Spline_initKnotsFromString (me.get(), degree, interiorKnots);
- return me;
- } catch (MelderError) {
- Melder_throw (U"MSpline not created from strings.");
- }
- }
- /******** ISplines ************************************************/
- double structISpline :: v_evaluate (double x) {
- if (x < our xmin || x > our xmax)
- return 0.0;
- double result = 0.0;
- for (integer i = 1; i <= numberOfCoefficients; i ++) {
- if (coefficients [i] != 0.0)
- result += coefficients [i] * NUMispline2 (knots.get(), degree, i, x);
- }
- return result;
- }
- void structISpline :: v_evaluateTerms (double x, double *terms) {
- for (integer i = 1; i <= numberOfCoefficients; i ++)
- terms [i] = NUMispline2 (knots.get(), degree, i, x);
- }
- integer structISpline :: v_getOrder () {
- return degree;
- }
- Thing_implement (ISpline, Spline, 0);
- autoISpline ISpline_create (double xmin, double xmax, integer degree, integer numberOfInteriorKnots) {
- try {
- autoISpline me = Thing_new (ISpline);
- integer numberOfCoefficients = numberOfInteriorKnots + degree;
- integer numberOfKnots = numberOfCoefficients + degree;
- Spline_init (me.get(), xmin, xmax, degree, numberOfCoefficients, numberOfKnots);
- return me;
- } catch (MelderError) {
- Melder_throw (U"ISpline not created.");
- }
- }
- autoISpline ISpline_createFromStrings (double xmin, double xmax, integer degree, conststring32 coef, conststring32 interiorKnots) {
- try {
- if (degree > Spline_MAXIMUM_DEGREE)
- Melder_throw (U"Degree should be <= 20.");
- autoISpline me = Thing_new (ISpline);
- FunctionTerms_initFromString (me.get(), xmin, xmax, coef, 1);
- Spline_initKnotsFromString (me.get(), degree, interiorKnots);
- return me;
- } catch (MelderError) {
- Melder_throw (U"ISpline not created from strings.");
- };
- }
- /*
- #define RationalFunction_members Function_members \
- Polynomial num, denum;
- #define RationalFunction_methods Function_methods
- class_create (RationalFunction, Function)
- RationalFunction RationalFunction_create (double xmin, double xmax,
- integer degree_num, integer degree_denum)
- {
- RationalFunction me = new (RationalFunction);
- if (! me || ! (my num = Polynomial_create (xmin, xmax, degree_num)) ||
- ! (my denum = Polynomial_create (xmin, xmax, degree_denum))) forget (me);
- return me;
- }
- RationalFunction RationalFunction_createFromString (I, double xmin, double xmax,
- char *num, char *denum)
- {
- RationalFunction me = new (RationalFunction); integer i;
- if (! (my num = Polynomial_createFromString (xmin, xmax, num)) ||
- ! (my denum = Polynomial_createFromString (xmin, xmax, denum))) forget (me);
- if (my denum -> v [1] != 1 && my denum -> v [1] != 0)
- {
- double q0 = my denum -> v [1];
- for (i=1; 1 <= my num ->numberOfCoefficients; i ++) my num -> v [i] /= q0;
- for (i=1; 1 <= my denum ->numberOfCoefficients; i ++) my denum -> v [i] /= q0;
- }
- return me;
- }
- // divide out common roots
- RationalFunction RationalFunction_simplify (RationalFunction me)
- {
- Roots num = nullptr, denum = nullptr; RationalFunction thee = nullptr;
- if (! (num = Polynomial_to_Roots (my num)) ||
- ! (denum = Polynomial_to_Roots (my denum))) goto end;
- }
- */
- /* end of file Polynomial.cpp */
|