optimize.py

""" Optimization (minimization) functions """
from __future__ import division, print_function, absolute_import, unicode_literals
#***************************************************************************************************
# Copyright 2015, 2019 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
# Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights
# in this software.
# Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except
# in compliance with the License.  You may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0 or in the LICENSE file in the root pyGSTi directory.
#***************************************************************************************************

import numpy as _np
import time as _time
import sys as _sys
import os as _os
import scipy.optimize as _spo

try:
    from scipy.optimize import Result as _optResult  # for earlier scipy versions
except:
    from scipy.optimize import OptimizeResult as _optResult  # for later scipy versions

from .customcg import fmax_cg


def minimize(fn, x0, method='cg', callback=None,
             tol=1e-10, maxiter=1000000, maxfev=None,
             stopval=None, jac=None):
    """
    Minimizes the function fn starting at x0.

    This is a gateway function to all other minimization routines within this
    module, providing a common interface to many different minimization methods
    (including and extending beyond those available from scipy.optimize).

    Parameters
    ----------
    fn : function
        The function to minimize.

    x0 : numpy array
        The starting point (argument to fn).

    method : string, optional
        Which minimization method to use.  Allowed values are:
        "simplex" : uses fmin_simplex
        "supersimplex" : uses fmin_supersimplex
        "customcg" : uses fmax_cg (custom CG method)
        "brute" : uses scipy.optimize.brute
        "basinhopping" : uses scipy.optimize.basinhopping with L-BFGS-B
        "swarm" : uses fmin_particle_swarm
        "evolve" : uses fmin_evolutionary (which uses DEAP)
        < methods available from scipy.optimize.minimize >

    callback : function, optional
        A callback function to be called in order to track optimizer progress.
        Should have signature: myCallback(x, f=None, accepted=None).  Note that
        create_obj_func_printer(...) function can be used to create a callback.

    tol : float, optional
        Tolerance value used for all types of tolerances available in a given method.

    maxiter : int, optional
        Maximum iterations.

    maxfev : int, optional
        Maximum function evaluations; used only when available, and defaults to maxiter.

    stopval : float, optional
        For basinhopping method only.  When f <= stopval then basinhopping outer loop
        will terminate.  Useful when a bound on the minimum is known.

    jac : function
        Jacobian function.

    Returns
    -------
    scipy.optimize.Result object
        Includes members 'x', 'fun', 'success', and 'message'.
    """

    if maxfev is None: maxfev = maxiter

    #Run Minimization Algorithm
    if method == 'simplex':
        solution = fmin_simplex(fn, x0, slide=1.0, tol=tol, maxiter=maxiter)

    elif method == 'supersimplex':
        solution = fmin_supersimplex(fn, x0, outer_tol=1.0, inner_tol=tol,
                                     max_outer_iter=100, min_inner_maxiter=100, max_inner_maxiter=maxiter)

    elif method == 'customcg':
        def fn_to_max(x):
            """ Function to maximize """
            f = fn(x); return -f if f is not None else None

        if jac is not None:
            def dfdx_and_bdflag(x):
                """ Returns derivative and boundary flag """
                j = -jac(x)
                bd = _np.zeros(len(j))  # never say fn is on boundary, since this is an analytic derivative
                return j, bd
        else:
            dfdx_and_bdflag = None

        # Note: even though we maximize, return value is negated to conform to min routines
        solution = fmax_cg(fn_to_max, x0, maxiter, tol, dfdx_and_bdflag, None)

    elif method == 'brute':
        ranges = [(0.0, 1.0)] * len(x0); Ns = 4  # params for 'brute' algorithm
        xmin, _ = _spo.brute(fn, ranges, (), Ns)  # jac=jac
        #print "DEBUG: Brute fmin = ",fmin
        solution = _spo.minimize(fn, xmin, method="Nelder-Mead", options={}, tol=tol, callback=callback, jac=jac)

    elif method == 'basinhopping':
        def _basin_callback(x, f, accept):
            if callback is not None: callback(x, f=f, accepted=accept)
            if stopval is not None and f <= stopval:
                return True  # signals basinhopping to stop
            return False
        solution = _spo.basinhopping(fn, x0, niter=maxiter, T=2.0, stepsize=1.0,
                                     callback=_basin_callback, minimizer_kwargs={'method': "L-BFGS-B", 'jac': jac})

        #DEBUG -- follow with Nelder Mead to make sure basinhopping found a minimum. (It seems to)
        #print "DEBUG: running Nelder-Mead:"
        #opts = { 'maxfev': maxiter, 'maxiter': maxiter }
        #solution = _spo.minimize(fn, solution.x, options=opts, method="Nelder-Mead", tol=1e-8, callback=callback)
        #print "DEBUG: done: best f = ",solution.fun

        solution.success = True  # basinhopping doesn't seem to set this...

    elif method == 'swarm':
        solution = fmin_particle_swarm(fn, x0, tol, maxiter, popsize=1000)  # , callback = callback)

    elif method == 'evolve':
        solution = fmin_evolutionary(fn, x0, num_generations=maxiter, num_individuals=500)

