vectorBinDisect32 problems

[bswift]

sorry this was blank at first...

I'm running some psLib testing code on rma and I'm getting a warning I don't
think I should be getting. First, here's the code:

========================================
#include <stdio.h>
#include <math.h>
#include "pslib_strict.h"
#include "psTest.h"

static psS32 testVectorFitPolynomial1D_01(void);

/* list of test functions in test suite */
testDescription tests[] = {

{testVectorFitPolynomial1D_01,2,"psVectorFitPolynomial1D_01",0,false},

{NULL}

};

/* basic main function: just need to define the test suite name */
int main(int argc, char* argv[])
{

psTraceSetLevel("VectorFitPolynomial1DOrd", 6);
psLogSetFormat("HLNM");
psLogSetLevel(PS_LOG_INFO);
return !runTestSuite(stderr,"psPolynomialExamples", tests, argc, argv);

}

# define X0 +5.0
# define X1 +3.0
# define X2 -1.0
# define X3 -2.0

static psPolynomial1D *makePolynomial1D(int Nord, psPolynomialType type) {

define the polynomial
psPolynomial1D *poly = psPolynomial1DAlloc(Nord, type);
switch (Nord) {

case 3:

poly->coeff[3] = X3;

case 2:

poly->coeff[2] = X2;

case 1:

poly->coeff[1] = X1;

case 0:

poly->coeff[0] = X0;
break;

default:

break;

}
return (poly);

}

static void fillVectors(psVector *x, psVector *y, int Nord, int Npts,
psPolynomialType type) {

psF32 X;, P1, P2, P3;
double c[] = {(double)X0,(double)X1,(double)X2,(double)X3};
double sv,d,dd;
double sum, Tn, Tn_1, Tn_2;

for (int i = 0; i < Npts; i++) {

X = (float)i/Npts - 0.5;
x->data.F32[i] = X;
switch (type) {

case PS_POLYNOMIAL_ORD:

switch (Nord) {

case 0:

y->data.F32[i] = X0;
break;

case 1:

y->data.F32[i] = X0 + X1*X;
break;

case 2:

y->data.F32[i] = X0 + X1*X + X2*X*X;
break;

case 3:

y->data.F32[i] = X0 + X1*X + X2*X*X + X3*X*X*X;
break;

default:

y->data.F32[i] = 0;

}
break;

case PS_POLYNOMIAL_CHEB:

switch (Nord) {

case 0:

y->data.F32[i] = c[0];
break;

case 1:

y->data.F32[i] =c[0] + c[1]*X;
break;

default:

Tn_1=X; T1
Tn_2=1; T0
sum=0.;
for (int j=2; j<=Nord; j++) {

Tn = 2.0*X*Tn_1 - Tn_2;
sum = sum + c[j]*Tn;
Tn_2 = Tn_1;
Tn_1 = Tn;

}
y->data.F32[i] = sum + c[1]*X + c[0] - 0.5*c[0];

}
break;

}

}
return;

}

# define NPTS 10

