InterPolynomial Class Reference

#include <IntPoly.h>

Inheritance diagram for InterPolynomial:

Collaboration diagram for InterPolynomial:

Public Member Functions
double *	NewtonCoefficient (double *)
void	ComputeLagrangeBasis (double *, unsigned nPtsTotal)
void	GenerateBasis (double rho, double rhosmall, Matrix data, ObjectiveFunction *of)
	InterPolynomial (unsigned _deg, double rho, Vector vBase, Matrix data, ObjectiveFunction *of)
int	findAGoodPointToReplace (int excludeFromT, double rho, Vector pointToAdd, double *modelStep=NULL)
void	replace (int t, Vector pointToAdd, double valueF)
int	maybeAdd (Vector pointToAdd, unsigned k, double rho, double valueF)
void	updateM (Vector newPoint, double valueF)
int	checkIfValidityIsInBound (Vector dd, unsigned k, double bound, double rho)
int	getGoodInterPolationSites (Matrix d, int k, double rho, Vector *v=NULL)
double	interpError (Vector Point)
void	translate (int k)
void	translate (Vector translation)
	~InterPolynomial ()
	InterPolynomial (const InterPolynomial &A)
InterPolynomial &	operator= (const InterPolynomial &A)
InterPolynomial	clone ()
void	copyFrom (InterPolynomial a)
void	copyFrom (Polynomial a)
	InterPolynomial (unsigned dim, unsigned deg)
Public Member Functions inherited from Polynomial
	Polynomial ()
	Polynomial (unsigned Dim, unsigned deg=0, double *data=0)
	Polynomial (unsigned Dim, double val)
	Polynomial (MultInd &)
	Polynomial (char *name)
unsigned	dim ()
unsigned	deg ()
unsigned	sz ()
	operator double * () const
	~Polynomial ()
	Polynomial (const Polynomial &A)
Polynomial &	operator= (const Polynomial &A)
Polynomial	clone ()
void	copyFrom (Polynomial a)
Polynomial	operator* (const double)
Polynomial	operator/ (const double)
Polynomial	operator+ (Polynomial)
Polynomial	operator- (Polynomial)
Polynomial	operator- (void)
Polynomial	operator+ (void)
Polynomial	operator+= (Polynomial)
Polynomial	operator-= (Polynomial)
Polynomial	operator*= (const double)
Polynomial	operator/= (const double)
double	shiftedEval (Vector Point, double minusVal)
double	operator() (Vector)
Polynomial	derivate (int i)
void	gradient (Vector P, Vector G)
void	gradientHessian (Vector P, Vector G, Matrix H)
void	translate (Vector translation)
int	operator== (const Polynomial q)
int	equals (Polynomial q)
void	print ()
void	save (char *name)
void	setFlag (unsigned int val)
void	unsetFlag (unsigned int val)
unsigned	queryFlag (unsigned int val)

Public Attributes
double	M
unsigned	nPtsUsed
unsigned	nUpdateOfM
Polynomial *	NewtonBasis
Vector *	NewtonPoints
Vector	vBase
double *	valuesF
int	kbest

Protected Member Functions
void	destroyCurrentBuffer ()
Protected Member Functions inherited from Polynomial
void	init (int _dim, int _deg, double *data=NULL)
void	destroyCurrentBuffer ()

Additional Inherited Members
Static Public Attributes inherited from Polynomial
static const unsigned int	NicePrint = 1
static const unsigned int	Warning = 2
static const unsigned int	Normalized = 4
static unsigned int	flags = Polynomial::Warning||Polynomial::NicePrint
static Polynomial	emptyPolynomial
Protected Types inherited from Polynomial
typedef struct Polynomial::PolynomialDataTag	PolynomialData
Protected Attributes inherited from Polynomial
PolynomialData *	d

Detailed Description

Definition at line 35 of file IntPoly.h.

Constructor & Destructor Documentation

InterPolynomial() [1/3]

InterPolynomial::InterPolynomial	(	unsigned	_deg,
		double	rho,
		Vector	vBase,
		Matrix	data,
		ObjectiveFunction *	of
	)

