MODULE module_DistriCorrection
! Modules for modIFication of a Ditribution

  USE module_definitions
  USE module_generic
  USE module_scientific

CONTAINS

!!!!!!! Subroutines
! randLowRemovePercent_HighEqui: Method to re-distribute from lower tail of thedistribution to the 
!   higher one following percentages preserving total amount. 

SUBROUTINE randLowRemovePercent_HighEqui3D(d1, d2, d3, dist, drng, rangeVals, dpercen, PercentRange,  &
  KindRemove, Nquants, Hquant, RepvalS, fillValue, dbg, newdist, redistributed, rngamount,            &
  Habovequant, aboveind)
! Method to re-distribute from lower tail of thedistribution to the higher one following percentages 
!   preserving total amount. This method, removes values from the lower end of the distribution and 
!   equi-distribute the total amount on a percentile on the high end

  IMPLICIT NONE

  INTEGER, INTENT(in)                                    :: d1, d2, d3, drng, dpercen, Nquants
  REAL(r_k), INTENT(in)                                  :: Hquant, fillValue
  REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in)             :: dist
  REAL(r_k), DIMENSION(drng), INTENT(in)                 :: rangeVals
  REAL(r_k), DIMENSION(dpercen), INTENT(in)              :: PercentRange
  CHARACTER(len=50), INTENT(in)                          :: KindRemove, RepvalS
  LOGICAL, INTENT(in)                                    :: dbg
  REAL(r_k), DIMENSION(d1,d2,d3), INTENT(out)            :: newdist
  REAL(r_k), DIMENSION(d1,d2,drng), INTENT(out)          :: rngamount
  INTEGER, DIMENSION(d1,d2,d3), INTENT(out)              :: redistributed, aboveind
  REAL(r_k), DIMENSION(d1,d2), INTENT(out)               :: Habovequant

! Local
  INTEGER                                                :: i, j, iv, ir, irm, ind, ia
  INTEGER                                                :: ierr
  INTEGER                                                :: irg, erg
  INTEGER                                                :: iHquant
  INTEGER                                                :: Nvals, Nrngvals, Nrmvals, Nredis, Nrmprev
  INTEGER                                                :: NaboveHquants
  INTEGER, DIMENSION(d3)                                 :: distind
  REAL(r_k), DIMENSION(d3)                               :: varvals
  REAL(r_k), DIMENSION(:),ALLOCATABLE                    :: varvalues, quants, newquants
  REAL(r_k)                                              :: one, zero
  REAL(r_k)                                              :: ndist, xdist, meandist
  REAL(r_k)                                              :: Repval, leakvals
  REAL(r_k)                                              :: iamount, initial_amount, eamount
  REAL(r_k)                                              :: TOTamount, dTOTamount
  INTEGER, DIMENSION(:), ALLOCATABLE                     :: rngind, rmind
  CHARACTER(len=4)                                       :: IS1, IS2
  CHARACTER(len=50)                                      :: spc, RS1
  CHARACTER(len=50)                                      :: fname
  CHARACTER(len=1000)                                    :: vectorS1, vectorS2, bspc
  LOGICAL, DIMENSION(d3)                                 :: search1, search2, search
  CHARACTER(len=256)                                     :: msg

