GridFire/build-config/cppad/include/cppad/core/jacobian.hpp

# ifndef CPPAD_CORE_JACOBIAN_HPP
# define CPPAD_CORE_JACOBIAN_HPP
/* --------------------------------------------------------------------------
CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-18 Bradley M. Bell

CppAD is distributed under the terms of the
             Eclipse Public License Version 2.0.

This Source Code may also be made available under the following
Secondary License when the conditions for such availability set forth
in the Eclipse Public License, Version 2.0 are satisfied:
      GNU General Public License, Version 2.0 or later.
---------------------------------------------------------------------------- */

/*
$begin Jacobian$$
$spell
    jac
    typename
    Taylor
    Jacobian
    DetLu
    const
$$


$section Jacobian: Driver Routine$$

$head Syntax$$
$icode%jac% = %f%.Jacobian(%x%)%$$


$head Purpose$$
We use $latex F : \B{R}^n \rightarrow \B{R}^m$$ to denote the
$cref/AD function/glossary/AD Function/$$ corresponding to $icode f$$.
The syntax above sets $icode jac$$ to the
Jacobian of $icode F$$ evaluated at $icode x$$; i.e.,
$latex \[
    jac = F^{(1)} (x)
\] $$

$head f$$
The object $icode f$$ has prototype
$codei%
    ADFun<%Base%> %f%
%$$
Note that the $cref ADFun$$ object $icode f$$ is not $code const$$
(see $cref/Forward or Reverse/Jacobian/Forward or Reverse/$$ below).

$head x$$
The argument $icode x$$ has prototype
$codei%
    const %Vector% &%x%
%$$
(see $cref/Vector/Jacobian/Vector/$$ below)
and its size
must be equal to $icode n$$, the dimension of the
$cref/domain/seq_property/Domain/$$ space for $icode f$$.
It specifies
that point at which to evaluate the Jacobian.

$head jac$$
The result $icode jac$$ has prototype
$codei%
    %Vector% %jac%
%$$
(see $cref/Vector/Jacobian/Vector/$$ below)
and its size is $latex m * n$$; i.e., the product of the
$cref/domain/seq_property/Domain/$$
and
$cref/range/seq_property/Range/$$
dimensions for $icode f$$.
For $latex i = 0 , \ldots , m - 1 $$
and $latex j = 0 , \ldots , n - 1$$
$latex \[.
    jac[ i * n + j ] = \D{ F_i }{ x_j } ( x )
\] $$


$head Vector$$
The type $icode Vector$$ must be a $cref SimpleVector$$ class with
$cref/elements of type/SimpleVector/Elements of Specified Type/$$
$icode Base$$.
The routine $cref CheckSimpleVector$$ will generate an error message
if this is not the case.

$head Forward or Reverse$$
This will use order zero Forward mode and either
order one Forward or order one Reverse to compute the Jacobian
(depending on which it estimates will require less work).
After each call to $cref Forward$$,
the object $icode f$$ contains the corresponding
$cref/Taylor coefficients/glossary/Taylor Coefficient/$$.
After a call to $code Jacobian$$,
the zero order Taylor coefficients correspond to
$icode%f%.Forward(0, %x%)%$$
and the other coefficients are unspecified.

$head Example$$
$children%
    example/general/jacobian.cpp
%$$
The routine
$cref/Jacobian/jacobian.cpp/$$ is both an example and test.
It returns $code true$$, if it succeeds and $code false$$ otherwise.

$end
-----------------------------------------------------------------------------
*/

//  BEGIN CppAD namespace
namespace CppAD {

template <class Base, class RecBase, class Vector>
void JacobianFor(ADFun<Base, RecBase> &f, const Vector &x, Vector &jac)
{   size_t i;
    size_t j;

    size_t n = f.Domain();
    size_t m = f.Range();

    // check Vector is Simple Vector class with Base type elements
    CheckSimpleVector<Base, Vector>();

    CPPAD_ASSERT_UNKNOWN( size_t(x.size())   == f.Domain() );
    CPPAD_ASSERT_UNKNOWN( size_t(jac.size()) == f.Range() * f.Domain() );

    // argument and result for forward mode calculations
    Vector u(n);
    Vector v(m);

    // initialize all the components
    for(j = 0; j < n; j++)
        u[j] = Base(0.0);

    // loop through the different coordinate directions
    for(j = 0; j < n; j++)
    {   // set u to the j-th coordinate direction
        u[j] = Base(1.0);

        // compute the partial of f w.r.t. this coordinate direction
        v = f.Forward(1, u);

        // reset u to vector of all zeros
        u[j] = Base(0.0);

        // return the result
        for(i = 0; i < m; i++)
            jac[ i * n + j ] = v[i];
    }
}
template <class Base, class RecBase, class Vector>
void JacobianRev(ADFun<Base, RecBase> &f, const Vector &x, Vector &jac)
{   size_t i;
    size_t j;

    size_t n = f.Domain();
    size_t m = f.Range();

    CPPAD_ASSERT_UNKNOWN( size_t(x.size())   == f.Domain() );
    CPPAD_ASSERT_UNKNOWN( size_t(jac.size()) == f.Range() * f.Domain() );

    // argument and result for reverse mode calculations
    Vector u(n);
    Vector v(m);

    // initialize all the components
    for(i = 0; i < m; i++)
        v[i] = Base(0.0);

    // loop through the different coordinate directions
    for(i = 0; i < m; i++)
    {   if( f.Parameter(i) )
        {   // return zero for this component of f
            for(j = 0; j < n; j++)
                jac[ i * n + j ] = Base(0.0);
        }
        else
        {
            // set v to the i-th coordinate direction
            v[i] = Base(1.0);

            // compute the derivative of this component of f
            u = f.Reverse(1, v);

            // reset v to vector of all zeros
            v[i] = Base(0.0);

            // return the result
            for(j = 0; j < n; j++)
                jac[ i * n + j ] = u[j];
        }
    }
}

template <class Base, class RecBase>
template <class Vector>
Vector ADFun<Base,RecBase>::Jacobian(const Vector &x)
{   size_t i;
    size_t n = Domain();
    size_t m = Range();

    CPPAD_ASSERT_KNOWN(
        size_t(x.size()) == n,
        "Jacobian: length of x not equal domain dimension for F"
    );

    // point at which we are evaluating the Jacobian
    Forward(0, x);

    // work factor for forward mode
    size_t workForward = n;

    // work factor for reverse mode
    size_t workReverse = 0;
    for(i = 0; i < m; i++)
    {   if( ! Parameter(i) )
            ++workReverse;
    }

    // choose the method with the least work
    Vector jac( n * m );
# ifdef CPPAD_FOR_TMB
    if( workForward < workReverse )
# else
    if( workForward <= workReverse )
# endif
        JacobianFor(*this, x, jac);
    else
        JacobianRev(*this, x, jac);

    return jac;
}

} // END CppAD namespace

# endif