#    elif method == 'homebrew':
#      solution = fmin_homebrew(fn, x0, maxiter)

    else:
        #Set options for different algorithms
        opts = {'maxiter': maxiter, 'disp': False}
        if method == "BFGS": opts['gtol'] = tol  # gradient norm tolerance
        elif method == "L-BFGS-B": opts['gtol'] = opts['ftol'] = tol  # gradient norm and fractional y-tolerance
        elif method == "Nelder-Mead": opts['maxfev'] = maxfev  # max fn evals (note: ftol and xtol can also be set)

        if method in ("BFGS", "CG", "Newton-CG", "L-BFGS-B", "TNC", "SLSQP", "dogleg", "trust-ncg"):  # use jacobian
            solution = _spo.minimize(fn, x0, options=opts, method=method, tol=tol, callback=callback, jac=jac)
        else:
            solution = _spo.minimize(fn, x0, options=opts, method=method, tol=tol, callback=callback)

    return solution


def fmin_supersimplex(fn, x0, outer_tol, inner_tol, max_outer_iter, min_inner_maxiter, max_inner_maxiter):
    """
    Minimize a function using repeated applications of the simplex algorithm.

    By varying the maximum number of iterations and repeatedly calling scipy's
    Nelder-Mead simplex optimization, this function performs as a robust (but
    slow) minimization.

    Parameters
    ----------
    fn : function
        The function to minimize.

    x0 : numpy array
        The starting point (argument to fn).

    outer_tol : float
        Tolerance of outer loop

    inner_tol : float
        Tolerance fo inner loop

    max_outer_iter : int
        Maximum number of outer-loop iterations

    min_inner_maxiter : int
        Minimum number of inner-loop iterations

    max_inner_maxiter : int
        Maxium number of outer-loop iterations

    Returns
    -------
    scipy.optimize.Result object
        Includes members 'x', 'fun', 'success', and 'message'.
    """
    f_init = fn(x0)
    f_final = f_init - 10 * outer_tol  # prime the loop
    x_start = x0

    i = 1
    cnt_at_same_maxiter = 1
    inner_maxiter = min_inner_maxiter

    while (f_init - f_final > outer_tol or inner_maxiter < max_inner_maxiter) and i < max_outer_iter:
        if f_init - f_final <= outer_tol and inner_maxiter < max_inner_maxiter:
            inner_maxiter *= 10; cnt_at_same_maxiter = 1
        if cnt_at_same_maxiter > 10 and inner_maxiter > min_inner_maxiter:
            inner_maxiter /= 10; cnt_at_same_maxiter = 1
        f_init = f_final

        print(">>> fmin_supersimplex: outer iteration %d (inner_maxiter = %d)" % (i, inner_maxiter))
        i += 1; cnt_at_same_maxiter += 1

        opts = {'maxiter': inner_maxiter, 'maxfev': inner_maxiter, 'disp': False}
        inner_solution = _spo.minimize(fn, x_start, options=opts, method='Nelder-Mead', callback=None, tol=inner_tol)

        if not inner_solution.success:
            print("WARNING: fmin_supersimplex inner loop failed (tol=%g, maxiter=%d): %s"
                  % (inner_tol, inner_maxiter, inner_solution.message))

        f_final = inner_solution.fun
        x_start = inner_solution.x
        print(">>> fmin_supersimplex: outer iteration %d gives min = %f" % (i, f_final))

    solution = _optResult()
    solution.x = inner_solution.x
    solution.fun = inner_solution.fun
    if i < max_outer_iter:
        solution.success = True
    else:
        solution.success = False
        solution.message = "Maximum iterations exceeded"
    return solution