int testVectorFitPolynomial1D_01(void) {

int retval=0;
int NORD, maxNORD=3;
float tolerance = 1.0E-5;
psPolynomialType type[] = {PS_POLYNOMIAL_ORD, PS_POLYNOMIAL_CHEB};

for (int k = 0; k<=1; k++) { loop over the 2 polynomial types

for (int j = 1; j <= maxNORD; j++) { loop over the polynomial orders

int nGood = 0, nBad = 0;
NORD = j;
printf("Testing NORD=%d\n", NORD);
psPolynomial1D *poly = makePolynomial1D(NORD, type[k]);
psPolynomial1D *polycopy = makePolynomial1D(NORD, type[k]);
psVector *f = psVectorAlloc(NPTS, PS_TYPE_F32);
psVector *fErr = psVectorAlloc(NPTS, PS_TYPE_F32);
psVector *x = psVectorAlloc(NPTS, PS_TYPE_F32);

fillVectors(x, f, NORD, NPTS, type[k]);
psTime *now = psTimeGetNow(PS_TIME_UTC);
const psRandom *r = psRandomAlloc(PS_RANDOM_TAUS, now->sec);
for (int i = 0; i < NPTS; i++) {

double errgen = psRandomGaussian(r);
double amp = 1.0E-2;
fErr->data.F32[i] = amp*errgen;
fErr->data.F32[i] = 0.;

}
psFree(now);
psFree(r);

f->data.F32[NPTS-1] = X0 + 10.*X1*x->data.F32[NPTS-1];
fErr->data.F32[NPTS-1] = 0.0000001;

psVectorFitPolynomial1D(poly, NULL, 0, f, fErr, x);
for (int i = 0; i <= NORD; i++) {

if (poly->coeff[i] <= polycopy->coeff[i]+tolerance

&&

poly->coeff[i] >= polycopy->coeff[i]-tolerance) {
nGood++;

} else {

nBad++;
printf("polycopy->coeff[%d]=%.20f

poly->coeff[%d]=%.20f\n",i,polycopy->coeff[i],i,poly->coeff[i]);

}
if (polycopy->coeffErr[i] != poly->coeffErr[i]) {

printf("polycopy->coeffErr[%d]=%.20f

poly->coeffErr[%d]=%.20f\n",i,polycopy->coeffErr[i],i,poly->coeffErr[i]);

}

}

psFree(poly);
psFree(polycopy);
psFree(f);
psFree(fErr);
psFree(x);

if (nBad > 0) {

psError(PS_ERR_UNKNOWN,true,"%d incorrect polynomial evaluations for

NORD=%d\n", nBad, NORD);

retval = retval + pow(2,NORD);

}
printf("NORD: %d retval: %d\n",NORD,retval);

}

}
return retval;

}

===============================================

here's one of the warnings it's giving me:

rma |W|vectorBinDisectF32

vectorBinDisectF32(): ordinate 0.987688 is outside vector range (-0.500000 -

0.400000).

This only happens in the second loop over k where the polynomial type is
PS_POLYNOMIAL_CHEB.

[gusciora@…]

The function psVectorFitPolynomial1D() is failing for this reason: when asked
to fit a Chebyshev polynomial to the x and y vectors, it assumes that the values
of the x vector span the range [-1.0:1.0]. In your test code, the x vector
spans the range [-0.5:0.4].

At present, I'm looking into how the current algorithm for deriving the
Chebyshev coefficients can be modified so that those x values need not fully
span the [-1.0:1.0] interval.

In the meantime, is it possible for you to modify your application so that the x
vectors are defined over the [-1.0:1.0] range?

[Paul Price]

Can't we treat the data over the range [-1,-0.5] and [0.5,1] as bad? Why is it
essential to have data over the entire range [-1,1] to fit a polynomial?

[bswift]

I used Chebyshev polynomials before spanning the same range [-0.5:0.4] and had
no problems there. I've used them to test psPolynomial1DEval(),
psPolynomial1DEvalVector(), and psVectorFitPolynomial1D() but left those tests
out of the code before. I can post those as well if you'd like.

[Paul Price]

No, the problem is in psVectorFitPolynomial.

[gusciora@…]

It's certainly possible to use Chebyshev polynomials over a range other than
[-1.0;1.0]. However, it's a tad difficult with our current framework for psLib.

For example, it's easy to scale the x-vector to [-1.0:1.0] upon entry to the
Chebyshev fitting routines. I've tried this, and it works. However, doing so
requires that the polynomial evaluation routines to scale the x-value as well.
And in order to do so, that routine would need to know the original range of the
x-vector: the problem is that, currently, there's no way to pass the original
range of the x-vector to the polynomial eval routine.

Here's one solution: create 2 private members of the polynomial struct called
->min and ->max. They refer to the original range of the x-vector when the
Chebyshev coefficients were derived, and they are set by the VectorFitPolynomial
routines as well. Then, when the polynomial eval routines are called, they
simply scale the x-value to [-1.0:1.0] based on those ->min and ->max members.
How does that sound?

Disclaimer: I am by no means and expert on this and I could easily be missing a
much simpler and obvious solution.

[Paul Price]

I think you want to use the standard, linear fitting method: generate the
matrix, invert, multiply. Otherwise, you're not using all the data.

[gusciora@…]

I'm not aware of a linear fitting method for Chebyshev polynomials (although I
know what that is for regular polynomials). Do you have a reference?

The current algorithm exploits the fact that Chebyshev polynomials are
orthogonal and at specific points (the Chebyshev nodes) in the [-1:1]
interval, they are either 0, 1, or -1. We calculate the coefficients based on
the function values at those points.

[Paul Price]

It's good to use the orthogonal property, but that means throwing away some of
the data --- by interpolating to pre-defined points, you're not using all the
data. This is not acceptable.

Please use the same algorithm as for the ordinary polynomials.

[gusciora@…]

This is the exact algorithm used in Numerical Recipes. The interpolation is
necessary since we're fitting the polynomials to predetermined data points,
not approximating a known function.

I can first fit an ordinary polynomial to the data, and then derive the
coefficients of the Chebyshev polynomials. This seems very odd to me.
However, if you make that a documented requirement, I can do it.

[Paul Price]

No, I mean use the same algorithm as for the ordinary polynomials (i.e.,
generate the least squares matrix and vector, invert the matrix, multiply by the
vector), but solve for the Chebyshev polynomial coefficients. That is, instead
of constructing the matrix using ordinary polynomials, you do it using the
Chebyshev polynomials.

Interpolating is not an option, since it throws away data.

[gusciora@…]

I implemented the least squares matrix and inverted the matrix. This worked.
However, it's CPU-intensive.

A faster algorithm is this. Simply fit an ordinary polynomial to the data,
then "reverse-engineer" those coefficients to determine the Chebyshev
coefficients: basically, for each c_0x^{i term in the ordinary polynomial, we
derive a linear set of equations based on all the x}i terms in all the
component Cheby polynomials, and then solve this set of equations to determine
the Chebyshev coefficients. This set of equations is already in
upper-triangular form so it's inexpensive to solve. This method is in the
current CVS tree.

One last note: I'd still like to keep the older method for deriving the Cheby
coefficients in the code in a commented-out form. The reason is that I
anticipate someone, later on, beginning a discussion that starts like this:
"Why bother using Chebyshev polys if in order to calculate their coefficients
we first have to derive the ordinary polynomial? Shouldn't we simply use the
ordinary polynomial instead?" The answer, I believe, is that you might use
the Chebyshev polys when you want a much faster algorithm for deriving them,
though you'll sacrifice a little accuracy. The older method provides that.

[Paul Price]

vectorFitPolynomial1DChebyFast should never be used, unless the input data is
regularly distributed and not masked (rare, if ever).

The "slow" implementation is the correct one. Could you please optimise it by:

• Using psPolynomial1DEvalVector within a loop on the polynomials instead of

evaluating inside a double loop.

• Calculating A and B in the same loop
• Using the fact that the least squares matrix is symmetric, so that only ~ 1/2

the calculations are required.

[Paul Price]

Bugs have been fixed... closing.