# Preliminaries: Teach Julia that functions form a vector space
import Base.+,Base.*,Base.zero
+(f::Function,g::Function) = x -> f(x)+g(x)
*(c::Number,f::Function) = x -> c*f(x)
*(f::Function,c::Number) = x -> c*f(x)
zero(Function) = x -> 0

function T(f::Function)
    return x->f(x+1)
end
    

[T(sin)(5) sin(6)]

# An example check that shifting is linear
# we check at x=5 that 2*T(sin)+3*T(cos) = T(2*sin+3*cos), where T denotes shift by one
[( 2*T(sin) + 3*T(cos) )(5) T( 2*sin + 3*cos )(5)]

using SymPy

x,a,b = Sym.(["x","a","b"]);

diff( a*cos(x) + b*sin(x) ,x)

f=([sin cos]*[0 1;-1 0]*[5,2])[1]
x=rand()
[f(x) 2sin(x)-5cos(x)] ## Check that it gives the right function

x = Sym("x")
f = a*sin(x) + b*cos(x)
Tf = subs(f,x,x+Sym("theta"))

expand_trig(Tf)

A = Rational.([
 1   1   2    6  24
 0  -1  -4  -18 -96
 0   0   1    9  72
 0   0   0   -1 -16
 0   0   0    0   1])./[1 1 2 6 24]

[1 x x^2 x^3 x^4]*A

Int.(inv(A))

D=[0 1 0 0 0
   0 0 2 0 0
   0 0 0 3 0
   0 0 0 0 4
   0 0 0 0 0]

[1 x x^2 x^3 x^4]

[1 x x^2 x^3 x^4]*D

## Now the derivative in a Laguerre Basis
Int.(A\D*A)

# Not working yet
# SymPy.mpmath[:laguerre](4,0,Sym("x")) 

# Derivative in standard basis
D

# Derivative in Laguerre basis
Int.(A\D*A)

A = rand(1:9,4,4)

B = rand(1:9,4,4)

kron(A,B)

X = rand(1:9,4,4)

# The vec operator strings a matrix column wise into one long column
[ kron(A,B)*vec(X)  vec(B*X*A')]

for j=1:4, i=1:4
    V=zeros(Int,4,4)
    V[i,j]=1
    display(V)
end