Possible bug in `derivative_cb_post` of `ReducedFunctional`

[Charles Zhang]

Hi there,

When printing out the current control values through derivative_cb_post, it seems all control values m are accidentally set to be m[-1]. eval_cb_post does not have this problem.

from dolfin import *
from dolfin_adjoint import *

set_log_level(ERROR)


# Create mesh
nx = 9
mesh = UnitSquareMesh(nx, nx)

# Define function space
V = FunctionSpace(mesh, "CG", 1)

# Define a source dependent on three controls (c, d and e)
# 'f = c - d**2 - e**2'
class SourceExpression(Expression):
    def __init__(self, c, d, e, derivative=None):
        self.c = c
        self.d = d
        self.e = e
        self.derivative = derivative

    def eval(self, value, x):
        if self.derivative is None:
            value[0] = float(self.c) - float(self.d)**2  - float(self.e)**2
        elif self.derivative == self.c:
            value[0] = 1.0
        elif self.derivative == self.d:
            value[0] = -2.0*float(self.d)
        elif self.derivative == self.e:
            value[0] = -2.0*float(self.e)

def forward(f):
    """
    forward model
    """
    u = Function(V, name='State')
    v = TestFunction(V)
    F = inner(grad(u), grad(v)) * dx - f * v *dx
    bc = DirichletBC(V, 0.0, "on_boundary")
    solve(F == 0, u, bc)
    return u


def derivative_cb(j, dj, m):
    """
    output after each derivative evaluation
    """
    print "<<<< grad"
    print 'j  = ', j
    print 'dj = ', [float(va) for va in dj]
    print 'm  = ', [float(va) for va in m]
    print ">>>>"

def eval_cb(j, m):
    """
    output after each function evaluation
    """
    print "<<<< eval"
    print 'j  = ', j
    print 'm  = ', [float(va) for va in m]
    print ">>>>"

c = Constant(2.0)
d = Constant(3.0)
e = Constant(4.0)

f = SourceExpression(c, d, e)
f.dependencies = c, d, e

# Provide the derivative coefficients
f.user_defined_derivatives = {c: SourceExpression(c, d, e, derivative=c),
                              d: SourceExpression(c, d, e, derivative=d),
                              e: SourceExpression(c, d, e, derivative=e)}

# run forward model once
u = forward(f)
# define target functional
J = Functional(-inner(u, u) * dx)
# define controls
controls = [Control(j) for j in f.dependencies]
# create ReducedFunctional object
rf = ReducedFunctional(J, controls,
                       derivative_cb_post=derivative_cb, eval_cb_post=eval_cb)

# minimize
f_opt = minimize(rf, bounds=((0.0, 0.0, 0.0), (10.0, 10.0, 10.0)), tol=1.0e-10)
print float(f_opt[0]), float(f_opt[1]), float(f_opt[2])

[Simon Funke]

It might be an issue with an old dolfin-adjoint version?

I am running dolfin and dolfin-adjoint master and get the following output. Is this what you expect?

<<<< eval
j  =  -0.84575501663
m  =  [2.0, 3.0, 4.0]
>>>>
Solving for State:0:1:Adjoint
This took 0.002199588s.
Solving adjoint solution took 0.034512125s
Solving for State:0:0:Adjoint
Solving adjoint solution took 0.002178631s
<<<< grad
j  =  -0.84575501663
dj =  [0.0735439144895308, -0.44126348693718465, 5.1068568001257e-311]
m  =  [2.0, 3.0, 4.0]
>>>>
<<<< eval
j  =  -1.0737901939
m  =  [1.926456085510469, 3.4412634869371845, 4.0]
>>>>
Solving for State:0:1:Adjoint
This took 0.006321137s.
Solving adjoint solution took 0.029120697s
Solving for State:0:0:Adjoint
Solving adjoint solution took 0.001853134s
<<<< grad
j  =  -1.0737901939
dj =  [0.08286748677975964, -0.5703377130188737, 5.754281524285e-311]
m  =  [1.926456085510469, 3.4412634869371845, 4.0]
>>>>
<<<< eval
j  =  -1.46292197903
m  =  [1.8435885987307095, 4.011601199956059, 4.0]
>>>>
Solving for State:0:1:Adjoint
This took 0.002523567s.
Solving adjoint solution took 0.128733089s
Solving for State:0:0:Adjoint
Solving adjoint solution took 0.001618369s
<<<< grad
j  =  -1.46292197903
dj =  [0.09672417481665403, -0.7760376315184977, 6.7164837951223e-311]
m  =  [1.8435885987307095, 4.011601199956059, 4.0]
>>>>
<<<< eval
j  =  -4.67744675471
m  =  [1.512118651611671, 6.292952052031555, 4.0]
>>>>
Solving for State:0:1:Adjoint
This took 0.001926787s.
Solving adjoint solution took 0.01477909s
Solving for State:0:0:Adjoint
Solving adjoint solution took 0.001712308s
<<<< grad
j  =  -4.67744675471
dj =  [0.17295330964766933, -2.1767737697058998, 1.2009801104673e-310]
m  =  [1.512118651611671, 6.292952052031555, 4.0]
>>>>
<<<< eval
j  =  -21.1536206021
m  =  [0.9735013730975209, 10.0, 4.0]
>>>>
Solving for State:0:1:Adjoint
This took 0.001918611s.
Solving adjoint solution took 0.030719964s
Solving for State:0:0:Adjoint
Solving adjoint solution took 0.001564958s
<<<< grad
j  =  -21.1536206021
dj =  [0.36780430343682824, -7.356086068736564, 2.5540167682933e-310]
m  =  [0.9735013730975209, 10.0, 4.0]
>>>>
<<<< eval
j  =  -21.2891168908
m  =  [0.6056970696606927, 10.0, 4.0]
>>>>
Solving for State:0:1:Adjoint
This took 0.002262619s.
Solving adjoint solution took 0.036086942s
Solving for State:0:0:Adjoint
Solving adjoint solution took 0.001755452s
<<<< grad
j  =  -21.2891168908
dj =  [0.36898038031686464, -7.379607606337294, 2.56218339398087e-310]
m  =  [0.6056970696606927, 10.0, 4.0]
>>>>
<<<< eval
j  =  -21.5131937689
m  =  [0.0, 10.0, 4.0]
>>>>
Solving for State:0:1:Adjoint
This took 0.002134776s.
Solving adjoint solution took 0.015530506s
Solving for State:0:0:Adjoint
Solving adjoint solution took 0.001565134s
<<<< grad
j  =  -21.5131937689
dj =  [0.370917133947199, -7.418342678943977, 2.57563212528055e-310]
m  =  [0.0, 10.0, 4.0]
>>>>
0.0 10.0 4.0
RUNNING THE L-BFGS-B CODE

           * * *

