Merge lp:~libadjoint/dolfin-adjoint/improved-optimisation into lp:dolfin-adjoint

[Simon Funke]

This branch inplements several improvements for the optimisation interface.

Solving an optimisation problem is now as simple as:

>> minimize(reduced_functional)

where reduced_functional is an implementation of the reduced functional (i.e. reduced_functional(m) must return the functional evaluation with parameter choice m).

I also wrote a new class ReducedFunctional that implements such a reduced functional based on the current annotation of the forward model (and using the new functional implementation, cheers for that).
An example how to use it is:
>> reduced_functional = ReducedFunctional(J, InitialConditionParameter(m))

A functional value for a new parameter set can then be easily evaluated with (make sure that the forward model is annotated before calling this):
>> func_value = reduced_functional(m)

I also interfaced new optimisation algorithms: conjugate gradient, BFGS, and the truncated Newton method

All tests pass, so it might be worth giving it a try.

[lp:~libadjoint/dolfin-adjoint/improved-optimisation]

356. By Simon Funke on 2012-08-17

- optimisation -> optimization (follow the standard in programming)
- bugfix for scipy.slsqp method, where the parameter full_output = True caused the optimisation modele to crash
- new "scale" parameter which can be used to rescale the problem

357. By Simon Funke on 2012-08-17

- Add a maximization function
- Minor fixes

358. By Simon Funke on 2012-08-19

Minimize/maximize now return the optimal control

359. By David Ham on 2012-08-29

merge from trunk

=== modified file 'dolfin_adjoint/__init__.py'
--- dolfin_adjoint/__init__.py	2012-07-13 18:23:30 +0000
+++ dolfin_adjoint/__init__.py	2012-08-29 10:56:40 +0000
@@ -26,5 +26,6 @@
 import gst
 import function
 import genericmatrix
-import optimisation
+import optimization
+import reduced_functional 
 from ui import *

=== renamed file 'dolfin_adjoint/optimisation.py' => 'dolfin_adjoint/optimization.py'
--- dolfin_adjoint/optimisation.py	2012-07-03 15:16:10 +0000
+++ dolfin_adjoint/optimization.py	2012-08-29 10:56:40 +0000
@@ -1,57 +1,86 @@
-from dolfin import *
-from dolfin_adjoint import * 
+import dolfin
+from dolfin import cpp, MPI, info_red, info
+from dolfin_adjoint import constant, utils
+from reduced_functional import ReducedFunctional
 import numpy
 import sys
 
-def get_global(m):
-    ''' Takes a distributed object and returns a numpy array that contains all global values '''
-    if type(m) == float:
-        return numpy.array(m)
-    if type(m) == constant.Constant:
-        a = numpy.zeros(m.value_size())
-        p = numpy.zeros(m.value_size())
-        m.eval(a, p)
-        return a
-    elif type(m) in (function.Function, functions.function.Function):
-        m_v = m.vector()
-        m_a = cpp.DoubleArray(m.vector().size())
-        try:
-            m.vector().gather(m_a, numpy.arange(m_v.size(), dtype='I'))
-            return numpy.array(m_a.array())
-        except TypeError:
-            m_a = m.vector().gather(numpy.arange(m_v.size(), dtype='I'))
-            return m_a 
-    else:
-        raise TypeError, 'Unknown parameter type %s.' % str(type(m)) 
-
-def set_local(m, m_global_array):
-    ''' Sets the local values of the distrbuted object m to the values contained in the global array m_global_array '''
-    if type(m) == constant.Constant:
-        if m.rank() == 0:
-            m.assign(m_global_array[0])
-        else:
-            m.assign(Constant(tuple(m_global_array)))
-    elif type(m) in (function.Function, functions.function.Function):
-        range_begin, range_end = m.vector().local_range()
-        m_a_local = m_global_array[range_begin:range_end]
-        m.vector().set_local(m_a_local)
-        m.vector().apply('insert')
-    else:
-        raise TypeError, 'Unknown parameter type'
+def get_global(m_list):
+    ''' Takes a (optional a list of) distributed object(s) and returns one numpy array containing their global values '''
+    if not isinstance(m_list, (list, tuple)):
+        m_list = [m_list]
+
+    m_global = []
+    for m in m_list:
+        # Parameters of type float
+        if m == None or type(m) == float:
+            m_global.append(m)
+        # Parameters of type Constant 
+        elif type(m) == constant.Constant:
+            a = numpy.zeros(m.value_size())
+            p = numpy.zeros(m.value_size())
+            m.eval(a, p)
+            m_global += a.tolist()
+        # Function parameters of type Function 
+        elif hasattr(m, "vector"): 
+            m_v = m.vector()
+            m_a = cpp.DoubleArray(m.vector().size())
+            try:
+                m.vector().gather(m_a, numpy.arange(m_v.size(), dtype='I'))
+                m_global += m_a.array().tolist()
+            except TypeError:
+                m_a = m.vector().gather(numpy.arange(m_v.size(), dtype='I'))
+                m_global += m_a.tolist()
+        else:
+            raise TypeError, 'Unknown parameter type %s.' % str(type(m)) 
+
+    return numpy.array(m_global, dtype='d')
+
+def set_local(m_list, m_global_array):
+    ''' Sets the local values of a (or optionally  a list of) distributed object(s) to the values contained in the global array m_global_array '''
+
+    if not isinstance(m_list, (list, tuple)):
+        m_list = [m_list]
+
+    offset = 0
+    for m in m_list:
+        # Parameters of type dolfin.Constant 
+        if type(m) == constant.Constant:
+            m.assign(constant.Constant(numpy.reshape(m_global_array[offset:offset+m.value_size()], m.shape())))
+            offset += m.value_size()    
+        # Function parameters of type dolfin.Function 
+        elif hasattr(m, "vector"): 
+            range_begin, range_end = m.vector().local_range()
+            m_a_local = m_global_array[offset + range_begin:offset + range_end]
+            m.vector().set_local(m_a_local)
+            m.vector().apply('insert')
+            offset += m.vector().size() 
+        else:
+            raise TypeError, 'Unknown parameter type'
 
 def serialise_bounds(bounds, m):
