# Demo to show how to use leasqr
#  Written by Doug
# June 2 2015
# First we will generate some data from a given formula.
# The formula is chosen as:
# y= 2*x(1)^2 +3*x(1) -4*x(2) +1*x(3)-2*x(3)^2-.7*x(4)*x(2)+5*x(5)+33
# Now setup the points where we measure the system 31 times
# this is a (31,5) matrix called mp (measurement parameters)


# This demo helped me to see how important the choice 
# of measurement points is.
# What we have is a surface in 5d space, and if
# we slice that surface with a 4d plain, or and other d plain
# then we get a line, and can curve fit to that line. Now if we 
# slice with any other plain we get a different line. We can get
# exact fitts to each of these lines, but each fit will have a
# different set of coefficients. What is needed is a system of
# orthogonal plain in 5d space and take the points on enough plains
# to span the 5d space. This is what the Taguchi method gives.
#
# For more on the Taguchi system see:
# https://controls.engin.umich.edu/wiki/index.php/Design_of_experiments_via_taguchi_methods:_orthogonal_arrays
# http://reliawiki.org/index.php/Taguchi_Orthogonal_Arrays#Three_Level_Designs
#
# I have used part of the L27 matrix for mp_taguchi
# and as can be seen this is the only one of my attempts
# that gives the correct answer.
# Thus with only 27 well chosen points and the leasqr function
# we can deduce the coefficients of the generating function.

 




mp0=[0:.1:3 ; 0:.5:15; 0:.2:6; 0:.3:9; 0:1:30]';

mp1=[0:.1:3 ; 15:-.5:0; 0:.2:6; 0:.3:9; 0:1:30]';
mp2=[0:.1:3 ; 0:.5:15; 6:-.2:0; 0:.3:9; 0:1:30]';
mp3=[0:.1:3 ; 0:.5:15; 0:.2:6; 9:-.3:0; 0:1:30]';

mp4=[0:.1:3 ; 0:.5:15; 0:.2:6; 0:.3:9; 30:-1:0]';
mp5=[0:.1:3 ; 15:-.5:0; 0:.2:6; 9:-.3:0; 0:1:30]';
mp6=[3:-.1:0 ; 0:.5:15; 0:.2:6; 0:.3:9; 0:1:30]';



mp_taguchi =[0,0,0,0,0
0,0,0,0,15
0,0,0,0,30
0,7,3,4,0
0,7,3,4,15
0,7,3,4,30
0,15,6,9,0
0,15,6,9,15
0,15,6,9,30
1.5,0,3,9,0
1.5,0,3,9,15
1.5,0,3,9,30
1.5,7,6,0,0
1.5,7,6,0,15
1.5,7,6,0,30
1.5,15,0,4,0
1.5,15,0,4,15
1.5,15,0,4,30
3,0,6,4,0
3,0,6,4,15
3,0,6,4,30
3,7,0,9,0
3,7,0,9,15
3,7,0,9,30
3,15,3,0,0
3,15,3,0,15
3,15,3,0,30 ];

# Here we can choose which set of experiment points to use 
#  You can try other combinations.

# mp=[mp0];
# mp=[mp0;mp1;mp2];
# mp=[mp0;mp1;mp2;mp3;mp4;mp5;mp6];

mp=mp_taguchi


# now make (create) the measured data using the chosen experiment points.

md= 2*mp(:,1).^2 +3*mp(:,1) -4*mp(:,2) +1*mp(:,3)-2*mp(:,3).^2-.7*mp(:,4).*mp(:,2)+5*mp(:,5)+33;

plot(md)

# now see if we can find the coefficients of the original function using leasqr
# I use the names used in the help section of leasqr

x=mp;
y=md;
pin=[1,1,1,1,1,.1,1]'

F=@(x, p)  p(1)*x(:,1).^2 +p(2)*x(:,1) +p(3)*x(:,2) +p(4)*mp(:,3) ...
+p(5)*mp(:,3).^2+p(6)*mp(:,4).*mp(:,2)+p(7)*mp(:,5)+33;

 [f1, p1, kvg1, iter1, corp1, covp1, covr1, stdresid1, Z1, r21] = ...
    leasqr (x,y, pin, F);

# p1 are the coefficients calculated by leasqr
p1
 
 
#  display the original coefficients for comparing
p_original=[2,3,-4,1,-2,-.7,5]'

# plot the calculated point on the original plot --- to see how well they fit.
hold on
plot(f1,'x')  
hold off