#! /usr/bin/env python3
#
def counterfeit_detection_compressed ( n, coin, correct, k ):

#*****************************************************************************80
#
## counterfeit_detection_compressed detects counterfeit coins.
#
#  Discussion:
#
#    We are given the weights of N coins, and the correct weight of a single 
#    coin.  We are asked to identify the counterfeit coins, that is, those
#    with incorrect weight.  We don't know how many such coins there are.
#
#    We use techniques from compressed sensing.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    23 March 2019
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Kurt Bryan, Tanya Leise,
#    Making do with less: an introduction to compressed sensing,
#    SIAM Review,
#    Volume 55, Number 3, September 2013, pages 547-566.
#
#  Parameters:
#
#    Input, integer n, the number of coins.
#
#    Input, real coin(n), the weights of the coins.
#
#    Input, real correct, the correct weight of a single coin.
#
#    Input, integer k, the number of rows to use in the sensing matrix.
#
#    Output, integer suspect, the index of the suspected counterfeit coin.
#
  import numpy as np
  from scipy.optimize import linprog
  from scipy.optimize import OptimizeWarning
  import warnings

  phi = np.random.randint ( 0, 2, size = ( k, n ) )

  b1 = np.matmul ( phi, coin )
  b2 = correct * np.sum ( phi, axis = 1 )
  b2 = np.reshape ( b2, ( k ) )
  b = b1 - b2
#
#  Find x which minimizes ||x||_1 subject to phi*x=b.
#
  f = np.ones ( n )
#
#  Without this, linprog will warn that the Phi matrix does not have
#  full row rank.  We don't care!
#
  with warnings.catch_warnings():
    warnings.filterwarnings ( 'ignore', category = OptimizeWarning )
    result = linprog ( f, A_ub=None, b_ub=None, A_eq=phi, b_eq=b, bounds = (0,None) )
#   result = linprog ( f, A_ub=None, b_ub=None, A_eq=phi, b_eq=b, bounds = (0,None), options={'tol':1.0e-8} )

  suspect = np.nonzero ( result.x )
  suspect = np.ravel ( suspect )

  return suspect

def counterfeit_detection_compressed_test ( k ):

#*****************************************************************************80
#
## counterfeit_detection_compressed_test tests counterfeit_detection_compressed.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    10 March 2019
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Kurt Bryan, Tanya Leise,
#    Making do with less: an introduction to compressed sensing,
#    SIAM Review,
#    Volume 55, Number 3, September 2013, pages 547-566.
#
#  Parameters:
#
#    Input, integer k, the number of random subsets to test.
#    If n = 100, k = 20 is reasonable.
#
  import numpy as np
  import platform

  print ( '' )
  print ( 'counterfeit_detection_compressed_test' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  Test counterfeit_detection_compressed,' )
  print ( '  which seeks to identify multiple counterfeit coins among n coins' )
  print ( '  using compressed sensing techniques.' )

  problem = 2

  if ( problem == 1 ):

    n = 7
    correct = 17.0
    coin = correct * np.ones ( n )
    fake = np.array ( [ 2 ] )
    fake_num = len ( fake )
    coin[fake[0]] = correct + 0.5 * np.random.random ( fake_num )

  else:
#
#  Very strange!  Setting all counterfeits to 17.5 improves chances of
#  getting a solution...
#
    n = 100
    correct = 17.0
    coin = correct * np.ones ( n )

    fake_num = 3
    fake = np.random.choice ( n, fake_num, replace = False )
    fake.sort( )

    for i in range ( 0, fake_num ):
#     coin[fake[i]] = correct + 0.5 * np.random.random ( 1 )
      coin[fake[i]] = correct + 0.5

  print ( '  There are %d coins to check.' % ( n ) )
  print ( '  The correct coin weight is %g' % ( correct ) )
#
#  Report the fakes.
#
  print ( '' )
  print ( '  There were %d fakes' % ( fake_num ) )
  print ( '' )
  print ( '     Fake  Index  Weight::' )
  print ( '' )
  for i in range ( 0, fake_num ):
    print ( '  %7d:  %4d  %g' % ( i, fake[i], coin[fake[i]] ) )
#
#  Detect and report the suspected fakes.
#
  suspect = counterfeit_detection_compressed ( n, coin, correct, k )
  suspect_num = len ( suspect )

  print ( '' )
  print ( '  %d random subsets were weighed.' % ( k ) )
  print ( '  The function found %d suspects.' % ( suspect_num ) )
  print ( '' )
  print ( '  Suspect  Index  Weight:' )
  print ( '' )
  for i in range ( 0, suspect_num ):
    print ( '  %7d:  %4d  %g' % ( i, suspect[i], coin[suspect[i]] ) )
#
#  Declare success or failure.
#
  if ( suspect_num != fake_num ):
    success = False
  elif ( np.array_equal ( fake, suspect ) ):
    success = True
  else:
    success = False

  print ( '' )
  if ( success ):
    print ( '  The test succeeded' )
  else:
    print ( '  The test failed.' )

  return

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( )
  k = 20
  counterfeit_detection_compressed_test ( k )
  timestamp ( )