Definition at line 228 of file IntPoly.cpp.

    : Polynomial(data.nColumn()-2, _deg), M(0.0), nUpdateOfM(0), NewtonBasis(NULL), NewtonPoints(NULL), kbest(-1)
{
    nPtsUsed=choose( _deg+d->dim,d->dim );
    if (data.nLine()==0)
    {
        printf( "InterPolynomial::InterPolynomial( double *) : No data\n");
        getchar(); exit(-1);
    }

    // Generate basis 
    if (_vBase==Vector::emptyVector)
    {
        double rhosmall=rho*ROOT2*.5;
        GenerateBasis(rho,rhosmall, data,NULL);
        if (!NewtonBasis) 
        {
            GenerateBasis(rho,rho*.5, data,of);
        }
    } else
    {
        int i,ii,j=0,k,mdim=dim();
        double norm, *p, *pb=_vBase;
        KeepBests kb(nPtsUsed*2,mdim);
//        fprintf(stderr,"Value Objective=%e\n", vBase);
        
        ii=data.lineIndex(_vBase);
        if (ii!=-1)
        {
            p=data[ii];
            kb.add(0.0,p[mdim],p); j++; 
        }

        i=data.nLine();
        while (i--)
        {
            if (i==ii) continue;
            p=data[i];
            if (p[mdim+1]) continue;
            norm=0; k=mdim;
            while (k--) norm+=sqr(p[k]-pb[k]);
            if ((rho==0.0)||
                ((norm<=4.00*rho*rho))//&&(norm>.25*rhomin*rhomin))
               ) 
            {
                kb.add(norm,p[mdim], p); j++;
            }
        }
        if (j<(int)nPtsUsed)
        {
            printf("Not Enough points in DB to build interpolation polynomial\n");
            return;
        }
        // we have retained only the 2*nPtsUsed best points:
        j=mmin(j,(int)(2*nPtsUsed));
        NewtonPoints=new Vector[j];
        valuesF=(double*)malloc(j*sizeof(double));
        for (i=0; i<j; i++)
        {
            valuesF[i]=kb.getValue(i);
            NewtonPoints[i]=Vector(mdim, kb.getOptValue(i));
        }
        //kbest=findK(valuesF,j,NULL,NewtonPoints);
        //vBase=NewtonPoints[kbest].clone();
        vBase=_vBase;
        for (i=0; i<j; i++) NewtonPoints[i]-=vBase;
        ComputeLagrangeBasis(valuesF, j);
    }
    if (NewtonBasis==NULL) return;

//    test();

    // Compute Interpolant
//    double *NewtonCoefPoly=NewtonCoefficient(_Yp);
    double *NewtonCoefPoly=valuesF;

    double *coc= NewtonCoefPoly+nPtsUsed-1;
    Polynomial *ppc= NewtonBasis+nPtsUsed-1;
    this->copyFrom((Polynomial)(*coc * *ppc));        //take highest degree

    int i=nPtsUsed-1;
    if (i)
      while(i--)
        (*this) += *(--coc) * *(--ppc);    
                // No reallocation here because of the order of
                // the summation
    //free(NewtonCoefPoly);
}

~InterPolynomial()

InterPolynomial::~InterPolynomial

(

)

Definition at line 819 of file IntPoly.cpp.

{
    destroyCurrentBuffer();
}

InterPolynomial() [2/3]

InterPolynomial::InterPolynomial

(

const InterPolynomial &

A

)

Definition at line 824 of file IntPoly.cpp.

{
    // shallow copy for inter poly.
    d=A.d;
    NewtonBasis=A.NewtonBasis;
    NewtonPoints=A.NewtonPoints;
//    ValuesF=A.ValuesF;
    M=A.M;
    nPtsUsed=A.nPtsUsed; 
    nUpdateOfM=A.nUpdateOfM;
    (d->ref_count)++;
    
}

InterPolynomial() [3/3]

InterPolynomial::InterPolynomial	(	unsigned	dim,
		unsigned	deg
	)

Definition at line 861 of file IntPoly.cpp.

                                                            : Polynomial(_dim, _deg), M(0.0), nUpdateOfM(0)
{   
    nPtsUsed=choose( _deg+_dim,_dim );
    NewtonBasis= new Polynomial[nPtsUsed];
    NewtonPoints= new Vector[nPtsUsed];
}

Member Function Documentation

checkIfValidityIsInBound()

int InterPolynomial::checkIfValidityIsInBound	(	Vector	dd,
		unsigned	k,
		double	bound,
		double	rho
	)

Definition at line 684 of file IntPoly.cpp.

