minor code cleaning, improved documentation and invented .txt export in...

minor code cleaning, improved documentation and invented .txt export in SplineTest.cpp as well as an .m to plot the created splines in matlab

include/ITABase/Math/LinearAlgebra.h

@@ -26,8 +26,13 @@ namespace ITABase
 {
 	namespace Math
 	{
-		//! Solves a Tridiagonal Bandmatrix with Thomas algorithm. Expects the 3 diagonals (lower,middel and upper) and the excitation vector.
-		ITA_BASE_API std::vector<float> BandmatrixSolver(const std::vector<float> vdLower, const std::vector<float> vdMiddle, const std::vector<float> vdUpper, const std::vector<float> vdExcitation);
+		//! Solves a Tridiagonal n x n Bandmatrix Problem with Thomas algorithm. Expects the 3 diagonals (lower,middle and upper) and the excitation vector.
+		/**
+		* Problems following the structure: e = M * x will be solved with this function.
+		* with M is a Matrix with only the 3 main diagonals nonequal zero, e the known excitation vector and x the vector of unknown variables.
+		* Since the matrix dimension is n x n the excitatio, the middle and the returned vector are n elements big. Therefore the lower and the upper diagonale needs to have a length of n-1.
+		*/
+		ITA_BASE_API std::vector<float> BandmatrixSolver(const std::vector<float>& vdLower, const std::vector<float>& vdMiddle, const std::vector<float>& vdUpper, const std::vector<float>& vdExcitation);
 	}
 }
 

include/ITABase/Math/PiecewisePolynomial.h

@@ -51,13 +51,15 @@ namespace ITABase
 		inline const std::vector<float>& BreakPoints() const { return vdBreakPoints; };
 
 		//! Evaluate the piecewise polynomial at a given x value.
-		float Value(const float xValue) const;
+		float Value(const float & xValue) const;
 		//! Evaluates the piecewise polynomial at all given x values.
-		std::vector<float> Value(const std::vector<float>& vdValues) const;
+		std::vector<float> Value(const std::vector<float>& vdXValues) const;
 		//! Finds the index of Interval at a given x value
 		int IntervalIndex(const float xValue) const;
 		//! Return a piecewise polynomial containing the derivations of the polynomials of this one
 		CPiecewisePolynomial Derivation() const;
+		//! Debug only?!
+		std::vector<float>::const_iterator getCoefficients() const { return vdCoefficients.begin(); }
 
 	private:
 		//! Returns the iterator to the first coefficient of the interval with given index

src/ITABase/Math/LinearAlgebra.cpp

@@ -5,30 +5,30 @@
 
 using namespace ITABase;
 
