#! /usr/bin/env python
#
def i4_to_halton ( dim_num, step, seed, leap, base ):

#*****************************************************************************80
#
## I4_TO_HALTON computes one element of a leaped Halton subsequence.
#
#  Discussion:
#
#    The DIM_NUM-dimensional Halton sequence is really DIM_NUM separate
#    sequences, each generated by a particular base.
#
#    This routine selects elements of a "leaped" subsequence of the
#    Halton sequence.  The subsequence elements are indexed by a
#    quantity called STEP, which starts at 0.  The STEP-th subsequence
#    element is simply element
#
#      SEED(1:DIM_NUM) + STEP * LEAP(1:DIM_NUM)
#
#    of the original Halton sequence.
#
#    An I4 is an integer value.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    08 November 2014
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    John Halton,
#    On the efficiency of certain quasi-random sequences of points
#    in evaluating multi-dimensional integrals,
#    Numerische Mathematik,
#    Volume 2, Number 1, December 1960, pages 84-90
#
#    John Halton, GB Smith,
#    Algorithm 247:
#    Radical-Inverse Quasi-Random Point Sequence,
#    Communications of the ACM,
#    Volume 7, Number 12, December 1964, pages 701-702.
#
#    Ladislav Kocis, William Whiten,
#    Computational Investigations of Low-Discrepancy Sequences,
#    ACM Transactions on Mathematical Software,
#    Volume 23, Number 2, June 1997, pages 266-294.
#
#  Parameters:
#
#    Input, integer DIM_NUM, the spatial dimension.
#    1 <= DIM_NUM is required.
#
#    Input, integer STEP, the index of the subsequence element.
#    0 <= STEP is required.
#
#    Input, integer SEED(DIM_NUM), the Halton sequence index
#    corresponding to STEP = 0.
#    0 <= SEED(1:DIM_NUM) is required.
#
#    Input, integer LEAP(DIM_NUM), the successive jumps in the
#    Halton sequence.  1 <= LEAP(1:DIM_NUM) is required.
#
#    Input, integer BASE(DIM_NUM), the Halton bases.
#    1 < BASE(1:DIM_NUM) is required.
#
#    Output, real R(DIM_NUM), the STEP-th element of the leaped
#    Halton subsequence.
#
  import numpy as np
  from sys import exit
#
#  Check the input.
#
  if ( dim_num < 1 ):
    print ( '' )
    print ( 'I4_TO_HALTON - Fatal error!' )
    print ( '  DIM_NUM < 1.' )
    exit ( 'I4_TO_HALTON - Fatal error!' )

  if ( step < 0 ):
    print ( '' )
    print ( 'I4_TO_HALTON - Fatal error!' )
    print ( ' STEP < 0.' )
    exit ( 'I4_TO_HALTON - Fatal error!' )

  if ( any ( seed < 0 ) ):
    print ( '' )
    print ( 'I4_TO_HALTON - Fatal error!' )
    print ( '  Some SEED(*) < 0.' )
    exit ( 'I4_TO_HALTON - Fatal error!' )

  if ( any ( leap < 1 ) ):
    print ( '' )
    print ( 'I4_TO_HALTON - Fatal error!' )
    print ( '  Some LEAP < 1.' )
    exit ( 'I4_TO_HALTON - Fatal error!' )

  if ( any ( base <= 1 ) ):
    print ( '' )
    print ( 'I4_TO_HALTON - Fatal error!' )
    print ( '  Some BASE <= 1.' )
    exit ( 'I4_TO_HALTON - Fatal error!' )
#
#  Calculate the data.
#
  r = np.zeros ( dim_num )

  for i in range ( 0, dim_num ):
 
    seed2 = seed[i] + step * leap[i]

    base_inv = 1.0 / base[i]

    while ( seed2 != 0 ):
      digit = ( seed2 % base[i] )
      r[i] = r[i] + digit * base_inv
      base_inv = base_inv / base[i]
      seed2 = seed2 / base[i]
 
  return r

def i4_to_halton_test ( ):

#*****************************************************************************80
#
## I4_TO_HALTON_TEST tests I4_TO_HALTON. 
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    08 November 2014
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform

  print ( '' )
  print ( 'I4_TO_HALTON_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  I4_TO_HALTON computes a Halton sequence.' )
  print ( '  The user specifies all data explicitly.' )
  print ( '' )
  print ( '  In this test, we call I4_TO_HALTON repeatedly.' )
  print ( '  We use distinct primes as bases.' )

  n = 11

  dim_num = 3
  seed = np.array ( [ 0, 0, 0 ] )
  leap = np.array ( [ 1, 1, 1 ] )
  base = np.array ( [ 2, 3, 5 ] )

  print ( '' )
  print ( '   I    R(0)      R(1)      R(2)' )
  print ( '' )
  for step in range ( 0, n ):
    r = i4_to_halton ( dim_num, step, seed, leap, base )
    print ( '  %2d  %8.4f  %8.4f  %8.4f' % ( step, r[0], r[1], r[2] ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'I4_TO_HALTON_TEST' )
  print ( '  Normal end of execution.' )
  return

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( )
  i4_to_halton_test ( )
  timestamp ( )