'''
Created on Jan 5, 2016

@author: anon
'''
import numpy as np
import scipy.linalg as la
import tensorflow as tf


def which_nodes(xmin, xmax, quadratures_nodes):
    # Fix the integration interval to [-1,1].
    nodes2 = 0.5 * (xmax - xmin) * quadratures_nodes + 0.5 * (xmax + xmin)
    return nodes2


def GaussLegendreQuadrature(N=30):
    """ Returns the quadratures_nodes anb weights of the Gaussian-Lenegendre quadrature """
    beta = np.array([(n + 1.0) / np.sqrt((2.0 * n + 1.0) * (2.0 * n + 3.0))
                     for n in range(N - 1)], dtype=np.float32)
    M = np.diag(beta, -1) + np.diag(beta, 1)
    nodes, V = la.eigh(M, overwrite_a=True, overwrite_b=True)
    weight = 2 * V[0, :] ** 2
    return nodes, weight


nodes, weight = GaussLegendreQuadrature()


def quadratures(funct, a=None, b=None, num_nodes=30, **kwargs):
    quadrature_window = b - a
    q_p_normalized, q_w = normalize_quadrature_points(a, b, num_nodes)
    intensity_q = funct(q_p_normalized, **kwargs)
    integrand_q = tf.multiply(q_w, intensity_q)

    integral = tf.multiply(tf.cast(quadrature_window * .5, dtype=tf.float64),
                           tf.reduce_sum(integrand_q, 2))

    return integral


def normalize_quadrature_points(a=None, b=None, num_nodes=30):
    quadratures_nodes, quadratures_weights = GaussLegendreQuadrature(num_nodes)
    x_shape = tf.shape(a)
    B, T = x_shape[0], x_shape[1]
    q_p = tf.reshape(quadratures_nodes, (1, 1, num_nodes))
    q_p = tf.tile(q_p, [B, T, 1])
    q_w = tf.cast(tf.reshape(quadratures_weights, (1, 1, num_nodes)), tf.float64)
    min_p_max = tf.multiply(tf.cast(.5, tf.float64), (a + b))
    max_m_min = tf.multiply(tf.cast(.5, tf.float64), (b - a))

    q_p_t = tf.cast(q_p, tf.float64) * max_m_min[:, :, None] + min_p_max[:, :, None]
    return q_p_t, q_w


"""
xmin=0.
xmax=1.
nodes2 = which_nodes(xmin,xmax,quadratures_nodes)
f = np.exp(nodes2)
(np.dot(f, weight))*((xmax-xmin)/2.)
"""
