! ΣΥΝΑΡΤΗΣΕΙΣ ΥΠΟΡΟΥΤΙΝΕΣ στην FORTRAN 95 :-)

! Άσκηση
! Δίνεται ο ακέραιος Ν (Ν<=100) και στην συνέχεια Ν ζεύγη (χ,ψ).
! Ζητείται να εκτελεστεί παλινδρόμηση / μοντελοποίηση με την μέθοδο των
! ελαχίστων τετραγώνων στα παραπάνω δεδομένα. Να υπολογιστούν η 
! γραμμική, εκθετική, λογαριθμική και δυναμική συνάρτηση παλινδρόμησης
! και τα αποτελέσματα να ταξινομηθούν σε φθίνουσα σειρά ρ^2.

! Λύση
!1. Διαβάζουμε καλά την άσκηση

!2. Δεδομένα: Ν, Α(100,2)
!   Ζητούμενα: Β(4,3) σε φθίνουσα σειρά ρ^2

!3. Θεωρία:

!3.1.
! Γραμμική  παλινδρόμηση y = A + B * x

! A = ( SUM(y) - B * SUM(x) ) / N

! B = ( N * SUM(x*y) - SUM(x) * SUM(y) ) / ( N * SUM(x^2) - SUM(x)^2 )

! R = ( N * SUM(x*y) - SUM(x) * SUM(y) ) / 
!      SQRT( ( N * SUM(x^2) - SUM(x)^2 ) * ( N * SUM(y^2) - SUM(y)^2 ) )

!3.2.
! Λογαριθμική παλινδρόμηση y = A + B * LN(x)

!3.3.
! Εκθετική παλινδρόμηση y = A * EXP(B*x)

!3.4.
! Δυναμική παλινδρόμηση y = A * x ^ B

!0. Άδεια χρήσης / License:

!    Regression analysis: Linear, Logarithnic, Exponential and Power
!    Copyright (C) May 2017 Ch Iossif - chiossif@yahoo.com

!    This program is free software: you can redistribute it and/or modify
!    it under the terms of the GNU General Public License as published by
!    the Free Software Foundation, either version 3 of the License, or
!    (at your option) any later version.

!    This program is distributed in the hope that it will be useful,
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!    GNU General Public License for more details.

!    You should have received a copy of the GNU General Public License
!    along with this program.  If not, see <http://www.gnu.org/licenses/>.

! ( η συνέχεια επί της οθόνης ;-) )

program regression

implicit none
real(16):: A(100,2) , B(4,3)
integer(2):: N , i, j, k, PIV(4)

! Read data into N and A
!read*,N
!read*,((A(i,j),j=1,2),i=1,N)

!DATA N/5/, A/10.0, 15.0, 20.0, 25.0, 30.0,95*0.0,1003.0,1005.0,1010.0,1008.0,1014.0,95*0.0/
!Results in Linear A=998.0 B=0.5 R=0.919018277

!DATA N/5/, A/29.0, 50.0, 74.0, 103.0, 118.0, 95*0.0, 1.6, 23.5, 38.0, 46.4, 48.9, 95*0.0/
!Results in Logarithmic A=-111.1283963 B=34.02014719 R=0.994013942

!DATA N/5/, A/6.9, 12.9, 19.8, 26.7, 35.1, 95*0.0, 21.4, 15.7, 12.1, 8.5, 5.2, 95*0.0/
!Results in Exponential A=30.49758743 B=-0.049203708 R=-0.997247351

DATA N/5/, A/28.0, 30.0, 33.0, 35.0, 38.0, 95*0.0, 2410.0, 3033.0, 3895.0, 4491.0, 5717.0, 95*0.0/
!Results in Power A=0.238801299 B=2.771865947 R=0.998906243

do i=1,N
  print*,i,A(i,1),A(i,2)
end do

! Call subroutines
call      linear_regression( A, N, B(1,1), B(1,2), B(1,3) )
call logarithmic_regression( A, N, B(2,1), B(2,2), B(2,3) )
call exponential_regression( A, N, B(3,1), B(3,2), B(3,3) )
call       power_regression( A, N, B(4,1), B(4,2), B(4,3) )

! Bubble sort B array by third column and keep track in PIV ;-)
! Initialize PIV
do i=1,4
  PIV(i)=i
end do

do i=1,3 !(3 = 4 - 1)
  do j=i+1,4
    if (ABS(B(i,3)).LE.ABS(B(j,3))) then
      do k=1,3
        call swap_real(B(i,k),B(j,k))
      end do
      call swap_int(PIV(i),PIV(j))
    endif
  end do
end do

!Display results
do i=1,4
  select case ( PIV(i) )
    case (1)
      print*,"Linear regression      : y = A + B * x"
    case (2)
      print*,"Logarithmic regression : y = A + B * LN(x)"
    case (3)
      print*,"Exponential regression : y = A * EXP(B*x)"
    case (4)
      print*,"Power regression       : y = A * x ^ B"
  end select
  print*,"A=",B(i,1)," B=",B(i,2)," R=",B(i,3)
