! regress.f90
! written by Tim Mitchell
! module contains routines to perform various regression operations
! contains: RelSig, RegressVectors, LinearLSRZeroVec, LinearLSRVec

module Regress

use FileNames
use PlainPerm

implicit none

contains

!*******************************************************************************
! calcs the relationship between two vectors (y=a+bx), r=Pearson correl coeff
!	n=no. valid data pairs, p=prob that r has been achieved by chance
! you may specify:
!	number of permutations to calc (lower=quicker), 
!	threshold of 'r' at which to calc 'p' (higher=quicker),
!	threshold of 'p' at which to discard 'p' information (higher=quicker)

subroutine RelSig (Exe,Wye,Aye,Bee,Enn,Are,Pea,SpecBeg,SpecEnd, &
		   SpecNPerm,SpecAreThresh,SpecPeaThresh)

real, dimension (:), pointer 		:: Exe,Wye,PermWye
integer, dimension (:,:), pointer 	:: Perm
integer, dimension (:), pointer 	:: Pair

real, intent(out) 		:: Aye,Bee,Enn,Are,Pea
integer, intent(in), optional	:: SpecNPerm		! default = 1000
real, intent(in), optional	:: SpecAreThresh	! default = 0.80
real, intent(in), optional	:: SpecPeaThresh	! default = 0.95
integer, intent(in), optional	:: SpecBeg,SpecEnd	! default = 1,NVal

integer, dimension(1) 	:: Seed
real, parameter 	:: MissVal = -999.0

real 	:: AreThresh,PeaThresh,Rand,PermAye,PermBee,PermEnn,PermAre
real    :: Sample,Larger,LargerThresh
integer :: AllocStat, SubLen,Beg,End
integer :: NVal,XVal,NPerm,XPerm,NPair,XPair,NSub,XSub

character (len=10) :: TimeText
character (len= 8) :: DateText

!write (99,*), "call RelSig"		! @@@@@@@@@@@@@@@@@@@
NVal = size(Exe)
if (size(Wye).NE.NVal) print*, "  > ##### RelSig: bad-len #####"
AreThresh = 0.80
if (present(SpecAreThresh)) AreThresh = SpecAreThresh
PeaThresh = 0.95
if (present(SpecPeaThresh)) PeaThresh = SpecPeaThresh
Beg = 1
if (present(SpecBeg)) Beg=SpecBeg 
End = NVal
if (present(SpecEnd)) End=SpecEnd

NPerm = 1000
if (present(SpecNPerm)) NPerm = SpecNPerm
NSub = 10 ; SubLen = nint(real(NPerm)/real(NSub)) ; NPerm=SubLen*NSub

!write (99,*), "call RegressVectors 1"		! @@@@@@@@@@@@@@@@@@@
call RegressVectors (Exe,Wye,Aye,Bee,Are,Enn,5,Beg=Beg,End=End)
!write (99,*), "done RegressVectors 1"		! @@@@@@@@@@@@@@@@@@@

Pea=MissVal
if (Are.GE.AreThresh) then		! 'r' is sufficiently high to bother with 
  NPair = Enn
  allocate (Perm(SubLen,NPair),Pair(NPair),PermWye(NVal),stat=AllocStat)	
  if (AllocStat.NE.0) print*, "  > ##### RelSig: Perm: alloc #####"
  PermWye=MissVal ; Perm=0
  
  XPair = 0
  do XVal = Beg,End			! link series of pairs to full data series
    if (Exe(XVal).NE.MissVal.AND.Wye(XVal).NE.MissVal) then
      XPair=XPair+1 ; Pair(XPair)=XVal
    end if
  end do
  
  Larger=0 ; Sample=NPerm ; LargerThresh=(1.0-PeaThresh)*Sample ; XSub=0
  do
    XSub=XSub+1 ; XPerm=0
!    write (99,*), "call FlatNoReplace"		! @@@@@@@@@@@@@@@@@@@
    call FlatNoReplace (Perm)		! plainperm.f90 - generate permutations