-    ''' Converts bounds to an array of tuples and serialises it in a parallel environment. '''
-    
-    bounds_arr = []
+    ''' Converts bounds to an array of (min, max) tuples and serialises it in a parallel environment. '''
+
+    # Convert the bounds into the canoncial array form [ [lower_bound1, lower_bound2, ... ], [upper_bound1, upper_bound2, ...] ]
+    if len(numpy.array(bounds).shape) == 1:
+        bounds = numpy.array([[b] for b in bounds])
+
+    if len(bounds) != 2:
+        raise ValueError, "The 'bounds' parameter must be of the form [lower_bound, upper_bound] for one parameter or [ [lower_bound1, lower_bound2, ...], [upper_bound1, upper_bound2, ...] ] for multiple parameters."
+
+    bounds_arr = [[], []]
     for i in range(2):
-        if type(bounds[i]) == int or type(bounds[i]) == float:
-            bounds_arr.append(bounds[i]*numpy.ones(m.vector().size()))
-        else:
-            bounds_arr.append(get_global(bounds[i]))
-            
+        for j in range(len(bounds[i])):
+            if type(bounds[i][j]) in [int,  float, numpy.int32, numpy.int64]:
+                bounds_arr[i] += (bounds[i][j]*numpy.ones(m[j].vector().size())).tolist()
+            else:
+                bounds_arr[i] += get_global(bounds[i][j]).tolist()
+
+    # Transpose and return the array to get the form [ [lower_bound1, upper_bound1], [lower_bound2, upper_bound2], ... ] 
     return numpy.array(bounds_arr).T
 
-def minimise_scipy_slsqp(J, dJ, m, bounds = None, **kwargs):
+def minimize_scipy_slsqp(J, dJ, m, bounds = None, **kwargs):
+    ''' Interface to the SQP algorithm in scipy '''
     from scipy.optimize import fmin_slsqp
 
     m_global = get_global(m)
@@ -65,9 +94,13 @@
         mopt = fmin_slsqp(J, m_global, fprime = dJ, bounds = bounds, **kwargs)
     else:
         mopt = fmin_slsqp(J, m_global, fprime = dJ, **kwargs)
-    set_local(m, mopt)
+    if type(mopt) == list:
+        mopt = mopt[0]
+    set_local(m, numpy.array(mopt))
+    return m
 
-def minimise_scipy_fmin_l_bfgs_b(J, dJ, m, bounds = None, **kwargs):
+def minimize_scipy_fmin_l_bfgs_b(J, dJ, m, bounds = None, **kwargs):
+    ''' Interface to the L-BFGS-B algorithm in scipy '''
     from scipy.optimize import fmin_l_bfgs_b
     
     m_global = get_global(m)
@@ -81,72 +114,144 @@
 
     mopt, f, d = fmin_l_bfgs_b(J, m_global, fprime = dJ, bounds = bounds, **kwargs)
     set_local(m, mopt)
