From c45e3cc40c179e96a727bf068b8eb18d9b70ea99 Mon Sep 17 00:00:00 2001
From: Jiri Pittner <jiri@pittnerovi.com>
Date: Mon, 22 Nov 2021 14:22:19 +0100
Subject: [PATCH] simple_linfit implemented

---
 simple.h | 50 +++++++++++++++++++++++++++++++++++++++++++++++---
 t.cc     | 26 +++++++++++++++++++++++++-
 2 files changed, 72 insertions(+), 4 deletions(-)

diff --git a/simple.h b/simple.h
index ef01db5..46dccb9 100644
--- a/simple.h
+++ b/simple.h
@@ -17,6 +17,7 @@
 */
 
 //this header defines some simple algorithms independent of external libraries
+//using small runtime-constant size matrices and vectors
 //particularly intended to embedded computers
 //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;}
 
 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 i,icol,irow,j,k,l,ll;
@@ -93,14 +94,57 @@ return det;
 
 template<typename T, int n>
 class simple_linfit {
-T mata[n][n];
-T matb[1][n];
 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
 #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
 
 
diff --git a/t.cc b/t.cc
index 7b0d529..7a86beb 100644
--- a/t.cc
+++ b/t.cc
@@ -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)
 }
 
-if(1)
+if(0)
 {
 //prepare random  mat3
 int seed;
@@ -2386,7 +2386,31 @@ cout<<"linear solve det="<<d<<endl;
 cout <<r;
 cout <<"det error="<<d-dd<<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;
 }
 
 }