{

    // check validity around x_k
    // bound is epsilon in the paper
    // return index of the worst point of J
    // if (j==-1) then everything OK : next : trust region step
    // else model step: replace x_j by x_k+d where d
    // is calculated with LAGMAX

    unsigned i,N=nPtsUsed, n=dim();
    int j;
    Polynomial *pp=NewtonBasis;
    Vector *xx=NewtonPoints, xk=xx[k],vd;
    Vector Distance(N);
    double *dist=Distance, *dd=dist, distMax, vmax, tmp;
    Matrix H(n,n);
    Vector GXk(n); //,D(n);

    for (i=0; i<N; i++) *(dd++)=xk.euclidianDistance(xx[i]);

    while (true)
    {
        dd=dist; j=-1; distMax=2*rho;
        for (i=0; i<N; i++)
        {
            if (*dd>distMax) { j=i; distMax=*dd; };
            dd++;
        }
        if (j<0) return -1;

        // to prevent to choose the same point once again:
        dist[j]=0;
        
        pp[j].gradientHessian(xk,GXk,H);
//        d=H.multiply(xk);
//        d.add(G);

        tmp=M*distMax*distMax*distMax;

        if (tmp*rho*(GXk.euclidianNorm()+0.5*rho*H.frobeniusNorm())>=bound)
        {
/*          vd=L2NormMinimizer(pp[j], xk, rho);
            vd+=xk;
            vmax=condorAbs(pp[j](vd));

            Vector vd2=L2NormMinimizer(pp[j], xk, rho);
            vd2+=xk;
            double vmax2=condorAbs(pp[j](vd));
            
            if (vmax<vmax2) { vmax=vmax2; vd=vd2; }
*/
            vd=LAGMAXModified(GXk,H,rho,vmax);
//            tmp=vd.euclidianNorm();
            vd+=xk;
            vmax=condorAbs(pp[j](vd));
            
            if (tmp*vmax>=bound) break;
        }
    }
    if (j>=0) ddv.copyFrom(vd);
    return j;
}

clone()

InterPolynomial InterPolynomial::clone

(

)

Definition at line 868 of file IntPoly.cpp.

{
    // a deep copy
    InterPolynomial m(d->dim, d->deg);
    m.copyFrom(*this);
    return m;
}

ComputeLagrangeBasis()

void InterPolynomial::ComputeLagrangeBasis	(	double *	yy,
		unsigned	nPtsTotal
	)

Definition at line 127 of file IntPoly.cpp.

{
#ifdef NOOBJECTIVEFUNCTION
   const double eps  =  1e-6;
#else
   const double eps  =  1e-10;
#endif
   const double good =  1 ;
   unsigned dim=d->dim,i;
  
   Vector *xx=NewtonPoints, xtemp;
   Polynomial *pp= new Polynomial[nPtsUsed], *qq=pp;
   NewtonBasis=pp;

   if (!pp)
   {
      printf("ComputeNewtonBasis( ... ) : Alloc for polynomials failed\n");
      getchar(); exit(251);
   }
  
   MultInd I(dim);
   for (i=0; i<nPtsUsed; i++)
   {
      *(qq++)=Polynomial(I);
      I++;
   }
    
   unsigned k, kmax;
   double v, vmax, vabs;
#ifdef VERBOSE
   printf("Constructing first quadratic ... (N=%i)\n",nPtsUsed);
#endif
   for (i=0; i<nPtsUsed; i++)
   {
#ifdef VERBOSE
      printf(".");
#endif
      vmax = vabs = 0;
      kmax = i;
      if (i==0)
      {
          // to be sure point 0 is always taken:
          vmax=(pp[0])( xx[0] );
      } else
      for (k=i; k<nPtsTotal; k++)    // Pivoting
      {
          v=(pp[i])( xx[k] );
          if (fabs(v) > vabs)
          {
              vmax = v;
              vabs = condorAbs(v);
              kmax = k;
          }
          if (fabs(v) > good ) break;
      }

      // Now, check ...
      if (fabs(vmax) < eps)
      {
          printf("Cannot construct Lagrange basis");
          delete[] NewtonBasis;
          NewtonBasis=NULL;
          return;
      }

      // exchange component i and k of NewtonPoints
      // fast because of shallow copy
      xtemp   =xx[kmax];
      xx[kmax]=xx[i];
      xx[i]   =xtemp;

      // exchange component i and k of newtonData
      v       =yy[kmax];
      yy[kmax]=yy[i];
      yy[i]   =v;

      pp[i]/=vmax;
      for (k=0;  k<i;        k++) pp[k] -= (pp[k])( xx[i] ) * pp[i];
      for (k=i+1; k<nPtsUsed; k++) pp[k] -= (pp[k])( xx[i] ) * pp[i];

      // Next polynomial, break if necessary
   }
#ifdef VERBOSE
   printf("\n");
#endif
}

