# Discrete objective functional in dolfin-adjoint

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 *

# 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)
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)*sin(pi*x) / (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

def __str__(self):
return "MisfitBlock"

# d/dval_i(0.5 * (val_i - obs_i) ** 2.0) = val_i - obs_i
return np.sum([inputs(*pt) - obs_vals[i] for i, pt in enumerate(obs_coords)]) * adj_inputs

def recompute_component(self, inputs, block_variable, idx, prepared):
return backend_get_misfit(inputs)

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__
File "/home/dlilien/.pyenv/versions/3.8.1/lib/python3.8/site-packages/scipy/optimize/_differentiable_functions.py", line 171, in _update_grad
File "/home/dlilien/.pyenv/versions/3.8.1/lib/python3.8/site-packages/scipy/optimize/_differentiable_functions.py", line 85, in update_grad
File "/home/dlilien/.pyenv/versions/3.8.1/lib/python3.8/site-packages/scipy/optimize/_differentiable_functions.py", line 82, in grad_wrapped
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
File "/home/dlilien/.local/lib/python3.8/site-packages/pyadjoint/tape.py", line 46, in wrapper
return function(*args, **kwargs)