with GENERIC_REAL_LINEAR_EQUATIONS ; generic type REAL is digits <> ; package GENERIC_REAL_LEAST_SQUARE_FIT is package REAL_LINEAR_EQUATIONS is new GENERIC_REAL_LINEAR_EQUATIONS ( REAL ) ; use REAL_LINEAR_EQUATIONS ; use REAL_POLYNOMIAL_ARITHMETIC ; use REAL_MATRIX_ARITHMETIC ; -- The purpose of this package is to provide a reliable and convenient -- means for fitting existing data by a few coefficients. The companion -- package check_fit provides the means to use the coefficients for -- interpolation and limited extrapolation. -- -- This package implements the least square fit. -- -- The problem is stated as follows : -- Given measured data for values of Y based on values of X1,X2 and X3. e.g. -- -- Y_actual X1 X2 X3 -- -------- ----- ----- ----- -- 32.5 1.0 2.5 3.7 -- 7.2 2.0 2.5 3.6 -- 6.9 3.0 2.7 3.5 -- 22.4 2.2 2.1 3.1 -- 10.4 1.5 2.0 2.6 -- 11.3 1.6 2.0 3.1 -- -- Find a, b and c such that Y_approximate = a * X1 + b * X2 + c * X3 -- and such that the sum of (Y_actual - Y_approximate) squared is minimized. -- -- The method for determining the coefficients a, b and c follows directly -- form the problem definition and mathematical analysis. (See more below) -- -- Y is called the dependent variable and X1 .. Xn the independent variables. -- The procedures below implements a few special cases and the general case. -- The number of independent variables can vary. -- The approximation equation may use powers of the independent variables -- The user may create additional independent variables e.g. X2 = SIN(X1) -- with the restriction that the independent variables are linearly -- independent. e.g. Xi not equal p Xj + q for all i,j,p,q -- -- -- -- The mathematical derivation of the least square fit is as follows : -- -- Given data for the independent variable Y in terms of the dependent -- variables S,T,U and V consider that there exists a function F -- such that Y = F(S,T,U,V) -- The problem is to find coefficients a,b,c and d such that -- Y_approximate = a * S + b * T + c * U + d * V -- and such that the sum of ( Y - Y_approximate ) squared is minimized. -- -- Note: a, b, c, d are scalars. S, T, U, V, Y, Y_approximate are vectors. -- -- To find the minimum of SUM( Y - Y_approximate ) ** 2 -- the derivatives must be taken with respect to a,b,c and d and -- all must equal zero simultaneously. The steps follow : -- -- SUM( Y - Y_approximate ) ** 2 = SUM( Y - a*S - b*T - c*U - d*V ) ** 2 -- -- d/da = -2 * S * SUM( Y - A*S - B*T - C*U - D*V ) -- d/db = -2 * T * SUM( Y - A*S - B*T - C*U - D*V ) -- d/dc = -2 * U * SUM( Y - A*S - B*T - C*U - D*V ) -- d/dd = -2 * V * SUM( Y - A*S - B*T - C*U - D*V ) -- -- Setting each of the above equal to zero and putting constant term on left -- the -2 is factored out, -- the independent variable is moved inside the summation -- -- SUM( a*S*S + b*S*T + c*S*U + d*S*V = S*Y ) -- SUM( a*T*S + b*T*T + c*T*U + d*T*V = T*Y ) -- SUM( a*U*S + b*U*T + c*U*U + d*U*V = U*Y ) -- SUM( a*V*S + b*V*T + c*V*U + d*V*V = V*Y ) -- -- Distributing the SUM inside yields -- -- a * SUM(S*S) + b * SUM(S*T) + c * SUM(S*U) + d * SUM(S*V) = SUM(S*Y) -- a * SUM(T*S) + b * SUM(T*T) + c * SUM(T*U) + d * SUM(T*V) = SUM(T*Y) -- a * SUM(U*S) + b * SUM(U*T) + c * SUM(U*U) + d * SUM(U*V) = SUM(U*Y) -- a * SUM(V*S) + b * SUM(V*T) + c * SUM(V*U) + d * SUM(V*V) = SUM(V*Y) -- -- To find the coefficients a,b,c and d solve the linear system of equations -- -- | SUM(S*S) SUM(S*T) SUM(S*U) SUM(S*V) | | a | | SUM(S*Y) | -- | SUM(T*S) SUM(T*T) SUM(T*U) SUM(T*V) | x | b | = | SUM(T*Y) | -- | SUM(U*S) SUM(U*T) SUM(U*U) SUM(U*V) | | c | | SUM(U*Y) | -- | SUM(V*S) SUM(V*T) SUM(V*U) SUM(V*V) | | d | | SUM(V*Y) | -- -- Some observations : -- S,T,U and V must be linearly independent. -- There must be more data sets (Y, S, T, U, V) than variables. -- The analysis did not depend on the number of independent variables -- A polynomial fit results from the substitutions S=1, T=X, U=X**2, V=X**3 -- The general case for any order polynomial follows, fit_pn. -- Any substitution such as three variables to various powers may be used. -- -- Y = C * X procedure FIT ( Y , X : REAL_VECTOR ; C : out REAL ) ; -- Y = C1 * X1 + C2 * X2 procedure FIT ( Y , X1 , X2 : REAL_VECTOR ; C1 , C2 : out REAL ) ; -- Y = C1 * X1 + C2 * X2 + C3 * X3 procedure FIT ( Y , X1 , X2 , X3 : REAL_VECTOR ; C1 , C2 , C3 : out REAL ) ; -- Y = C1 * X1 + C2 * X2 + C3 * X3 + C4 * X4 procedure FIT ( Y , X1 , X2 , X3 , X4 : REAL_VECTOR ; C1 , C2 , C3 , C4 : out REAL ) ; -- Y = EVALUATE ( C , X ) -- Y = C(0) + C(1)*X + C(2)*X**2 + C(3)*X**3 + ... -- to C(C'LAST)*X**C'LAST procedure FIT ( Y , X : REAL_VECTOR ; C : out REAL_POLYNOMIAL ) ; -- Y = C11 * X1 + C12 * X1**2 + C21 * X2 + C22 * X2**2 procedure FIT ( Y , X1 , X2 : REAL_VECTOR ; C11 , C12 , C21 , C22 : out REAL ) ; -- Y = C11 * X1 + C12 * X1**2 + -- C21 * X2 + C22 * X2**2 + -- C31 * X3 + C32 * X3**2 procedure FIT ( Y , X1 , X2 , X3 : REAL_VECTOR ; C11 , C12 , C21 , C22 , C31 , C32 : out REAL ) ; -- Y = C11 * X1 + C12 * X1**2 + C13 * X1**3 + -- C21 * X2 + C22 * X2**2 + C23 * X2**3 procedure FIT ( Y , X1 , X2 : REAL_VECTOR ; C11 , C12 , C13 , C21 , C22 , C23 : out REAL ) ; -- Y = C(1,1)*X1 + C(1,2)*X1**2 + C(1,3)*X1**3 + ... -- to C(1,C'LAST)*X1**C'LAST + ... -- C(2,1)*X2 + C(2,2)*X2**2 + C(2,3)*X2**3 + ... -- to C(2,C'LAST)*X2**C'LAST procedure FIT ( Y , X1 , X2 : REAL_VECTOR ; C : in out REAL_MATRIX ) ; -- Y = C(1,1)*X1 + C(1,2)*X1**2 + C(1,3)*X1**3 + ... -- to C(1,C'LAST)*X1**C'LAST + ... -- C(2,1)*X2 + C(2,2)*X2**2 + C(2,3)*X2**3 + ... -- to C(2,C'LAST)*X2**C'LAST + ... -- C(3,1)*X3 + C(3,2)*X3**2 + C(3,3)*X3**3 + ... -- to C(3,C'LAST)*X3**C'LAST procedure FIT ( Y , X1 , X2 , X3 : REAL_VECTOR ; C : in out REAL_MATRIX ) ; -- Y = C(1,1)*X1 + C(1,2)*X1**2 + C(1,3)*X1**3 + ... -- to C(1,C'LAST)*X1**C'LAST + ... -- C(2,1)*X2 + C(2,2)*X2**2 + C(2,3)*X2**3 + ... -- to C(2,C'LAST)*X2**C'LAST + ... -- C(3,1)*X3 + C(3,2)*X3**2 + C(3,3)*X3**3 + ... -- to C(3,C'LAST)*X3**C'LAST + ... -- C(4,1)*X4 + C(4,2)*X4**2 + C(4,3)*X3**3 + ... -- to C(4,C'LAST)*X4**C'LAST procedure FIT ( Y , X1 , X2 , X3 , X4 : REAL_VECTOR ; C : in out REAL_MATRIX ) ; -- Y(K) = C(I,J) * X(I,K) ** J summed over I in C'RANGE(1) -- over J in C'RANGE(2) procedure FIT ( Y : REAL_VECTOR ; X : REAL_MATRIX ; C : in out REAL_MATRIX ) ; end GENERIC_REAL_LEAST_SQUARE_FIT ; -- with TEXT_IO ; use TEXT_IO ; -- only needed for testing package body GENERIC_REAL_LEAST_SQUARE_FIT is -- -- The mathmatical derivation of the least square fit is as follows : -- -- Given data for the independent variable Y in terms of the dependent -- variables S,T,U and V consider that there exists a function F -- such that Y = F(S,T,U,V) -- The problem is to find coefficients A,B,C and D such that -- Y_approximate = A * S + B * T + C * U + D * V -- and such that the sum of ( Y - Y_approximate ) squared is minimized. -- -- To find the minimum of SUM( Y - Y_approximate ) ** 2 -- the derivitives must be taken with respect to S,T,U and V and -- all must equal zero simultaneously. The steps follow : -- -- SUM( Y - Y_approximate ) ** 2 = SUM( Y - A*S - B*T - C*U - D*V ) ** 2 -- -- d/dS = -2 * S * SUM( Y - A*S - B*T - C*U - D*V ) -- d/dT = -2 * T * SUM( Y - A*S - B*T - C*U - D*V ) -- d/dU = -2 * U * SUM( Y - A*S - B*T - C*U - D*V ) -- d/dV = -2 * V * SUM( Y - A*S - B*T - C*U - D*V ) -- -- Setting each of the above equal to zero and putting constant term on left -- the -2 is factored out, -- the independent variable is moved inside the summation -- -- SUM( A*S*S + B*S*T + C*S*U + D*S*V = S*Y ) -- SUM( A*T*S + B*T*T + C*T*U + D*T*V = T*Y ) -- SUM( A*U*S + B*U*T + C*U*U + D*U*V = U*Y ) -- SUM( A*V*S + B*V*T + C*V*U + D*V*V = V*Y ) -- -- Distributing the SUM inside yeilds -- -- A * SUM(S*S) + B * SUM(S*T) + C * SUM(S*U) + D * SUM(S*V) = SUM(S*Y) -- A * SUM(T*S) + B * SUM(T*T) + C * SUM(T*U) + D * SUM(T*V) = SUM(T*Y) -- A * SUM(U*S) + B * SUM(U*T) + C * SUM(U*U) + D * SUM(U*V) = SUM(U*Y) -- A * SUM(V*S) + B * SUM(V*T) + C * SUM(V*U) + D * SUM(V*V) = SUM(V*Y) -- -- To find the coefficients A,B,C and D solve the linear system of equations -- -- | SUM(S*S) SUM(S*T) SUM(S*U) SUM(S*V) | | A | | SUM(S*Y) | -- | SUM(T*S) SUM(T*T) SUM(T*U) SUM(T*V) | x | B | = | SUM(T*Y) | -- | SUM(U*S) SUM(U*T) SUM(U*U) SUM(U*V) | | C | | SUM(U*Y) | -- | SUM(V*S) SUM(V*T) SUM(V*U) SUM(V*V) | | D | | SUM(V*Y) | -- -- Some observations : -- S,T,U and V must be linearly independent -- The analysis did not depend on the number of independent variables -- A polynomial fit results from the substitutions S=1, T=X, U=X**2, V=X**3 -- The general case for any order polynomial follows. -- Any substitiution such as just using odd powers of X follow. -- The "Y" vector can be expanded to a matrix yeilding a matrix -- of coefficients -- -- Y = C * X procedure FIT ( Y , X : REAL_VECTOR ; C : out REAL ) is SUM_X_X : REAL := 0.0 ; SUM_X_Y : REAL := 0.0 ; begin if X'LENGTH /= Y'LENGTH then raise CONSTRAINT_ERROR ; end if ; if X'LENGTH < 1 then raise CONSTRAINT_ERROR ; end if ; for I in X'RANGE loop SUM_X_X := SUM_X_X + X ( I ) * X ( I ) ; SUM_X_Y := SUM_X_Y + X ( I ) * Y ( I - X'FIRST + Y'FIRST ) ; end loop ; -- No need to test abs(SUM_X_X) = 0.0 and raise numeric error, it happens C := SUM_X_Y / SUM_X_X ; end FIT ; -- Y = C1 * X1 + C2 * X2 procedure FIT ( Y , X1 , X2 : REAL_VECTOR ; C1 , C2 : out REAL ) is M : REAL_MATRIX ( 0 .. 1 , 0 .. 1 ) := ( 0 .. 1 =>( 0 .. 1 => 0.0 )) ; V : REAL_VECTOR ( 0 .. 1 ) := ( 0 .. 1 => 0.0 ) ; C : REAL_VECTOR ( 0 .. 1 ) ; BIAS_X1 : constant INTEGER := Y'FIRST - X1'FIRST ; BIAS_X2 : constant INTEGER := Y'FIRST - X2'FIRST ; begin if X1'LENGTH /= X2'LENGTH or X2'LENGTH /= Y'LENGTH then raise CONSTRAINT_ERROR ; end if ; if Y'LENGTH < 2 then raise CONSTRAINT_ERROR ; end if ; for I in Y'RANGE loop M ( 0 , 0 ) := M ( 0 , 0 ) + X1 ( I - BIAS_X1 ) * X1 ( I - BIAS_X1 ) ; M ( 0 , 1 ) := M ( 0 , 1 ) + X1 ( I - BIAS_X1 ) * X2 ( I - BIAS_X2 ) ; M ( 1 , 0 ) := M ( 1 , 0 ) + X2 ( I - BIAS_X2 ) * X1 ( I - BIAS_X1 ) ; M ( 1 , 1 ) := M ( 1 , 1 ) + X2 ( I - BIAS_X2 ) * X2 ( I - BIAS_X2 ) ; V ( 0 ) := V ( 0 ) + X1 ( I - BIAS_X1 ) * Y ( I ) ; V ( 1 ) := V ( 1 ) + X2 ( I - BIAS_X2 ) * Y ( I ) ; end loop ; C := LINEAR_EQUATIONS ( M , V ) ; C1 := C ( 0 ) ; C2 := C ( 1 ) ; end FIT ; -- Y = C1 * X1 + C2 * X2 + C3 * X3 procedure FIT ( Y , X1 , X2 , X3 : REAL_VECTOR ; C1 , C2 , C3 : out REAL ) is M : REAL_MATRIX ( 0 .. 2 , 0 .. 2 ) := ( 0 .. 2 =>( 0 .. 2 => 0.0 )) ; V : REAL_VECTOR ( 0 .. 2 ) := ( 0 .. 2 => 0.0 ) ; C : REAL_VECTOR ( 0 .. 2 ) ; BIAS_X1 : constant INTEGER := Y'FIRST - X1'FIRST ; BIAS_X2 : constant INTEGER := Y'FIRST - X2'FIRST ; BIAS_X3 : constant INTEGER := Y'FIRST - X3'FIRST ; begin if X1'LENGTH /= X2'LENGTH or X2'LENGTH /= X3'LENGTH or X3'LENGTH /= Y'LENGTH then raise CONSTRAINT_ERROR ; end if ; if Y'LENGTH < 3 then raise CONSTRAINT_ERROR ; end if ; for I in Y'RANGE loop M ( 0 , 0 ) := M ( 0 , 0 ) + X1 ( I - BIAS_X1 ) * X1 ( I - BIAS_X1 ) ; M ( 0 , 1 ) := M ( 0 , 1 ) + X1 ( I - BIAS_X1 ) * X2 ( I - BIAS_X2 ) ; M ( 0 , 2 ) := M ( 0 , 2 ) + X1 ( I - BIAS_X1 ) * X3 ( I - BIAS_X3 ) ; M ( 1 , 0 ) := M ( 1 , 0 ) + X2 ( I - BIAS_X2 ) * X1 ( I - BIAS_X1 ) ; M ( 1 , 1 ) := M ( 1 , 1 ) + X2 ( I - BIAS_X2 ) * X2 ( I - BIAS_X2 ) ; M ( 1 , 2 ) := M ( 1 , 2 ) + X2 ( I - BIAS_X2 ) * X3 ( I - BIAS_X3 ) ; M ( 2 , 0 ) := M ( 2 , 0 ) + X3 ( I - BIAS_X3 ) * X1 ( I - BIAS_X1 ) ; M ( 2 , 1 ) := M ( 2 , 1 ) + X3 ( I - BIAS_X3 ) * X2 ( I - BIAS_X2 ) ; M ( 2 , 2 ) := M ( 2 , 2 ) + X3 ( I - BIAS_X3 ) * X3 ( I - BIAS_X3 ) ; V ( 0 ) := V ( 0 ) + X1 ( I - BIAS_X1 ) * Y ( I ) ; V ( 1 ) := V ( 1 ) + X2 ( I - BIAS_X2 ) * Y ( I ) ; V ( 2 ) := V ( 2 ) + X3 ( I - BIAS_X3 ) * Y ( I ) ; end loop ; C := LINEAR_EQUATIONS ( M , V ) ; C1 := C ( 0 ) ; C2 := C ( 1 ) ; C3 := C ( 2 ) ; end FIT ; -- Y = C1 * X1 + C2 * X2 + C3 * X3 + C4 * X4 procedure FIT ( Y , X1 , X2 , X3 , X4 : REAL_VECTOR ; C1 , C2 , C3 , C4 : out REAL ) is N : constant INTEGER := 3 ; M : REAL_MATRIX ( 0 .. N , 0 .. N ) := ( 0 .. N =>( 0 .. N => 0.0 )) ; V : REAL_VECTOR ( 0 .. N ) := ( 0 .. N => 0.0 ) ; C : REAL_VECTOR ( 0 .. N ) ; X_POWER : REAL_VECTOR ( 0 .. N ) ; BIAS_X1 : constant INTEGER := Y'FIRST - X1'FIRST ; BIAS_X2 : constant INTEGER := Y'FIRST - X2'FIRST ; BIAS_X3 : constant INTEGER := Y'FIRST - X3'FIRST ; BIAS_X4 : constant INTEGER := Y'FIRST - X4'FIRST ; begin if X1'LENGTH /= X2'LENGTH or X2'LENGTH /= X3'LENGTH or X3'LENGTH /= X4'LENGTH or X4'LENGTH /= Y'LENGTH then raise CONSTRAINT_ERROR ; end if ; if Y'LENGTH < 4 then raise CONSTRAINT_ERROR ; end if ; for K in Y'RANGE loop X_POWER ( 0 ) := X1 ( K - BIAS_X1 ) ; X_POWER ( 1 ) := X2 ( K - BIAS_X2 ) ; X_POWER ( 2 ) := X3 ( K - BIAS_X3 ) ; X_POWER ( 3 ) := X4 ( K - BIAS_X4 ) ; for I in 0 .. N loop for J in 0 .. N loop M ( I , J ) := M ( I , J ) + X_POWER ( I ) * X_POWER ( J ) ; end loop ; V ( I ) := V ( I ) + Y ( K ) * X_POWER ( I ) ; end loop ; end loop ; C := LINEAR_EQUATIONS ( M , V ) ; C1 := C ( 0 ) ; C2 := C ( 1 ) ; C3 := C ( 2 ) ; C4 := C ( 3 ) ; end FIT ; -- Y = EVALUATE ( C , X ) -- Y = C(0) + C(1)*X + C(2)*X**2 + C(3)*X**3 + ... -- to C(C'LAST)*X**C'LAST procedure FIT ( Y , X : REAL_VECTOR ; C : out REAL_POLYNOMIAL ) is N : constant INTEGER := C'LENGTH - 1 ; M : REAL_MATRIX ( 0 .. N , 0 .. N ) := ( 0 .. N =>( 0 .. N => 0.0 )) ; V : REAL_VECTOR ( 0 .. N ) := ( 0 .. N => 0.0 ) ; X_POWER : REAL_VECTOR ( 0 .. 2 * N ) ; BIAS : constant INTEGER := X'FIRST - Y'FIRST ; begin if X'LENGTH /= Y'LENGTH then raise CONSTRAINT_ERROR ; end if ; if X'LENGTH < C'LENGTH then raise CONSTRAINT_ERROR ; end if ; X_POWER ( 0 ) := 1.0 ; for K in X'RANGE loop for I in 1 .. 2 * N loop X_POWER ( I ) := X_POWER ( I - 1 ) * X ( K ) ; end loop ; for I in 0 .. N loop for J in 0 .. N loop M ( I , J ) := M ( I , J ) + X_POWER ( I + J ) ; end loop ; V ( I ) := V ( I ) + Y ( K - BIAS ) * X_POWER ( I ) ; end loop ; end loop ; V := LINEAR_EQUATIONS ( M , V ) ; for I in V'RANGE loop C ( I ) := V ( I ) ; end loop ; end FIT ; -- Y = C11 * X1 + C12 * X1**2 + C21 * X2 + C22 * X2**2 procedure FIT ( Y , X1 , X2 : REAL_VECTOR ; C11 , C12 , C21 , C22 : out REAL ) is M : REAL_MATRIX ( 0 .. 3 , 0 .. 3 ) := ( 0 .. 3 =>( 0 .. 3 => 0.0 )) ; V : REAL_VECTOR ( 0 .. 3 ) := ( 0 .. 3 => 0.0 ) ; X_POWER : REAL_VECTOR ( 0 .. 3 ) ; C : REAL_VECTOR ( 0 .. 3 ) ; BIAS_X1 : constant INTEGER := Y'FIRST - X1'FIRST ; BIAS_X2 : constant INTEGER := Y'FIRST - X2'FIRST ; begin if X1'LENGTH /= X2'LENGTH or X2'LENGTH /= Y'LENGTH then raise CONSTRAINT_ERROR ; end if ; if Y'LENGTH < 4 then raise CONSTRAINT_ERROR ; end if ; for K in Y'RANGE loop X_POWER ( 0 ) := X1 ( K - BIAS_X1 ) ; X_POWER ( 1 ) := X_POWER ( 0 ) ** 2 ; X_POWER ( 2 ) := X2 ( K - BIAS_X2 ) ; X_POWER ( 3 ) := X_POWER ( 2 ) ** 2 ; for I in 0 .. 3 loop for J in 0 .. 3 loop M ( I , J ) := M ( I , J ) + X_POWER ( I ) * X_POWER ( J ) ; end loop ; V ( I ) := V ( I ) + X_POWER ( I ) * Y ( K ) ; end loop ; end loop ; C := LINEAR_EQUATIONS ( M , V ) ; C11 := C ( 0 ) ; C12 := C ( 1 ) ; C21 := C ( 2 ) ; C22 := C ( 3 ) ; end FIT ; -- Y = C11 * X1 + C12 * X1**2 + -- C21 * X2 + C22 * X2**2 + -- C31 * X3 + C32 * X3**2 procedure FIT ( Y , X1 , X2 , X3 : REAL_VECTOR ; C11 , C12 , C21 , C22 , C31 , C32 : out REAL ) is M : REAL_MATRIX ( 0 .. 5 , 0 .. 5 ) := ( 0 .. 5 =>( 0 .. 5 => 0.0 )) ; V : REAL_VECTOR ( 0 .. 5 ) := ( 0 .. 5 => 0.0 ) ; X_POWER : REAL_VECTOR ( 0 .. 5 ) ; C : REAL_VECTOR ( 0 .. 5 ) ; BIAS_X1 : constant INTEGER := Y'FIRST - X1'FIRST ; BIAS_X2 : constant INTEGER := Y'FIRST - X2'FIRST ; BIAS_X3 : constant INTEGER := Y'FIRST - X3'FIRST ; begin if X1'LENGTH /= X2'LENGTH or X2'LENGTH /= X3'LENGTH or X3'LENGTH /= Y'LENGTH then raise CONSTRAINT_ERROR ; end if ; if Y'LENGTH < 6 then raise CONSTRAINT_ERROR ; end if ; for K in Y'RANGE loop X_POWER ( 0 ) := X1 ( K - BIAS_X1 ) ; X_POWER ( 1 ) := X_POWER ( 0 ) ** 2 ; X_POWER ( 2 ) := X2 ( K - BIAS_X2 ) ; X_POWER ( 3 ) := X_POWER ( 2 ) ** 2 ; X_POWER ( 4 ) := X3 ( K - BIAS_X3 ) ; X_POWER ( 5 ) := X_POWER ( 4 ) ** 2 ; for I in M'RANGE(1) loop for J in M'RANGE(2) loop M ( I , J ) := M ( I , J ) + X_POWER ( I ) * X_POWER ( J ) ; end loop ; V ( I ) := V ( I ) + X_POWER ( I ) * Y ( K ) ; end loop ; end loop ; C := LINEAR_EQUATIONS ( M , V ) ; C11 := C ( 0 ) ; C12 := C ( 1 ) ; C21 := C ( 2 ) ; C22 := C ( 3 ) ; C31 := C ( 4 ) ; C32 := C ( 5 ) ; end FIT ; -- Y = C11 * X1 + C12 * X1**2 + C13 * X1**3 + -- C21 * X2 + C22 * X2**2 + C23 * X2**3 procedure FIT ( Y , X1 , X2 : REAL_VECTOR ; C11 , C12 , C13 , C21 , C22 , C23 : out REAL ) is M : REAL_MATRIX ( 0 .. 5 , 0 .. 5 ) := ( 0 .. 5 =>( 0 .. 5 => 0.0 )) ; V : REAL_VECTOR ( 0 .. 5 ) := ( 0 .. 5 => 0.0 ) ; X_POWER : REAL_VECTOR ( 0 .. 5 ) ; C : REAL_VECTOR ( 0 .. 5 ) ; BIAS_X1 : constant INTEGER := Y'FIRST - X1'FIRST ; BIAS_X2 : constant INTEGER := Y'FIRST - X2'FIRST ; begin if X1'LENGTH /= X2'LENGTH or X2'LENGTH /= Y'LENGTH then raise CONSTRAINT_ERROR ; end if ; if Y'LENGTH < 6 then raise CONSTRAINT_ERROR ; end if ; for K in Y'RANGE loop X_POWER ( 0 ) := X1 ( K - BIAS_X1 ) ; X_POWER ( 1 ) := X_POWER ( 0 ) ** 2 ; X_POWER ( 2 ) := X_POWER ( 0 ) ** 3 ; X_POWER ( 3 ) := X2 ( K - BIAS_X2 ) ; X_POWER ( 4 ) := X_POWER ( 3 ) ** 2 ; X_POWER ( 5 ) := X_POWER ( 3 ) ** 3 ; for I in M'RANGE(1) loop for J in M'RANGE(2) loop M ( I , J ) := M ( I , J ) + X_POWER ( I ) * X_POWER ( J ) ; end loop ; V ( I ) := V ( I ) + X_POWER ( I ) * Y ( K ) ; end loop ; end loop ; C := LINEAR_EQUATIONS ( M , V ) ; C11 := C ( 0 ) ; C12 := C ( 1 ) ; C13 := C ( 2 ) ; C21 := C ( 3 ) ; C22 := C ( 4 ) ; C23 := C ( 5 ) ; end FIT ; -- Y = C(1,1)*X1 + C(1,2)*X1**2 + C(1,3)*X1**3 + ... -- to C(1,C'LAST)*X1**C'LAST + ... -- C(2,1)*X2 + C(2,2)*X2**2 + C(2,3)*X2**3 + ... -- to C(2,C'LAST)*X2**C'LAST procedure FIT ( Y , X1 , X2 : REAL_VECTOR ; C : in out REAL_MATRIX ) is N : constant INTEGER := C'LENGTH ( 2 ) - 1 ; N2 : constant INTEGER := 2 * N + 1 ; M : REAL_MATRIX ( 0 .. N2 , 0 .. N2 ) := ( 0 .. N2 =>( 0 .. N2 => 0.0 )) ; V : REAL_VECTOR ( 0 .. N2 ) := ( 0 .. N2 => 0.0 ) ; C_VEC : REAL_VECTOR ( 0 .. N2 ) ; X_POWER : REAL_VECTOR ( 0 .. N2 ) ; BIAS_X1 : constant INTEGER := Y'FIRST - X1'FIRST ; BIAS_X2 : constant INTEGER := Y'FIRST - X2'FIRST ; begin if X1'LENGTH /= X2'LENGTH or X2'LENGTH /= Y'LENGTH then raise CONSTRAINT_ERROR ; end if ; if C'LENGTH ( 1 ) /= 2 then raise CONSTRAINT_ERROR ; end if ; if Y'LENGTH < 2 * C'LENGTH(2) then raise CONSTRAINT_ERROR ; end if ; for K in Y'RANGE loop X_POWER ( 0 ) := X1 ( K - BIAS_X1 ) ; X_POWER ( N + 1 ) := X2 ( K - BIAS_X2 ) ; for I in 1 .. N loop X_POWER ( I ) := X_POWER ( I - 1 ) * X_POWER ( 0 ) ; X_POWER ( I + N + 1 ) := X_POWER ( I + N ) * X_POWER ( N + 1 ) ; end loop ; for I in M'RANGE(1) loop for J in M'RANGE(2) loop M ( I , J ) := M ( I , J ) + X_POWER ( I ) * X_POWER ( J ) ; end loop ; V ( I ) := V ( I ) + X_POWER ( I ) * Y ( K ) ; end loop ; end loop ; C_VEC := LINEAR_EQUATIONS ( M , V ) ; VECTOR_TO_ROW ( C_VEC( 0 .. N ) , C'FIRST ( 1 ) , C) ; VECTOR_TO_ROW ( C_VEC( N + 1 .. N2 ) , C'LAST ( 1 ) , C) ; end FIT ; -- Y = C(1,1)*X1 + C(1,2)*X1**2 + C(1,3)*X1**3 + ... -- to C(1,C'LAST)*X1**C'LAST + ... -- C(2,1)*X2 + C(2,2)*X2**2 + C(2,3)*X2**3 + ... -- to C(2,C'LAST)*X2**C'LAST + ... -- C(3,1)*X3 + C(3,2)*X3**2 + C(3,3)*X3**3 + ... -- to C(3,C'LAST)*X3**C'LAST procedure FIT ( Y , X1 , X2 , X3 : REAL_VECTOR ; C : in out REAL_MATRIX ) is N : constant INTEGER := C'LENGTH ( 2 ) - 1 ; N2 : constant INTEGER := 3 * N + 2 ; M : REAL_MATRIX ( 0 .. N2 , 0 .. N2 ) := ( 0 .. N2 =>( 0 .. N2 => 0.0 )) ; V : REAL_VECTOR ( 0 .. N2 ) := ( 0 .. N2 => 0.0 ) ; C_VEC : REAL_VECTOR ( 0 .. N2 ) ; X_POWER : REAL_VECTOR ( 0 .. N2 ) ; BIAS_X1 : constant INTEGER := Y'FIRST - X1'FIRST ; BIAS_X2 : constant INTEGER := Y'FIRST - X2'FIRST ; BIAS_X3 : constant INTEGER := Y'FIRST - X3'FIRST ; begin if X1'LENGTH /= X2'LENGTH or X2'LENGTH /= X3'LENGTH or X3'LENGTH /= Y'LENGTH then raise CONSTRAINT_ERROR ; end if ; if C'LENGTH ( 1 ) /= 3 then raise CONSTRAINT_ERROR ; end if ; if Y'LENGTH < 3 * C'LENGTH(2) then raise CONSTRAINT_ERROR ; end if ; for K in Y'RANGE loop X_POWER ( 0 ) := X1 ( K - BIAS_X1 ) ; X_POWER ( N + 1 ) := X2 ( K - BIAS_X2 ) ; X_POWER ( 2 * N + 2 ) := X3 ( K - BIAS_X3 ) ; for I in 1 .. N loop X_POWER ( I ) := X_POWER ( I - 1 ) * X_POWER ( 0 ) ; X_POWER ( I + N + 1 ) := X_POWER ( I + N ) * X_POWER ( N + 1 ) ; X_POWER ( I + 2 * N + 2 ) := X_POWER ( I + 2 * N + 1 ) * X_POWER ( 2 * N + 2 ) ; end loop ; for I in 0 .. N2 loop for J in 0 .. N2 loop M ( I , J ) := M ( I , J ) + X_POWER ( I ) * X_POWER ( J ) ; end loop ; V ( I ) := V ( I ) + X_POWER ( I ) * Y ( K ) ; end loop ; end loop ; C_VEC := LINEAR_EQUATIONS ( M , V ) ; VECTOR_TO_ROW ( C_VEC( 0 .. N ) , C'FIRST ( 1 ) , C) ; VECTOR_TO_ROW ( C_VEC( N + 1 .. 2 * N + 1 ) , C'FIRST ( 1 ) + 1 , C) ; VECTOR_TO_ROW ( C_VEC( 2 * N + 2 .. N2 ) , C'FIRST ( 1 ) + 2 , C) ; end FIT ; -- Y = C(1,1)*X1 + C(1,2)*X1**2 + C(1,3)*X1**3 + ... -- to C(1,C'LAST)*X1**C'LAST + ... -- C(2,1)*X2 + C(2,2)*X2**2 + C(2,3)*X2**3 + ... -- to C(2,C'LAST)*X2**C'LAST + ... -- C(3,1)*X3 + C(3,2)*X3**2 + C(3,3)*X3**3 + ... -- to C(3,C'LAST)*X3**C'LAST + ... -- C(4,1)*X4 + C(4,2)*X4**2 + C(4,3)*X3**3 + ... -- to C(4,C'LAST)*X4**C'LAST procedure FIT ( Y , X1 , X2 , X3 , X4 : REAL_VECTOR ; C : in out REAL_MATRIX ) is N : constant INTEGER := C'LENGTH ( 2 ) - 1 ; N2 : constant INTEGER := 4 * N + 3 ; M : REAL_MATRIX ( 0 .. N2 , 0 .. N2 ) := ( 0 .. N2 =>( 0 .. N2 => 0.0 )) ; V : REAL_VECTOR ( 0 .. N2 ) := ( 0 .. N2 => 0.0 ) ; C_VEC : REAL_VECTOR ( 0 .. N2 ) ; X_POWER : REAL_VECTOR ( 0 .. N2 ) ; BIAS_X1 : constant INTEGER := Y'FIRST - X1'FIRST ; BIAS_X2 : constant INTEGER := Y'FIRST - X2'FIRST ; BIAS_X3 : constant INTEGER := Y'FIRST - X3'FIRST ; BIAS_X4 : constant INTEGER := Y'FIRST - X4'FIRST ; begin if X1'LENGTH /= X2'LENGTH or X2'LENGTH /= X3'LENGTH or X3'LENGTH /= X4'LENGTH or X4'LENGTH /= Y'LENGTH then raise CONSTRAINT_ERROR ; end if ; if C'LENGTH ( 1 ) /= 4 then raise CONSTRAINT_ERROR ; end if ; if Y'LENGTH < 4 * C'LENGTH(2) then raise CONSTRAINT_ERROR ; end if ; for K in Y'RANGE loop X_POWER ( 0 ) := X1 ( K - BIAS_X1 ) ; X_POWER ( N + 1 ) := X2 ( K - BIAS_X2 ) ; X_POWER ( 2 * N + 2 ) := X3 ( K - BIAS_X3 ) ; X_POWER ( 3 * N + 3 ) := X4 ( K - BIAS_X4 ) ; for I in 1 .. N loop X_POWER ( I ) := X_POWER ( I - 1 ) * X_POWER ( 0 ) ; X_POWER ( I + N + 1 ) := X_POWER ( I + N ) * X_POWER ( N + 1 ) ; X_POWER ( I + 2 * N + 2 ) := X_POWER ( I + 2 * N + 1 ) * X_POWER ( 2 * N + 2 ) ; X_POWER ( I + 3 * N + 3 ) := X_POWER ( I + 3 * N + 2 ) * X_POWER ( 3 * N + 3 ) ; end loop ; for I in 0 .. N2 loop for J in 0 .. N2 loop M ( I , J ) := M ( I , J ) + X_POWER ( I ) * X_POWER ( J ) ; end loop ; V ( I ) := V ( I ) + X_POWER ( I ) * Y ( K ) ; end loop ; end loop ; C_VEC := LINEAR_EQUATIONS ( M , V ) ; VECTOR_TO_ROW ( C_VEC( 0 .. N ) , C'FIRST ( 1 ) , C) ; VECTOR_TO_ROW ( C_VEC( N + 1 .. 2 * N + 1 ) , C'FIRST ( 1 ) + 1 , C) ; VECTOR_TO_ROW ( C_VEC( 2 * N + 2 .. 3 * N + 2 ) , C'FIRST ( 1 ) + 2 , C) ; VECTOR_TO_ROW ( C_VEC( 3 * N + 3 .. N2 ) , C'FIRST ( 1 ) + 3 , C) ; end FIT ; -- Y(K) = C(I,J) * X(I,K) ** J summed over I in C'RANGE(1) -- over J in C'RANGE(2) procedure FIT ( Y : REAL_VECTOR ; X : REAL_MATRIX ; C : in out REAL_MATRIX ) is N : constant INTEGER := C'LENGTH ( 1 ) * C'LENGTH ( 2 ) - 1 ; M : REAL_MATRIX ( 0 .. N , 0 .. N ) := ( 0 .. N =>( 0 .. N => 0.0 )) ; V : REAL_VECTOR ( 0 .. N ) := ( 0 .. N => 0.0 ) ; C_V : REAL_VECTOR ( 0 .. N ) ; X_POWER : REAL_VECTOR ( 0 .. N ) ; BIAS : constant INTEGER := X'FIRST ( 2 ) - Y'FIRST ; K1 : INTEGER := 0 ; begin if X'LENGTH ( 2 ) /= Y'LENGTH then raise CONSTRAINT_ERROR ; -- X and Y must have same number of data points end if ; if C'LENGTH ( 1 ) /= X'LENGTH ( 1 ) then raise CONSTRAINT_ERROR ; -- C and X must have same number of rows end if ; if Y'LENGTH < C'LENGTH ( 1 ) * C'LENGTH ( 2 ) then raise NUMERIC_ERROR ; -- insufficient data end if ; for K in X'RANGE ( 2 ) loop for I in INTEGER(0) .. C'LENGTH ( 1 ) - 1 loop X_POWER ( I * C'LENGTH( 2 )) := X ( I + X'FIRST , K ) ; for J in INTEGER(1) .. C'LENGTH ( 2 ) - 1 loop X_POWER ( I * C'LENGTH( 2 ) + J) := X_POWER ( I * C'LENGTH( 2 ) + J - 1) * X ( I + X'FIRST , K ) ; end loop ; end loop ; for I in M'RANGE(1) loop for J in M'RANGE(2) loop M ( I , J ) := M ( I , J ) + X_POWER ( I ) * X_POWER ( J ) ; end loop ; V ( I ) := V ( I ) + Y ( K - BIAS ) * X_POWER ( I ) ; end loop ; end loop ; C_V := LINEAR_EQUATIONS ( M , V ) ; for I in C'RANGE ( 1 ) loop VECTOR_TO_ROW ( C_V( K1 .. K1 + C'LENGTH( 2 ) - 1) , I , C) ; K1 := K1 + C'LENGTH ( 2 ) ; end loop ; end FIT ; -- end GENERIC_REAL_LEAST_SQUARE_FIT ;