<commit_before>// ///////////////////////////////////////////////////////////////////////////
// tenh/utility/optimization.hpp by <NAME>, created 2013/12/09
// Copyright Leap Motion Inc.
// ///////////////////////////////////////////////////////////////////////////

#ifndef TENH_UTILITY_OPTIMIZATION_HPP_
#define TENH_UTILITY_OPTIMIZATION_HPP_

#include "tenh/core.hpp"

#include <cmath>
#include <iostream>

#include "tenh/implementation/innerproduct.hpp"
#include "tenh/implementation/vector.hpp"
#include "tenh/implementation/vee.hpp"
#include "tenh/interop/eigen_invert.hpp"

namespace Tenh {

template <typename T>
bool isNaN(T const &x)
{
    return false;
}

template<>
bool isNaN(float const &x)
{
    return std::isnan(x);
}

template<>
bool isNaN(double const &x)
{
    return std::isnan(x);
}

template<>
bool isNaN(long double const &x)
{
    return std::isnan(x);
}

// adaptive minimization which uses, in order of availability/preference,
// 1. Newton's method, 2. conjugate gradient, and 3. gradient descent.
template <typename InnerProductId_, typename ObjectiveFunction_, typename BasedVectorSpace_, typename Scalar_, typename GuessUseArrayType_, typename Derived_>
ImplementationOf_t<BasedVectorSpace_,Scalar_> minimize (ObjectiveFunction_ const &func,
                                                        ImplementationOf_t<BasedVectorSpace_,Scalar_,GuessUseArrayType_,Derived_> const &guess,
                                                        Scalar_ tolerance,
                                                        Scalar_ *minimum = NULL)
{
    typedef ImplementationOf_t<BasedVectorSpace_,Scalar_> VectorType;
    typedef typename InnerProduct_f<BasedVectorSpace_,InnerProductId_,Scalar_>::T VectorInnerProductType;
    typedef ImplementationOf_t<typename DualOf_f<BasedVectorSpace_>::T,Scalar_> CoVectorType;
    typedef typename InnerProduct_f<typename DualOf_f<BasedVectorSpace_>::T,InnerProductId_,Scalar_>::T CoVectorInnerProductType;
    typedef SymmetricPowerOfBasedVectorSpace_c<2,BasedVectorSpace_> Sym2;
    typedef typename ObjectiveFunction_::D2 HessianType;
    typedef ImplementationOf_t<SymmetricPowerOfBasedVectorSpace_c<2,BasedVectorSpace_>,Scalar_> HessianInverseType;
    STATIC_ASSERT_TYPES_ARE_EQUAL(VectorType, typename ObjectiveFunction_::V);
    STATIC_ASSERT_TYPES_ARE_EQUAL(Scalar_, typename ObjectiveFunction_::Out);
    static int const LINE_SEARCH_SAMPLE_COUNT = 50;
    static Scalar_ const STEP_SCALE = Scalar_(1);
    static Scalar_ const MAX_STEP_SIZE = Scalar_(-1);
    // static Scalar_ const MAX_STEP_SIZE = Scalar_(1);
    static int const MAX_ITERATION_COUNT = 20;
    static Scalar_ const GRADIENT_DESCENT_STEP_SIZE = Scalar_(0.25);
    static bool const PRINT_DEBUG_OUTPUT = false;
    static Scalar_ const EPSILON = 1e-10;

    VectorInnerProductType vector_innerproduct;
    CoVectorInnerProductType covector_innerproduct;

    AbstractIndex_c<'i'> i;
    AbstractIndex_c<'j'> j;

    VectorType current_approximation = guess;
    int iteration_count = 0;
    int gradient_descent = 0;
    int conjugate_gradient = 0;
    int newtons_method = 0;
    bool tolerance_was_attained = false;

    while (iteration_count < MAX_ITERATION_COUNT)
    {
        Scalar_ current_value = func.function(current_approximation);
        CoVectorType g = func.D_function(current_approximation);
        Scalar_ g_squared_norm = covector_innerproduct(g, g);//g(i)*covector_innerproduct.split(i*j)*g(j);
        Scalar_ g_norm = std::sqrt(g_squared_norm);

        if (PRINT_DEBUG_OUTPUT)
        {
            std::cout << current_approximation << "\n" << g << " " << g_norm <<  " " << current_value << std::endl;
        }

        if (g_norm <= tolerance)
        {
            tolerance_was_attained = true;
            if (minimum != NULL)
                *minimum = current_value;
            break;
        }

        HessianType h(func.D2_function(current_approximation));
        HessianInverseType hinv(Static<WithoutInitialization>::SINGLETON);
        VectorType step(Static<WithoutInitialization>::SINGLETON);
        Scalar_ d;

        if (invert_2tensor(h, hinv))
        {
            step(i).no_alias() = hinv.split(i*j) * -g(j);
            d = h(step, step);//h(a)*(step(i)*step(j)).bundle(i*j,Sym2(),a);
        }
        else
        {
            if (PRINT_DEBUG_OUTPUT)
            {
                std::cout << "Hessian wasn't invertible." << std::endl;
            }
            d = 0; // less than EPSILON.
        }

        if (isNaN(d) || d < EPSILON) // h isn't postive definite along step so fall back to conjugate gradient
        {
            VectorType v(Static<WithoutInitialization>::SINGLETON);
            v(j).no_alias() = g(i) * covector_innerproduct.split(i*j);
            // d = g(i)*v(i);
            d = h(v, v);//v(i) * h.split(i*j) * v(j);

            if (isNaN(d) || d < EPSILON) // h isn't positive definite along g either, gradient descent
            {
                step(i).no_alias() = -GRADIENT_DESCENT_STEP_SIZE * g(i) / g_norm;
                if (PRINT_DEBUG_OUTPUT)
                {
                    std::cout << "Gradient descent. With step = " << step << std::endl;<commit_msg>Replace Hessian inversion with SVD solve.<commit_after>// ///////////////////////////////////////////////////////////////////////////
// tenh/utility/optimization.hpp by <NAME>, created 2013/12/09
// Copyright Leap Motion Inc.
// ///////////////////////////////////////////////////////////////////////////

#ifndef TENH_UTILITY_OPTIMIZATION_HPP_
#define TENH_UTILITY_OPTIMIZATION_HPP_

#include "tenh/core.hpp"

#include <cmath>
#include <iostream>

#include "tenh/implementation/innerproduct.hpp"
#include "tenh/implementation/vector.hpp"
#include "tenh/implementation/vee.hpp"
#include "tenh/interop/eigen_invert.hpp"
#include "tenh/interop/eigen_svd.hpp"

namespace Tenh {

template <typename T>
bool isNaN(T const &x)
{
    return false;
}

template<>
bool isNaN(float const &x)
{
    return std::isnan(x);
}

template<>
bool isNaN(double const &x)
{
    return std::isnan(x);
}

template<>
bool isNaN(long double const &x)
{
    return std::isnan(x);
}

// adaptive minimization which uses, in order of availability/preference,
// 1. Newton's method, 2. conjugate gradient, and 3. gradient descent.
template <typename InnerProductId_, typename ObjectiveFunction_, typename BasedVectorSpace_, typename Scalar_, typename GuessUseArrayType_, typename Derived_>
ImplementationOf_t<BasedVectorSpace_,Scalar_> minimize (ObjectiveFunction_ const &func,
                                                        ImplementationOf_t<BasedVectorSpace_,Scalar_,GuessUseArrayType_,Derived_> const &guess,
                                                        Scalar_ tolerance,
                                                        Scalar_ *minimum = NULL)
{
    typedef ImplementationOf_t<BasedVectorSpace_,Scalar_> VectorType;
    typedef typename InnerProduct_f<BasedVectorSpace_,InnerProductId_,Scalar_>::T VectorInnerProductType;
    typedef ImplementationOf_t<typename DualOf_f<BasedVectorSpace_>::T,Scalar_> CoVectorType;
    typedef typename InnerProduct_f<typename DualOf_f<BasedVectorSpace_>::T,InnerProductId_,Scalar_>::T CoVectorInnerProductType;
    typedef SymmetricPowerOfBasedVectorSpace_c<2,BasedVectorSpace_> Sym2;
    typedef typename ObjectiveFunction_::D2 HessianType;
    typedef ImplementationOf_t<TensorProductOfBasedVectorSpaces_c<TypeList_t<typename DualOf_f<BasedVectorSpace_>::T, TypeList_t<typename DualOf_f<BasedVectorSpace_>::T> > >,Scalar_> Hessian2TensorType;
    STATIC_ASSERT_TYPES_ARE_EQUAL(VectorType, typename ObjectiveFunction_::V);
    STATIC_ASSERT_TYPES_ARE_EQUAL(Scalar_, typename ObjectiveFunction_::Out);
    static int const LINE_SEARCH_SAMPLE_COUNT = 50;
    static Scalar_ const STEP_SCALE = Scalar_(1);
    static Scalar_ const MAX_STEP_SIZE = Scalar_(-1);
    // static Scalar_ const MAX_STEP_SIZE = Scalar_(1);
    static int const MAX_ITERATION_COUNT = 20;
    static Scalar_ const GRADIENT_DESCENT_STEP_SIZE = Scalar_(0.25);
    static bool const PRINT_DEBUG_OUTPUT = false;
    static Scalar_ const EPSILON = 1e-10;

    VectorInnerProductType vector_innerproduct;
    CoVectorInnerProductType covector_innerproduct;

    AbstractIndex_c<'i'> i;
    AbstractIndex_c<'j'> j;

    VectorType current_approximation = guess;
    int iteration_count = 0;
    int gradient_descent = 0;
    int conjugate_gradient = 0;
    int newtons_method = 0;
    bool tolerance_was_attained = false;

    while (iteration_count < MAX_ITERATION_COUNT)
    {
        Scalar_ current_value = func.function(current_approximation);
        CoVectorType g = func.D_function(current_approximation);
        CoVectorType minus_g(Static<WithoutInitialization>::SINGLETON);
        Scalar_ g_squared_norm = covector_innerproduct(g, g);//g(i)*covector_innerproduct.split(i*j)*g(j);
        Scalar_ g_norm = std::sqrt(g_squared_norm);
        minus_g(i) = -g(i);

        if (PRINT_DEBUG_OUTPUT)
        {
            std::cout << current_approximation << "\n" << g << " " << g_norm <<  " " << current_value << std::endl;
        }

        if (g_norm <= tolerance)
        {
            tolerance_was_attained = true;
            if (minimum != NULL)
                *minimum = current_value;
            break;
        }

        HessianType h(func.D2_function(current_approximation));
        Hessian2TensorType h2(Static<WithoutInitialization>::SINGLETON);
        VectorType step(Static<WithoutInitialization>::SINGLETON);
        Scalar_ d;
        h2(i*j) = h.split(i*j);

        SVD_solve(h2,step,minus_g);
        d = h(step, step);//h(a)*(step(i)*step(j)).bundle(i*j,Sym2(),a);

        if (isNaN(d) || d < EPSILON) // h isn't postive definite along step so fall back to conjugate gradient
        {
            VectorType v(Static<WithoutInitialization>::SINGLETON);
            v(j).no_alias() = g(i) * covector_innerproduct.split(i*j);
            // d = g(i)*v(i);
            d = h(v, v);//v(i) * h.split(i*j) * v(j);

            if (isNaN(d) || d < EPSILON) // h isn't positive definite along g either, gradient descent
            {
                step(i).no_alias() = -GRADIENT_DESCENT_STEP_SIZE * g(i) / g_norm;
                if (PRINT_DEBUG_OUTPUT)
                {
                    std::cout << "Gradient descent. With step = " << step << std::endl;
{'fraction_non_alphanumeric': 0.13084022303403192, 'fraction_numerical': 0.0071140165352816765, 'mean_word_length': 2.7746734397677795, 'pattern_counts': {'":': 0, '<': 112, '<?xml version=': 0, '>': 56, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2204312', 'n_tokens_mistral': 3297, 'n_tokens_neox': 3078, 'n_words': 644}