Previous topic

Hessian of Two Particle Coulomb Potential

Next topic

Negative Binomial Regression

This Page

Minimal Surface ProblemΒΆ

In this example it is shown how the gradient of the function

def O_tilde(u):
    """ this is the objective function"""
    M = numpy.shape(u)[0]
    h = 1./(M-1)
    return M**2*h**2 + numpy.sum(0.25*( (u[1:,1:] - u[0:-1,0:-1])**2 + (u[1:,0:-1] - u[0:-1, 1:])**2))

can be computed conveniently using ALGOPY. It is the objective function of a minimal surface problem. With boundary constraints and an additional constraint in the domain in form of a cylinder one obtains a plot like

../_images/mayavi_3D_plot.png

At first the initial values are specified and the function evaluation is traced

# INITIAL VALUES
M = 30
h = 1./M
u = numpy.zeros((M,M),dtype=float)
u[0,:]=  [numpy.sin(numpy.pi*j*h/2.) for j in range(M)]
u[-1,:] = [ numpy.exp(numpy.pi/2) * numpy.sin(numpy.pi * j * h / 2.) for j in range(M)]
u[:,0]= 0
u[:,-1]= [ numpy.exp(i*h*numpy.pi/2.) for i in range(M)]

# trace the function evaluation and store it in cg
cg = algopy.CGraph()
Fu = algopy.Function(u)
Fy = O_tilde(Fu)
cg.trace_off()
cg.independentFunctionList = [Fu]
cg.dependentFunctionList = [Fy]

One can now compute the gradient of the objective function as follows

def dO_tilde(u):
    # use ALGOPY to compute the gradient
    g = cg.gradient([u])[0]

    # on the edge the analytical solution is fixed
    # so search direction must be zero on the boundary

    g[:,0]  = 0
    g[0,:]  = 0
    g[:,-1] = 0
    g[-1,:] = 0
    return g

The complete code is

import numpy
import algopy

def O_tilde(u):
    """ this is the objective function"""
    M = numpy.shape(u)[0]
    h = 1./(M-1)
    return M**2*h**2 + numpy.sum(0.25*( (u[1:,1:] - u[0:-1,0:-1])**2 + (u[1:,0:-1] - u[0:-1, 1:])**2))



# INITIAL VALUES
M = 30
h = 1./M
u = numpy.zeros((M,M),dtype=float)
u[0,:]=  [numpy.sin(numpy.pi*j*h/2.) for j in range(M)]
u[-1,:] = [ numpy.exp(numpy.pi/2) * numpy.sin(numpy.pi * j * h / 2.) for j in range(M)]
u[:,0]= 0
u[:,-1]= [ numpy.exp(i*h*numpy.pi/2.) for i in range(M)]

# trace the function evaluation and store it in cg
cg = algopy.CGraph()
Fu = algopy.Function(u)
Fy = O_tilde(Fu)
cg.trace_off()
cg.independentFunctionList = [Fu]
cg.dependentFunctionList = [Fy]


def dO_tilde(u):
    # use ALGOPY to compute the gradient
    g = cg.gradient([u])[0]

    # on the edge the analytical solution is fixed
    # so search direction must be zero on the boundary

    g[:,0]  = 0
    g[0,:]  = 0
    g[:,-1] = 0
    g[-1,:] = 0
    return g


def projected_gradients(x0, ffcn,dffcn, box_constraints, beta = 0.5, delta = 10**-3, epsilon = 10**-2, max_iter = 1000, line_search_max_iter = 100):
    """
    INPUT:	box_constraints		[L,U], where L (resp. U) vector or matrix with the lower (resp. upper) bounds
    """
    x = x0.copy()
    L = numpy.array(box_constraints[0])
    U = numpy.array(box_constraints[1])
    def pgn(s):
        a = 1.* (x>L)
        b = 1.*(abs(x-L) <0.00001)
        c = 1.*(s>0)
        d = numpy.where( a + (b*c))
        return numpy.sum(s[d])

    def P(x, s, alpha):
        x_alpha = x + alpha * s
        a = x_alpha-L
        b = U - x_alpha
        return x_alpha - 1.*(a<0) * a + b * 1. * (b<0)


    s = - dffcn(x)
    k = 0
    while pgn(s)>epsilon and k<= max_iter:
        k +=1
        s = - dffcn(x)
        for m in range(line_search_max_iter):
            #print 'm=',m
            alpha = beta**m
            x_alpha = P(x,s,alpha)
            if ffcn( x_alpha ) - ffcn(x) <= - delta * numpy.sum(s* (x_alpha - x)):
                break
        x_old = x.copy()
        x = x_alpha

    return x_old,s


# Setup of the optimization

# X AND Y PARTITION
x_grid = numpy.linspace(0,1,M)
y_grid = numpy.linspace(0,1,M)

# BOX CONSTRAINTS
lo = 2.5
L = numpy.zeros((M,M),dtype=float)

for n in range(M):
    for m in range(M):
        L[n,m] = 2.5 * ( (x_grid[n]-0.5)**2 + (y_grid[m]-0.5)**2 <= 1./16)

U = 100*numpy.ones((M,M),dtype=float)

Z,s = projected_gradients(u,O_tilde,dO_tilde,[L,U])


# # Plot with MAYAVI
x = y = list(range(numpy.shape(Z)[0]))

try:
    import enthought.mayavi.mlab as mlab
except:
    import mayavi.mlab as mlab
mlab.figure()
mlab.view(azimuth=130)
s = mlab.surf(x, y, Z, representation='wireframe', warp_scale='auto', line_width=1.)


mlab.savefig('./mayavi_3D_plot.png')
# mlab.show()