end do
end program regression

!
! Linear regression      : y = A + B * x
!
subroutine linear_regression( Z, N, A, B, R)

implicit none
real(16),intent(IN)::Z(100,2)
integer(2),intent(IN)::N
real(16),intent(OUT)::A, B, R
real(16)::x(100), y(100), SUMXY, SUMM, SUMSQR
integer(2)::i

do i=1,N
  x(i)=Z(i,1)
  y(i)=Z(i,2)
end do

B = ( N * SUMXY(x,y,N) - SUMM(x,N) * SUMM(y,N) ) / ( N * SUMSQR(x,N) - SUMM(x,N)**2 )
A = ( SUMM(y,N) - B * SUMM(x,N) ) / N
R = ( N * SUMXY(x,y,N) - SUMM(x,N) * SUMM(y,N) ) / &
SQRT( ( N * SUMSQR(x,N) - SUMM(x,N)**2 ) * ( N * SUMSQR(y,N) - SUMM(y,N)**2 ) )

end subroutine linear_regression

!
! Logarithmic regression : y = A + B * LN(x)
!
subroutine logarithmic_regression( Z, N, A, B, R)

implicit none
real(16),intent(IN)::Z(100,2)
integer(2),intent(IN)::N
real(16),intent(OUT)::A, B, R
real(16)::NZ(100,2)
integer(2)::i

do i=1,N
  NZ(i,1) = LOG(Z(i,1))
  NZ(i,2) = Z(i,2)
end do

call linear_regression( NZ, N, A, B, R)

end subroutine logarithmic_regression

!
! Exponential regression : y = A * EXP(B*x)
!
subroutine exponential_regression( Z, N, A, B, R)

implicit none
real(16),intent(IN)::Z(100,2)
integer(2),intent(IN)::N
real(16),intent(OUT)::A, B, R
real(16)::NZ(100,2)
integer(2)::i

do i=1,N
  NZ(i,1) = Z(i,1)
  NZ(i,2) = LOG(Z(i,2))
end do

call linear_regression( NZ, N, A, B, R)
A=EXP(A)

end subroutine exponential_regression

!
! Power regression       : y = A * x ^ B
!
subroutine power_regression( Z, N, A, B, R)

implicit none
real(16),intent(IN)::Z(100,2)
integer(2),intent(IN)::N
real(16),intent(OUT)::A, B, R
real(16)::NZ(100,2)
integer(2)::i

do i=1,N
  NZ(i,1) = LOG(Z(i,1))
  NZ(i,2) = LOG(Z(i,2))
end do

call linear_regression( NZ, N, A, B, R)
A=EXP(A)

end subroutine power_regression

!
! Swap integers
!
subroutine swap_int( i , j )

implicit none
integer(2),intent(INOUT)::i,j
integer(2)::k

k=i
i=j
j=k

end subroutine swap_int

!
! Swap reals
!
subroutine swap_real( x , y )

implicit none
real(16),intent(INOUT)::x,y
real(16)::z

z=x
x=y
y=z

end subroutine swap_real

!
! Summation of array SUM(x)
!
real(16) function summ(A,N)
implicit none
real(16),intent(IN)::A(100)
integer(2),intent(IN)::N
real(16)::s
integer(2)::i

s=0.0
do i=1,n
  s=s+A(i)
end do
summ=s

end function summ

!
! Summation of product of two arrays SUM(x*y)
!
real(16) function sumxy(A,B,N)
implicit none
real(16),intent(IN)::A(100),B(100)
integer(2),intent(IN)::N
real(16)::s
integer(2)::i

s=0.0
do i=1,n
  s=s+A(i)*B(i)
end do
sumxy=s

end function sumxy

!
! Summation of squared array SUM(x*x)
!
real(16) function sumsqr(A,N)
implicit none
real(16),intent(IN)::A(100)
integer(2),intent(IN)::N
real(16)::s
integer(2)::i

s=0.0
do i=1,n
  s=s+A(i)**2
end do
sumsqr=s

end function sumsqr


! INPUT (first DATA uncommented)

!DATA N/5/, A/10.0, 15.0, 20.0, 25.0, 30.0,95*0.0,1003.0,1005.0,1010.0,1008.0,1014.0,95*0.0/
!Results in Linear A=998.0 B=0.5 R=0.919018277 OK !!!

!DATA N/5/, A/29.0, 50.0, 74.0, 103.0, 118.0, 95*0.0, 1.6, 23.5, 38.0, 46.4, 48.9, 95*0.0/
!Results in Logarithmic A=-111.1283963 B=34.02014719 R=0.994013942 OK !!!

!DATA N/5/, A/6.9, 12.9, 19.8, 26.7, 35.1, 95*0.0, 21.4, 15.7, 12.1, 8.5, 5.2, 95*0.0/
!Results in Exponential A=30.49758743 B=-0.049203708 R=-0.997247351