copyFrom() [1/2]

void InterPolynomial::copyFrom

(

InterPolynomial

a

)

Definition at line 876 of file IntPoly.cpp.

{
    if (m.d->dim!=d->dim)
    {
        printf("poly: copyFrom: dim do not agree\n");
        getchar(); exit(254);
    }
    if (m.d->deg!=d->deg)
    {
        printf("poly: copyFrom: degree do not agree\n");
        getchar(); exit(254);
    }
    Polynomial::copyFrom(m);
    M=m.M;
//    nPtsUsed=m.nPtsUsed; // not useful because dim and degree already agree.
    nUpdateOfM=m.nUpdateOfM;
// ValuesF        

    int i=nPtsUsed;
    while (i--) 
    { 
//      NewtonBasis[i]=m.NewtonBasis[i];
//      NewtonPoints[i]=m.NewtonPoints[i]; 
        NewtonBasis[i]=m.NewtonBasis[i].clone();
        NewtonPoints[i]=m.NewtonPoints[i].clone(); 
    } 
}

copyFrom() [2/2]

void InterPolynomial::copyFrom

(

Polynomial

a

)

Definition at line 904 of file IntPoly.cpp.

{
    Polynomial::copyFrom(m);
}

destroyCurrentBuffer()

void InterPolynomial::destroyCurrentBuffer ( )

protected

Definition at line 805 of file IntPoly.cpp.

{
    if (!d) return;
    if (d->ref_count==1)
    {
        if (NewtonBasis)
        {
            delete[] NewtonBasis;
            delete[] NewtonPoints;
            free(valuesF);
        }
    }
}

findAGoodPointToReplace()

int InterPolynomial::findAGoodPointToReplace	(	int	excludeFromT,
		double	rho,
		Vector	pointToAdd,
		double *	modelStep = NULL
	)

Definition at line 551 of file IntPoly.cpp.

{
    //not tested

    // excludeFromT is set to k if not success from optim and we want
    //        to be sure that we keep the best point
    // excludeFromT is set to -1 if we can replace the point x_k by
    //        pointToAdd(=x_k+d) because of the success of optim.
    
    // chosen t: the index of the point inside the newtonPoints
    //            which will be replaced.

    Vector *xx=NewtonPoints;
    Vector XkHat;
    if (excludeFromT>=0) XkHat=xx[excludeFromT];
    else XkHat=pointToAdd;

    int t=-1, i, N=nPtsUsed;
    double a, aa, maxa=-1.0, maxdd=0;
    Polynomial *pp=NewtonBasis;

  //  if (excludeFromT>=0) maxa=1.0;

    for (i=0; i<N; i++)
    {
        if (i==excludeFromT) continue;
        aa=XkHat.euclidianDistance(xx[i]);
        if (aa==0.0) return -1;
        a=aa/rho;
        // because of the next line, rho is important:
        a=mmax(a*a*a,1.0);
        a*=condorAbs(pp[i] (pointToAdd));

        if (a>maxa)
        {
            t=i; maxa=a; maxdd=aa;
        }
    }
    if (maxd) *maxd=maxdd;
    return t;
}

GenerateBasis()

void InterPolynomial::GenerateBasis	(	double	rho,
		double	rhosmall,
		Matrix	data,
		ObjectiveFunction *	of
	)

Definition at line 324 of file IntPoly.cpp.

