! module PermMod
! deals with permutation testing
! contains: PValue

module PermMod

use SortMod
use ArrayMod

implicit none

contains

!*******************************************************************************
! all subscripts are supplied in the call
! ConCall and SusCall are two (not necessarily equivalent vectors), 
!    containing control and suspected diff vals, and sprinkled with MissVal
! PermMin and PermIdeal must be filled
! requires arraymod.f90 to be among main program Use statements

function PValue (ConCall,SusCall,PermMin,PermIdeal)

real PValue

real, pointer, dimension (:) 		:: ConCall, SusCall, ConValid, SusValid
real, pointer, dimension (:) 		:: Sums

integer, pointer, dimension (:) 	:: Ranks
integer, pointer, dimension (:,:) 	:: SampleKey

integer, intent (in) 	:: PermMin, PermIdeal

real, parameter :: MissVal = -999.0

integer :: AllocStat
integer :: ConN, SusN, PoolN, PossN, MissN, SampleN, HalfN
integer :: XSus, XCon, XMem, XPoss, XSample
integer :: QCalcAll
integer :: Integ, NextInteg, FinalPlacing, CurrentRank

!*******************************************
! perform initial checks and preparations

PValue = MissVal

call UnBlankReal (ConCall,ConValid)	! weed out missing values and reform arrays
call UnBlankReal (SusCall,SusValid)
ConN = size (ConValid,1)
SusN = size (SusValid,1)
if (ConN.LE.0) return			! check that the sample sizes are non-zero
if (SusN.LE.0) return

PoolN = ConN + SusN			! pooled sample

PossN = 0 				! total permutations in pooled sample 
XMem  = PoolN + 1			! only calc fully if < PermIdeal
do
  XMem = XMem - 1
  
  if (PossN.EQ.0) then
  	PossN = XMem
  else
  	PossN = PossN * XMem
  end if
  
  if (PossN.GT.PermIdeal) exit
  if (XMem.EQ.(PoolN-SusN+1)) exit
end do

QCalcAll = 1				! either calc all possible
if (PossN.GT.PermIdeal) QCalcAll = 2	! or     calc a random sample
					! or     return a MissVal
if (PermMin.EQ.MissVal.OR.PermMin.LE.0) QCalcAll = -1
if (PermIdeal.EQ.MissVal.OR.PermIdeal.LE.0) QCalcAll = -1
if (PossN.LT.PermMin) QCalcAll = -1

!*******************************************
! calc p-value

if (QCalcAll.EQ.1) then						! using all possible perms
  SampleN = PossN
  call KeyAllPoss (SampleN,PoolN,SusN,SampleKey)		! this mod sub: get key   
else if (QCalcAll.EQ.2) then					! random sampling from perms
  SampleN = PermIdeal
  call KeySample (SampleN,PoolN,SusN,SampleKey)			! this mod sub: get key  

  do XSus = 1, SusN						! replace first sample member with suspected
    SampleKey(1,XSus) = XSus
  end do
  do XCon = 1, ConN
    SampleKey(1,(SusN+XCon)) = MissVal
  end do
end if