-
-optimisation_algorithms_dict = {'scipy.l_bfgs_b': ('The L-BFGS-B implementation in scipy.', minimise_scipy_fmin_l_bfgs_b),
-                                'scipy.slsqp': ('The SLSQP implementation in scipy.', minimise_scipy_slsqp) }
-
-def print_optimisation_algorithms():
-    ''' Prints the available optimisation algorithms '''
-
-    print 'Available optimisation algorithms:'
-    for function_name, (description, func) in optimisation_algorithms_dict.iteritems():
+    return m
+
+def minimize_scipy_tnc(J, dJ, m, bounds = None, **kwargs):
+    from scipy.optimize import fmin_tnc
+    
+    m_global = get_global(m)
+
+    # Shut up all processors but the first one.
+    if MPI.process_number() != 0:
+        kwargs['iprint'] = -1
+
+    if bounds:
+        bounds = serialise_bounds(bounds, m)
+
+    mopt, nfeval, rc = fmin_tnc(J, m_global, fprime = dJ, bounds = bounds, **kwargs)
+    set_local(m, mopt)
+    return m
+
+def minimize_scipy_cg(J, dJ, m, **kwargs):
+    from scipy.optimize import fmin_cg
+    
+    m_global = get_global(m)
+
+    # Shut up all processors but the first one.
+    if MPI.process_number() != 0:
+        kwargs['iprint'] = -1
+
+    mopt, fopt, func_calls, grad_calls, warnflag, allvecs = fmin_cg(J, m_global, fprime = dJ, **kwargs)
+    set_local(m, mopt)
+    return m
+
+def minimize_scipy_bfgs(J, dJ, m, **kwargs):
+    from scipy.optimize import fmin_bfgs
+    
+    m_global = get_global(m)
+
+    # Shut up all processors but the first one.
+    if MPI.process_number() != 0:
+        kwargs['iprint'] = -1
+
+    mopt, fopt, gopt, Bopt, func_calls, grad_calls, warnflag, allvecs = fmin_bfgs(J, m_global, fprime = dJ, **kwargs)
+    set_local(m, mopt)
+    return m
+
+optimization_algorithms_dict = {'scipy.l_bfgs_b': ('The L-BFGS-B implementation in scipy.', minimize_scipy_fmin_l_bfgs_b),
+                                'scipy.slsqp': ('The SLSQP implementation in scipy.', minimize_scipy_slsqp),
+                                'scipy.tnc': ('The truncated Newton algorithm implemented in scipy.', minimize_scipy_tnc), 
+                                'scipy.cg': ('The nonlinear conjugate gradient algorithm implemented in scipy.', minimize_scipy_cg), 
+                                'scipy.bfgs': ('The BFGS implementation in scipy.', minimize_scipy_bfgs), 
+                                }
+
+def print_optimization_algorithms():
+    ''' Prints the available optimization algorithms '''
+
+    print 'Available optimization algorithms:'
+    for function_name, (description, func) in optimization_algorithms_dict.iteritems():
         print function_name, ': ', description
 
-def minimise(reduced_functional, functional, parameter, m, algorithm, **kwargs):
+def minimize(reduced_func, algorithm = 'scipy.l_bfgs_b', scale = 1.0, **kwargs):
     ''' Solves the minimisation problem with PDE constraint:
 
-           min_m functional(u, m) 
+           min_m func(u, m) 
              s.t. 
            e(u, m) = 0
            lb <= m <= ub
            g(m) <= u
            
-        where m is the control variable, u is the solution of the PDE system e(u, m) = 0, functional is the functional of interest and lb, ub and g(m) constraints the control variables. 
-        The optimisation problem is solved using a gradient based optimisation algorithm and the functional gradients are computed by solving the associated adjoint system.
+        where m is the control variable, u is the solution of the PDE system e(u, m) = 0, func is the functional of interest and lb, ub and g(m) constraints the control variables. 
+        The optimization problem is solved using a gradient based optimization algorithm and the functional gradients are computed by solving the associated adjoint system.
 
-        The functional arguments are as follows:
-        * 'reduced_functional' must be a python function the implements the reduced functional (i.e. functional(u(m), m)). That is, it takes m as a parameter, solves the model and returns the functional value. 
-        * 'functional' must be a dolfin_adjoint.functional object describing the functional of interest
-        * 'parameter' must be a dolfin_adjoint.parameter that is to be minimised
-        * 'm' must contain the control values. The optimisation algorithm uses these values as a initial guess and updates them after each optimisation iteration. The optimal control values can be accessed by reading 'm' after calling minimise.
-        * 'bounds' is an optional keyword parameter to support control constraints: bounds = (lb, ub). lb and ub can either be floats to enforce a global bound or a dolfin.Function to define a varying bound.
-        * 'algorithm' specifies the optimistation algorithm to be used to solve the problem. The available algorithms can be listed with the print_optimisation_algorithms function.
+        The function arguments are as follows:
+        * 'reduced_func' must be a ReducedFunctional object. 
+        * 'algorithm' specifies the optimization algorithm to be used to solve the problem. The available algorithms can be listed with the print_optimization_algorithms function.
+        * 'scale' is a factor to scale to problem. Use a negative number to solve a maximisation problem.
+        * 'bounds' is an optional keyword parameter to support control constraints: bounds = (lb, ub). lb and ub must be of the same type than the parameters m. 
         
-        Additional arguments specific for the optimisation algorithms can be added to the minimise functions (e.g. iprint = 2). These arguments will be passed to the underlying optimisation algorithm. For detailed information about which arguments are supported for each optimisation algorithm, please refer to the documentaton of the optimisation algorithm.
+        Additional arguments specific for the optimization algorithms can be added to the minimize functions (e.g. iprint = 2). These arguments will be passed to the underlying optimization algorithm. For detailed information about which arguments are supported for each optimization algorithm, please refer to the documentaton of the optimization algorithm.
         '''
-    def dJ_array(m_array):
+
+    def reduced_func_deriv_array(m_array):
+        ''' An implementation of the reduced functional derivative that accepts the parameter as an array ''' 
 
         # In the case that the parameter values have changed since the last forward run, 
-        # we need to rerun the forward model with the new parameters
+        # we first need to rerun the forward model with the new parameters to have the 
+        # correct forward solutions
+        m = [p.data() for p in reduced_func.parameter]
         if (m_array != get_global(m)).any():
-            reduced_functional_array(m_array) 
+            reduced_func_array(m_array) 
 