{
    if (deg()!=2)
    {
        printf("The GenerateBasis function only works for interpolation polynomials of degree 2.\n");
        exit(244);
    }
    int i,j,k,l,dim=data.nColumn()-2, N=(dim+1)*(dim+2)/2,
        *failed=(int*)malloc(N*sizeof(int)),nToCompute, nl=data.nLine(),best;
    Vector Base=data.getLine(0,dim);
    VectorInt vi(dim);
    valuesF=(double*)malloc((N+dim)*sizeof(double));
    double dBase=((double**)data)[0][dim],*vrho=valuesF+N,
           *base=Base,*dcur,*p,*pb, norm, normbest, r; // value objective function
    
    NewtonPoints=new Vector[N];
    Vector *cp, cur; // cp: current Point
    NewtonPoints[N-1]=Base.clone(); valuesF[N-1]=dBase; 

    // generate perfect sampling site
    for (i=0; i<dim; i++) 
    { 
        cur=Base.clone();
        cur[i]+=rho;
        NewtonPoints[i]=cur;
        vi[i]=i;
    }

    nToCompute=0; cp=NewtonPoints;
    // try to find good points in DB, which are near perfect sampling sites
    for (i=0; i<dim; i++)
    {
        normbest=INF; pb=cp[i];
        for (j=1; j<nl; j++)
        {
            p=data[j];
            if (p[dim+1]) continue;
            norm=0; k=dim;
            while (k--) norm+=sqr(p[k]-pb[k]);
            if (normbest>norm) { normbest=norm; best=j; }
        }
        if (normbest<rhosmall*rhosmall)
        {
            cp[i]=data.getLine(best,dim);
            valuesF[i]=((double**)data)[best][dim];
        } else
        {
            cur=cp[nToCompute]; cp[nToCompute]=cp[i]; cp[i]=cur;
            l=vi[nToCompute];   vi[nToCompute]=vi[i]; vi[i]=l;
            valuesF[i]=valuesF[nToCompute];
            nToCompute++;
        }
    }

    if ((!of)&&(nToCompute)) { free(failed); free(valuesF); delete[] NewtonPoints; return; }

    // if some points have not been found in DB, start evaluation
    while (nToCompute>0)
    {
        // parallelize here !
        calculateNParallelJob(nToCompute,valuesF,NewtonPoints,of,failed);

        k=0;
        for (j=0; j<nToCompute; j++)
        {
            of->saveValue(cp[j],valuesF[j],failed[j]);
            if (failed[j])
            {
                dcur=cp[j];
                for (i=0; i<dim; i++)
                {
                    dcur[i]=base[i]+(base[i]-dcur[i])*.7;
                }

                // group all missing values at the beginning (around NewtonPoints+dim)
                cur=cp[k]; cp[k]=cp[j]; cp[j]=cur;
                l=vi[k];   vi[k]=vi[j]; vi[j]=l;
                valuesF[j]=valuesF[k];
                k++;
            }
        }
        nToCompute=k;
    }

    // re-order vectors (based on vi)
    for (j=0; j<dim-1; j++)
    {
        k=j; while (vi[k]!=j) k++;
        if (k==j) continue;
        cur=cp[k];    cp[k]=cp[j];           cp[j]=cur;
        l=vi[k];      vi[k]=vi[j];           vi[j]=l;
        r=valuesF[k]; valuesF[k]=valuesF[j]; valuesF[j]=r;
    }

    // select again some good sampling sites
    cp=NewtonPoints+dim;
    for (j=0; j<dim; j++)
    {
        cur=Base.clone();
        dcur=NewtonPoints[j];
        r=dcur[j]-base[j];
        if (valuesF[j]<dBase) { cur[j]+=r*2.0; vrho[j]=r; }
        else { cur[j]-=r; vrho[j]=-r; }
        *(cp++)=cur;
    }
    for (j=0; j<dim; j++)
    {
        for (k=0; k<j; k++)
        {
            cur=Base.clone();
            cur[j]+=vrho[j];
            cur[k]+=vrho[k];
            *(cp++)=cur;
        }
    }

    nToCompute=0;
    cp=NewtonPoints+dim;
    // try to find good points in DB, which are near perfect sampling sites
    for (i=0; i<N-dim-1; i++)
    {
        normbest=INF; pb=cp[i];
        for (j=1; j<nl; j++)
        {
            p=data[j];
            if (p[dim+1]) continue;
            norm=0; k=dim;
            while (k--) norm+=sqr(p[k]-pb[k]);
            if (normbest>norm) { normbest=norm; best=j; }
        }
        if (normbest<rhosmall*rhosmall)
        {
            cp[i]=data.getLine(best,dim);
            valuesF[i+dim]=((double**)data)[best][dim];
        } else
        {
            cur=cp[nToCompute];
            cp[nToCompute]=cp[i];
            cp[i]=cur;
            valuesF[i+dim]=valuesF[nToCompute+dim];
            nToCompute++;
        }
    }

    if ((!of)&&(nToCompute)) { free(failed); free(valuesF); delete[] NewtonPoints; return; }

    while (nToCompute>0)
    {
        // parallelize here !
        calculateNParallelJob(nToCompute,valuesF+dim,cp,of,failed);

        k=0;
        for (j=0; j<nToCompute; j++)
        {
            of->saveValue(cp[j],valuesF[j+dim],failed[j]);
            if (failed[j])
            {
                dcur=cp[j];
                for (i=0; i<dim; i++)
                {
                    dcur[i]=base[i]+(base[i]-dcur[i])*.65;
                }

                // group all missing values at the beginning (around NewtonPoints+dim)
                cur=cp[k];
                cp[k]=cp[j];
                cp[j]=cur;
                valuesF[j+dim]=valuesF[k+dim];
                k++;
            }
        }
        nToCompute=k;

    }
    free(failed); 

    // "insertion sort" to always place best points first
    for (i=0; i<N-1; i++)
    {
        kbest=findK(valuesF+i,N-i,of,NewtonPoints+i)+i; 
        
        // to accept infeasible points for i>=1 :
        of=NULL;

        cur=NewtonPoints[i]; NewtonPoints[i]=NewtonPoints[kbest]; NewtonPoints[kbest]=cur;
        r=valuesF[i];        valuesF[i]=valuesF[kbest];           valuesF[kbest]=r;
    }
    kbest=0;
    vBase=NewtonPoints[0].clone();
    for (i=0; i<N; i++) NewtonPoints[i]-=vBase;
    ComputeLagrangeBasis(valuesF, N);
}

