# -*- coding: utf-8 -*-
"""
functions to generate spikes from the fLIF model with or without adaptation variable

cf description of the model in the paper "An integrate-and-fire model to generate spike trains with long memory"
by A.Richard, P.Orio and E.Tanré

last modif: 5/09/17
"""

import numpy as np
import numpy.random as npr

def autocov(k,H):
    return 0.5*(np.absolute(k-1)**(2*H) + np.absolute(k+1)**(2*H) - 2*np.absolute(k)**(2*H))


def fracNoise(n,H):
    # generates 2^n points of a fractional noise (ie increments of fractional Brownian motion) of Hurst parameter H.
    # by Davies-Harte method 

    n=min(25,n)
    N=2**n
    C =  np.zeros([2*N])
    k = np.arange(0,N)
    C = np.concatenate((autocov(k,H),np.zeros([1]),autocov(k[:0:-1],H)),axis=0)
    EigV = np.fft.fft(C)
    EigV = EigV.real
    V1 = npr.normal(size=N)   
    V2 = npr.normal(size=N)
    V=(V1[1:N]+V2[1:N]*1j)/np.sqrt(2)
    EigV = np.sqrt(EigV)/np.sqrt(2*N)
    W =  np.concatenate((EigV[0:1]*V1[0:1],EigV[1:N]*V,EigV[N:N+1]*V2[0:1],EigV[N+1:2*N]*V[::-1].conjugate()),axis=0)  
    ZZ = np.fft.fft(W)
    Z = ZZ[0:N].real
    return Z

def fracNoiseN(N,H):
    # generates N points of a fractional noise (ie increments of fractional Brownian motion) of Hurst parameter H.
    # by Davies-Harte method 
    
    N=int(N)
    C =  np.zeros([2*N])
    k = np.arange(0,N)
    C = np.concatenate((autocov(k,H),np.zeros([1]),autocov(k[:0:-1],H)),axis=0)
    EigV = np.fft.fft(C)
    EigV = EigV.real
    V1 = npr.normal(size=N)   
    V2 = npr.normal(size=N)
    V=(V1[1:N]+V2[1:N]*1j)/np.sqrt(2)
    EigV = np.sqrt(EigV)/np.sqrt(2*N)
    W =  np.concatenate((EigV[0:1]*V1[0:1],EigV[1:N]*V,EigV[N:N+1]*V2[0:1],EigV[N+1:2*N]*V[::-1].conjugate()),axis=0)  
    ZZ = np.fft.fft(W)
    Z = ZZ[0:N].real
    return Z