-        dJdm = utils.compute_gradient(functional, parameter)
+        dJdm = utils.compute_gradient(reduced_func.functional, reduced_func.parameter)
         dJdm_global = get_global(dJdm)
 
-        if dolfin.parameters["optimisation"]["test_gradient"]:
-            minconv = utils.test_gradient_array(reduced_functional_array, dJdm_global, m_array, 
-                                                seed = dolfin.parameters["optimisation"]["test_gradient_seed"])
+        # Perform the gradient test
+        if dolfin.parameters["optimization"]["test_gradient"]:
+            minconv = utils.test_gradient_array(reduced_func_array, scale * dJdm_global, m_array, 
+                                                seed = dolfin.parameters["optimization"]["test_gradient_seed"])
             if minconv < 1.9:
                 raise RuntimeWarning, "A gradient test failed during execution."
             else:
                 info("Gradient test succesfull.")
-            reduced_functional_array(m_array) 
-
-        return dJdm_global 
-
-    def reduced_functional_array(m_array):
-
-        # Reset any prior annotation of the adjointer as we are about to rerun the forward model.
-        solving.adj_reset()
-        # If functional is a FinalFunctinal, we need to set the activated flag to False
-        if hasattr(functional, 'activated'):
-            functional.activated = False
-
+            reduced_func_array(m_array) 
+
+        return scale * dJdm_global 
+
+    def reduced_func_array(m_array):
+        ''' An implementation of the reduced functional that accepts the parameter as an array '''
+        # In case the annotation is not reused, we need to reset any prior annotation of the adjointer before reruning the forward model.
+        if not reduced_func.replays_annotation:
+            solving.adj_reset()
+
+        # Set the parameter values and execute the reduced functional
+        m = [p.data() for p in reduced_func.parameter]
         set_local(m, m_array)
-        return reduced_functional(m)
-
-    if algorithm not in optimisation_algorithms_dict.keys():
-        raise ValueError, 'Unknown optimisation algorithm ' + algorithm + '. Use the print_optimisation_algorithms to get a list of the available algorithms.'
-
-    optimisation_algorithms_dict[algorithm][1](reduced_functional_array, dJ_array, m, **kwargs)
+        return scale * reduced_func(m)
+
+    if algorithm not in optimization_algorithms_dict.keys():
+        raise ValueError, 'Unknown optimization algorithm ' + algorithm + '. Use the print_optimization_algorithms to get a list of the available algorithms.'
+
+    return optimization_algorithms_dict[algorithm][1](reduced_func_array, reduced_func_deriv_array, [p.data() for p in reduced_func.parameter], **kwargs)
+
+def maximize(reduced_func, algorithm = 'scipy.l_bfgs_b', scale = 1.0, **kwargs):
+    ''' Solves the maximisation problem with PDE constraint:
+
+           max_m func(u, m) 
+             s.t. 
+           e(u, m) = 0
+           lb <= m <= ub
+           g(m) <= u
+           
+        where m is the control variable, u is the solution of the PDE system e(u, m) = 0, func is the functional of interest and lb, ub and g(m) constraints the control variables. 
+        The optimization problem is solved using a gradient based optimization algorithm and the functional gradients are computed by solving the associated adjoint system.
+
+        The function arguments are as follows:
+        * 'reduced_func' must be a ReducedFunctional object. 
+        * 'algorithm' specifies the optimization algorithm to be used to solve the problem. The available algorithms can be listed with the print_optimization_algorithms function.
+        * 'scale' is a factor to scale to problem. Use a negative number to solve a maximisation problem.
+        * 'bounds' is an optional keyword parameter to support control constraints: bounds = (lb, ub). lb and ub must be of the same type than the parameters m. 
+        
+        Additional arguments specific for the optimization algorithms can be added to the minimize functions (e.g. iprint = 2). These arguments will be passed to the underlying optimization algorithm. For detailed information about which arguments are supported for each optimization algorithm, please refer to the documentaton of the optimization algorithm.
+        '''
+    return minimize(reduced_func, algorithm, scale = -scale, **kwargs)

=== modified file 'dolfin_adjoint/options.py'
--- dolfin_adjoint/options.py	2012-06-11 11:13:00 +0000
+++ dolfin_adjoint/options.py	2012-08-29 10:56:40 +0000
@@ -7,7 +7,7 @@
 adj_params.add("fussy_replay", True)
 adj_params.add("stop_annotating", False)
 
-opt_params = Parameters("optimisation")
+opt_params = Parameters("optimization")
 opt_params.add("test_gradient", False)
 opt_params.add("test_gradient_seed", 0.0001)
 

=== modified file 'dolfin_adjoint/parameter.py'
--- dolfin_adjoint/parameter.py	2012-07-30 20:44:01 +0000
+++ dolfin_adjoint/parameter.py	2012-08-29 10:56:40 +0000
@@ -23,6 +23,7 @@
     '''coeff: the variable whose initial condition you wish to perturb.
        perturbation: the perturbation direction in which you wish to compute the gradient. Must be a Function.'''
 
+    self.coeff = coeff
     self.var = None 
     # Find the first occurance of the coeffcient
     for t in range(adjglobals.adjointer.timestep_count):