getGoodInterPolationSites()

int InterPolynomial::getGoodInterPolationSites	(	Matrix	d,
		int	k,
		double	rho,
		Vector *	v = NULL
	)

Definition at line 750 of file IntPoly.cpp.

{
    //not tested

    unsigned i, N=nPtsUsed, n=dim();
    int ii, j,r=0;
    Polynomial *pp=NewtonBasis;
    Vector *xx=NewtonPoints, xk,vd;
    Vector Distance(N);
    double *dist=Distance, *dd=dist, distMax, vmax;
    Matrix H(n,n);
    Vector GXk(n);
    if (k>=0) xk=xx[k]; else xk=*v;

    for (i=0; i<N; i++) *(dd++)=xk.euclidianDistance(xx[i]);

    for (ii=0; ii<d.nLine(); ii++)
    {
        dd=dist; j=-1; 
        distMax=-1.0;
        for (i=0; i<N; i++)
        {
            if (*dd>distMax) { j=i; distMax=*dd; };
            dd++;
        }
        // to prevent to choose the same point once again:
        dist[j]=-1.0;
        
        if (distMax>2*rho) r++;
        pp[j].gradientHessian(xk,GXk,H);
        vd=LAGMAXModified(GXk,H,rho,vmax);
        vd+=xk;

        d.setLine(ii,vd);
    }
    return r;
}

interpError()

double InterPolynomial::interpError

(

Vector

Point

)

Definition at line 536 of file IntPoly.cpp.

{
    unsigned i=nPtsUsed;
    double sum=0,a;
    Polynomial *pp=NewtonBasis;
    Vector *xx=NewtonPoints;

    while (i--)
    {
        a=Point.euclidianDistance(xx[i]);
        sum+=condorAbs(pp[i]( Point ))*a*a*a;
    }
    return M*sum;
}

maybeAdd()

int InterPolynomial::maybeAdd	(	Vector	pointToAdd,
		unsigned	k,
		double	rho,
		double	valueF
	)

Definition at line 649 of file IntPoly.cpp.

{

    unsigned i,N=nPtsUsed;
    int j = 0;
    Vector *xx=NewtonPoints, xk=xx[k];
    double distMax=-1.0,dd;
 /*
    Polynomial *pp=NewtonBasis;
    Vector *xx=NewtonPoints, xk=xx[k],vd;
    Matrix H(n,n);
    Vector GXk(n); //,D(n);
*/
    // find worst point/newton poly

    for (i=0; i<N; i++)
    {
        dd=xk.euclidianDistance(xx[i]);
        if (dd>distMax) { j=i; distMax=dd; };
    }
    dd=xk.euclidianDistance(pointToAdd);

    // no tested:

    if (condorAbs(NewtonBasis[j](pointToAdd))*distMax*distMax*distMax/(dd*dd*dd)>1.0) 
    {
        printf("good point found.\n");
        replace(j, pointToAdd, valueF);
        return 1;
    }
    return 0;
}