def fmin_simplex(fn, x0, slide=1.0, tol=1e-8, maxiter=1000):
    """
    Minimizes a function using a custom simplex implmentation.

    This was used primarily to check scipy's Nelder-Mead method
    and runs much slower, so there's not much reason for using
    this method.

    Parameters
    ----------
    fn : function
        The function to minimize.

    x0 : numpy array
        The starting point (argument to fn).

    slide : float, optional
        Affects initial simplex point locations

    tol : float, optional
        Relative tolerance as a convergence criterion.

    maxiter : int, optional
        Maximum iterations.

    Returns
    -------
    scipy.optimize.Result object
        Includes members 'x', 'fun', 'success', and 'message'.
    """

    # Setup intial values
    n = len(x0)
    f = _np.zeros(n + 1)
    x = _np.zeros((n + 1, n))

    x[0] = x0

    # Setup intial X range
    for i in range(1, n + 1):
        x[i] = x0
        x[i, i - 1] = x0[i - 1] + slide

    # Setup intial functions based on x's just defined
    for i in range(n + 1):
        f[i] = fn(x[i])

    # Main Loop operation, loops infinitly until break condition
    counter = 0
    while True:
        low = _np.argmin(f)
        high = _np.argmax(f)
        counter += 1

        # Compute Migration
        d = (-(n + 1) * x[high] + sum(x)) / n

        # Break if value is close
        if _np.sqrt(_np.dot(d, d) / n) < tol or counter == maxiter:
            solution = _optResult()
            solution.x = x[low]; solution.fun = f[low]
            if counter < maxiter:
                solution.success = True
            else:
                solution.success = False
                solution.message = "Maximum iterations exceeded"
            return solution

        newX = x[high] + 2.0 * d
        newF = fn(newX)

        if newF <= f[low]:
            # Bad news, new value is lower than any other point => replace high point with new values
            x[high] = newX
            f[high] = newF
            newX = x[high] + d
            newF = fn(newX)

            # Check if need to expand
            if newF <= f[low]:
                x[high] = newX
                f[high] = newF

        else:
            # Good news, new value is higher than lowest point

            # Check if need to contract
            if newF <= f[high]:
                x[high] = newX
                f[high] = newF
            else:
                # Contraction
                newX = x[high] + 0.5 * d
                newF = fn(newX)
                if newF <= f[high]:
                    x[high] = newX
                    f[high] = newF
                else:
                    for i in range(len(x)):
                        if i != low:
                            x[i] = (x[i] - x[low])
                            f[i] = fn(x[i])