@@ -58,6 +59,9 @@
     else:
       return None
 
+  def data(self):
+    return self.coeff
+
 class ScalarParameter(DolfinAdjointParameter):
   '''This Parameter is used as input to the tangent linear model (TLM)
   when one wishes to compute dJ/da, where a is a single scalar parameter.'''
@@ -136,6 +140,9 @@
     else:
       return None
 
+  def data(self):
+    return self.a
+
 class ScalarParameters(DolfinAdjointParameter):
   '''This Parameter is used as input to the tangent linear model (TLM)
   when one wishes to compute dJ/dv . delta v, where v is a vector of scalar parameters.'''
@@ -194,3 +201,6 @@
       dJdv[i] = out
 
     return dJdv
+
+  def data(self):
+    return self.v

=== added file 'dolfin_adjoint/reduced_functional.py'
--- dolfin_adjoint/reduced_functional.py	1970-01-01 00:00:00 +0000
+++ dolfin_adjoint/reduced_functional.py	2012-08-29 10:56:40 +0000
@@ -0,0 +1,79 @@
+import libadjoint
+from dolfin_adjoint import adjlinalg, adjrhs, constant 
+from dolfin_adjoint.adjglobals import adjointer
+
+class DummyEquation(object):
+    pass
+
+class ReducedFunctional(object):
+    def __init__(self, functional, parameter):
+        ''' Creates a reduced functional object, that evaluates the functional value for a given parameter value.
+            The arguments are as follows:
+            * 'functional' must be a dolfin_adjoint.Functional object. 
+            * 'parameter' must be a single or a list of dolfin_adjoint.DolfinAdjointParameter objects.
+            '''
+        self.functional = functional
+        if not isinstance(parameter, (list, tuple)):
+            parameter = [parameter]
+        self.parameter = parameter
+        # This flag indicates if the functional evaluation is based on replaying the forward annotation. 
+        self.replays_annotation = True
+        self.eqns = []
+
+    def eval_callback(self, value):
+        ''' This function is called before the reduced functional is evaluated.
+            It is intended to be overwritten by the user, for example to plot the control values 
+            that are passed into the callback as "value". ''' 
+        pass
+
+    def __call__(self, value):
+        ''' Evaluates the reduced functional for the given parameter value, by replaying the forward model.
+            Note: before using this evaluation, make sure that the forward model has been annotated. '''
+
+        self.eval_callback(value)
+        if not isinstance(value, (list, tuple)):
+            value = [value]
+        if len(value) != len(self.parameter):
+            raise ValueError, "The number of parameters must equal the number of parameter values."
+
+        # Update the parameter values
+        for i in range(len(value)):
+            if type(value[i]) == constant.Constant:
+                # Constants are not duplicated in the annotation. That is, changing a constant that occurs
+                # in the forward model will also change the forward replay with libadjoint.
+                # However, this is not the case for functions...
+                pass
+            elif hasattr(value[i], 'vector'):
+                # ... since these are duplicated and then occur as rhs in the annotation. 
+                # Therefore, we need to update the right hand side callbacks for
+                # the equation that targets the associated variable.
+
+                # Create a RHS object with the new control values
+                init_rhs = adjlinalg.Vector(value[i]).duplicate()
+                init_rhs.axpy(1.0, adjlinalg.Vector(value[i]))
+                rhs = adjrhs.RHS(init_rhs)
+                # Register the new rhs in the annotation
+                eqn = DummyEquation() 
+                eqn_nb = self.parameter[i].var.equation_nb(adjointer)
+                eqn.equation = adjointer.adjointer.equations[eqn_nb]
+                # Store the equation as a class variable in order to keep a python reference in the memory
+                self.eqns.append(eqn)
+                rhs.register(self.eqns[-1])
+            else:
+                raise NotImplementedError, "The ReducedFunctional class currently only works for parameters that are Functions"
+
+
+        # Replay the annotation and evaluate the functional
+        func_value = 0.
+        for i in range(adjointer.equation_count):
+            (fwd_var, output) = adjointer.get_forward_solution(i)
+
+            storage = libadjoint.MemoryStorage(output)
+            storage.set_overwrite(True)
+            adjointer.record_variable(fwd_var, storage)
+            if i == adjointer.timestep_end_equation(fwd_var.timestep):
+                func_value += adjointer.evaluate_functional(self.functional, fwd_var.timestep)
+
+            #adjglobals.adjointer.forget_forward_equation(i)
+        return func_value
+

=== modified file 'dolfin_adjoint/solving.py'
--- dolfin_adjoint/solving.py	2012-06-29 10:25:47 +0000
+++ dolfin_adjoint/solving.py	2012-08-29 10:56:40 +0000
@@ -118,6 +118,8 @@
   # /before/ we map the coefficients -> dependencies,
   # so that libadjoint records the dependencies with the right timestep number.
   if not linear:
+    # Register the initial condition before the first nonlinear solve
+    register_initial_conditions([[u, adjglobals.adj_variables[u]],], linear=False)
     var = adjglobals.adj_variables.next(u)
   else:
     var = None

