diff_gradient.h

Go to the documentation of this file.

/* =========================================================================== *
   |
   |  Copyright (c) 1994-2011 by Kobus Barnard (author)
   |
   |  Personal and educational use of this code is granted, provided that this
   |  header is kept intact, and that the authorship is not misrepresented, that
   |  its use is acknowledged in publications, and relevant papers are cited.
   |
   |  For other use contact the author (kobus AT cs DOT arizona DOT edu).
   |
   |  Please note that the code in this file has not necessarily been adequately
   |  tested. Naturally, there is no guarantee of performance, support, or fitness
   |  for any particular task. Nonetheless, I am interested in hearing about
   |  problems that you encounter.
   |
   |  Author:  Ernesto Brau, Zewei Jiang
 * =========================================================================== */

/* $Id$ */

#ifndef DIFF_GRADIENT_H
#define DIFF_GRADIENT_H

#include <m_cpp/m_vector.h>
#include <diff_cpp/diff_util.h>
#include <vector>

namespace kjb {

template<class Func, class Model, class Adapter>
Vector gradient_cfd
(
    const Func& f,
    const Model& x,
    const std::vector<double>& dx,
    const Adapter& adapter
)
{
    const size_t D = adapter.size(&x);

    Vector G(D);
    Model y = x;

    for(size_t i = 0; i < D; i++)
    {
        double yi = adapter.get(&y, i);

        move_param(y, i, dx[i], adapter);
        double fxp = f(y);

        move_param(y, i, -2*dx[i], adapter);
        double fxm = f(y);

        G[i] = (fxp - fxm) / (2*dx[i]);

        // return to original spot
        adapter.set(&y, i, yi);
    }

    return G;
}

template<class Func, class Vec>
inline
Vector gradient_cfd
(
    const Func& f,
    const Vec& x,
    const std::vector<double>& dx
)
{
    return gradient_cfd(f, x, dx, Vector_adapter<Vec>());
}

template<class Func, class Model, class Adapter>
Vector gradient_ffd
(
    const Func& f,
    const Model& x,
    const std::vector<double>& dx,
    const Adapter& adapter
)
{
    const size_t D = adapter.size(&x);

    Vector G(D);
    double fx = f(x);
    Model y = x;

    for(size_t i = 0; i < D; i++)
    {
        double yi = adapter.get(&y, i);
        //y[i] += dx[i];
        move_param(y, i, dx[i], adapter);
        double fxp = f(y);

        G[i] = (fxp - fx) / dx[i];

        // return to original spot
        //y[i] -= dx[i];
        //move_param(y, i, -dx[i], adapter);
        adapter.set(&y, i, yi);
    }

    return G;
}

template<class Func, class Vec>
inline
Vector gradient_ffd
(
    const Func& f,
    const Vec& x,
    const std::vector<double>& dx
)
{
    return gradient_ffd(f, x, dx, Vector_adapter<Vec>());
}

/* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */
/*                        INDEPENDENT VERSIONS OF GRADIENT                    */
/*                                                                            */
/* These versions use "smart" functions that may exploit indepdendence prop-  */
/* erties to avoid full function evaluation.  For any dimension i, f(x) may   */
/* be decomposed the sum of terms involing x_i and those not.  The latter     */
/* terms cancel in the finite-differences equation, and can be ignored.  Thus */
/* the functions below receive f(x,i) which omits terms not involving x_i,    */
/* often at great performance advantage.                                      */
/* \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ */

template<class Func, class Model, class Adapter>
Vector gradient_ind_cfd
(
    const Func& f,
    const Model& x,
    const std::vector<double>& dx,
    const Adapter& adapter
)
{
    const size_t D = adapter.size(&x);

    Vector G(D);
    Model y = x;

    for(size_t i = 0; i < D; i++)
    {
        double yi = adapter.get(&y, i);

        move_param(y, i, dx[i], adapter);
        double fxp = f(y, i);

        move_param(y, i, -2*dx[i], adapter);
        double fxm = f(y, i);

        G[i] = (fxp - fxm) / (2*dx[i]);

        // return to original spot
        adapter.set(&y, i, yi);
    }

    return G;
}

template<class Func, class Vec>
inline
Vector gradient_ind_cfd
(
    const Func& f,
    const Vec& x,
    const std::vector<double>& dx
)
{
    return gradient_ind_cfd(f, x, dx, Vector_adapter<Vec>());
}

} //namespace kjb

#endif /*DIFF_GRADIENT_H */