#TODO err_crit is never used?
def fmin_particle_swarm(f, x0, err_crit, iter_max, popsize=100, c1=2, c2=2):
    """
    A simple implementation of the Particle Swarm Optimization Algorithm.
    Pradeep Gowda 2009-03-16

    Parameters
    ----------
    f : function
        The function to minimize.

    x0 : numpy array
        The starting point (argument to fn).

    err_crit : float
        Critical error (i.e. tolerance).  Stops when error < err_crit.

    iter_max : int
        Maximum iterations.

    popsize : int, optional
        Population size.  Larger populations are better at finding the global
        optimum but make the algorithm take longer to run.

    c1 : float, optional
        Coefficient describing a particle's affinity for it's (local) maximum.

    c2 : float, optional
        Coefficient describing a particle's affinity for the best maximum any
        particle has seen (the current global max).

    Returns
    -------
    scipy.optimize.Result object
        Includes members 'x', 'fun', 'success', and 'message'.
    """
    dimensions = len(x0)
    LARGE = 1e10

    class Particle:
        """ Particle "container" class """
        pass

    #initialize the particles
    particles = []
    for i in range(popsize):
        p = Particle()
        p.params = x0 + 2 * (_np.random.random(dimensions) - 0.5)
        p.best = p.params[:]
        p.fitness = LARGE  # large == bad fitness
        p.v = _np.zeros(dimensions)
        particles.append(p)

    # let the first particle be the global best
    gbest = particles[0]; ibest = 0
    # bDoLocalFitnessOpt = False

    #DEBUG
    #if False:
    #    import pickle as _pickle
    #    bestGaugeMx = _pickle.load(open("bestGaugeMx.debug"))
    #    lbfgsbGaugeMx = _pickle.load(open("lbfgsbGaugeMx.debug"))
    #    cgGaugeMx = _pickle.load(open("cgGaugeMx.debug"))
    #    initialGaugeMx = x0.reshape( (4,4) )
    #
    #    #DEBUG: dump line cut to plot
    #    nPts = 100
    #    print "DEBUG: best offsets = \n", bestGaugeMx - initialGaugeMx
    #    print "DEBUG: lbfgs offsets = \n", lbfgsbGaugeMx - initialGaugeMx
    #    print "DEBUG: cg offsets = \n", cgGaugeMx - initialGaugeMx
    #
    #    print "# DEBUG plot"
    #    #fDebug = open("x0ToBest.dat","w")
    #    #fDebug = open("x0ToLBFGS.dat","w")
    #    fDebug = open("x0ToCG.dat","w")
    #    #fDebug = open("LBFGSToBest.dat","w")
    #    #fDebug = open("CGToBest.dat","w")
    #    #fDebug = open("CGToLBFGS.dat","w")
    #
    #    for i in range(nPts+1):
    #        alpha = float(i) / nPts
    #        #matM = (1.0-alpha) * initialGaugeMx + alpha*bestGaugeMx
    #        #matM = (1.0-alpha) * initialGaugeMx + alpha*lbfgsbGaugeMx
    #        matM = (1.0-alpha) * initialGaugeMx + alpha*cgGaugeMx
    #        #matM = (1.0-alpha) * lbfgsbGaugeMx + alpha*bestGaugeMx
    #        #matM = (1.0-alpha) * cgGaugeMx + alpha*bestGaugeMx
    #        #matM = (1.0-alpha) * cgGaugeMx + alpha*lbfgsbGaugeMx
    #        print >> fDebug, "%g %g" % (alpha, f(matM.flatten()))
    #    exit()
    #
    #
    #    fDebug = open("lineDataFromX0.dat","w")
    #    min_offset = -1; max_offset = 1
    #    for i in range(nPts+1):
    #        offset = min_offset + float(i)/nPts * (max_offset-min_offset)
    #        print >> fDebug, "%g" % offset,
    #
    #        for k in range(len(x0)):
    #            x = x0.copy(); x[k] += offset
    #            try:
    #                print >> fDebug, " %g" % f(x),
    #            except:
    #                print >> fDebug, " nan",
    #        print >> fDebug, ""
    #
    #    print >> fDebug, "#END DEBUG plot"
    #    exit()
    #END DEBUG

    #err = 1e10
    for iter_num in range(iter_max):
        w = 1.0  # - i/iter_max

        #bDoLocalFitnessOpt = bool(iter_num > 20 and abs(lastBest-gbest.fitness) < 0.001 and iter_num % 10 == 0)
        # lastBest = gbest.fitness
        # minDistToBest = 1e10; minV = 1e10; maxV = 0 #DEBUG

        for (ip, p) in enumerate(particles):
            fitness = f(p.params)

            #if bDoLocalFitnessOpt:
            #    opts = {'maxiter': 100, 'maxfev': 100, 'disp': False }
            #    local_soln = _spo.minimize(f,p.params,options=opts, method='L-BFGS-B',callback=None, tol=1e-2)
            #    p.params = local_soln.x
            #    fitness = local_soln.fun

            if fitness < p.fitness:  # low 'fitness' is good b/c we're minimizing
                p.fitness = fitness
                p.best = p.params

            if fitness < gbest.fitness:
                gbest = p; ibest = ip

            v = w * p.v + c1 * _np.random.random() * (p.best - p.params) \
                + c2 * _np.random.random() * (gbest.params - p.params)
            p.params = p.params + v
            for (i, pv) in enumerate(p.params):
                p.params[i] = ((pv + 1) % 2) - 1  # periodic b/c on box between -1 and 1

            #from .. import tools as tools_
            #matM = p.params.reshape( (4,4) )  #DEBUG
            #minDistToBest = min(minDistToBest, _tools.frobeniusdist(
            #                                    bestGaugeMx,matM)) #DEBUG
            #minV = min( _np.linalg.norm(v), minV)
            #maxV = max( _np.linalg.norm(v), maxV)

        #print "DB: min diff from best = ", minDistToBest #DEBUG
        #print "DB: min,max v = ", (minV,maxV)

        #if False: #bDoLocalFitnessOpt:
        #    opts = {'maxiter': 100, 'maxfev': 100, 'disp': False }
        #    print "initial fun = ",gbest.fitness,
        #    local_soln = _spo.minimize(f,gbest.params,options=opts, method='L-BFGS-B',callback=None, tol=1e-5)
        #    gbest.params = local_soln.x
        #    gbest.fitness = local_soln.fun
        #    print "  final fun = ",gbest.fitness

        print("Iter %d: global best = %g (index %d)" % (iter_num, gbest.fitness, ibest))

        #if err < err_crit:  break  #TODO: stopping condition

    ## Uncomment to print particles
    #for p in particles:
    #    print 'params: %s, fitness: %s, best: %s' % (p.params, p.fitness, p.best)

    solution = _optResult()
    solution.x = gbest.params; solution.fun = gbest.fitness
    solution.success = True