- std::vector<float> Math::BandmatrixSolver(const std::vector<float> vdLower,const std::vector<float> vdMiddle, const std::vector<float> vdUpper, const std::vector<float> vdExcitation)
+ std::vector<float> Math::BandmatrixSolver(const std::vector<float>& vdLower, const std::vector<float>& vdMiddle, const std::vector<float>& vdUpper, const std::vector<float>& vdExcitation)
 {
-	//solving bandmatrix with "numerical recipees in C" page 51
+	//solving bandmatrix with Thomas Algorithm. Implementation according to "numerical recipees in C", William H. Press -2nd edition 1992, ISBN 0-521-43108-5, page 51
 
 	if (vdMiddle[0] == 0)
 		ITA_EXCEPT_INVALID_PARAMETER("Matrix needs to be semi positiv definit.");
 	if (vdMiddle.size() != vdExcitation.size() || vdLower.size() != vdUpper.size() || vdExcitation.size() != vdLower.size() + 1)
 		ITA_EXCEPT_INVALID_PARAMETER("Dimensions of arguments missmatch.");
 
-	int n = vdMiddle.size();
+	const int iNum = vdMiddle.size();
 
-	std::vector<float> vdGamma(n);
+	std::vector<float> vdGamma(iNum);
 	float fBeta = vdMiddle[0];
-	std::vector<float> vdB(n); //solution of bandmatrix
+	std::vector<float> vdSolution(iNum); //
 
-	vdB[0] = vdExcitation[0] / vdMiddle[0];
-	for (int j = 1; j < n; j++) {
-		vdGamma[j] = vdUpper[j - 1] / fBeta;
-		fBeta = vdMiddle[j] - vdLower[j - 1] * vdGamma[j];
-		vdB[j] = (vdExcitation[j] - vdLower[j - 1] * vdB[j-1]) / fBeta;
+	vdSolution[0] = vdExcitation[0] / vdMiddle[0];
+	for (int idx = 1; idx < iNum; idx++) {
+		vdGamma[idx] = vdUpper[idx - 1] / fBeta;
+		fBeta = vdMiddle[idx] - vdLower[idx - 1] * vdGamma[idx];
+		vdSolution[idx] = (vdExcitation[idx] - vdLower[idx - 1] * vdSolution[idx-1]) / fBeta;
 	}
-	for (int j = n-2; j >= 0; j--) {//backsubstitution
-		vdB[j] -= vdGamma[j+1] * vdB[j + 1];
+	for (int idx = iNum-2; idx >= 0; idx--) {//backsubstitution
+		vdSolution[idx] -= vdGamma[idx+1] * vdSolution[idx + 1];
 	}
 	
-return vdB;
+return vdSolution;
 }

src/ITABase/Math/PiecewisePolynomial.cpp

@@ -16,7 +16,7 @@ CPiecewisePolynomial::CPiecewisePolynomial(const std::vector<float>& vdBreakPoin
 	if (iOrder < 0)
 		ITA_EXCEPT_INVALID_PARAMETER("Polynom order must be greater or equal zero.");
 
-	const int iMinimumNBreakPoints = std::max(2, iOrder + 1); // At least one interval
+	const int iMinimumNBreakPoints = 2; // At least one interval
 	if ( vdBreakPoints.size() < iMinimumNBreakPoints )
 		ITA_EXCEPT_INVALID_PARAMETER("Number of break points not sufficient for polynom order: Expected " + std::to_string(iMinimumNBreakPoints) + " but got " + std::to_string(vdBreakPoints.size()) + ".");
 
@@ -26,26 +26,25 @@ CPiecewisePolynomial::CPiecewisePolynomial(const std::vector<float>& vdBreakPoin
 
 }
 
-float CPiecewisePolynomial::Value(const float xValue) const
+float CPiecewisePolynomial::Value(const float & xValue) const
 {
-	// Check if x-value is inside boundaries
-	// Find index of interval
+	//a(x-x0)^(n) + b(x-x0)^(n-1) + c(x-x0)^(n-2) +...
 	int iInterval = IntervalIndex(xValue);
-	//TODO: get corresponding coefficients and calculate values
 	std::vector<float>::const_iterator itCoef = CoefficientsOfInterval(iInterval);
-	//f(x) = a1(x−x1)^3 + b1(x−x1)^2 + c1(x−x1) + d1
-	float x0 = vdBreakPoints[iInterval];
-	if (iOrder == 3)
-		return itCoef[0] * pow(xValue - x0, 3) + itCoef[1] * pow(xValue - x0, 2) + itCoef[2] * (xValue - x0) + itCoef[3];
-	else
-		ITA_EXCEPT_INVALID_PARAMETER("currently only splines of order 3");
+	const float x0 = vdBreakPoints[iInterval];
+	float fResult = 0;
+
+	for (int idx = 0; idx <= iOrder;  idx++) {
+		fResult += itCoef[idx] * pow(xValue - x0, iOrder - idx);
+	}
+	return fResult;
 }
 
-std::vector<float> ITABase::CPiecewisePolynomial::Value(const std::vector<float>& vdValues) const
+std::vector<float> ITABase::CPiecewisePolynomial::Value(const std::vector<float>& vdXValues) const
 {
-	std::vector<float> result(vdValues.size());
-	for (int i = 0; i < vdValues.size(); i++)
-		result[i] = Value(vdValues[i]);
+	std::vector<float> result(vdXValues.size());
+	for (int i = 0; i < vdXValues.size(); i++)
+		result[i] = Value(vdXValues[i]);
 	return result;
 }
 

tests/SplineTest.cpp

@@ -22,6 +22,7 @@
 #include <ITABase/Math/Spline.h>
 #include <ITAStopWatch.h>
 #include <iostream>
+#include <fstream>
 
 //using namespace std;
 bool MyCompare(std::vector<float> a, std::vector<float> b, float epsilon) { 
@@ -38,16 +39,17 @@ bool MyCompare(std::vector<float> a, std::vector<float> b, float epsilon) {
 
 int main(int iNumInArgs, char* pcInArgs[])
 {
-	/*beispiel von TM interaktiv
+	/*
+	//beispiel von TM interaktiv, look if it's working in general
 	std::vector<float> vdSupportingPoints{1,7,12,15,19};
 	std::vector<float> vdDataPoints = {8,10,7,8,7};
 	int iPolynomialOrder = 3;
 	IW_ITA_BASE_SPLINE::ITABase::CPiecewisePolynomial myPoly1 = IW_ITA_BASE_SPLINE::ITABase::Spline(vdSupportingPoints,vdDataPoints,iPolynomialOrder);
 	if MyCompare(std::cout << "Quick and dirty test: "<< myPoly1.Value(5)<<" -hier soll grob 10 rauskommen;)";
-	
+
 	float epsilon = 0.1;
 	//Matlab spline. DataPoints sind aus randi
-	
+
 	std::vector<float> vdDataPoints = { 13,14,2,14,10,2,5,9,15,15,3,15,11,4,0,-1,-5,-15,-11,-3,2 };
 	std::vector<float> vdX = { -10,  float(1.5),  float(-7.898),3, float(3.3), float(9.9),-2 };
 	std::vector<float> vdResult = { float(13.0000),  float(17.4325),    float(2.4282),    float(8.0000),    float(9.4143),   float(11.4819),    float(15.0000) };
@@ -56,9 +58,11 @@ int main(int iNumInArgs, char* pcInArgs[])
 	std::cout << "MATLAB SPLINE TEST 1" << std::endl;
 	if (MyCompare(vdResult, myPoly1.Value(vdX),epsilon))
 		std::cout << "Test 1 bestanden" << std::endl;
-	
+
 	std::vector<float> vdResult2 = { float(77.0000),    float(7.9485),   float(18.7512),   float(50.0000),   float(70.7647),   float(36.3051),   float(28.0000)
-	*///Test 2
+	*/
+	//Test 2, Speedtest
+	/*
 	std::vector<float> vdSupportingPoints{ -10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 };
 	std::vector<float> vdDataPoints2{ 77,    80,    19,    49,    45,    65,    71,    76,    28,    68,    66,    17,12,    50,    96,    35,    59,    23,    76,    26,    51 };
 	std::vector<float> vdX = { -10,  float(1.5),  float(-7.898),float(3.3), float(3.3), float(9.9),float(-2.5),float(8.9),float(0.01),float(-4.3) };
@@ -66,7 +70,7 @@ int main(int iNumInArgs, char* pcInArgs[])
 	std::vector<float> result(10);
 	ITAStopWatch myWatch = ITAStopWatch();
 	std::cout << "Rueckgabe von piecwisePolynomial Objekten" << std::endl;
-	
+
 	for (int k = 0; k < 100; k++) {
 		myWatch.start();
 		myPoly = ITABase::CubicSpline(vdSupportingPoints, vdDataPoints2);
@@ -75,12 +79,53 @@ int main(int iNumInArgs, char* pcInArgs[])
 	std::cout << myWatch << std::endl;
 	myWatch.reset();
 	std::cout << "Rueckgabe von vector Objekten" << std::endl;
-	
+
 	for (int k = 0; k < 100; k++) {
 		myWatch.start();
 		result = ITABase::CubicSpline(vdSupportingPoints, vdDataPoints2, vdX);
 		myWatch.stop();
 	}
 	std::cout << myWatch << std::endl;
+	*/
+
+	//Test 3, export Coefficients to .txt to plot them in Matlab
+	std::vector<float> vdSupportingPoints{ -10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 };
+	std::vector<float> vdDataPoints{ 77,    80,    19,    49,    45,    65,    71,    76,    28,    68,    66,    17,12,    50,    96,    35,    59,    23,    76,    26,    51 };
+	std::vector<float> vdX = { -10,  float(1.5),  float(-7.898),float(3.3), float(3.3), float(9.9),float(-2.5),float(8.9),float(0.01),float(-4.3) };
+	IW_ITA_BASE_SPLINE::ITABase::CPiecewisePolynomial myPoly = ITABase::CubicSpline(vdSupportingPoints, vdDataPoints);
+
+	std::vector<float>::const_iterator itCoeff = myPoly.getCoefficients();//Access to private property...
+	std::ofstream myFile;
+	myFile.open("points.txt");
+	for (int k = 0; k < vdDataPoints.size(); k++) {
+		myFile << vdSupportingPoints[k] << ", " << vdDataPoints[k] << std::endl;
+	}
+	myFile.close();
+
+	myFile.open("coeff.txt");
+	for (int k = 0; k < vdDataPoints.size() - 1; k++) {
+		myFile << itCoeff[4 * k] << ", " << itCoeff[4 * k + 1] << ", " << itCoeff[4 * k + 2] << ", " << itCoeff[4 * k + 3] << std::endl;
+	}
+	myFile.close();
+	//Test 4 last check up
+
+	//should equal the SupportingPoint
+	std::cout << "At supportion Point 0, the calculated value should be 66: " << myPoly.Value(0) << std::endl;
+	//should get specific Value
+	std::cout << "At random Point x=0.5 the spline should interpolate ~40.9: " << myPoly.Value(0.5) << std::endl;
+	//deviation test
+	std::vector<float> vdBreakPoints = { 0, 1 };
+	std::vector<float> vdCoeff = { 1,1,1};
+	int iOrder = 2;
+	IW_ITA_BASE_SPLINE::ITABase::CPiecewisePolynomial myPoly2 = ITABase::CPiecewisePolynomial(vdBreakPoints,vdCoeff,iOrder);
+	myPoly2 = myPoly2.Derivation();
+	std::cout << "The deviation of x^2+x+1 should be: " << myPoly2.getCoefficients()[0] << " * x + " << myPoly2.getCoefficients()[1] << std::endl;
+	std::cout << "so the new polynoial only got " << myPoly2.NumCoefficients() << " coefficients." << std::endl;
+	//should trhow exception
+	myPoly2.Value(2);
+	
+
+	
+	
 	return 0;
 }

tests/coeff.txt

+-24.4868, 0, 27.4868, 77
+58.434, -73.4604, -45.9736, 80
+-54.2494, 101.842, -17.5923, 19
+33.5635, -60.9064, 23.343, 49
+-22.0045, 39.784, 2.22051, 45
+16.4545, -26.2295, 15.775, 65
+-30.8134, 23.1339, 12.6795, 71
+54.7992, -69.3063, -33.4929, 76
+-47.3834, 95.0913, -7.70791, 28
+4.7343, -47.0588, 40.3245, 68
+23.4462, -32.8559, -39.5902, 66
+-7.51898, 37.4826, -34.9636, 17
+5.62973, 14.9256, 17.4446, 12
+-50, 31.8148, 64.1851, 50
+79.3701, -118.185, -22.1851, 96
+-75.4804, 119.925, -20.4449, 35
+77.5513, -106.516, -7.03549, 59
+-85.725, 126.138, 12.5868, 23
+73.3487, -131.037, 7.68815, 76
+-29.6697, 89.0092, -34.3395, 26

tests/points.txt

+-10, 77
+-9, 80
+-8, 19
+-7, 49
+-6, 45
+-5, 65
+-4, 71
+-3, 76
+-2, 28
+-1, 68
+0, 66
+1, 17
+2, 12
+3, 50
+4, 96
+5, 35
+6, 59
+7, 23
+8, 76
+9, 26
+10, 51

tests/splineTest.m

+%Plots the desired Spline interpolation for the DataPoint Values at the
+%Supporting Point location indicated in the provided .txt file and compares
+%the MATLAB spline against the spline created with C++ by reading another
+%.txt file
+%
+%Expected file Formats:
+%   coeff.txt
+%       a1, b1, c1, d1
+%       a2, b2, c2, d2
+%       a3, b3, c3, d3
+%       ...,...,...,...
+%   points.txt
+%       sup1, data1
+%       sup2, data2
+%       sup3, data3
+%       ... , ...
+% The SupportingPoints need to be sorted ascending!
+
+%% read Data from files
+FID = fopen('points.txt');
+if FID == -1
+    disp('unable to read indicated file');
+    return
+end
+C = textscan(FID,'%f %f','Delimiter',',');
+supPoints = C{1};
+dataPoints = C{2};
+fclose(FID);
+
+%Coeffs
+FID = fopen('coeff.txt');
+if FID == -1
+    disp('unable to read indicated file');
+    return
+end
+
+coeff = zeros(size(supPoints,1)-1,4);%alocate
+tline = fgetl(FID);
+k=1;
+while ischar(tline)
+    C = textscan(tline,'%f','Delimiter',',');
+    coeff(k,:) = C{1,1};
+    
+    tline = fgetl(FID);
+    k = k+1;
+end
+fclose(FID);
+
+%% plot
+myPoly = mkpp(supPoints,coeff);
+
+grain = (supPoints(end)-supPoints(1))/(10*size(supPoints,1)); %how dense shall the plotted points be?
+xx = [supPoints(1) : grain : supPoints(end)];
+
+plot(supPoints,dataPoints,'o',xx,spline(supPoints,dataPoints,xx));
+hold on
+grid on
+title 'Splines: Matlab vs. ITA'
+plot(xx,ppval(myPoly,xx));
+plot(xx,spline(supPoints,dataPoints,xx)-ppval(myPoly,xx));
+legend('Data','Matlab','ITA','Diff');
+hold off