! usesort.f90
! written by Tim Mitchell
! module contains assorted routines that require SortMod
! contains: RobustLSR, DeTrendCol
! requires: SortMod higher up in the use statements

module UseSort

implicit none

contains

!*******************************************************************************
! performs robust regression between x and y according to
! Emerson, J. D., and D. C. Hoaglin. 1983.  Resistant lines 
!  for y versus x.  pp. 129-65. In D. C. Hoaglin, F. Mosteller,
!  and J. W. Tukey, eds., Understanding Robust and               
!  Exploratory Data Analysis.  John Wiley & Sons.
! received in the form of an IDL program from Mark New on 26.9.00
! pro robust_reg,xx,yy,a,b,yfit,miss=miss
! translated into an f90 program by Tim Mitchell on 26.9.00
! last modified 27.09.00
! y=a+bx

subroutine RobustLSR (ValueN,ExeSeries,WyeSeries,WyeFit,Aye,Bee,Iter)

use SortMod

real, pointer, dimension (:) 		:: ExeSeries, WyeSeries, WyeFit, ExeOp, WyeOp, Result
real, pointer, dimension (:)		:: ExeLeft, ExeMid, ExeRight
real, pointer, dimension (:)		:: WyeLeft, WyeMid, WyeRight
real, allocatable, dimension (:) 	:: ExeSorted, WyeSorted

integer, pointer, dimension (:)		:: Ranks
integer, pointer, dimension (:)		:: OrderLeft, OrderMid, OrderRight

real, intent (out) 	:: Aye, Bee

integer, intent (in) 	:: ValueN
integer, intent (out) 	:: Iter

real, parameter :: MissVal = -999.0
real, parameter :: VeryLarge = 1000000.0

real :: ExeMedLeft, ExeMedMid, ExeMedRight, WyeMedLeft, WyeMedMid, WyeMedRight
real :: Aye0, Aye1, Bee0, Bee1, Bee2, Are0, Are1, DeeAreBee0, DeeAreBee1, DeeBee1
real :: OldAye, OldBee, ChangeAye, ChangeBee

integer :: XValue, XMedian, XIter
integer :: AllocStat, ExitAll
integer :: ValidN, MissN
integer :: XLeftMin,XLeftMax,XMidMin,XMidMax,XRightMin,XRightMax,LeftN,MidN,RightN

!****************************************
! initialise