#    if iter_num < maxiter:
#        solution.success = True
#    else:
#        solution.success = False
#        solution.message = "Maximum iterations exceeded"
    return solution


def fmin_evolutionary(f, x0, num_generations, num_individuals):
    """
    Minimize a function using an evolutionary algorithm.

    Uses python's deap package to perform an evolutionary
    algorithm to find a function's global minimum.

    Parameters
    ----------
    fn : function
        The function to minimize.

    x0 : numpy array
        The starting point (argument to fn).

    num_generations : int
        The number of generations to carry out. (similar to the number
        of iterations)

    num_individuals : int
        The number of individuals in each generation.  More individuals
        make finding the global optimum more likely, but take longer
        to run.


    Returns
    -------
    scipy.optimize.Result object
        Includes members 'x', 'fun', 'success', and 'message'.
    """

    import deap.creator as _creator
    import deap.base as _base
    import deap.tools as _tools
    numParams = len(x0)

    # Create the individual class
    _creator.create("FitnessMin", _base.Fitness, weights=(-1.0,))
    _creator.create("Individual", list, fitness=_creator.FitnessMin)

    # Create initialization functions
    toolbox = _base.Toolbox()
    toolbox.register("random", _np.random.random)
    toolbox.register("individual", _tools.initRepeat, _creator.Individual,
                     toolbox.random, n=numParams)  # fn to init an individual from a list of numParams random numbers
    # fn to create a population (still need to specify n)
    toolbox.register("population", _tools.initRepeat, list, toolbox.individual)

    # Create operation functions
    def _evaluate(individual):
        return f(_np.array(individual)),  # note: must return a tuple

    toolbox.register("mate", _tools.cxTwoPoint)
    toolbox.register("mutate", _tools.mutGaussian, mu=0, sigma=0.5, indpb=0.1)
    toolbox.register("select", _tools.selTournament, tournsize=3)
    toolbox.register("evaluate", _evaluate)

    # Create the population
    pop = toolbox.population(n=num_individuals)

    # Evaluate the entire population
    fitnesses = list(map(toolbox.evaluate, pop))
    for ind, fit in zip(pop, fitnesses):
        ind.fitness.values = fit

    PROB_TO_CROSS = 0.5
    PROB_TO_MUTATE = 0.2

    # Initialize statistics
    stats = _tools.Statistics(key=lambda ind: ind.fitness.values)
    stats.register("avg", _np.mean)
    stats.register("std", _np.std)
    stats.register("min", _np.min)
    stats.register("max", _np.max)
    logbook = _tools.Logbook()

    #Run algorithm
    for g in range(num_generations):
        record = stats.compile(pop)
        logbook.record(gen=g, **record)
        print("Gen %d: %s" % (g, record))

        # Select the next generation individuals
        offspring = toolbox.select(pop, len(pop))
        # Clone the selected individuals
        offspring = list(map(toolbox.clone, offspring))

        # Apply crossover on the offspring
        for child1, child2 in zip(offspring[::2], offspring[1::2]):
            if _np.random.random() < PROB_TO_CROSS:
                toolbox.mate(child1, child2)
                del child1.fitness.values
                del child2.fitness.values

        # Apply mutation on the offspring
        for mutant in offspring:
            if _np.random.random() < PROB_TO_MUTATE:
                toolbox.mutate(mutant)
                del mutant.fitness.values

        # Evaluate the individuals with an invalid fitness
        invalid_ind = [ind for ind in offspring if not ind.fitness.valid]
        fitnesses = toolbox.map(toolbox.evaluate, invalid_ind)
        for ind, fit in zip(invalid_ind, fitnesses):
            ind.fitness.values = fit

        # The population is entirely replaced by the offspring
        pop[:] = offspring

    #get best individual and return params
    indx_min_fitness = _np.argmin([ind.fitness.values[0] for ind in pop])
    best_params = _np.array(pop[indx_min_fitness])

    solution = _optResult()
    solution.x = best_params; solution.fun = pop[indx_min_fitness].fitness.values[0]
    solution.success = True
    return solution


