import numpy as np
from numpy import dtype
import scipy
from numpy import transpose
from scipy import matmul
from scipy.linalg import expm,bandwidth,norm,solve,cosm,sinm,tanm,eigh,eig_banded,eigh_tridiagonal,cholesky_banded,cho_solve_banded,signm,convolution_matrix
from numpy.lib.scimath import sqrt as csqrt
from itertools import product
from numpy import asarray, Inf, dot, floor, eye, diag, exp, \
     product, logical_not, ravel, transpose, conjugate, \
     cast, log, ogrid, imag, real, absolute, amax, sign, \
     isfinite, sqrt, identity, single, ceil, log2
import math

def expmm(A, q=False):
    
    A = asarray(A)
    A_L1 = norm(A,1)
    # print("a_l1=",A_L1)
    n_squarings = 0

    if A.dtype == 'float64' or A.dtype == 'complex128':
        if A_L1 < 1.495585217958292e-002:
            U,V = _pade3(A)
        elif A_L1 < 2.539398330063230e-001:
            U,V = _pade5(A)
        elif A_L1 < 9.504178996162932e-001:
            U,V = _pade7(A)
        elif A_L1 < 2.097847961257068e+000:
            U,V = _pade9(A)
        else:
            print("ssss")
            maxnorm = 5.371920351148152
            n_squarings = max(0, int(ceil(log2(A_L1 / maxnorm))))
            print("n_squariings",n_squarings)
            A = A / 2**n_squarings
            print("a=",A)
            U,V = _pade13(A)
    elif A.dtype == 'float32' or A.dtype == 'complex64':
        if A_L1 < 4.258730016922831e-001:
            U,V = _pade3(A)
        elif A_L1 < 1.880152677804762e+000:
            U,V = _pade5(A)
        else:
            maxnorm = 3.925724783138660
            n_squarings = max(0, int(ceil(log2(A_L1 / maxnorm))))
            A = A / 2**n_squarings
            U,V = _pade7(A)
    else:
        raise ValueError("invalid type: "+str(A.dtype))

    P = U + V  # p_m(A) : numerator
    # print("p=", P)
    Q = -U + V # q_m(A) : denominator
    # print("q=",q)
    R = solve(Q,P)
    # print("r=",R)
    # squaring step to undo scaling
    for i in range(n_squarings):
        R = dot(R,R)

    return R

# implementation of Pade approximations of various degree using the algorithm presented in [Higham 2005]

def _pade3(A):
    b = (120., 60., 12., 1.)
    ident = eye(A.shape[0], A.shape[1], dtype=A.dtype)
    A2 = dot(A,A)
    U = dot(A , (b[3]*A2 + b[1]*ident))
    V = b[2]*A2 + b[0]*ident
    return U,V

def _pade5(A):
    b = (30240., 15120., 3360., 420., 30., 1.)
    ident = eye(A.shape[0], A.shape[1], dtype=A.dtype)
    A2 = dot(A,A)
    A4 = dot(A2,A2)
    U = dot(A, b[5]*A4 + b[3]*A2 + b[1]*ident)
    V = b[4]*A4 + b[2]*A2 + b[0]*ident
    return U,V

def _pade7(A):
    b = (17297280., 8648640., 1995840., 277200., 25200., 1512., 56., 1.)
    ident = eye(A.shape[0], A.shape[1], dtype=A.dtype)
    A2 = dot(A,A)
    A4 = dot(A2,A2)
    A6 = dot(A4,A2)
    U = dot(A, b[7]*A6 + b[5]*A4 + b[3]*A2 + b[1]*ident)
    V = b[6]*A6 + b[4]*A4 + b[2]*A2 + b[0]*ident
    return U,V

def _pade9(A):
    b = (17643225600., 8821612800., 2075673600., 302702400., 30270240.,
                2162160., 110880., 3960., 90., 1.)
    ident = eye(A.shape[0], A.shape[1], dtype=A.dtype)
    A2 = dot(A,A)
    A4 = dot(A2,A2)
    A6 = dot(A4,A2)
    A8 = dot(A6,A2)
    U = dot(A, b[9]*A8 + b[7]*A6 + b[5]*A4 + b[3]*A2 + b[1]*ident)
    V = b[8]*A8 + b[6]*A6 + b[4]*A4 + b[2]*A2 + b[0]*ident
    return U,V

