I am attempting a PDE-constrained minimization problem where the objective functional is defined only at discrete points (geophysical measurements of a continuous field, sampled at O(1000) points in x-y space). I have the problem working with toy data and a continuous objective functional, but need to switch to discrete to use real data. Mathematically, I need to switch from
J = \int_\Omega\frac{1}{2}\left(y_{model}-y_{obs}\right)^2dx
to
J = \sum_i\frac{1}{2}\left(y^i_{model} - y^i_{obs}\right)^2
I am getting various pyadjoint errors when I attempt to use an OverloadedFunctional. MWE, modified from the Poisson example, follows (works with continuous J, fails with the discrete version I attempted). I assume that I am returning the wrong type from backend_misfit or evaluate_adj_component, but I do not really understand what to change or why this manifests in the DirichletBCs. Any help is appreciated. MWE:
import numpy as np
from dolfin import *
from dolfin_adjoint import *
from pyadjoint import Block
from pyadjoint.overloaded_function import overload_function
# Simple forward Poisson problem on unit square
n = 100
mesh = UnitSquareMesh(n, n)
V = FunctionSpace(mesh, "CG", 1)
W = FunctionSpace(mesh, "DG", 0)
f = interpolate(Expression("0.0", degree=1), W)
u = Function(V, name='State')
v = TestFunction(V)
F = (inner(grad(u), grad(v)) - f*v)*dx
bc = DirichletBC(V, 0.0, "on_boundary")
solve(F == 0, u, bc)
# Define functional of interest and the reduced functional
d = Expression("sin(pi*x[0])*sin(pi*x[1]) / (2.0 * pi * pi)", degree=3)
# Now try a discrete version by sampling the continuous version at 100 random points in the square
# Note that in practice, these would be measurements and not be interpolated from a function/expression
obs_coords = np.random.rand(100, 2)
obs_vals = np.array([d(*coords) for coords in obs_coords])
def get_misfit(func):
return 0.5 * np.sum([(func(*pt) - obs_vals[i]) ** 2.0 for i, pt in enumerate(obs_coords)])
backend_get_misfit = get_misfit
class MisfitBlock(Block):
def __init__(self, func, **kwargs):
super(MisfitBlock, self).__init__()
self.kwargs = kwargs
self.add_dependency(func)
def __str__(self):
return "MisfitBlock"
def evaluate_adj_component(self, inputs, adj_inputs, block_variable, idx, prepared=None):
# d/dval_i(0.5 * (val_i - obs_i) ** 2.0) = val_i - obs_i
return np.sum([inputs[0](*pt) - obs_vals[i] for i, pt in enumerate(obs_coords)]) * adj_inputs[0]
def recompute_component(self, inputs, block_variable, idx, prepared):
return backend_get_misfit(inputs[0])
get_misfit = overload_function(get_misfit, MisfitBlock)
alpha = Constant(1e-6)
# The code will run if I use the following (integrated) functional
# J = assemble((0.5*inner(u-d, u-d))*dx)
# But it fails with my attempt at a discrete functional
J = get_misfit(u)
control = Control(f)
rf = ReducedFunctional(J, control)
f_opt = minimize(rf, bounds=(0.0, 0.8), tol=1.0e-10, options={"gtol": 1.0e-10, "factr": 0.0})
Running on 2019.1.0 the error is:
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 0.000e+00 (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton solver finished in 0 iterations and 0 linear solver iterations.
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 0.000e+00 (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton solver finished in 0 iterations and 0 linear solver iterations.
/home/dlilien/.local/lib/python3.8/site-packages/pyadjoint/optimization/optimization.py:133: OptimizeWarning: Unknown solver options: factr
res = scipy_minimize(J, m_global, method=method, bounds=bounds, **kwargs)
Traceback (most recent call last):
File "example_discrete.py", line 66, in <module>
f_opt = minimize(rf, bounds=(0.0, 0.8), tol=1.0e-10, options={"gtol": 1.0e-10, "factr": 0.0})
File "/home/dlilien/.local/lib/python3.8/site-packages/pyadjoint/tape.py", line 46, in wrapper
return function(*args, **kwargs)
File "/home/dlilien/.local/lib/python3.8/site-packages/pyadjoint/optimization/optimization.py", line 254, in minimize
opt = algorithm(rf_np, **kwargs)
File "/home/dlilien/.local/lib/python3.8/site-packages/pyadjoint/optimization/optimization.py", line 133, in minimize_scipy_generic
res = scipy_minimize(J, m_global, method=method, bounds=bounds, **kwargs)
File "/home/dlilien/.pyenv/versions/3.8.1/lib/python3.8/site-packages/scipy/optimize/_minimize.py", line 617, in minimize
return _minimize_lbfgsb(fun, x0, args, jac, bounds,
File "/home/dlilien/.pyenv/versions/3.8.1/lib/python3.8/site-packages/scipy/optimize/lbfgsb.py", line 306, in _minimize_lbfgsb
sf = _prepare_scalar_function(fun, x0, jac=jac, args=args, epsilon=eps,
File "/home/dlilien/.pyenv/versions/3.8.1/lib/python3.8/site-packages/scipy/optimize/optimize.py", line 261, in _prepare_scalar_function
sf = ScalarFunction(fun, x0, args, grad, hess,
File "/home/dlilien/.pyenv/versions/3.8.1/lib/python3.8/site-packages/scipy/optimize/_differentiable_functions.py", line 95, in __init__
self._update_grad()
File "/home/dlilien/.pyenv/versions/3.8.1/lib/python3.8/site-packages/scipy/optimize/_differentiable_functions.py", line 171, in _update_grad
self._update_grad_impl()
File "/home/dlilien/.pyenv/versions/3.8.1/lib/python3.8/site-packages/scipy/optimize/_differentiable_functions.py", line 85, in update_grad
self.g = grad_wrapped(self.x)
File "/home/dlilien/.pyenv/versions/3.8.1/lib/python3.8/site-packages/scipy/optimize/_differentiable_functions.py", line 82, in grad_wrapped
return np.atleast_1d(grad(x, *args))
File "/home/dlilien/.local/lib/python3.8/site-packages/pyadjoint/optimization/optimization.py", line 66, in <lambda>
dJ = lambda m: rf_np.derivative(m, forget=forget, project=project)
File "/home/dlilien/.local/lib/python3.8/site-packages/pyadjoint/tape.py", line 46, in wrapper
return function(*args, **kwargs)
File "/home/dlilien/.local/lib/python3.8/site-packages/pyadjoint/reduced_functional_numpy.py", line 70, in derivative
dJdm = self.rf.derivative()
File "/home/dlilien/.local/lib/python3.8/site-packages/pyadjoint/reduced_functional.py", line 61, in derivative
derivatives = compute_gradient(self.functional,
File "/home/dlilien/.local/lib/python3.8/site-packages/pyadjoint/drivers.py", line 29, in compute_gradient
tape.evaluate_adj(markings=True)
File "/home/dlilien/.local/lib/python3.8/site-packages/pyadjoint/tape.py", line 140, in evaluate_adj
self._blocks[i].evaluate_adj(markings=markings)
File "/home/dlilien/.local/lib/python3.8/site-packages/pyadjoint/tape.py", line 46, in wrapper
return function(*args, **kwargs)
File "/home/dlilien/.local/lib/python3.8/site-packages/pyadjoint/block.py", line 126, in evaluate_adj
prepared = self.prepare_evaluate_adj(inputs, adj_inputs, relevant_dependencies)
File "/home/dlilien/.local/lib/python3.8/site-packages/fenics_adjoint/solving.py", line 182, in prepare_evaluate_adj
adj_sol, adj_sol_bdy = self._assemble_and_solve_adj_eq(dFdu_form, dJdu)
File "/home/dlilien/.local/lib/python3.8/site-packages/fenics_adjoint/solving.py", line 220, in _assemble_and_solve_adj_eq
bc.apply(dJdu)
TypeError: apply(): incompatible function arguments. The following argument types are supported:
1. (self: dolfin.cpp.fem.DirichletBC, arg0: dolfin::GenericVector) -> None
2. (self: dolfin.cpp.fem.DirichletBC, arg0: dolfin::GenericMatrix) -> None
3. (self: dolfin.cpp.fem.DirichletBC, arg0: dolfin::GenericMatrix, arg1: dolfin::GenericVector) -> None
4. (self: dolfin.cpp.fem.DirichletBC, arg0: dolfin::GenericVector, arg1: dolfin::GenericVector) -> None
5. (self: dolfin.cpp.fem.DirichletBC, arg0: dolfin::GenericMatrix, arg1: dolfin::GenericVector, arg2: dolfin::GenericVector) -> None
Invoked with: <dolfin.fem.dirichletbc.DirichletBC object at 0x7ff4a7560680>, -2.226135444562327