NewtonCoefficient()

double* InterPolynomial::NewtonCoefficient

(

double *

)

operator=()

InterPolynomial & InterPolynomial::operator=

(

const InterPolynomial &

A

)

Definition at line 838 of file IntPoly.cpp.

{
    // shallow copy
    if (this != &A)
    {
        destroyCurrentBuffer();

        d=A.d;
        NewtonBasis=A.NewtonBasis;
        NewtonPoints=A.NewtonPoints;
//        ValuesF=A.ValuesF;
        M=A.M;
        nPtsUsed=A.nPtsUsed; 
        nUpdateOfM=A.nUpdateOfM;
        kbest=A.kbest;
        vBase=A.vBase;
        valuesF=A.valuesF;

        (d->ref_count) ++ ;
    }
    return *this;
}

replace()

void InterPolynomial::replace	(	int	t,
		Vector	pointToAdd,
		double	valueF
	)

Definition at line 622 of file IntPoly.cpp.

{
    //not tested
    updateM(pointToAdd, valueF);
    if (t<0) return;

    Vector *xx=NewtonPoints;
    Polynomial *pp=NewtonBasis, t1;
    int i, N=nPtsUsed;
    double t2=(pp[t]( pointToAdd ));

    if (t2==0) return;

    t1=pp[t]/=t2;
  
    for (i=0; i<t; i++)   pp[i]-= pp[i]( pointToAdd )*t1;
    for (i=t+1; i<N; i++) pp[i]-= pp[i]( pointToAdd )*t1;
    xx[t].copyFrom(pointToAdd);

    // update the coefficents of general poly.

    valueF-=(*this)(pointToAdd);
    if (condorAbs(valueF)>1e-11) (*this)+=valueF*pp[t];
    
//    test();
}

translate() [1/2]

void InterPolynomial::translate

(

int

k

)

translate() [2/2]

void InterPolynomial::translate

(

Vector

translation

)

Definition at line 790 of file IntPoly.cpp.

{
    if (translation==Vector::emptyVector) return;
    vBase+=translation;
    Polynomial::translate(translation);
    int i=nPtsUsed;
    while (i--) NewtonBasis[i].translate(translation);
    i=nPtsUsed;
    while (i--) if (NewtonPoints[i]==translation) NewtonPoints[i].zero();
                else NewtonPoints[i]-=translation;
}

updateM()

void InterPolynomial::updateM	(	Vector	newPoint,
		double	valueF
	)

Definition at line 519 of file IntPoly.cpp.

{
    //not tested
    unsigned i=nPtsUsed;
    double sum=0,a;
    Polynomial *pp=NewtonBasis;
    Vector *xx=NewtonPoints;

    while (i--)
    {
        a=newPoint.euclidianDistance(xx[i]);
        sum+=condorAbs(pp[i]( newPoint ))*a*a*a;
    }
    M=mmax(M, condorAbs((*this)(newPoint)-valueF)/sum);
    nUpdateOfM++;
}

Member Data Documentation

kbest

int InterPolynomial::kbest

Definition at line 49 of file IntPoly.h.

M

double InterPolynomial::M

Definition at line 39 of file IntPoly.h.

NewtonBasis

Polynomial* InterPolynomial::NewtonBasis

Definition at line 43 of file IntPoly.h.

NewtonPoints

Vector* InterPolynomial::NewtonPoints

Definition at line 47 of file IntPoly.h.

nPtsUsed

unsigned InterPolynomial::nPtsUsed

Definition at line 40 of file IntPoly.h.

nUpdateOfM

unsigned InterPolynomial::nUpdateOfM

Definition at line 40 of file IntPoly.h.

valuesF

double* InterPolynomial::valuesF

Definition at line 48 of file IntPoly.h.

vBase

Vector InterPolynomial::vBase

Definition at line 47 of file IntPoly.h.

---

The documentation for this class was generated from the following files:

• xmipp/external/condor/IntPoly.h
• xmipp/external/condor/IntPoly.cpp