#!/usr/bin/env python
# * coding: utf8 *
'''
You can set the values in the simulation of dynamic neural fields of the type:
∂ ∂ ⌠ xy
(η + γ + 1) V(x,t) = I(x,t) + ⎮ K(xy) S( V(y, t  ) ) d²y
∂t ∂t ⌡Ω c
where # V(x,t) is the potential of a neural population at position x and time t
# Ω is the domain of integration of size lxl (mm²)
# c is the velocity of an action potential (mm/s)
# γ is the first order derivative value
# η is the second order derivative value
# I(x,t) is the input at position x and time t
# K(x) is the synaptic neighborhood function from [0,√2l] > ℝ
# S(x) is the firing rate of a single neuron from ℝ⁺ > ℝ
The integration is made over the finite 2d domain [l/2,+l/2]x[l/2,+l/2] discretized
into n x n elements considered as a toric surface, during a period of t seconds.
'''
'''Which data to show in the graph.'''
# 1 = show V, potential matrix updates
# 2 = show V0, potential matrix at time=0 (does not update V)
# 3 = show I, input matrix
# 4 = show K, kernel matrix
showData = 1
'''Temporal values.'''
endTime = 1 # simulation duration (a float in seconds or 1:infinity)
dt = 0.001 # temporal discretization (delta t in seconds)
'''Derivative values.'''
gamma = 1.0 # γ first order
eta = 0.0 # η second order
'''Axonal transmission speed value.'''
c = 500.0 # mm/s
'''Field space values  applies to length and also width of square field.'''
l = 30.0 # field size
n = 512 # number of field discretized units
# This sets up the square field, x. Do not change the next 3 lines!!!! **********
import numpy as np
a,b= np.meshgrid(np.arange(l/2.0,l/2.0,l/float(n)),np.arange(l/2.0,l/2.0,l/float(n)))
x = np.sqrt(a**2+b**2)
# Do not change the previous 3 lines! *******************************************
'''This is the field voltage at time=0, V0.
You can delete/add/change variables but you must initialize a V0 that is a numpy array of size n*n.'''
V0 = np.zeros( (n,n) ) # our V at t=0 that will be used in the simulation
'''Noise applied to the voltage at t>=0, noiseVcont.
This variable is multiplied by a matrix of random numbers reset every epoc.
The noiseVcont variable can be None (for no continuous noise) or a numpy array of size n*n.'''
noiseVcont = np.exp((a**2/32.0+b**2/32.0))/(np.pi*32.0) * 0.1 * np.sqrt(dt)
'''This is data for the second order calculation, Uexcite.
You can delete/add/change variables but if eta is not 0.0,
you must initialize a Uexcite that is a numpy array of size n*n.
If eta (above) == 0.0, then Uexcite can be None, but this is not neccesary.'''
Uexcite = np.zeros((n,n)) # # AXEL: set Uexcite to zero
'''This is the input from external source, I.
You can delete/add/change variables but you must initialize an I that uses x.'''
Gamma = 20.0 # Γ value
sigma = 5.65685425 # σ value, cannot be 0.0 (division by 0)
# Gaussian value for I of the form I_ + Γ *exp(x² /σ²) / (π.σ²)
I = Gamma * (np.exp(1 * x**2 / sigma**2) / (sigma**2 * np.pi))
'''This is the synaptic connectivity kernel, K.
You can delete/add/change variables but you must initialize a K that uses x.'''
K = 4*np.exp(x/3) / (18*np.pi)
'''This is the firing rate, S.
You can delete/add/change variables but you must keep the function name S and return the firing rate.'''
def updateS(V): # V is the passed in field voltage
S0 = 1.0 # S: maximum frequency
alpha = 10000.0 # α: steepness at the threshold
theta = 0.005 # θ: firing threshold
return S0 / (1.0 + np.exp(1*alpha*(Vtheta)))