allocate (Sums(SampleN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: PValue: Allocation failure #####"
Sums = 0.0

if (QCalcAll.GT.0) then
  call CalcSampleSums						! internal sub: calc sums

  deallocate (SampleKey, stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: PValue: Deallocation failure #####"
  
  call RankVector (SampleN,MissN,Sums,Ranks)			! sortmod.f90: ranks sums
  
  Integ = Ranks(1)
  HalfN = SampleN / 2
  
  if      (Integ.LT.SampleN.AND.Integ.GT.1) then  
    XSample = XSample + 1

    if      (Integ.GT.HalfN) then  
	NextInteg = Integ + 1
    else if (Integ.LT.HalfN) then
	NextInteg = Integ - 1
    end if
  
    FinalPlacing = 0
    CurrentRank  = 1
    XSample      = 1

    do
      if (Ranks(XSample).EQ.NextInteg) then
        if (Sums(XSample).EQ.Sums(CurrentRank)) then
        	Integ = Ranks(XSample)
        	CurrentRank = XSample
        	XSample = 1
        end if
      else
        FinalPlacing = 1
      end if
      
      if (Integ.EQ.SampleN.OR.Integ.EQ.1) exit
      if (FinalPlacing.EQ.1) exit
    end do
  end if
  
  if (Integ.EQ.SampleN) Integ = SampleN - 1			! prevents p=+1.0
  if (Integ.EQ.1)       Integ = 2				! prevents p=-1.0
  if (Integ.LT.HalfN)   Integ = Integ - SampleN - 1		! converts to negative equivalent
  
  PValue = real(Integ) / real(SampleN)
end if

deallocate (Sums,ConValid,SusValid, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: PValue: Deallocation failure #####"
  
contains

!*******************************************
! calc p-value using a random sample

subroutine CalcSampleSums

do XSample = 1, SampleN
    do XSus = 1, SusN
      if (SampleKey(XSample,XSus).NE.MissVal) then
	Sums(XSample) = Sums(XSample) + SusValid(XSus)
      end if
    end do

    do XCon = 1, ConN
      if (SampleKey(XSample,(SusN+XCon)).NE.MissVal) then
	Sums(XSample) = Sums(XSample) + ConValid(XCon)
      end if
    end do
end do
  
end subroutine CalcSampleSums

end function PValue

!*******************************************************************************
! give it PossN,PoolN,SusN and it gives back a key for the sample (PossN,PoolN)
!   MissVal = not part of permutation
!   positive integer = position in permutation (noramlly only .NE.MissVal is required)

subroutine KeyAllPoss (PossN,PoolN,SusN,SampleKey)

integer, pointer, dimension (:,:) 	:: SampleKey

real, parameter :: MissVal = -999.0

integer :: AllocStat
integer :: SusN, PoolN, PossN
integer :: XSus, XPoss, XPool, XMem
integer :: MemMove, MemMovePos, NextBlank, FinishedPerm

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

allocate (SampleKey(PossN,PoolN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: KeyAllPoss: Allocation failure #####"
SampleKey = MissVal

do XSus = 1, SusN
  SampleKey (1,XSus) = XSus
end do

XPoss = 1

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

do
  XPoss        = XPoss + 1			! loop once for each permutation
  MemMove      = SusN				! initialise member to move to the last member
  FinishedPerm = 0				! flags when perm has been made
  
  do XPool = 1, PoolN				! initialise current perm to the last perm
    SampleKey (XPoss,XPool) = SampleKey ((XPoss-1),XPool)
  end do

  do
    MemMovePos = 0				! find position of member to move (index)
    do
      MemMovePos = MemMovePos + 1
      if (SampleKey(XPoss,MemMovePos).EQ.MemMove) exit
    end do
    
    NextBlank = MemMovePos			! find next blank (MissVal or index)
    do
      NextBlank = NextBlank + 1
      if (NextBlank.GT.PoolN)                    NextBlank = MissVal
      if (NextBlank.EQ.MissVal)                  exit
      if (SampleKey(XPoss,NextBlank).EQ.MissVal) exit
    end do
    
    if (NextBlank.EQ.MissVal) then		! cannot move member down
      SampleKey (XPoss,MemMovePos) = MissVal
      MemMove = MemMove - 1
    else					! can move member down
      SampleKey (XPoss,MemMovePos) = MissVal
      SampleKey (XPoss,NextBlank)  = MemMove	! ...move member down
      
      if (MemMove.LT.SusN) then			! ...add in rest of members to first blanks, in order
        NextBlank = 0
        
        do XMem = (MemMove+1), SusN
          do
            NextBlank = NextBlank + 1
            if (SampleKey(XPoss,NextBlank).EQ.MissVal) exit
          end do
    
          SampleKey (XPoss,NextBlank) = XMem
        end do
      end if
      
      FinishedPerm = 1
    end if
    
    if (FinishedPerm.EQ.1) exit
  end do

  if (XPoss.EQ.PossN) exit			! exit when order is inverted  
end do

end subroutine KeyAllPoss

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

subroutine KeySample (SampleN,PoolN,SusN,SampleKey)

integer, pointer, dimension (:,:) 	:: SampleKey

real, parameter :: MissVal = -999.0

real :: Rand

integer :: AllocStat
integer :: SusN, PoolN, SampleN
integer :: XSample, XMem
integer :: Index

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

allocate (SampleKey(SampleN,PoolN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: KeySample: Allocation failure #####"
SampleKey = MissVal

call RandomiseSeed				! arraymod.f90

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

do XSample = 1, SampleN
  do XMem = 1, SusN
    do
      do
        call random_number (Rand)
        Rand  = Rand * (PoolN+1)		! avoids reduction of 1,PoolN probabilities
        Index = nint (Rand)
        if (Index.GE.1.AND.Index.LE.PoolN) exit
      end do
    
      if (SampleKey(XSample,Index).EQ.MissVal) SampleKey(XSample,Index) = XMem
      if (SampleKey(XSample,Index).EQ.XMem)    exit
    end do
  end do
end do

end subroutine KeySample

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

end module PermMod