#
#	Test code for nmm.py (ver.1.0)
#
# 2020-06-15
# 2020-10-15
# by Kenar
#


import numpy as np
import nmm

pi = np.pi
sqrt = np.sqrt
sin = np.sin
exp = np.exp
abs = np.abs
arctan = np.arctan

names={}


# Test functions for optimization
# https://en.wikipedia.org/wiki/Test_functions_for_optimization
#
# Function Minimization
# http://distfiles.gentoo.org/distfiles/mntutorial.pdf
#
# JORGE J. MORE
# https://www.caam.rice.edu/~yzhang/caam454/nls/MGH.pdf
#
# R. Fletcher and M. J. D. Powell
# http://bioinfo.ict.ac.cn/~dbu/AlgorithmCourses/Lectures/Fletcher-Powell.pdf
#


a = 100
b = 0.01

#
#	t is a tuple of len(t) > 1
#


def f1(t):
    x,y = t
    #return a*(x+y-1)**2 + b*(x-y-1/2)**4
    #return (x-1/2)**2 + (y-1/2)**4
    #return (x-1/8)**2 + (y-3/4)**2

# minimum point of f2,g1,g2 is (3/4,1/4)

def f2(t):
  x,y = t
  return a*(x+y-1)**2 + b*(x-y-1/2)**4

def g1(t):
    x,y = t
    return a*abs(x+y-1) + b*(x-y-1/2)**4

def g2(t):
    x,y = t
    return a*sqrt(abs(x+y-1)) + b*(x-y-1/2)**4

# extended rosenbrock # ref: More
def rosenbrock(x):
	s = 0
	for k in range(0,len(x)-1):
		s += 100*(x[k+1]-x[k]**2)**2 + (x[k]-1)**2
	return s
  # f = 0 at (1,1)	# More for dim=2
  # f = 0 at (1,1,...,1)	# More

names[f1]="f1"
names[f2]="f2"
names[g1]="g1"
names[g2]="g2"


# Helical Valley # ref: Fletcher-Powell, or More
def helical_valley(t):
  x1,x2,x3 = t
  # -1/4 < th < 3/4
  # -2.4 < x3 < 7.5
  """- NOTE: the th(x1,x2) below fails with an error
  def th(x1,x2):
    if x1 > 0: return arctan(x2/x1)/(2*pi)
    else: return 1/2 + arctan(x2/x1)/(2*pi)
  -"""
  th = lambda x1,x2: np.where(x1>0, arctan(x2/x1)/(2*pi), 1/2 + arctan(x2/x1)/(2*pi))
  r = lambda x1,x2: sqrt(x1**2+x2**2)
  return 100*((x3 - 10*th(x1, x2))**2 + (r(x1,x2) - 1)**2) + x3**2
  # f = 0 at (1,0,0)	# More

# Powell badly scaled function # ref: More
def powell(t):
  x1,x2 = t
  return (10**4*x1*x2 - 1)**2 + (exp(-x1)+exp(-x2) - 1.0001)**2
  # f = 0 at (x1,x2)=(1.098 e-5, 9.106)	# More
  # comments by Kenar:
  # note that this function is symmetric with respect to  (x1,x2)
  # two global minimums at (1.09815936e-05, 9.10614648e+00) and (9.10614648e+00, 1.09815936e-05)
  # and one local minimum at (-0.00995081, -0.00994396)

# Biggs EXP6 function	# ref: More
def biggs(t):
  x1,x2,x3,x4,x5,x6 = t
  m = 13; s = 0
  for k in range(1,m+1):
    a = 0.1*k
    b = x3*exp(-a*x1) - x4*exp(-a*x2) + x6*exp(-a*x5)
    c = exp(-a) - 5*exp(-10*a) + 3*exp(-4*a)
    s += (b-c)**2
  return s
  # f = 5.65565 e-3 if m = 13	--- local minimum # More
  # f = 0 at (1,10,1,5,4,3)		--- real minimum # More
  
# https://en.wikipedia.org/wiki/Test_functions_for_optimization
def beale(t):
  x,y = t
  a = 1.5 - x + x*y
  b = 2.25 - x + x*y**2
  c = 2.625 - x + x*y**3
  return a**2 + b**2 + c**2
  # should have min 0 at (3,0.5)
    
