-- ******************************************************************************************************************************
--
--                                                          L I N R E G
--
--  Program:      LINREG
--
--  Programmer:   Dr. David G. Simpson
--                Department of Physical Science
--                Prince George's Community College
--                Largo, Maryland  20774
--
--  Date:         January 22, 2002
--
--  Language:     Ada-83
--
--  Description:  This program performs a linear regression analysis for a set of data given as (x,y) pairs.  The output from
--                the program is the slope and y-intercept of the least-squares best fit straight line through the data points.
--
-- ******************************************************************************************************************************

with Text_IO; use Text_IO;                                                          -- use Text_IO package

procedure LinReg is

   package Real_IO is new Float_IO (Float); use Real_IO;                            -- instantiate generic Float_IO package




-- *****************************************************************************************************************************
--
--  Variable declarations
--
-- *****************************************************************************************************************************

   b     : Float;                                                                   -- y-intercpt of least-squares best fit line
   m     : Float;                                                                   -- slope of least-squares best fit line
   n     : Float := 0.0;                                                            -- number of data points
   r     : Float;                                                                   -- squared correlation coefficient
   sumx  : Float := 0.0;                                                            -- sum of x
   sumx2 : Float := 0.0;                                                            -- sum of x**2
   sumxy : Float := 0.0;                                                            -- sum of x * y
   sumy  : Float := 0.0;                                                            -- sum of y
   sumy2 : Float := 0.0;                                                            -- sum of y**2
   x     : Float;                                                                   -- input x data
   y     : Float;                                                                   -- input y data





-- ******************************************************************************************************************************
--
--  Sqrt()
--
--  Calculate the square root of a number using Newton's method.  (Ada-83 has no standard math package!)
--
-- ******************************************************************************************************************************

   function Sqrt (X : in Float) return Float is                                     -- square root function

      Xn     : Float;                                                               -- n-th interation esimate of square root
      Del    : Float;                                                               -- difference from previous iteration
      ConvLim : Float := 1.0e-10;                                                   -- convergence limit

   begin

      if X < 0.0 then                                                               -- if argument is less than zero..
         Put_Line ("Error in Sqrt.");                                               -- ..then imaginary result; print error msg
         return 0.0;                                                                -- ..and return a value of 0.0.
      end if;

      Xn := X/2.0;                                                                  -- initial estimate for square root is X/2
      for n in 1..100 loop                                                          -- loop a maximum of 100 times
         Del := (Xn*Xn - X)/(2.0*Xn);                                               -- find difference for this iteration
         if abs(Del) < ConvLim then                                                 -- if less than convergence limit..
            exit;                                                                   -- ..then converged; exit loop
         end if;
         Xn := Xn - Del;                                                            -- else update Xn for this iteration
      end loop;

      return Xn;                                                                    -- return final result

   end Sqrt;





-- ******************************************************************************************************************************
--
--  Main program code
--
-- ******************************************************************************************************************************

begin

   Put_Line ("LINREG - Perform linear regression");                                 -- print introductory message
   Put_Line (" (Enter X=-9999 to stop data entry and compute linear regression.)");

   loop                                                                             -- loop for all data points
      New_Line;
      Put ("Enter x:  ");                                                           -- prompt for x
      Get (x);                                                                      -- read x
      if x = -9999.0 then                                                           -- if no more data..
         exit;                                                                      -- ..then exit loop
      end if;                                                                       -- else prompt for y
      Put ("Enter y:  ");                                                           -- prompt for y
      Get (y);                                                                      -- read y

      n := n + 1.0;                                                                 -- increment number of data points by 1
      sumx  := sumx + x;                                                            -- compute sum of x
      sumx2 := sumx2 + x * x;                                                       -- compute sum of x**2
      sumxy := sumxy + x * y;                                                       -- compute sum of x * y
      sumy  := sumy + y;                                                            -- compute sum of y
      sumy2 := sumy2 + y * y;                                                       -- compute sum of y**2
   end loop;

   m := (n * sumxy  -  sumx * sumy) / (n * sumx2 - sumx**2);                        -- compute slope
   b := (sumy * sumx2  -  sumx * sumxy) / (n * sumx2  -  sumx**2);                  -- compute y-intercept
   r := (sumxy - sumx * sumy / n) /                                                 -- compute correlation coefficient
            Sqrt((sumx2 - sumx**2/n) * (sumy2 - sumy**2/n));

   New_Line;

   Put ("Slope        m = ");                                                       -- print results: slope
   Put (m, Fore => 2, Aft => 6, Exp => 3);

   New_Line;

   Put ("y-intercept  b = ");                                                       -- print y-intercept
   Put (b, Fore => 2, Aft => 6, Exp => 3);

   New_Line;

   Put ("Correlation  r = ");                                                       -- print correlation coefficient
   Put (r, Fore => 2, Aft => 6, Exp => 3);

   New_Line;

end LinReg;