Spektre - 2 months ago 16

C++ Question

**[Prologue]**

This **Q&A** is meant to explain more clearly the inner working of my approximations search class which I first published here

I was requested for more detailed info about this few times already (for various reasons) so I decided to write

How to approximate values/parameters in Real domain (

`double`

- real domain (precision)
`double`

**C++**language- configurable precision of approximation
- known interval for search
- fitted value/parameter is not strictly monotonic or not function at all

Answer

**Approximation search**

This is analogy to binary search but without its restrictions that searched function/value/parameter must be strictly monotonic function while sharing the `O(n.log(n))`

complexity.

**For example Let assume following problem**

We have known function `y=f(x)`

and want to find `x0`

such that `y0=f(x0)`

. This can be basically done by inverse function to `f`

but there are many functions that we do not know how to compute inverse to it. So how to compute this in such case?

**knowns**

`y=f(x)`

- input function`y0`

- wanted point`y`

value`a0,a1`

- solution`x`

interval range

**Unknowns**

`x0`

- wanted point`x`

value must be in range`x0=<a0,a1>`

**Algorithm**

**probe some points**`x(i)=<a0,a1>`

evenly dispersed along the range with some step`da`

So for example

`x(i)=a0+i*da`

where`i={ 0,1,2,3... }`

**for each**`x(i)`

compute the distance/error`ee`

of the`y=f(x(i))`

This can be computed for example like this:

`ee=fabs(f(x(i))-y0)`

but any other metrics can be used too.**remember point**`aa=x(i)`

with minimal distance/error`ee`

**stop when**`x(i)>a1`

**recursively increase accuracy**so first restrict the range to search only around found solution for example:

`a0'=aa-da; a1'=aa+da;`

then increase precision of search by lowering search step:

`da'=0.1*da;`

if

`da'`

is not too small or if max recursions count is not reached then go to**#1****found solution is in**`aa`

This is what I have in mind:

On the left side is the initial search illustrated (bullets **#1,#2,#3,#4**). On the right side next recursive search (bullet **#5**). This will recursively loop until desired accuracy is reached (number of recursions). Each recursion increase the accuracy `10`

times (`0.1*da`

). The gray vertical lines represent probed `x(i)`

points.

**Here the C++ source code for this:**

```
//---------------------------------------------------------------------------
//--- approx ver: 1.01 ------------------------------------------------------
//---------------------------------------------------------------------------
#ifndef _approx_h
#define _approx_h
#include <math.h>
//---------------------------------------------------------------------------
class approx
{
public:
double a,aa,a0,a1,da,*e,e0;
int i,n;
bool done,stop;
approx() { a=0.0; aa=0.0; a0=0.0; a1=1.0; da=0.1; e=NULL; e0=NULL; i=0; n=5; done=true; }
approx(approx& a) { *this=a; }
~approx() {}
approx* operator = (const approx *a) { *this=*a; return this; }
//approx* operator = (const approx &a) { ...copy... return this; }
void init(double _a0,double _a1,double _da,int _n,double *_e)
{
if (_a0<=_a1) { a0=_a0; a1=_a1; }
else { a0=_a1; a1=_a0; }
da=fabs(_da);
n =_n ;
e =_e ;
e0=-1.0;
i=0; a=a0; aa=a0;
done=false; stop=false;
}
void step()
{
if ((e0<0.0)||(e0>*e)) { e0=*e; aa=a; } // better solution
if (stop) // increase accuracy
{
i++; if (i>=n) { done=true; a=aa; return; } // final solution
a0=aa-fabs(da);
a1=aa+fabs(da);
a=a0; da*=0.1;
a0+=da; a1-=da;
stop=false;
}
else{
a+=da; if (a>a1) { a=a1; stop=true; } // next point
}
}
};
//---------------------------------------------------------------------------
#endif
//---------------------------------------------------------------------------
```

**This is how to use it:**

```
approx aa;
double ee,x,y,x0,y0=here_your_known_value;
// a0, a1, da,n, ee
for (aa.init(0.0,10.0,0.1,6,&ee); !aa.done; aa.step())
{
x = aa.a; // this is x(i)
y = f(x) // here compute the y value for whatever you want to fit
ee = fabs(y-y0); // compute error of solution for the approximation search
}
```

in the rem above `for (aa.init(...`

are the operand named. The `a0,a1`

is the interval on which the `x(i)`

is probed, `da`

is initial step between `x(i)`

and `n`

is the number of recursions. so if `n=6`

and `da=0.1`

the final max error of `x`

fit will be `~0.1/10^6=0.0000001`

. The `&ee`

is pointer to variable where the actual error will be computed. I choose pointer so there are not collisions when nesting this.

**[notes]**

This approximation search can be nested to any dimensionality (but of coarse you need to be careful about the speed) see some examples

- Approximation of n points to the curve with the best fit
- Curve fitting with y points on repeated x positions (Galaxy Spiral arms)
- Increasing accuracy of solution of transcendental equation
- Find Minimum area ellipse enclosing a set of points in c++

In case of non-function fit and the need of getting "all" the solutions you can use recursive subdivision of search interval after solution found to check for another solution. See example:

**What you should be aware of?**

you have to carefully choose the search interval `<a0,a1>`

so it contains the solution but is not too wide (or it would be slow). Also initial step `da`

is very important if it is too big you can miss local min/max solutions or if too small the thing will got too slow (especially for nested multidimensional fits).