#ifndef __monotonic_function_space_h__
#define __monotonic_function_space_h__

namespace toposmooth
{
/*
    We are given a tree $T$ with fixed scalar values at the root node and at the leaf nodes.
    At every node $i$ other than the root and leaf nodes, we store a degree of freedom $theta_i$.
    Let the function $f$ map from a node in $T$ to a scalar value.
    In other words, $f: i \in T \rightarrow \mathbb{R}$.
    Let
                  { fixed scalar value if $i$ is the root node or one of the leaf nodes
        $f(i)$ =  { lerp( max( { f( $j$ ) for $j$ in leafs of $i$ } ), f( parent of $i$ ), 0 < t_i < 1 ) = max( leafs of $i$ ) + t_i * ( f( parent of $i$ ) - max( { f( $j$ ) for $j$ in leafs of $i$ } ) ) otherwise
    
    So for the root and leaf nodes, we store a fixed scalar value.
    For the other nodes, we store a $t_i$.  To constrain $t_i$ to lie between 0 and 1, we actually store a $theta_i$ and use the identity $t_i = .5 + .5 * cos( theta_i )$.
    
    NOTE: This assumes that values are monotonically descending.  If values
          should be monotonically ascending instead, use min() instead of max().
*/

// The TNC optimizer goes crazy if we use float's instead of double's.
typedef double real_t;

class tree_ops_t
{
public:
    typedef long node_t;
    
    virtual ~tree_ops_t() {}
    
    // This function should return true if the node corresponding to the 'node' parameter
    // is a fixed node (the root or a leaf node) and false otherwise.
    virtual bool is_fixed( node_t node ) const = 0;
    
    // This function should return the fixed value corresponding to the 'node' parameter;
    // the 'node' parameter passed to this function will always be fixed (the root or a leaf node).
    virtual real_t fixed_value( node_t node ) const = 0;
    
    // This function should return the current theta value corresponding to the 'node' parameter;
    // the 'node' parameter passed to this function will never be fixed (the root or a leaf node).
    virtual real_t theta( node_t node ) const = 0;
    
    // This function should return the maximum value among all leaf nodes having 'node' as
    // an ancestor.  The 'node' parameter passed to this function will never be a leaf node.
    // NOTE: This assumes that values are monotonically descending.  If values
    //       should be monotonically ascending instead, return min() instead of max().
    virtual real_t max_leaf_value( node_t node ) const = 0;
    
    // This function should return the parent node_t of the 'node' parameter.
    // If the 'node' parameter is the root node, the function should return 'node'.
    // In other words, the root node should be considered to be its own parent.
    virtual node_t parent_node( node_t node ) const = 0;
    
    // This function should return true if the node corresponding to the 'old' parameter
    // is an ancestor of the node corresponding to the 'young' parameter and false otherwise.
    // If 'old' and 'young' are the same, false should be returned.
    virtual bool is_ancestor( node_t old, node_t young ) const = 0;
    
    // This function computes and returns the scalar value of the node given by 'node_i'.
    virtual real_t f_i( node_t node_i ) const;
    
    // This function computes and returns the scalar value of the node given by 'node_i'.
    // It does not use recursion, but instead uses dynamic memory allocation.
    // The optional parameter 'depth' is a hint as to how many nodes deep 'node_i' is;
    // If 'node_i' is the root, its depth is 0.
    real_t f_i_nonrecursive( node_t node_i, unsigned int depth = 0 ) const;
    
    // This function computes and returns the scalar value of the node given by 'node_i_and_lineage[0]',
    // where 'node_i_and_lineage' is an array of length 'num_nodes', which must be at least one,
    // and 'node_i_and_lineage[i+1]' is the parent of 'node_i_and_lineage[i]', all the way to the
    // root.  The node given by 'node_i_and_lineage[num_nodes-1]' must always be fixed.
    real_t f_i( node_t* node_i_and_lineage, int num_nodes ) const;
    
    // This function computes and returns the derivative of the scalar value of the node
    // given by 'node_i' with respect to the theta of 'node_j'.
    // Upon return, the optional output parameter 'constant_out' will be
    // set to true if the result will always be zero (and false otherwise).
    real_t f_i_d_theta_j( node_t node_i, node_t node_j, bool* constant_out = 0 ) const;
    // This version doesn't perform any constant checks, so the user must only call this version
    // with the same inputs for which a call to the usual function would have return
    // true in 'constant_out'.
    real_t f_i_d_theta_j_nonconstant( node_t node_i, node_t node_j ) const;
    // A version of f_i_d_theta_j_nonconstant() that doesn't make any recursive calls.
    // It takes an additional, optional parameter estimating the depth of 'j' in the tree.
    real_t f_i_d_theta_j_nonconstant_nonrecursive( node_t node_i, node_t node_j, unsigned int depth_j = 0 ) const;
    
    // This function computes and returns the second derivative of the scalar value of the node
    // given by 'node_i' with respect to both the theta of 'node_j' and the theta of 'node_k'.
    // Upon return, the optional output parameter 'constant_out' will be
    // set to true if the result will always be zero (and false otherwise).
    real_t f_i_d_theta_j_theta_k( node_t node_i, node_t node_j, node_t node_k, bool* constant_out = 0 ) const;
    // This version doesn't perform any constant checks, so the user must only call this version
    // with the same inputs for which a call to the usual function would have return
    // true in 'constant_out'.
    real_t f_i_d_theta_j_theta_k_nonconstant( node_t node_i, node_t node_j, node_t node_k ) const;
    // A version of f_i_d_theta_j_theta_k_nonconstant() that doesn't make any recursive calls.
    // It takes an additional, optional parameter estimating the depth of 'j' or 'k' in the tree,
    // whichever is smaller.
    real_t f_i_d_theta_j_theta_k_nonconstant_nonrecursive( node_t node_i, node_t node_j, node_t node_k, unsigned int depth_j = 0 ) const;
    
public:
    // This function computes and returns the interpolating 't' value between 0 and 1 based on the
    // theta value of the node given by 'node_i'.
    real_t t_i( node_t node_i ) const;
    
    // This function computes and returns the derivative of the
    // interpolating 't' value between 0 and 1 based on the
    // theta value of the node given by 'node_i'.
    real_t dt_i( node_t node_i ) const;
    
    // This function computes and returns the second derivative of the
    // interpolating 't' value between 0 and 1 based on the
    // theta value of the node given by 'node_i'.
    real_t d2t_i( node_t node_i ) const;
};

}
// ~namespace toposmooth

#endif /* __monotonic_function_space_h__ */