!!!!!!! Variables
! d[1/2/3]: size of the values to treat
! dist: values to re-distribute
! drng: number of ranges provided
! rangeVals: ',' values from which establish the ranges of removals
! dpercen: number of percentages provided
! PercentRange: ',' percentages of removing to apply to each range (one value less than [rangeVals])
! KindRemove: kind of methodology to select which values to remove
!   'rand': randomly selection within the correspondent range
! Nquants: quantity of percentiles to compute
! Hquant: percentile of the distribution of values above which total removed amount will be equi-distributed
! RepvalS: replacing value to be assigned at the removed values
!   'min': initial minimal value of the distribution
!   'mean': initial mean value of the distribution
!   '[value]': any given value otherwise
! fillValue: missing value
!
! Outcome:
! newdist: Values of the new distribution
! redistributed: List of positions changed for each range
! rngamount: Amount re-distributed from each range
! Habovequant: value from which the amount is distributed
! aboveind: Indices of the values which received the distributed amount

  fname = 'randLowRemovePercent_HighEqui'

  spc = ' '
  bspc = Nstrings(spc,4)
  one = 1.
  zero = 0.

  IF (dbg) PRINT *,"  Re-aranging distribution of ", d3," values following '"//TRIM(fname)//"' method"

  ! Checking for the right amount of ranges and percentages
  IF (dpercen /= drng - 1) THEN
    WRITE(IS1,'(I4)')drng
    WRITE(IS2,'(I4)')dpercen
    vectorS1 = vectorR_KS(dpercen, PercentRange)
    vectorS2 = vectorR_KS(drng, rangeVals)

    msg = 'provided number of percentages: ' // TRIM(vectorS1) // ' (' // TRIM(IS2) //                &
         ') does not fit for: ' // TRIM(IS1) // 'ranges:' // TRIM(vectorS2)
    CALL ErrMsg(msg, fname, -1)
  END IF

  IF (dbg) THEN
    PRINT *,TRIM(bspc) // 'Ranges modIFication _______'
    DO ir=1, drng-1
      WRITE(RS1,'(F5.2)')PercentRange(ir)*100
      PRINT *,TRIM(bspc) // '  ' // TRIM(RS1) // ' % of (', rangeVals(ir), ',', rangeVals(ir+1), ']'
    END DO
  END IF

  IF (ALLOCATED(quants)) DEALLOCATE(quants)
  ALLOCATE(quants(Nquants+1), STAT=ierr)
  msg = "Problem allocating 'quants'"
  CALL ErrMsg(msg, fname, ierr)

  ! Checking for the upper-quantile
  quants = RangeR_K(Nquants+1, zeroRK, 100.*oneRK)

  iHquant = Index1DArrayR_K(quants, Nquants+1, hquant*100.*oneRK)

  IF (iHquant == -1) THEN
    WRITE(RS1,'(F50.10)')hquant*100.
    WRITE(IS1,'(I4)')Nquants
    vectorS1 = vectorR_KS(Nquants+1, quants)
    msg = 'provided high percentile ' // TRIM(RS1) // ' does not fit with ' // TRIM(IS1) //           &
      ' quantiles: ' // TRIM(vectorS1)
    CALL ErrMsg(msg, fname, iHquant)
  END IF

  IF (dbg) PRINT *,TRIM(bspc) // 'To be added above:', quants(iHquant), ' [', iHquant, '] percentile' 

  ! getting the right replacing value
  SELECT CASE (TRIM(repvalS))
    CASE ('min')
      repval = ndist
    CASE ('mean')
      repval = meandist
    CASE DEFAULT
      READ(RepvalS,'(F50.25)', IOSTAT=ierr)Repval
      msg = "For not 'min' or 'mean', a value is required as RepValS" // CHAR(10) // " provided: '" //&
        TRIM(RepvalS) // "'"
      CALL ErrMsg(msg, fname, ierr)
  END SELECT

  ! Future new distribution
  newdist = dist

  ! Total amounts for each range
  rngamount = zero

  ! Distribution indices
  distind = rangeI(d3,1,d3)

  ! Initialization of indices
  redistributed = -1
  aboveind = -1

  !! Looping on all dimensions
  DO i=1, d1
    DO j=1, d2

      varvals = dist(i,j,:)
      search = (varvals == fillValue)
      Nvals = COUNT(.NOT.search)

      doingnothing: IF (Nvals == 0) THEN
        PRINT *,TRIM(infomsg)
        PRINT *,'  ' // TRIM(fname) // ': all missed values for point:', i,', ',j, '!!'
        newdist(i,j,:) = dist(i,j,:)
        CYCLE

      ELSE

        IF (ALLOCATED(varvalues)) DEALLOCATE(varvalues)
        ALLOCATE(varvalues(Nvals))

        varvalues = PACK(varvals, MASK=.NOT.search)

        ! Amount of previous indices removed
        Nrmprev = 0

        ndist = MINVAL(varvalues)
        xdist = MAXVAL(varvalues)
        meandist = SUM(varvalues)/(d3*one)

        IF (xdist < rangeVals(drng)) THEN
          PRINT *,TRIM(infomsg)
          PRINT *,'  ' // TRIM(fname) // ': all values for point:', i,', ',j, ' are below ' //        &
            'ranges upper-end:', rangeVals(drng), ' !!'
          newdist(i,j,:) = dist(i,j,:)
          CYCLE 
        END IF

        CALL quantilesR_K(d3, varvalues, Nquants+1, quants)

        iamount = SUM(varvalues)
        initial_amount = iamount + 0.
    
        IF (dbg) THEN
          PRINT *,TRIM(bspc) // 'values _______'
          PRINT *,TRIM(bspc) // '  initial total amount:', iamount
          PRINT *,TRIM(bspc) // '  min:', ndist, 'max:', xdist, 'mean:', meandist
          PRINT *,TRIM(bspc) // '  quantiles:'
          CALL PrintQuantilesR_K(d3, varvalues, Nquants+1, quants, Nstrings(spc,6))
        END IF

        IF (dbg) THEN
          PRINT *, TRIM(bspc) // 'Loooping ranges _______'
        END IF
        bspc = Nstrings(spc,6)

        ! Getting total amount to remove and remove it from the value
        TOTamount = zero

        ranges: DO ir=1, drng-1
          irg = ir
          erg = ir + 1

          IF (dbg) THEN
            WRITE(RS1, '(F6.2)')PercentRange(ir)*100
            PRINT *,TRIM(bspc) // '* re-aranging ' // TRIM(RS1) // ' % of interval (', rangeVals(irg), ', ',  &
              rangeVals(erg), ']'
          END IF
  
          ! Searching the indices
          WHERE (varvalues > rangeVals(irg)) 
            search1 = .TRUE.
          ELSEWHERE
            search1 = .FALSE.
          END WHERE
  
          WHERE (varvalues <= rangeVals(erg)) 
            search2 = .TRUE.
          ELSEWHERE
            search2 = .FALSE.
          END WHERE

          WHERE (search1 .AND. search2)
            search = .TRUE.
          ELSEWHERE
            search = .FALSE.
          END WHERE

          ! Values within the range
          Nrngvals = COUNT(search)          
          IF (ALLOCATED(rngind)) DEALLOCATE(rngind)
          ALLOCATE(rngind(Nrngvals), STAT=ierr)
          msg = "Problem allocating 'rngind'"
          CALL ErrMsg(msg, fname, ierr)

          irm = 1
          DO iv=1, d3
            IF (search(iv)) THEN
              rngind(irm) = iv
              irm = irm + 1
            END IF
          END DO

          ! Number of values to remove
          Nrmvals = INT(Nrngvals*PercentRange(ir))

          IF (Nrmvals < 2) THEN
            PRINT *,TRIM(warnmsg)
            PRINT *,'  ' // TRIM(fname) // ': very few values to remove:', Nrmvals, ' !!'
            PRINT *,'    for '//TRiM(RS1)//' % of interval (', rangeVals(irg), ',', rangeVals(erg), ']'
          END IF

          ! Indices to remove
          IF (ALLOCATED(rmind)) DEALLOCATE(rmind)
          ALLOCATE(rmind(Nrmvals), STAT=ierr)
          msg = "Problem allocating 'rmind'"
          CALL ErrMsg(msg, fname, ierr)

          ! Removing kind
          Kremove: SELECT CASE (TRIM(KindRemove))
  
            CASE ('rand')
              CALL rand_sample(Nrngvals, Nrmvals, rmind)

            CASE DEFAULT
              msg = "Removing kind '" // TRIM(KindRemove) // "' not ready" // CHAR(10) //               &
                "Available ones: 'rand'"
              CALL ErrMsg(msg, fname, -1)
  
          END SELECT Kremove

          leakvals = zero
          DO irm=1, Nrmvals
            ind = rngind(rmind(irm))
            ! keeping the indices of the redistributed data
            redistributed(i,j,Nrmprev+irm) = ind
            rngamount(i,j,ir) = rngamount(i,j,ir) + (varvalues(ind) - repval)
            ! Keeping track of what is lost or gain during the transformation (for repval /= 0.)
            leakvals = leakvals + (repval - varvalues(ind))

            newdist(i,j,ind) = repval
          END DO

          TOTamount = TOTamount + rngamount(i,j,ir)
          IF (dbg) THEN
            PRINT *,TRIM(bspc) // '  removing:', Nrmvals, 'values from:', Nrngvals
            PRINT *,TRIM(bspc) // '  from indices:', rngind, 're-distributing:', rngind(rmind(1:Nrmvals))
            PRINT *,TRIM(bspc) // '  total amount to remove:', TOTamount
          END IF
          Nrmprev = Nrmprev + Nrmvals

        END DO ranges

        NonullRM: IF (TOTamount /= zero) THEN
          ! Equi-distributing the TOTamount on values >= hquant
          Habovequant(i,j) = quants(iHquant)
  
          IF (Habovequant(i,j) < rangeVals(drng)) THEN
            PRINT *,TRIM(infomsg)
            PRINT *,'    point:', i,', ',j,' has a High-quantile:', quants(iHquant), ' smaller' //    &
              ' than upper limit of ranges:', rangeVals(drng), '!!'
            !PRINT *,'    Trying to re-distribute on all the values above:', rangeVals(drng)
            Habovequant(i,j) = oneRK

            ! Searching the indices
            !WHERE (varvalues >= rangeVals(drng)) 
            WHERE (varvalues > oneRK) 
              search = .TRUE.
            ELSEWHERE
              search = .FALSE.
            END WHERE
            PRINT *,'    Re-distribute ', TOTamount, ' on', COUNT(search),' values above:', oneRK

          ELSE
            ! Searching the indices
            WHERE (varvalues >= Habovequant(i,j)) 
              search = .TRUE.
            ELSEWHERE
              search = .FALSE.
            END WHERE
          END IF
      
          ! Quantity to add at each above quantile
          NaboveHquants = COUNT(search)
          dTOTamount = TOTamount / NaboveHquants
  
          IF (dbg) THEN 
            PRINT *,TRIM(bspc) // 'incrementing:', NaboveHquants, 'values above:', Habovequant(i,j),  &
              ' by', dTOTamount
            PRINT *,TRIM(bspc) // '  replacing value:', repval
          END IF

          irm = 1
          DO iv=1, d3
            IF (search(iv)) THEN
              aboveind(i,j,irm) = iv
              newdist(i,j,iv) = varvalues(iv) + dTOTamount
              irm = irm + 1
            END IF
          END DO

          ! Final total amount of the new distribution
          eamount = SUM(newdist(i,j,:))

          IF (dbg) THEN
            Nredis = COUNT(redistributed(i,j,:) /= -1)
            PRINT *,TRIM(bspc), Nredis, 'Values removed from ranges ________'
            DO ia=1, Nredis
              ind = redistributed(i,j,ia)
              PRINT *,TRIM(bspc) // '    ', ia, ' prev:', varvalues(ind), 'new:', newdist(i,j,ind)
            END DO
            PRINT *,TRIM(bspc), NaboveHquants, 'Values above ________'
            DO ia=1, NaboveHquants
              ind = aboveind(i,j,ia)
              PRINT *,TRIM(bspc) // '    ', ia, 'prev:', varvalues(ind), 'new:', newdist(i,j,ind)
            END DO
          END IF

          ! Check
          IF (ABS(iamount - eamount) > 1.e-5) THEN
            PRINT *, warnmsg
            PRINT *,'  ' // TRIM(fname) // ': final quantity:', eamount, ' &', iamount, ' does not' // &
              ' fit: abs(iamount - eamount) > 1.e-5:', abs(iamount-eamount), ' TOTamount:', TOTamount
            STOP
          END IF

          IF (dbg) THEN
            varvalues = newdist(i,j,:)
            ndist = MINVAL(varvalues)
            xdist = MAXVAL(varvalues)
            meandist = SUM(varvalues)/(d3*one)

            IF (ALLOCATED(newquants)) DEALLOCATE(newquants)
            ALLOCATE(newquants(Nquants+1), STAT=ierr)
            msg = "Problem allocating 'newquants'"
            CALL ErrMsg(msg, fname, ierr)
  
            CALL quantilesR_K(d3, varvalues, Nquants+1, newquants)

            PRINT *,TRIM(bspc) // '  new distribution values _______'
            PRINT *,TRIM(bspc) // '  final total amount:', eamount
            PRINT *,TRIM(bspc) // '  min:', ndist, 'max:', xdist, 'mean:', meandist
            PRINT *,TRIM(bspc) // '  initial quantiles:'
            CALL PrintQuantilesR_K(20, varvalues, Nquants+1, quants, Nstrings(spc,8))
            PRINT *,TRIM(bspc) // '  new quantiles:'
            CALL PrintQuantilesR_K(20, varvalues, Nquants+1, newquants, Nstrings(spc,8))
            !STOP
          END IF

        ELSE

          PRINT *,TRIM(warnmsg)
          PRINT *,'  ' // TRIM(fname) // ': Any value has been removed !!'
          Habovequant(i,j) = zero

        END IF NonullRM

      END IF doingnothing
    END DO
    PRINT *,'  ' // TRIM(fname) // ': done', i, ' of', d1
  END DO

  IF (ALLOCATED(varvalues)) DEALLOCATE(varvalues)
  IF (ALLOCATED(rngind)) DEALLOCATE(rngind)
  IF (ALLOCATED(rmind)) DEALLOCATE(rmind)
  IF (ALLOCATED(quants)) DEALLOCATE(quants)
  IF (ALLOCATED(newquants)) DEALLOCATE(newquants)

  RETURN

END SUBROUTINE randLowRemovePercent_HighEqui3D

END MODULE module_DistriCorrection