randomwalk

import numpy as np
np.random.seed( 101 )
def f( x,y ):
    return -((x-4)**2 + (y+2)**2)
# Random walk algorithm begins here. ###############################
X = np.linspace( -5,5,401 ) # set up a grid in x
Y = np.linspace( -5,5,401 ) # set up a grid in y
xy = np.random.randint( 401,size=(2,) )
u = np.array( ( 0,-1 ) )  # "up" looking at array in quadrant 4
d = np.array( ( 0,+1 ) )  # "down" looking at array in quadrant 4
l = np.array( ( -1,0 ) )  # "left" looking at array in quadrant 4
r = np.array( ( +1,0 ) )  # "right" looking at array in quadrant 4

num_steps = 500000
best_xy = xy
best_f = f( xy[ 0 ],xy[ 1 ] )
steps = np.empty( ( num_steps,3 ) )
for i in range( num_steps ):
    # Take a random step, 25% chance in each direction.
    trial_xy = xy.copy()
    chance = np.random.uniform()
    if chance < 0.25:
        trial_xy = ( xy + u ) % Y.shape[ 0 ]
    elif chance < 0.5:
        trial_xy = ( xy + d ) % Y.shape[ 0 ]
    elif chance < 0.75:
        trial_xy = ( xy + l ) % X.shape[ 0 ]
    else:
        trial_xy = ( xy + r ) % X.shape[ 0 ]

    if f( X[ trial_xy[ 0 ] ],Y[ trial_xy[ 1 ] ] ) > best_f:
        # If the solution improves, accept it.
        best_f = f( X[ trial_xy[ 0 ] ],Y[ trial_xy[ 1 ] ] )
        best_xy = trial_xy.copy()
        xy = trial_xy.copy()
    else:
        # If the solution does not improve, sometimes accept it.
        chance = np.random.uniform()
        if chance < 0.25:
            xy = trial_xy.copy()

    if i % 100 == 0:
        print( xy,f( X[ trial_xy[ 0 ] ],Y[ trial_xy[ 1 ] ] ),best_xy,best_f )

    # Log steps for plotting later.
    steps[ i,: ] = np.array( ( xy[ 0 ],xy[ 1 ],f( X[ trial_xy[ 0 ] ],Y[ trial_xy[ 1 ] ] ) ) )