=== modified file 'dolfin_adjoint/ui.py'
--- dolfin_adjoint/ui.py	2012-07-13 18:23:30 +0000
+++ dolfin_adjoint/ui.py	2012-08-29 10:56:40 +0000
@@ -17,3 +17,5 @@
 from constant import Constant
 from unimplemented import *
 from timeforms import dt, TimeMeasure, START_TIME, FINISH_TIME
+from reduced_functional import ReducedFunctional
+from optimization import minimize, maximize, print_optimization_algorithms

=== modified file 'dolfin_adjoint/utils.py'
--- dolfin_adjoint/utils.py	2012-07-03 16:25:55 +0000
+++ dolfin_adjoint/utils.py	2012-08-29 10:56:40 +0000
@@ -321,7 +321,14 @@
   return min(convergence_order(with_gradient))
 
 def compute_gradient(J, param, forget=True):
-  dJdparam = None
+  try:
+    scalar = False
+    dJdparam = [None for i in range(len(param))]
+    lparam = param
+  except TypeError:
+    scalar = True
+    dJdparam = [None]
+    lparam = [param]
   last_timestep = adjglobals.adjointer.timestep_count
 
   for i in range(adjglobals.adjointer.timestep_count):
@@ -334,14 +341,15 @@
     adjglobals.adjointer.record_variable(adj_var, storage)
     fwd_var = libadjoint.Variable(adj_var.name, adj_var.timestep, adj_var.iteration)
 
-    out = param.inner_adjoint(adjglobals.adjointer, output.data, i, fwd_var)
-    dJdparam = _add(dJdparam, out)
+    for j in range(len(lparam)):
+      out = lparam[j].inner_adjoint(adjglobals.adjointer, output.data, i, fwd_var)
+      dJdparam[j] = _add(dJdparam[j], out)
 
-    if last_timestep > adj_var.timestep:
-      # We have hit a new timestep, and need to compute this timesteps' \partial J/\partial m contribution
-      last_timestep = adj_var.timestep
-      out = param.partial_derivative(adjglobals.adjointer, J, adj_var.timestep)
-      dJdparam = _add(dJdparam, out)
+      if last_timestep > adj_var.timestep:
+        # We have hit a new timestep, and need to compute this timesteps' \partial J/\partial m contribution
+        last_timestep = adj_var.timestep
+        out = lparam[j].partial_derivative(adjglobals.adjointer, J, adj_var.timestep)
+        dJdparam[j] = _add(dJdparam[j], out)
 
     if forget is None:
       pass
@@ -350,7 +358,10 @@
     else:
       adjglobals.adjointer.forget_adjoint_values(i)
 
-  return dJdparam
+  if scalar:
+    return dJdparam[0]
+  else:
+    return dJdparam
 
 def test_scalar_parameter_adjoint(J, a, dJda, seed=None):
   info_blue("Running Taylor remainder convergence analysis for the adjoint model ... ")

=== modified file 'tests/burgers_oo/burgers_oo.py'
--- tests/burgers_oo/burgers_oo.py	2012-07-02 18:37:16 +0000
+++ tests/burgers_oo/burgers_oo.py	2012-08-29 10:56:40 +0000
@@ -76,15 +76,15 @@
     print "Running adjoint ... "
 
     J = Functional(forward*forward*dx*dt[FINISH_TIME])
-    for (adjoint, var) in compute_adjoint(J, forget=False):
-      pass
+    dJdm = compute_gradient(J, InitialConditionParameter(forward), forget = False)
 
     def Jfunc(ic):
       forward = main(ic, annotate=False)
       return assemble(forward*forward*dx)
 
-    minconv = test_initial_condition_adjoint(Jfunc, ic, adjoint, seed=1.0e-5)
+    minconv = test_initial_condition_adjoint(Jfunc, ic, dJdm, seed=1.0e-5)
     if minconv < 1.9:
+      info_red("Test failed. Convergence rate is %f < 1.9", minconv)
       sys.exit(1)
 
     dJ = assemble(derivative(forward_copy*forward_copy*dx, forward_copy))
@@ -93,4 +93,5 @@
     ic.vector()[:] = ic_copy.vector()
     minconv = test_initial_condition_tlm(Jfunc, dJ, ic, seed=1.0e-5)
     if minconv < 1.9:
+      info_red("Test failed. Convergence rate is %d < 1.9", minconv)
       sys.exit(1)

=== modified file 'tests/optimal_control_mms/optimal_control_mms.py'
--- tests/optimal_control_mms/optimal_control_mms.py	2012-07-14 09:16:08 +0000
+++ tests/optimal_control_mms/optimal_control_mms.py	2012-08-29 10:56:40 +0000
@@ -1,5 +1,6 @@
 """ Solves a MMS problem with smooth control """
 
+import sys
 from dolfin import *
 from dolfin_adjoint import *
 
@@ -24,14 +25,15 @@
     u_d = 1/(2*pi**2)*sin(pi*x[0])*sin(pi*x[1]) 
 
     J = Functional((inner(u-u_d, u-u_d))*dx*dt[FINISH_TIME])
-    def Jfunc(m):
-      solve_pde(u, V, m)
-      return assemble(inner(u-u_d, u-u_d)*dx)
+
+    # Run the forward model once to create the annotation
+    solve_pde(u, V, m)
 
     # Run the optimisation 