def inpoint(x,ep):
  # ep is list of coordinates
  n = len(ep)
  xx = np.matrix(x+[1])
  xx = np.transpose(xx)
  epp = [0]*n
  for k in range(0,n):
    epp[k] = ep[k]+[1]
  A = np.matrix(epp)
  A= np.transpose(A)
  #print(xx)
  #print(A)
  z = np.linalg.solve(A,xx)
  #print(z)
  for zz in z:
    if zz <=0: return False
  return True
  
def test1_inpoint():
  x = [1/4,1/4]
  ep = [[0,0],[1,0],[0,1]]
  print(inpoint(x,ep))
  
#test1_inpoint()

def test1_search():
  search = nmm.anmsearch
  #search = nmm.snmsearch
  errs =(1e-1,1e-2,1e-3,1e-4,1e-5,1e-6,1e-7,1e-8,1e-9,1e-10)
  x0 = (0,0)
  mx = (3/4,1/4)
  for f in (f2,g1,g2):
    print(names[f])
    for err in errs:
      mp = search(f,x0,err,scale=0.1)
      print(err,mp,f(mp))
      no = np.linalg.norm(mp-mx)
      print(no<err, no,mp,mx)

#test1_search()

def test4_search():
  search = nmm.anmsearch
  of = rosenbrock	# objective function
  err = 1e-5
  for dim in range(2,10):
    x0 = np.random.uniform(-.5, .5, dim)
    mp = search(of,x0,err); print(dim,err,mp)
  # f = 0 at (1,1)	# More for dim=2
  # f = 0 at (1,1,...,1)	# More

#test4_search()

def test6_search():
  search = nmm.anmsearch
  of = helical_valley	# objective function
  err = 1e-5
  dim = 3
  for k in range(0,5):
    x0 = (-1,0,0) + np.random.uniform(-.5, .5, dim)
    mp = search(of,x0,err); print(x0,mp,of(mp))
  # f = 0 at (1,0,0)	# More

#test6_search()

def test6_search():
  search = nmm.anmsearch
  of = helical_valley	# objective function
  err = 1e-5
  dim = 3
  for k in range(0,5):
    x0 = (-1,0,0) + np.random.uniform(-.5, .5, dim)
    mp = search(of,x0,err); print(x0,mp,of(mp))
  # f = 0 at (1,0,0)	# More

#test6_search()

def test7_search():
  search = nmm.anmsearch
  of = powell	# objective function
  err = 1e-10
  dim = 2
  for k in range(0,5):
    x0 = (0,1) + np.random.uniform(-.5, .5, dim)
    mp = search(of,x0,err); print(x0,mp,of(mp))
  for k in range(0,5):
    x0 = (1,0) + np.random.uniform(-.5, .5, dim)
    mp = search(of,x0,err); print(x0,mp,of(mp))
  # f = 0 at (x1,x2)=(1.098 e-5, 9.106)	# More
  # comments by Kenar:
  # note that this function is symmetric with respect to  (x1,x2)
  # two global minimums at (1.09815933e-05, 9.10614674e+00) and (9.10614674e+00, 1.09815933e-05)
  # and one local minimum at (-0.00994809, -0.00994809)

test7_search()

def test8_search():
  search = nmm.anmsearch
  of = biggs	# objective function
  err = 1e-5
  dim = 6
  x0 = (1,10,1,5,4,3); print(x0,of(x0))
  for k in range(0,10):
    x0 = (1,2,1,1,1,1) + np.random.uniform(-.2, .2, dim)
    mp = search(of,x0,err); print(x0,mp,of(mp))
  # f = 5.65565 e-3 if m = 13	--- local minimum # More
  # f = 0 at (1,10,1,5,4,3)	--- real minimum # More

#test8_search()

def test9_search():
  search = nmm.anmsearch
  of = beale	# objective function
  err = 1e-5
  dim = 2
  print("beale")
  for k in range(0,5):
    x0 = (3,0.5) + np.random.uniform(-1, 1, dim)
    mp = search(of,x0,err); print(x0,mp,of(mp))
  # should have min 0 at (3,0.5)


#test9_search()