allocate (ExeOp  (ValueN), &
          WyeOp  (ValueN), &
          WyeFit (ValueN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: RobustLSR: Allocation failure #####"

Aye    = MissVal
Bee    = MissVal
WyeFit = MissVal

OldAye    = MissVal
OldBee    = MissVal
ChangeAye = MissVal
ChangeBee = MissVal

Iter    = 0
ExitAll = 0

ValidN = 0
do XValue = 1, ValueN
  if (ExeSeries(XValue).EQ.MissVal.OR.WyeSeries(XValue).EQ.MissVal) then
    ExeOp(XValue) = MissVal
    WyeOp(XValue) = MissVal
  else
    ExeOp(XValue) = ExeSeries(XValue)
    WyeOp(XValue) = WyeSeries(XValue)
    ValidN = ValidN + 1
  end if
end do

!****************************************
! main routine

if (ValidN.LT.9) then
	ExitAll = 1
	Iter = 0 - ValidN
end if

if (ExitAll.EQ.0) then
  Iter = 1
  
  allocate (Ranks    (ValueN), &
            ExeSorted(ValueN), &
            WyeSorted(ValueN), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: RobustLSR: Allocation failure #####"
  ExeSorted = MissVal
  WyeSorted = MissVal
  Ranks = MissVal
  
  call QuickSort (Reals=ExeLeft,OrderValid=Ranks,NMiss=MissN)
  
!  call RankVector (ValueN,MissN,ExeOp,Ranks)		! sort data
  
  do XValue = 1, ValidN
    ExeSorted (Ranks(XValue)) = ExeOp (XValue)
    WyeSorted (Ranks(XValue)) = WyeOp (XValue)
  end do
  
  XLeftMin  = 1							! specify thirds
  XLeftMax  = floor ( ValidN / 3.0 )
  LeftN     = XLeftMax
  
  XMidMin   = XLeftMax + 1
  XMidMax   = floor ( 2.0 * ValidN / 3.0 )
  MidN      = XMidMax - LeftN
  
  XRightMin = XMidMax + 1
  XRightMax = ValidN
  RightN    = XRightMax - MidN - LeftN
    
  allocate (ExeLeft (LeftN), &					! calc medians
            WyeLeft (LeftN), &
            ExeMid  (MidN),  &
            WyeMid  (MidN),  &
            ExeRight(RightN),&
            WyeRight(RightN),stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: RobustLSR: Allocation failure #####"
  
  ExeLeft  = ExeSorted (XLeftMin :XLeftMax )
  ExeMid   = ExeSorted (XMidMin  :XMidMax  )
  ExeRight = ExeSorted (XRightMin:XRightMax)
  
  WyeLeft  = WyeSorted (XLeftMin :XLeftMax )
  WyeMid   = WyeSorted (XMidMin  :XMidMax  )
  WyeRight = WyeSorted (XRightMin:XRightMax)
  
  call QuickSort (Reals=ExeLeft,NMiss=MissN)
  ExeMedLeft  = ExeLeft(nint(real(LeftN-MissN)/2.0))
  call QuickSort (Reals=WyeLeft,NMiss=MissN)
  WyeMedLeft  = ExeLeft(nint(real(LeftN-MissN)/2.0))
  
  call QuickSort (Reals=ExeRight,NMiss=MissN)
  ExeMedRight  = ExeRight(nint(real(RightN-MissN)/2.0))
  call QuickSort (Reals=WyeRight,NMiss=MissN)
  WyeMedRight  = ExeRight(nint(real(RightN-MissN)/2.0))
  
  call QuickSort (Reals=ExeMid,NMiss=MissN)
  ExeMedMid  = ExeMid(nint(real(MidN-MissN)/2.0))
  call QuickSort (Reals=WyeMid,NMiss=MissN)
  WyeMedMid  = ExeMid(nint(real(MidN-MissN)/2.0))
    
!  ExeMedLeft  = MedianValue (LeftN,  ExeLeft)
!  WyeMedLeft  = MedianValue (LeftN,  WyeLeft)
  
!  ExeMedMid   = MedianValue (MidN,   ExeMid)
!  WyeMedMid   = MedianValue (MidN,   WyeMid)

!  ExeMedRight = MedianValue (RightN, ExeRight)  
!  WyeMedRight = MedianValue (RightN, WyeRight)
  
  deallocate (ExeMid,WyeLeft,WyeRight,WyeMid, stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: RobustLSR: Deallocation3 failure #####"
  
  if ((ExeMedRight - ExeMedLeft).NE.0) then
  	Bee0 = (WyeMedRight - WyeMedLeft) / (ExeMedRight - ExeMedLeft)
  else
  	Bee0 = VeryLarge
  end if
  
  allocate (Result(ValidN),stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: RobustLSR: Allocation failure #####"

  do XValue = 1, ValidN
    Result(XValue) = (WyeSorted(XValue) - Bee0) * (ExeSorted(XValue) - ExeMedMid)
  end do
  
!  Aye0 = MedianValue (ValidN,Result)

  call QuickSort (Reals=Result,NMiss=MissN)
  Aye0  = Result(nint(real(ValidN-MissN)/2.0))

  do XValue = 1, ValidN
    Result(XValue) = WyeSorted(XValue) - (Bee0 * (ExeSorted(XValue) - ExeMedMid))
  end do
  
  ExeLeft  = Result (XLeftMin :XLeftMax )
  ExeRight = Result (XRightMin:XRightMax)
  
!  DeeAreBee0 = MedianValue (RightN,ExeRight) - MedianValue (LeftN,ExeLeft)

  call QuickSort (Reals=ExeRight,NMiss=MissN)
  DeeAreBee0 = ExeRight(nint(real(RightN-MissN)/2.0))

  call QuickSort (Reals=ExeLeft,NMiss=MissN)
  DeeAreBee0 = DeeAreBee0 - ExeLeft(nint(real(LeftN-MissN)/2.0))

  XIter = 0							! calculate b1 and drb1 iteratively
  do
    XIter = XIter + 1
    Iter = Iter + 1
    
    if ((ExeMedRight-ExeMedLeft).NE.0) then
    	DeeBee1 = DeeAreBee0 / (ExeMedRight-ExeMedLeft)
    else
    	DeeBee1 = VeryLarge
    end if

    if (XIter.EQ.1) Bee1 = Bee0 + DeeBee1
    
    do XValue = 1, ValidN
      Result(XValue) = WyeSorted(XValue) - (Bee1 * (ExeSorted(XValue) - ExeMedMid))
    end do
    
    ExeLeft  = Result (XLeftMin :XLeftMax )
    ExeRight = Result (XRightMin:XRightMax)
  
!    DeeAreBee1 = MedianValue (RightN,ExeRight) - MedianValue (LeftN,ExeLeft)

    call QuickSort (Reals=ExeRight,NMiss=MissN)
    DeeAreBee1 = ExeRight(nint(real(RightN-MissN)/2.0))

    call QuickSort (Reals=ExeLeft,NMiss=MissN)
    DeeAreBee1 = DeeAreBee1 - ExeLeft(nint(real(LeftN-MissN)/2.0))

    if (DeeAreBee1.GT.-0.001.AND.DeeAreBee1.LT.0.001) then
      do XValue = 1, ValidN
        Result(XValue) = WyeSorted(XValue) - (Bee1 * ExeSorted(XValue))
      end do
      
      call QuickSort (Reals=Result,NMiss=MissN)
      Aye  = Result(nint(real(ValidN-MissN)/2.0))

!      Aye = MedianValue (ValidN,Result)

      Bee = Bee1
      
      do XValue = 1, ValueN
        if (ExeSeries(XValue).NE.MissVal) then
          WyeFit(XValue) = Aye + (Bee * ExeSeries(XValue))
        end if
      end do
      
      ExitAll = 1    
    end if
					! interpolate between first and second estimates    
    if ((DeeAreBee1-DeeAreBee0).NE.0) then
      Bee2       = Bee1 - (DeeAreBee1 * ( (Bee1-Bee0) / (DeeAreBee1-DeeAreBee0) ) )
    else
      Bee2       = Bee1 - (DeeAreBee1 * VeryLarge)
    end if
    
    Bee0       = Bee1
    Bee1       = Bee2
    DeeAreBee0 = DeeAreBee1
    
    if (ExitAll.EQ.1.OR.XIter.EQ.10) exit
  end do

  deallocate (ExeSorted, WyeSorted, Result, Ranks, ExeLeft, ExeRight, stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: RobustLSR: Deallocation1 failure #####"
  								
end if

deallocate (ExeOp,WyeOp, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: RobustLSR: Deallocation2 failure #####"
  								
end subroutine RobustLSR

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

subroutine DeTrendCol (ExeSeries,WyeSeries)

real, pointer, dimension (:) :: ExeSeries, WyeSeries, WyeFit

real, parameter :: MissVal = -999.0

real :: Aye, Bee

integer :: Len, Iter, AllocStat, XValue

Len = size (ExeSeries)
if (Len.EQ.size(WyeSeries)) then
  call RobustLSR (Len,ExeSeries,WyeSeries,WyeFit,Aye,Bee,Iter)
  
  deallocate (WyeFit,stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: DeTrendCol: Deallocation3 failure #####"
  
  if (Aye.NE.MissVal.AND.Bee.NE.MissVal) then
    do XValue = 1, Len
      if (ExeSeries(XValue).NE.MissVal.AND.WyeSeries(XValue).NE.MissVal) then
        WyeSeries (XValue) = WyeSeries (XValue) - Aye - (Bee * ExeSeries (XValue))
      else
        WyeSeries (XValue) = MissVal
      end if
    end do   
  else
    do XValue = 1, Len
        WyeSeries (XValue) = MissVal
    end do   
  end if
else
  print*, "  > DeTrendCol: Input vectors have incompatible sizes."
  WyeSeries = MissVal
end if

end subroutine DeTrendCol

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

end module UseSort