Machine precision = 2.220D-16
 N =            3     M =           10

At X0         0 variables are exactly at the bounds

At iterate    0    f= -8.45755D-01    |proj g|=  4.41263D-01

At iterate    1    f= -1.07379D+00    |proj g|=  5.70338D-01
  ys=-5.764E-02  -gs= 2.001E-01 BFGS update SKIPPED

At iterate    2    f= -2.11536D+01    |proj g|=  3.67804D-01
  ys=-4.478E+01  -gs= 3.820E+00 BFGS update SKIPPED

At iterate    3    f= -2.15132D+01    |proj g|=  2.57563-310

           * * *

Tit   = total number of iterations
Tnf   = total number of function evaluations
Tnint = total number of segments explored during Cauchy searches
Skip  = number of BFGS updates skipped
Nact  = number of active bounds at final generalized Cauchy point
Projg = norm of the final projected gradient
F     = final function value

           * * *

   N    Tit     Tnf  Tnint  Skip  Nact     Projg        F
    3      3      7      3     2     1   2.576-310  -2.151D+01
  F =  -21.513193768937537     

CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL            

 Cauchy                time 0.000E+00 seconds.
 Subspace minimization time 0.000E+00 seconds.
 Line search           time 0.000E+00 seconds.

 Total User time 0.000E+00 seconds.

Not that I had to change the script slightly to get it to work on dolfin master:

from dolfin import *
from dolfin_adjoint import *

set_log_level(ERROR)


# Create mesh
nx = 9
mesh = UnitSquareMesh(nx, nx)

# Define function space
V = FunctionSpace(mesh, "CG", 1)

# Define a source dependent on three controls (c, d and e)
# 'f = c - d**2 - e**2'
class SourceExpression(Expression):
    def __init__(self, **kwargs):
        self.c = kwargs["c"]
        self.d = kwargs["d"]
        self.e = kwargs["e"]
        self.derivative = kwargs.get("derivative", None)

    def eval(self, value, x):
        if self.derivative is None:
            value[0] = float(self.c) - float(self.d)**2  - float(self.e)**2
        elif self.derivative == self.c:
            value[0] = 1.0
        elif self.derivative == self.d:
            value[0] = -2.0*float(self.d)
        elif self.derivative == self.e:
            value[0] = -2.0*float(self.e)

def forward(f):
    """
    forward model
    """
    u = Function(V, name='State')
    v = TestFunction(V)
    F = inner(grad(u), grad(v)) * dx - f * v *dx
    bc = DirichletBC(V, 0.0, "on_boundary")
    solve(F == 0, u, bc)
    return u


def derivative_cb(j, dj, m):
    """
    output after each derivative evaluation
    """
    print "<<<< grad"
    print 'j  = ', j
    print 'dj = ', [float(va) for va in dj]
    print 'm  = ', [float(va) for va in m]
    print ">>>>"

def eval_cb(j, m):
    """
    output after each function evaluation
    """
    print "<<<< eval"
    print 'j  = ', j
    print 'm  = ', [float(va) for va in m]
    print ">>>>"

c = Constant(2.0)
d = Constant(3.0)
e = Constant(4.0)

f = SourceExpression(c=c, d=d, e=e, degree=2)
f.dependencies = c, d, e

# Provide the derivative coefficients
f.user_defined_derivatives = {c: SourceExpression(c=c, d=d, e=2, derivative=c, degree=2),
                              d: SourceExpression(c=c, d=d, e=2, derivative=d, degree=2),
                              e: SourceExpression(c=c, d=d, e=2, derivative=e, degree=2)}

# run forward model once
u = forward(f)
# define target functional
J = Functional(-inner(u, u) * dx)
# define controls
controls = [Control(j) for j in f.dependencies]
# create ReducedFunctional object
rf = ReducedFunctional(J, controls,
                       derivative_cb_post=derivative_cb, eval_cb_post=eval_cb)

# minimize
f_opt = minimize(rf, bounds=((0.0, 0.0, 0.0), (10.0, 10.0, 10.0)), tol=1.0e-10)
print float(f_opt[0]), float(f_opt[1]), float(f_opt[2])

[Charles Zhang]

Thanks, Simon. I just pulled the dev-dolfin-adjoint version at quay.io, and the code works fine now.