#def fmin_homebrew(f, x0, maxiter):
#    """
#    Cooked up by Erik, this algorithm is similar to basinhopping but with some tweaks.
#
#    Parameters
#    ----------
#    fn : function
#        The function to minimize.
#
#    x0 : numpy array
#        The starting point (argument to fn).
#
#    maxiter : int
#        The maximum number of iterations.
#
#    Returns
#    -------
#    scipy.optimize.Result object
#        Includes members 'x', 'fun', 'success', and 'message'.
#    """
#
#    STEP = 0.01
#    MAX_STEPS = int(2.0 / STEP) # allow a change of at most 2.0
#    MAX_DIR_TRIES = 1000
#    T = 1.0
#
#    global_best_params = cur_x0 = x0
#    global_best = cur_f = f(x0)
#    N = len(x0)
#    trial_x0 = x0.copy()
#
#    for it in range(maxiter):
#
#        #Minimize using L-BFGS-B
#        opts = {'maxiter': maxiter, 'maxfev': maxiter, 'disp': False }
#        soln = _spo.minimize(f,trial_x0,options=opts, method='L-BFGS-B',callback=None, tol=1e-8)
#
#        # Update global best
#        if soln.fun < global_best:
#            global_best_params = soln.x
#            global_best = soln.fun
#
#        #check if we accept the new minimum
#        if soln.fun < cur_f or _np.random.random() < _np.exp( -(soln.fun - cur_f)/T ):
#            cur_x0 = soln.x; cur_f = soln.fun
#            print "Iter %d: f=%g accepted -- global best = %g" % (it, cur_f, global_best)
#        else:
#            print "Iter %d: f=%g declined" % (it, cur_f)
#
#        trial_x0 = None; numTries = 0
#        while trial_x0 is None and numTries < MAX_DIR_TRIES:
#            #choose a random direction
#            direction = _np.random.random( N )
#            numTries += 1
#
#            #print "DB: test dir %d" % numTries #DEBUG
#
#            #kick solution along random direction until the value of f starts to get smaller again (if it ever does)
#            #  (this indicates we've gone over a maximum along this direction)
#            last_f = cur_f
#            for i in range(1,MAX_STEPS):
#                test_x = cur_x0 + i*STEP * direction
#                test_f = f(test_x)
#                #print "DB: test step=%f: f=%f" % (i*STEP, test_f)
#                if test_f < last_f:
#                    trial_x0 = test_x
#                    print "Found new direction in %d tries, new f(x0) = %g" % (numTries,test_f)
#                    break
#                last_f = test_f
#
#        if trial_x0 is None:
#            raise ValueError("Maximum number of direction tries exceeded")
#
#    solution = _optResult()
#    solution.x = global_best_params; solution.fun = global_best
#    solution.success = True
##    if it < maxiter:
##        solution.success = True
##    else:
##        solution.success = False
##        solution.message = "Maximum iterations exceeded"
#    return solution