-    optimisation.minimise(Jfunc, J, InitialConditionParameter(m), m, algorithm = 'scipy.l_bfgs_b', pgtol=1e-16, factr=1, bounds = (-1, 1), iprint = 1, maxfun = 20)
-    #optimisation.minimise(Jfunc, J, InitialConditionParameter(m), m, algorithm = 'scipy.slsqp', bounds = (-1, 1), iprint = 3, iter = 60)
-    Jfunc(m)
+    reduced_func = ReducedFunctional(J, InitialConditionParameter(m))
+    minimize(reduced_func, algorithm = 'scipy.l_bfgs_b', pgtol=1e-16, factr=1, bounds = (-1, 1), iprint = 1, maxfun = 20)
+    #minimize(reduced_func, algorithm = 'scipy.slsqp', bounds = (-1, 1), iprint = 3, iter = 60)
+    solve_pde(u, V, m)
 
     m_analytic = sin(pi*x[0])*sin(pi*x[1]) 
     u_analytic = 1/(2*pi**2)*sin(pi*x[0])*sin(pi*x[1])
@@ -49,6 +51,7 @@
     control_errors.append(control_error)
     state_errors.append(state_error)
     element_sizes.append(1./n)
+    adj_reset()
 
 info_green("Control errors: " + str(control_errors))
 info_green("Control convergence: " + str(convergence_order(control_errors, base = 2)))
@@ -56,7 +59,9 @@
 info_green("State convergence: " + str(convergence_order(state_errors, base = 2)))
 
 if min(convergence_order(control_errors)) < 2.0:
+    info_red("Convergence order below tolerance") 
     sys.exit(1)
 if min(convergence_order(state_errors)) < 4.0:
+    info_red("Convergence order below tolerance") 
     sys.exit(1)
 info_green("Test passed")    

=== renamed directory 'tests/optimisation' => 'tests/optimization'
=== renamed file 'tests/optimisation/optimisation.py' => 'tests/optimization/optimization.py'
--- tests/optimisation/optimisation.py	2012-07-03 16:25:55 +0000
+++ tests/optimization/optimization.py	2012-08-29 10:56:40 +0000
@@ -4,9 +4,10 @@
 
 from dolfin import *
 from dolfin_adjoint import *
+import libadjoint
 
 dolfin.set_log_level(ERROR)
-dolfin.parameters["optimisation"]["test_gradient"] = True 
+dolfin.parameters["optimization"]["test_gradient"] = True 
 
 n = 10
 mesh = UnitInterval(n)
@@ -30,12 +31,13 @@
 
     t = 0.0
     end = 0.2
+    adjointer.time.start(t)
     while (t <= end):
         solve(F == 0, u_next, bc, annotate=annotate)
         u.assign(u_next, annotate=annotate)
 
         t += float(timestep)
-        adj_inc_timestep()
+        adj_inc_timestep(time=t, finished = t>end)
 
 if __name__ == "__main__":
 
@@ -43,21 +45,29 @@
     u = Function(ic, name='Velocity')
 
     J = Functional(u*u*dx*dt[FINISH_TIME])
-    def Jfunc(ic):
-      u.assign(ic)
-      main(u, annotate=True)
-      return assemble(u*u*dx)
+
+    # Run the model once to create the annotation
+    u.assign(ic)
+    main(u, annotate=True)
 
     # Run the optimisation 
     lb = project(Expression("-1"),  V)
-    optimisation.minimise(Jfunc, J, InitialConditionParameter(u), ic, algorithm = 'scipy.l_bfgs_b', pgtol=1e-6, factr=1e5, bounds = (lb, 1), iprint = 1)
+    
+    # Define the reduced funtional
+    reduced_functional = ReducedFunctional(J, InitialConditionParameter(u))
+
+    # Run the optimisation problem with gradient tests and L-BFGS-B
+    u_opt = minimize(reduced_functional, algorithm = 'scipy.l_bfgs_b', pgtol=1e-6, factr=1e5, bounds = (lb, 1), iprint = 1)
     ic = project(Expression("sin(2*pi*x[0])"),  V)
 
-    # For performance reasons, switch the gradient test off
-    dolfin.parameters["optimisation"]["test_gradient"] = False 
-    optimisation.minimise(Jfunc, J, InitialConditionParameter(u), ic, algorithm = 'scipy.slsqp', bounds = (lb, 1), iprint = 2, acc = 1e-10)
+    # Run the problem again with SQP, this time for performance reasons with the gradient test switched off
+    dolfin.parameters["optimization"]["test_gradient"] = False 
+    u_opt = minimize(reduced_functional, algorithm = 'scipy.slsqp', bounds = (lb, 1), iprint = 2, acc = 1e-10)
 
     tol = 1e-9
-    if Jfunc(ic) > tol:
-        print 'Test failed: Optimised functional value exceeds tolerance: ' , Jfunc(ic), ' > ', tol, '.'
+    final_functional = reduced_functional(u_opt)
+    print "Final functional value: ", final_functional
+    if final_functional > tol:
+        print 'Test failed: Optimised functional value exceeds tolerance: ' , final_functional, ' > ', tol, '.'
         sys.exit(1)
+

