#! /usr/bin/env python
#
def condition_sample1 ( n, a, m ):

#*****************************************************************************80
#
## CONDITION_SAMPLE1 estimates the L1 condition number of a matrix.
#
#  Discussion:
#
#    A naive sampling method is used.
#
#    Only "forward" sampling is used, that is, we only look at results
#    of the form y=A*x.
#
#    Presumably, solving systems A*y=x would give us a better idea of
#    the inverse matrix.
#
#    Moreover, a power sequence y1 = A*x, y2 = A*y1, ... and the same for
#    the inverse might work better too.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    04 October 2012
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer N, the dimension of the matrix.
#
#    Input, real A(N,N), the matrix.
#
#    Input, integer M, the number of samples to use.
#
#    Output, real VALUE, an estimate of the L1 condition number.
#
  import numpy as np
  from r8vec_norm_l1 import r8vec_norm_l1
  from r8vec_uniform_unit import r8vec_uniform_unit

  a_norm = 0.0
  ainv_norm = 0.0
  seed = 123456789

  for i in range ( 0, m ):

    x, seed = r8vec_uniform_unit ( n, seed )
    x_norm = r8vec_norm_l1 ( n, x )
    ax = np.dot ( a, x )
    ax_norm = r8vec_norm_l1 ( n, ax )

    if ( ax_norm == 0.0 ):
      value = 0.0
      return value

    a_norm    = max ( a_norm,    ax_norm / x_norm  )
    ainv_norm = max ( ainv_norm, x_norm  / ax_norm )

  value = a_norm * ainv_norm

  return value

def condition_sample1_test ( ):

#*****************************************************************************80
#
## CONDITION_SAMPLE1_TEST tests CONDITION_SAMPLE1.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    06 July 2015
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform
  from combin import combin
  from combin import combin_inverse
  from conex1 import conex1
  from conex1 import conex1_inverse
  from conex2 import conex2
  from conex2 import conex2_inverse
  from conex3 import conex3
  from conex3 import conex3_inverse
  from conex4 import conex4
  from conex4 import conex4_inverse
  from kahan import kahan
  from kahan import kahan_inverse
  from r8mat_norm_l1 import r8mat_norm_l1
  from r8mat_uniform_01 import r8mat_uniform_01

  m_test = np.array ( [ 10, 1000, 100000 ] )

  print ( '' )
  print ( 'CONDITION_SAMPLE1_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  CONDITION_SAMPLE1 estimates the L1 condition number using sampling' )
  print ( '  for a matrix in general storage.' )
  print ( '' )
  print ( '  Matrix                 Samples Order   Condition        Estimate' )
#
#  Combinatorial matrix.
#
  name = 'Combinatorial'
  n = 4
  alpha = 2.0
  beta = 3.0
  a = combin ( alpha, beta, n )
  a_inverse = combin_inverse ( alpha, beta, n )
  a_norm_l1 = r8mat_norm_l1 ( n, n, a )
  a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse )
  cond_l1 = a_norm_l1 * a_inverse_norm_l1
  print ( '' )
  for i in range ( 0, 3 ):
    m = m_test[i]
    cond = condition_sample1 ( n, a, m )
    print ( '  %20s  %8d  %4d  %14.6g  %14.6g' % ( name, m, n, cond_l1, cond ) )
#
#  CONEX1
#
  name = 'CONEX1'
  n = 4
  alpha = 100.0
  a = conex1 ( alpha )
  a_inverse = conex1_inverse ( alpha )
  a_norm_l1 = r8mat_norm_l1 ( n, n, a )
  a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse )
  cond_l1 = a_norm_l1 * a_inverse_norm_l1
  print ( '' )
  for i in range ( 0, 3 ):
    m = m_test[i]
    cond = condition_sample1 ( n, a, m )
    print ( '  %20s  %8d  %4d  %14.6g  %14.6g' % ( name, m, n, cond_l1, cond ) )
#
#  CONEX2
#
  name = 'CONEX2'
  n = 3
  alpha = 100.0
  a = conex2 ( alpha )
  a_inverse = conex2_inverse ( alpha )
  a_norm_l1 = r8mat_norm_l1 ( n, n, a )
  a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse )
  cond_l1 = a_norm_l1 * a_inverse_norm_l1
  print ( '' )
  for i in range ( 0, 3 ):
    m = m_test[i]
    cond = condition_sample1 ( n, a, m )
    print ( '  %20s  %8d  %4d  %14.6g  %14.6g' % ( name, m, n, cond_l1, cond ) )
#
#  CONEX3
#
  name = 'CONEX3'
  n = 5
  a = conex3 ( n )
  a_inverse = conex3_inverse ( n )
  a_norm_l1 = r8mat_norm_l1 ( n, n, a )
  a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse )
  cond_l1 = a_norm_l1 * a_inverse_norm_l1
  print ( '' )
  for i in range ( 0, 3 ):
    m = m_test[i]
    cond = condition_sample1 ( n, a, m )
    print ( '  %20s  %8d  %4d  %14.6g  %14.6g' % ( name, m, n, cond_l1, cond ) )
#
#  CONEX4
#
  name = 'CONEX4'
  n = 4
  a = conex4 ( )
  a_inverse = conex4_inverse ( )
  a_norm_l1 = r8mat_norm_l1 ( n, n, a )
  a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse )
  cond_l1 = a_norm_l1 * a_inverse_norm_l1
  print ( '' )
  for i in range ( 0, 3 ):
    m = m_test[i]
    cond = condition_sample1 ( n, a, m )
    print ( '  %20s  %8d  %4d  %14.6g  %14.6g' % ( name, m, n, cond_l1, cond ) )
#
#  KAHAN
#
  name = 'KAHAN'
  n = 4
  alpha = 0.25
  a = kahan ( alpha, n, n )
  a_inverse = kahan_inverse ( alpha, n )
  a_norm_l1 = r8mat_norm_l1 ( n, n, a )
  a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse )
  cond_l1 = a_norm_l1 * a_inverse_norm_l1
  print ( '' )
  for i in range ( 0, 3 ):
    m = m_test[i]
    cond = condition_sample1 ( n, a, m )
    print ( '  %20s  %8d  %4d  %14.6g  %14.6g' % ( name, m, n, cond_l1, cond ) )
#
#  Random
#
  seed = 123456789

  for j in range ( 0, 5 ):
    name = 'RANDOM'
    n = 4
    a, seed = r8mat_uniform_01 ( n, n, seed )
    a_inverse = np.linalg.inv ( a )
    a_norm_l1 = r8mat_norm_l1 ( n, n, a )
    a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse )
    cond_l1 = a_norm_l1 * a_inverse_norm_l1
    print ( '' )
    for i in range ( 0, 3 ):
      m = m_test[i]
      cond = condition_sample1 ( n, a, m )
      print ( '  %20s  %8d  %4d  %14.6g  %14.6g' % ( name, m, n, cond_l1, cond ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'CONDITION_SAMPLE1_TEST' )
  print ( '  Normal end of execution.' )
  return

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