! module PlainPerm
! simple routines for permutation testing
! contains: PartitionSample,FlatNoReplace

module PlainPerm

use FileNames

implicit none

contains

!*******************************************************************************
! A and B are two samples; under H0, they are two samples from the same pop
! this routine returns randomly partitioned A and B from the combined sample
! supply with the four vectors and start/finish

subroutine PartitionSample (SampleA,SampleB,PermA,PermB,BegA,EndA,BegB,EndB)

real, pointer, dimension (:) 	:: SampleA,SampleB 	! (NValue)
real, pointer, dimension (:,:) 	:: PermA,PermB		! (NPerm,NValue)

integer,pointer,dimension(:,:)  :: PoolPerm		! (NPerm,NValue)
integer,pointer,dimension(:)    :: Map		 	! (NPool)

integer, intent(in) :: BegA,EndA,BegB,EndB		! XValue indices

real, parameter :: MissVal = -999.0

integer :: AllocStat
integer :: XValidAye,XValidBee,XAye,XBee,XPerm,XPool
integer :: NValidAye,NValidBee,NAye,NBee,NPerm,NPool

! write (99,*), "running PartitionSample..."		! @@@@@@@@@@@@@@@@@@@@

NPerm=size(PermA,1) ; NAye=size(PermA,2) ; NBee=size(PermB,2)
if (size(SampleA).NE.NAye.OR.size(SampleB).NE.NBee) &
		print*, "  > ##### PartitionSample: Sample* #####"
if (size(PermB,1).NE.NPerm) &
		print*, "  > ##### PartitionSample: Perm* #####"
PermA=MissVal ; PermB=MissVal

NValidAye=0 ; NValidBee=0		! find valid sizes
do XAye = BegA,EndA
  if (SampleA(XAye).NE.MissVal) NValidAye=NValidAye+1
end do
do XBee = BegB,EndB
  if (SampleB(XBee).NE.MissVal) NValidBee=NValidBee+1
end do
NPool=NValidAye+NValidBee
if (NValidAye.LT.2.OR.NValidBee.LT.2.OR.NPool.LT.8) NPool=0

if (NPool.GT.0) then		! sufficient in each sample
  allocate (PoolPerm(NPerm,NPool),Map(NPool), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### PartitionSample: Alloc #####"
  XPool=0
  do XAye = BegA,EndA		! store sample indices in Pool
    if (SampleA(XAye).NE.MissVal) then
      XPool=XPool+1 ; Map(XPool)=XAye
    end if
  end do
  do XBee = BegB,EndB
    if (SampleB(XBee).NE.MissVal) then
      XPool=XPool+1 ; Map(XPool)=XBee
    end if
  end do
  
  call FlatNoReplace (PoolPerm)	! flat random version of pooled sample
  
  do XPerm = 1, NPerm		! each permutation
    do XPool = 1, NValidAye		! each valid place in sample A
      if (PoolPerm(XPerm,XPool).LE.NValidAye) then	! datum from A
        PermA(XPerm,Map(XPool)) = SampleA(Map(PoolPerm(XPerm,XPool)))
      else						! datum from B
        PermA(XPerm,Map(XPool)) = SampleB(Map(PoolPerm(XPerm,XPool)))
      end if
    end do
    do XPool = (NValidAye+1), NPool	! each valid place in sample B 
      if (PoolPerm(XPerm,XPool).LE.NValidAye) then	! datum from A
        PermB(XPerm,Map(XPool)) = SampleA(Map(PoolPerm(XPerm,XPool)))
      else						! datum from B
        PermB(XPerm,Map(XPool)) = SampleB(Map(PoolPerm(XPerm,XPool)))
      end if
    end do
  end do
    
  deallocate (PoolPerm,Map, stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### PartitionSample: Dealloc #####"
end if

! write (99,*), "done PartitionSample..."		! @@@@@@@@@@@@@@@@@@@@

end subroutine PartitionSample

!*******************************************************************************
! supply with a suitably sized array (NPerm,NVal)
! returned, filled with flat-random indices
! each XPerm has vector (1...NVal) filled with each value (1...NVal) once only
! this is suitable for selection without replacement

subroutine FlatNoReplace (Perm)

integer, pointer, dimension (:,:) 	:: Perm
logical, pointer, dimension (:) 	:: Free
integer, dimension(1) 			:: Seed

real, parameter :: MissVal = -999.0

integer :: NPerm,NVal,XPerm,XVal,XPlace,Match,NFree,XFree,AllocStat
real :: Rand
character (len=20) :: TimeText20
character (len=10) :: TimeText
character (len= 8) :: DateText

! write (99,*), "running FlatNoReplace..."	! @@@@@@@@@@@@@@@@@@@@

NPerm = size(Perm,1) ; NVal = size(Perm,2) ; Perm=MissVal

allocate (Free(NVal),stat=AllocStat)		
if (AllocStat.NE.0) print*, "  > ##### FlatNoReplace: Free: alloc #####"

call date_and_time (DateText, TimeText)
TimeText20 = adjustr(TimeText(8:10))
Seed(1) = GetIntFromText(TimeText20)		! filenames.f90
! call random_seed (put=Seed)			! compaq fortran version
call random_seed 

do XPerm = 1, NPerm		! each perm
  Free=.TRUE. ; NFree=NVal
  
  do XPlace = 1, NVal		! each place
    call random_number (Rand)
    Match = ceiling((Rand*NFree))	! random index = 1...NFree
    
    XVal=0 ; XFree=0
    do					! iterate to find this free value
      XVal=XVal+1
      if (Free(XVal)) XFree=XFree+1
      if (XFree.EQ.Match) then
        Perm(XPerm,XPlace) = XVal
        Free(XVal) = .FALSE.
        NFree=NFree-1
      end if
      if (Perm(XPerm,XPlace).GT.0) exit
    end do
  end do
end do

deallocate (Free,stat=AllocStat)		
if (AllocStat.NE.0) print*, "  > ##### FlatNoReplace: Free: dealloc #####"

! write (99,*), "done FlatNoReplace..."	! @@@@@@@@@@@@@@@@@@@@

end subroutine FlatNoReplace

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

end module PlainPerm