/* ----------------------------------------------------------------------- *//** 
 * 
 * @file newton.hpp
 * 
 * Generic implementaion of Newton's method, in the fashion of
 * user-definied aggregates. They should be called by actually database
 * functions, after arguments are properly parsed.
 *
 *//* ----------------------------------------------------------------------- */

#include <dbconnector/dbconnector.hpp>

#ifndef MADLIB_MODULES_CONVEX_ALGO_NEWTON_HPP_
#define MADLIB_MODULES_CONVEX_ALGO_NEWTON_HPP_

namespace madlib {

namespace modules {

namespace convex {

// use Eigen
using namespace madlib::dbal::eigen_integration;
    
// The reason for using ConstState instead of const State to reduce the
// template type list: flexibility to high-level for mutability control
// More: cast<ConstState>(MutableState) may not always work
template <class State, class ConstState, class Task>
class Newton {
public:
    typedef State state_type;
    typedef ConstState const_state_type;
    typedef typename Task::tuple_type tuple_type;
    typedef typename Task::model_type model_type;

    static void transition(state_type &state, const tuple_type &tuple);
    static void merge(state_type &state, const_state_type &otherState);
    static void final(state_type &state);
};

template <class State, class ConstState, class Task>
void
Newton<State, ConstState, Task>::transition(state_type &state,
        const tuple_type &tuple) {
    // accumulating the gradient & hessian
    Task::gradient(
            state.task.model,
            tuple.indVar,
            tuple.depVar,
            state.algo.gradient);
    Task::hessian(
            state.task.model,
            tuple.indVar,
            tuple.depVar,
            state.algo.hessian);
}

template <class State, class ConstState, class Task>
void
Newton<State, ConstState, Task>::merge(state_type &state,
        const_state_type &otherState) {
    // merging accumulated gradient & hessian
    state.algo.gradient += otherState.algo.gradient;
    state.algo.hessian += otherState.algo.hessian;
}

template <class State, class ConstState, class Task>
void
Newton<State, ConstState, Task>::final(state_type &state) {
    // w_{k+1} = w_k - H_k^{-1} g_k
    // instead of computing the inverse explicitly,
    // we solve H_k p_k = g_k, and w_{k+1} = w_k - 1 * p_k

    // solve in place, subject to change if gradient is somehow needed
    // ldlt decomposition is chosen for its numerical stability, see 
    // http://eigen.tuxfamily.org/dox-3.0/TutorialLinearAlgebra.html
    state.algo.gradient = state.algo.hessian.ldlt().solve(state.algo.gradient);
    state.task.model -= state.algo.gradient;
}

} // namespace convex

} // namespace modules

} // namespace madlib

#endif