!    write (99,*), "done FlatNoReplace"		! @@@@@@@@@@@@@@@@@@@
    
    do
      XPerm=XPerm+1
      
      do XPair = 1, NPair		! fill permed wye vector
        PermWye(Pair(XPair)) = Wye(Pair(Perm(XPerm,XPair)))
      end do
    
!      write (99,*), "call RegressVectors 2"		! @@@@@@@@@@@@@@@@@@@
      call RegressVectors (Exe,PermWye,PermAye,PermBee,PermAre,PermEnn,5, &
      			   Beg=Beg,End=End)
!      write (99,*), "done RegressVectors 2"		! @@@@@@@@@@@@@@@@@@@
      
      if (PermEnn.LT.5) then
        Sample=Sample-1
        LargerThresh=(1.0-PeaThresh)*Sample
      else if (PermAre.GT.Are) then
        Larger=Larger+1
      end if
    
      if (XPerm.EQ.SubLen.OR.Larger.GT.LargerThresh) exit
    end do
    
    if (XSub.EQ.NSub.OR.Larger.GT.LargerThresh) exit
  end do
  
  if (Larger.LE.LargerThresh) Pea=(Sample-Larger)/Sample
!  write (99,"(2f6.2,4i5)"), Are,Pea,XSub,XPerm,nint(Sample),nint(Larger) ! @@@@@@@@@@@@@@@@@@@@@@@@@
  
  deallocate (Perm,Pair,PermWye,stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### RelSig: Perm: dealloc #####"
end if
!write (99,*), "done RelSig"		! @@@@@@@@@@@@@@@@@@@

end subroutine RelSig

!*******************************************************************************
! produces a,b,r from y=a+bx, given x,y paired vectors

subroutine RegressVectors (Exe,Wye,Aye,Bee,Are,Enn,Min,RSS,RMSE,Beg,End)

real, dimension (:), pointer :: Exe,Wye
integer, intent(in) :: Min			! min pairs required for calc, .GE.2
real, intent(out) :: Aye,Bee,Are,Enn
real, intent(out), optional :: RSS	! residual sum of squares
real, intent(out), optional :: RMSE	! root-mean-squared-error (untested 28.4.03)
integer, intent(in), optional :: Beg,End	! only treat this range

real, parameter :: MissVal = -999.0
real    :: RealX,RealY
real    :: SumX,SumY,SumXX,SumYY,SumXY,SumN
real	:: Num,Den,One,Two
integer :: NVal,XVal,Val0,Val1

NVal = size(Exe) ; if (NVal.NE.size(Wye)) print*, "  > @@@@@ RegressVectors size"
Val0=1 ; Val1=NVal
if (present(Beg)) Val0=Beg
if (present(End)) Val1=End
if (Val0.GT.Val1.OR.Val0.LT.1.OR.Val1.GT.NVal) print*, "  > @@@@@ RegressVectors Beg,End"

SumX=0 ; SumY=0 ; SumXX=0 ; SumYY=0 ; SumXY=0 ; SumN=0
Aye=0 ; Bee=0 ; Are=0 ; Enn=0

do XVal = Val0,Val1
  if (Exe(XVal).NE.MissVal.AND.Wye(XVal).NE.MissVal) then
    SumX = SumX + Exe(XVal) 
    SumY = SumY + Wye(XVal)
    SumXX = SumXX + (Exe(XVal) * Exe(XVal)) 
    SumYY = SumYY + (Wye(XVal) * Wye(XVal)) 
    SumXY = SumXY + (Exe(XVal) * Wye(XVal))
    SumN  = SumN  + 1
  end if
end do

if (SumN.GE.real(Min)) then
  Den = (SumN*SumXX)-(SumX*SumX)
  if (Den.NE.0) then
    Num = (SumN*SumXY)-(SumX*SumY)
    Bee = Num / Den
    Aye = (SumY-(Bee*SumX))/SumN
    One = (SumXX/SumN)-((SumX/SumN)*(SumX/SumN))
    Two = (SumYY/SumN)-((SumY/SumN)*(SumY/SumN))
    
    if (One.GT.0.AND.Two.GT.0) then
      Den = sqrt(One) * sqrt(Two)
    else
      Den = 0.0
    end if
    
    if (Den.NE.0) then
      Num = (SumXY/SumN)-((SumX/SumN)*(SumY/SumN))
      Are = Num / Den
      Enn = SumN
    else
      Are = 1.0
      Enn = SumN
    end if
  end if
end if

if ((present(RSS).OR.present(RMSE))) then
 if (Enn.GT.0) then
  SumYY=0 ; SumN=0
  do XVal = Val0,Val1
    if (Exe(XVal).NE.MissVal.AND.Wye(XVal).NE.MissVal) then
      SumYY=SumYY+((Wye(XVal)-(Aye+(Bee*Exe(XVal))))**2)
      SumN=SumN+1
    end if
  end do
  if (present(RSS))  RSS  = SumYY
  if (present(RMSE)) RMSE = sqrt (SumYY / SumN)
 else
  if (present(RSS))  RSS  = MissVal
  if (present(RMSE)) RMSE = MissVal
 end if
end if

end subroutine RegressVectors

!*******************************************************************************

subroutine LinearLSRZeroVec (ExeVec,WyeVec,Aye)

real, dimension (:), pointer 	:: ExeVec, WyeVec
real, intent (out) 		:: Aye

real, parameter :: MissVal = -999.0

integer :: XYear, YearN

real :: SumXX, SumXY, SumEn

Aye = MissVal

YearN = size (ExeVec)
if (size(WyeVec).NE.YearN) print*, "  > ##### ERROR: LinearLSRZeroVec: vector mismatch #####"

SumXX = 0.0
SumXY = 0.0
SumEn = 0.0

do XYear = 1, YearN
    if (ExeVec(XYear).NE.MissVal.AND.WyeVec(XYear).NE.MissVal) then
      SumXX = SumXX + (ExeVec(XYear) * ExeVec(XYear))
      SumXY = SumXY + (ExeVec(XYear) * WyeVec(XYear))
      SumEn = SumEn + 1
    end if
end do

if (SumEn.GE.2.AND.SumXX.GT.0) Aye = SumXY / SumXX

end subroutine LinearLSRZeroVec

!*******************************************************************************

subroutine LinearLSRVec (Exe,Wye,Aye,Bee,Pea)

real, dimension (:), pointer 	:: Exe, Wye
real, intent (out) 		:: Bee, Aye, Pea	! y=a+bx, Pearson c/c

real, parameter :: MissVal = -999.0

integer :: XDatum, DataN

real :: SumX, SumY, SumXX, SumYY, SumXY, En, Num, Den

DataN = size (Exe)

Bee = MissVal
Aye = MissVal
Pea = MissVal

Num = 0.0
Den = 0.0

SumX  = 0.0
SumY  = 0.0
SumXX = 0.0
SumYY = 0.0
SumXY = 0.0
En    = 0.0

do XDatum = 1, DataN
    if (Exe(XDatum).NE.MissVal.AND.Wye(XDatum).NE.MissVal) then
      SumX  = SumX  + Exe(XDatum)
      SumY  = SumY  + Wye(XDatum)
      SumXX = SumXX + Exe(XDatum) * Exe(XDatum)
      SumYY = SumYY + Wye(XDatum) * Wye(XDatum)
      SumXY = SumXY + Exe(XDatum) * Wye(XDatum)
      En    = En    + 1.0
    end if
end do

if (En.GT.1) then
  Num = (En*SumXY)-(SumX*SumY)
  Den = (En*SumXX)-(SumX*SumX)
  if (Den.NE.0) Bee = Num / Den

  if (Bee.NE.MissVal) Aye = (SumY-(Bee*SumX))/En
    
  Num = (SumXY/En)-((SumX/En)*(SumY/En)) 
  Den =	sqrt((SumXX/En)-((SumX/En)*(SumX/En)))*sqrt((SumYY/En)-((SumY/En)*(SumY/En)))  
  if (Den.NE.0) Pea = Num / Den
end if

end subroutine LinearLSRVec

end module Regress