def create_obj_func_printer(objFunc, startTime=None):
    """
    Create a callback function that prints the value of an objective function.

    Parameters
    ----------
    objFunc : function
        The objective function to print.

    startTime : float (optional)
        A reference starting time to use when printing elapsed times. If None,
        then the system time when this function is called is used (which is
        often what you want).

    Returns
    -------
    function
        A callback function which prints objFunc.
    """
    if startTime is None:
        startTime = _time.time()  # for reference point of obj func printer

    def print_obj_func(x, f=None, accepted=None):
        """Just print the objective function value (used to monitor convergence in a callback) """
        if f is not None and accepted is not None:
            print("%5ds %22.10f %s" % (_time.time() - startTime, f, 'accepted' if accepted else 'not accepted'))
        else:
            result = objFunc(x)
            duration = _time.time() - startTime
            try:
                print("%5ds %22.10f" % (duration, result))
            except TypeError:  # Objfun returns vector, not scalar
                print('%5ds %s' % (duration, result))
    return print_obj_func


def _fwd_diff_jacobian(f, x0, eps=1e-10):
    y0 = f(x0).copy()
    M = len(y0)
    N = len(x0)
    jac = _np.empty((M, N), 'd')

    for j in range(N):
        #print('Adding eps to {}'.format(j))
        xj = x0.copy(); xj[j] += eps
        yj = f(xj).copy()
        #print('y0, yj')
        #print(y0[48:52])
        #print(yj[48:52])
        df = (yj - y0) / eps  # df_dxj
        jac[:, j] = df
        #print(df[48:52])

    return jac


def check_jac(f, x0, jacToCheck, eps=1e-10, tol=1e-6, errType='rel',
              verbosity=1):
    """
    Checks a jacobian function using finite differences.

    Parameters
    ----------
    f : function
        The function to check.

    x0 : numpy array
        The point at which to check the jacobian.

    jacToCheck : function
        A function which should compute the jacobian of f at x0.

    eps : float, optional
        Epsilon to use in finite difference calculations of jacobian.

    tol : float, optional
        The allowd tolerance on the relative differene between the
        values of the finite difference and jacToCheck jacobians
        if errType == 'rel' or the absolute difference if errType == 'abs'.

    verbosity : int, optional
        Controls how much detail is printed to stdout.

    Returns
    -------
    errSum : float
        The total error between the jacobians.
    errs : list
        List of (row,col,err) tuples giving the error for each row and column.
    ffd_jac : numpy array
        The computed forward-finite-difference jacobian.
    """
    orig_stdout = _sys.stdout
    devnull = open(_os.devnull, 'w')

    try:
        _sys.stdout = devnull  # redirect stdout to null during the many f(x) calls
        fd_jac = _fwd_diff_jacobian(f, x0, eps)
    finally:
        _sys.stdout = orig_stdout
        devnull.close()

    assert(jacToCheck.shape == fd_jac.shape)
    M, N = jacToCheck.shape

    errSum = 0; errs = []
    if errType == 'rel':
        for i in range(M):
            for j in range(N):
                err = _np.abs(fd_jac[i, j] - jacToCheck[i, j]) / (_np.abs(fd_jac[i, j]) + 1e-10)
                if err > tol: errs.append((i, j, err))
                errSum += err

    elif errType == 'abs':
        for i in range(M):
            for j in range(N):
                err = _np.abs(fd_jac[i, j] - jacToCheck[i, j])
                if err > tol:
                    errs.append((i, j, err))
                    if verbosity > 1:
                        print("JAC CHECK (%d,%d): %g vs %g (diff = %g)" %
                              (i, j, fd_jac[i, j], jacToCheck[i, j], fd_jac[i, j] - jacToCheck[i, j]))
                errSum += err

    errs.sort(key=lambda x: -x[2])

    if len(errs) > 0:
        maxabs = _np.max(_np.abs(jacToCheck))
        max_err_ratio = _np.max([x[2] / maxabs for x in errs])
        if verbosity > 0:
            if max_err_ratio > 0.01:
                print("Warning: jacobian_check has max err/jac_max = %g (jac_max = %g)" % (max_err_ratio, maxabs))

    if verbosity > 0:
        print("check_jac %s [err = %g]" % (("ERROR" if len(errs) else "OK"), errSum))

    return errSum, errs, fd_jac