=== renamed directory 'tests/optimisation_scalar' => 'tests/optimization_scalar'
=== renamed file 'tests/optimisation_scalar/optimisation_scalar.py' => 'tests/optimization_scalar/optimization_scalar.py'
--- tests/optimisation_scalar/optimisation_scalar.py	2012-07-02 18:37:16 +0000
+++ tests/optimization_scalar/optimization_scalar.py	2012-08-29 10:56:40 +0000
@@ -3,7 +3,7 @@
 import sys
 
 dolfin.set_log_level(ERROR)
-dolfin.parameters["optimisation"]["test_gradient"] = True
+dolfin.parameters["optimization"]["test_gradient"] = True
 
 n = 10
 mesh = UnitInterval(n)
@@ -33,19 +33,16 @@
 
 if __name__ == "__main__":
   nu = Constant(0.0001)
+  # Run the forward model once to have the annotation
   main(nu)
 
   J = Functional(inner(u, u)*dx*dt[FINISH_TIME])
 
-  def Jhat(nu): 
-    u.assign(ic)
-    main(nu)
-    return assemble(inner(u, u)*dx)
-
   # Run the optimisation 
-  optimisation.minimise(Jhat, J, ScalarParameter(nu), nu, 'scipy.slsqp', iprint = 2)
+  reduced_functional = ReducedFunctional(J, ScalarParameter(nu))
+  nu_opt = minimize(reduced_functional, 'scipy.slsqp', iprint = 2)
 
   tol = 1e-4
-  if Jhat(nu) > tol:
-    print 'Test failed: Optimised functional value exceeds tolerance: ' , Jhat(nu), ' > ', tol, '.'
+  if reduced_functional(nu_opt) > tol:
+    print 'Test failed: Optimised functional value exceeds tolerance: ', reduced_functional(nu_opt), ' > ', tol, '.'
     sys.exit(1)

=== added directory 'tests/reduced_functional_evaluation'
=== added file 'tests/reduced_functional_evaluation/reduced_functional_evaluation.py'
--- tests/reduced_functional_evaluation/reduced_functional_evaluation.py	1970-01-01 00:00:00 +0000
+++ tests/reduced_functional_evaluation/reduced_functional_evaluation.py	2012-08-29 10:56:40 +0000
@@ -0,0 +1,54 @@
+''' A simple test that compares the functional value computed manually and with libadjoints functional_evaluation.
+    Writting this test was motivated by the bug described on https://bugs.launchpad.net/dolfin-adjoint/+bug/1032291 '''
+from dolfin import *
+from dolfin_adjoint import *
+import sys
+
+# Global settings
+set_log_level(ERROR)
+
+mesh = UnitSquare(10, 10) 
+V = FunctionSpace(mesh, "DG", 1)
+
+u_new = Function(V, name = "u_new")
+u_old = Function(V, name = "u_old")
+u_test = TestFunction(V)
+
+T = 2.
+t = 0.
+dlt = 1.
+F1 = ( inner((u_new - u_old)/dlt, u_test)*dx - inner(Constant(1.), u_test)*dx )
+#solve(inner(u_new, u_test)*dx == 0, u_new)
+
+adjointer.time.start(t)
+man_func_value = 0.
+print "+++++++++++++ INITIAL RUN +++++++++" 
+man_func_value_contr = 0.5*assemble(inner(u_new, u_new)*dx)
+while t < T:
+
+  solve(F1 == 0, u_new)
+  u_old.assign(u_new)
+
+  t += dlt 
+  man_func_value_contr = assemble(inner(u_new, u_new)*dx)
+  if t>=T:
+    man_func_value += 0.5*man_func_value_contr
+  else:
+    man_func_value += man_func_value_contr
+  adj_inc_timestep(time=t, finished = t>=T)
+
+info_green("Manually computed functional value: %f", man_func_value)
+adj_html("forward.html", "forward")
+print
+print "+++++++++++++ REPLAY +++++++++" 
+u_new.vector()[:] = 0.
+u_old.vector()[:] = 0.
+J = Functional(inner(u_new, u_new)*dx*dt) 
+reduced_functional = ReducedFunctional(J, InitialConditionParameter(u_old))
+reduced_functional_value = reduced_functional(u_old)
+info_green("Functional value from reduced functional: %f", reduced_functional_value)
+
+if abs(reduced_functional_value - man_func_value) > 1e-13:
+  info_red("Test failed. Error: %f", abs(reduced_functional_value - man_func_value))
+else:
+  info_green("Test passed")

=== modified file 'tests/test.py'
--- tests/test.py	2012-07-14 17:27:55 +0000
+++ tests/test.py	2012-08-29 10:56:40 +0000
@@ -13,7 +13,7 @@
              'navier_stokes': 'mpirun -n 2 python navier_stokes.py',
              'svd_simple': 'mpirun -n 2 python svd_simple.py',
              'gst_mass': 'mpirun -n 2 python gst_mass.py',
-             'optimisation': 'mpirun -n 2 python optimisation.py',
+             'optimization': 'mpirun -n 2 python optimization.py',
              'optimal_control_mms': 'mpirun -n 2 python optimal_control_mms.py',
              'differentiability-dg-upwind': None,
              'differentiability-stokes': None,