simple_linfit implemented

This commit is contained in:
Jiri Pittner 2021-11-22 14:22:19 +01:00
parent 060163d4c4
commit c45e3cc40c
2 changed files with 72 additions and 4 deletions

View File

@ -17,6 +17,7 @@
*/ */
//this header defines some simple algorithms independent of external libraries //this header defines some simple algorithms independent of external libraries
//using small runtime-constant size matrices and vectors
//particularly intended to embedded computers //particularly intended to embedded computers
//it should be compilable separately from LA as well as being a part of LA //it should be compilable separately from LA as well as being a part of LA
@ -37,7 +38,7 @@ namespace LA_Simple {
#define SWAP(a,b) {T temp=(a);(a)=(b);(b)=temp;} #define SWAP(a,b) {T temp=(a);(a)=(b);(b)=temp;}
template<typename T, int n, int m> template<typename T, int n, int m>
T simple_gaussj(T (&a)[n][n],T (&b)[m][n]) //returns determinant, m is number of rhs to solve T simple_gaussj(T (&a)[n][n],T (&b)[m][n]) //returns determinant, m is number of rhs to solve, inverse in a, solution in b
{ {
int indxc[n],indxr[n],ipiv[n]; int indxc[n],indxr[n],ipiv[n];
int i,icol,irow,j,k,l,ll; int i,icol,irow,j,k,l,ll;
@ -93,14 +94,57 @@ return det;
template<typename T, int n> template<typename T, int n>
class simple_linfit { class simple_linfit {
T mata[n][n];
T matb[1][n];
public: public:
T fitmat[n][n];
T rhsmat[1][n];
T fitcoef[n];
int npoints;
void clear() {npoints=0; memset(&fitmat[0][0],0,n*n*sizeof(T)); memset(&rhsmat[0][0],0,1*n*sizeof(T)); memset(&fitcoef[0],0,n*sizeof(T));};
simple_linfit() {clear();}
void input(const T (&funcs)[n], const T y)
{
++npoints;
for(int i=0; i<n; ++i)
{
for(int j=0; j<=i; ++j)
{
T tmp=funcs[i]*funcs[j];
fitmat[i][j] += tmp;
if(i!=j) fitmat[j][i] += tmp;
}
rhsmat[0][i] += funcs[i]*y;
}
}
T solve(bool preserve=false)
{
//for(int i=0; i<n; ++i) {for(int j=0; j<n; ++j) std::cout <<fitmat[i][j]<<" "; std::cout<<std::endl;}
//for(int j=0; j<n; ++j) std::cout <<rhsmat[0][j]<<" "; std::cout<<std::endl;
if(npoints<n) return 0;
if(preserve)
{
T fitwork[n][n];memcpy(fitwork,fitmat,n*n*sizeof(T));
T rhswork[1][n];memcpy(rhswork,rhsmat,1*n*sizeof(T));
T det = simple_gaussj(fitwork,rhswork);
memcpy(&fitcoef[0],&rhswork[0][0],n*sizeof(T));
return det;
}
T det = simple_gaussj(fitmat,rhsmat);
memcpy(&fitcoef[0],&rhsmat[0][0],n*sizeof(T));
return det;
}
}; };
//stream I/O //stream I/O
#ifndef AVOID_STDSTREAM #ifndef AVOID_STDSTREAM
template <typename T, int n>
std::ostream& operator<<(std::ostream &o, const simple_linfit<T,n> &f)
{
for(int i=0; i<n; ++i) o<<f.fitcoef[i]<<" ";
return o;
}
#endif #endif

26
t.cc
View File

@ -2354,7 +2354,7 @@ cout<<"eival error = "<<(w-www).norm()<<endl;
cout<<"eivec error = "<<(m.diffabs(vvv)).norm()<<endl; //just ignore signs due to arb. phases (not full check) cout<<"eivec error = "<<(m.diffabs(vvv)).norm()<<endl; //just ignore signs due to arb. phases (not full check)
} }
if(1) if(0)
{ {
//prepare random mat3 //prepare random mat3
int seed; int seed;
@ -2386,7 +2386,31 @@ cout<<"linear solve det="<<d<<endl;
cout <<r; cout <<r;
cout <<"det error="<<d-dd<<endl; cout <<"det error="<<d-dd<<endl;
cout <<"solution error="<<(r-rr).norm()<<endl; cout <<"solution error="<<(r-rr).norm()<<endl;
}
if(1)
{
//prepare random mat3
int seed;
int f=open("/dev/random",O_RDONLY);
if(sizeof(int)!=read(f,&seed,sizeof(int))) laerror("cannot read /dev/random");
close(f);
srand(seed);
simple_linfit<double,3> fit;
double funcs[3];
for(int i=0; i<100; ++i)
{
double x = 10*(2.*random()/(1. + RAND_MAX) - 1.);
funcs[0]=1;
funcs[1]=x;
funcs[2]=x*x;
double y = 2*funcs[2] -3*funcs[1] + funcs[0];
//cout <<"test "<<x<<" "<<y<<endl;
fit.input(funcs,y);
}
double det = fit.solve();
cout <<"det= "<<det<<" fit "<<fit<<endl;
} }
} }