def _pade13(A):
    b = (64764752532480000., 32382376266240000., 7771770303897600.,
    1187353796428800., 129060195264000., 10559470521600., 670442572800.,
    33522128640., 1323241920., 40840800., 960960., 16380., 182., 1.)
    ident = eye(A.shape[0], A.shape[1], dtype=A.dtype)
    # print("ident =",ident)
    A2 = dot(A,A)
    # print("a2",A2)
    A4 = dot(A2,A2)
    # print("a4",A4)
    A6 = dot(A4,A2)
    # print("a6",A6)
    U = dot(A,dot(A6, b[13]*A6 + b[11]*A4 + b[9]*A2) + b[7]*A6 + b[5]*A4 + b[3]*A2 + b[1]*ident)
    print("u=",U)
    V = dot(A6, b[12]*A6 + b[10]*A4 + b[8]*A2) + b[6]*A6 + b[4]*A4 + b[2]*A2 + b[0]*ident
    print("v=",V)
    return U,V

def dft(n, scale=None):
    temp = -2j * np.pi * np.arange(n) / n
    print("temp=",temp)
    # omegas2 = np.exp(-2j * np.pi * np.arange(n) / n)
    # print("omegas2=",omegas2)
    omegas = np.exp(-2j * np.pi * np.arange(n) / n).reshape(-1, 1)
    print("omegas=",omegas)
    m = omegas ** np.arange(n)
    # print("m=",m)



    # tt =  -0.5-0.866025404j
    # print("tt1=",tt ** 1)
    # print("tt2=",tt ** 2)
    # print("tt3=",tt ** 3)
    # print("tt4=",tt ** 4)
    # print("tt5=",tt ** 5)
    # print("tt6=",tt ** 6)
    if scale == 'sqrtn':
        m /= math.sqrt(n)
    elif scale == 'n':
        m /= n
    return m

def hadamard(n, dtype=int):
    """
    Construct an Hadamard matrix.
    Constructs an n-by-n Hadamard matrix, using Sylvester's
    construction. `n` must be a power of 2.
    Parameters
    ----------
    n : int
        The order of the matrix. `n` must be a power of 2.
    dtype : dtype, optional
        The data type of the array to be constructed.
    Returns
    -------
    H : (n, n) ndarray
        The Hadamard matrix.
    Notes
    -----
    .. versionadded:: 0.8.0
    Examples
    --------
    >>> from scipy.linalg import hadamard
    >>> hadamard(2, dtype=complex)
    array([[ 1.+0.j,  1.+0.j],
           [ 1.+0.j, -1.-0.j]])
    >>> hadamard(4)
    array([[ 1,  1,  1,  1],
           [ 1, -1,  1, -1],
           [ 1,  1, -1, -1],
           [ 1, -1, -1,  1]])
    """

    # This function is a slightly modified version of the
    # function contributed by Ivo in ticket #675.

    if n < 1:
        lg2 = 0
    else:
        lg2 = int(math.log(n, 2))
    if 2 ** lg2 != n:
        raise ValueError("n must be an positive integer, and n must be "
                         "a power of 2")

    H = np.array([[1]], dtype=dtype)
    print(lg2)
    # Sylvester's construction
    for i in range(0, lg2):
        print(i)
        print(np.hstack((H, H)))
        print(np.hstack((H, -H)))
        H = np.vstack((np.hstack((H, H)), np.hstack((H, -H))))
        print("H=",H)

    return H

def helmert(n, full=False):
    """
    Create an Helmert matrix of order `n`.
    This has applications in statistics, compositional or simplicial analysis,
    and in Aitchison geometry.
    Parameters
    ----------
    n : int
        The size of the array to create.
    full : bool, optional
        If True the (n, n) ndarray will be returned.
        Otherwise the submatrix that does not include the first
        row will be returned.
        Default: False.
    Returns
    -------
    M : ndarray
        The Helmert matrix.
        The shape is (n, n) or (n-1, n) depending on the `full` argument.
    Examples
    --------
    >>> from scipy.linalg import helmert
    >>> helmert(5, full=True)
    array([[ 0.4472136 ,  0.4472136 ,  0.4472136 ,  0.4472136 ,  0.4472136 ],
           [ 0.70710678, -0.70710678,  0.        ,  0.        ,  0.        ],
           [ 0.40824829,  0.40824829, -0.81649658,  0.        ,  0.        ],
           [ 0.28867513,  0.28867513,  0.28867513, -0.8660254 ,  0.        ],
           [ 0.2236068 ,  0.2236068 ,  0.2236068 ,  0.2236068 , -0.89442719]])
    """
    H = np.tril(np.ones((n, n)), -1) - np.diag(np.arange(n))
    print("H=",H)
    d = np.arange(n) * np.arange(1, n+1)
    print("d=",d)
    H[0] = 1
    d[0] = n
    print("H2=",H)
    print("d2=",d)
    H_full = H / np.sqrt(d)[:, np.newaxis]
    if full:
        return H_full
    else:
        return H_full[1:]


# A = convolution_matrix([-1, 4, -2, 3, 5], 5, mode='same')
# print(A)

B =helmert(5, full=True)
print(B)
# a = [[1,2,3], [1,2,1], [1,1,1]]
# print(signm(a))