!DATA N/5/, A/28.0, 30.0, 33.0, 35.0, 38.0, 95*0.0, 2410.0, 3033.0, 3895.0, 4491.0, 5717.0, 95*0.0/
!Results in Power A=0.238801299 B=2.771865947 R=0.998906243 OK !!!

! OUTPUT

!      1   10.0000000000000000000000000000000000         1003.00000000000000000000000000000000      
!      2   15.0000000000000000000000000000000000         1005.00000000000000000000000000000000      
!      3   20.0000000000000000000000000000000000         1010.00000000000000000000000000000000      
!      4   25.0000000000000000000000000000000000         1008.00000000000000000000000000000000      
!      5   30.0000000000000000000000000000000000         1014.00000000000000000000000000000000      
! Exponential regression : y = A * EXP(B*x)
! A=   998.044637243313850405833714503806847        B=   4.95908410330473680128592728910886569E-0004  R=  0.919223898019449099581531227081690047      
! Linear regression      : y = A + B * x
! A=   998.000000000000000000000000000000000        B=  0.500000000000000000000000000000000000        R=  0.919018277617259685474995768382255906      
! Power regression       : y = A * x ^ B
! A=   982.050359107448159656957515008676855        B=   8.91316670635555559602100975564975024E-0003  R=  0.906490241716705324055658392081061782      
! Logarithmic regression : y = A + B * LN(x)
! A=   981.718277334977443710839502528224186        B=   8.98431924443863937001400424109889427        R=  0.906046716991343096115115792511713688      

!      1   29.0000000000000000000000000000000000         1.60000002384185791015625000000000000      
!      2   50.0000000000000000000000000000000000         23.5000000000000000000000000000000000      
!      3   74.0000000000000000000000000000000000         38.0000000000000000000000000000000000      
!      4   103.000000000000000000000000000000000         46.4000015258789062500000000000000000      
!      5   118.000000000000000000000000000000000         48.9000015258789062500000000000000000      
! Logarithmic regression : y = A + B * LN(x)
! A=  -111.128401937398062149857337388595604        B=   34.0201486701216915637000264462552176        R=  0.994013949456385815438002526587515196      
! Linear regression      : y = A + B * x
! A=  -6.37520368413094068253121046364516149        B=  0.508759415765386026264204016893651892        R=  0.955162965364315071875493316227293135      
! Power regression       : y = A * x ^ B
! A=   1.26304310200376217329882522148123265E-0003  B=   2.30422079193350905930954184203261347        R=  0.911190974363028861692676542459145221      
! Exponential regression : y = A * EXP(B*x)
! A=   1.79227449750179045865073540218806479        B=   3.22840634590068408619378640126738471E-0002  R=  0.820318023634031915222416330221578424      

!      1   6.90000009536743164062500000000000000         21.3999996185302734375000000000000000      
!      2   12.8999996185302734375000000000000000         15.6999998092651367187500000000000000      
!      3   19.7999992370605468750000000000000000         12.1000003814697265625000000000000000      
!      4   26.7000007629394531250000000000000000         8.50000000000000000000000000000000000      
!      5   35.0999984741210937500000000000000000         5.19999980926513671875000000000000000      
! Exponential regression : y = A * EXP(B*x)
! A=   30.4975878497284387508890819654099962        B=  -4.92037102186781262308280607496918960E-0002  R= -0.997247346357093354078303188491690199      
! Logarithmic regression : y = A + B * LN(x)
! A=   40.7012435153554213062984612984103209        B=  -9.82072487408214791343201339846950676        R= -0.996988102672733615553019211629376269      
! Linear regression      : y = A + B * x
! A=   23.9127674735158623446747174591356302        B= -0.558814977925160485399809626026184647        R= -0.985968008914535460909884351522129665      
! Power regression       : y = A * x ^ B
! A=   118.496995611507495157096496138366048        B= -0.822463737495838408803314127447302297        R= -0.959121091839313438593678761282510011      

!      1   28.0000000000000000000000000000000000         2410.00000000000000000000000000000000      
!      2   30.0000000000000000000000000000000000         3033.00000000000000000000000000000000      
!      3   33.0000000000000000000000000000000000         3895.00000000000000000000000000000000      
!      4   35.0000000000000000000000000000000000         4491.00000000000000000000000000000000      
!      5   38.0000000000000000000000000000000000         5717.00000000000000000000000000000000      
! Power regression       : y = A * x ^ B
! A=  0.238801068533772185093093444064770990        B=   2.77186615763770614776217725423968703        R=  0.998906255123514751386921209482980014      
! Exponential regression : y = A * EXP(B*x)
! A=   233.096123400178868850906410568698228        B=   8.46092320759881172024444207033381846E-0002  R=  0.997038300218968929838975925801104019      
! Linear regression      : y = A + B * x
! A=  -6707.55414012738853503184713375796190        B=   323.681528662420382165605095541401290        R=  0.996776699183564239957061509542566832      
! Logarithmic regression : y = A + B * LN(x)
! A=  -32822.4026632438505187203414979876959        B=   10541.2191458053039859444363956575950        R=  0.992725531080602923006148729585692377