using PyPlot
matplotlib[:rc]("text", usetex = true)
matplotlib[:rc]("figure", figsize = (4, 3), dpi = 300)

function rhs(x, h, t, τs, τc, σ = tanh)
    return -x - σ(h(t - τs) - h(t - τs - τc))
end

using ApproxFun
using NLsolve

function residual!(res, xs, basis, τs, τc, σ)
    f = Fun(basis, xs)
    df = f'
    resid = t -> df(t) - rhs(f(t), f, t, τs, τc, σ)
    res .= coefficients(Fun(resid, basis, length(xs)))
    return nothing
end

τs = 0.2
τc = 1.0
κ = 5.0
σ = x -> tanh(κ * x)

residuals = Array{Float64}(0)
amplitudes = similar(residuals)
n = 20
periods = linspace(0.6, 1.6, 2000)
for period in periods
    s = Fourier(0..period)
    initial_f = Fun(x -> 0.25sin(2pi * x / period), s, n)
    res = nlsolve(
        (res, xs) -> residual!(res, xs, s, τs, τc, σ),
        coefficients(initial_f),
        iterations = 10)
    push!(residuals, res.residual_norm)
    push!(amplitudes, maximum(abs.(res.zero)))
end

plot(periods, residuals, label="residual")
plot(periods, amplitudes, label="amplitude")
yscale("log")
ylim(ymin=1e-8)
legend()
grid()
xlabel("periodicity")

# signal-to-noise ratio (significance)
plot(periods, amplitudes ./ residuals, label="amplitude / residual", "-")
yscale("log")
#ylim(1e-2, 1e2)
#xlim(0.4, 1.6)
legend()
grid()
xlabel("periodicity")
axhline(1, color="black", ls="--")

using Optim

function prominence(period, τs, τc, σ, n=20)
    s = Fourier(0..period)
    initial_f = Fun(x -> 0.25sin(2pi * x / period), s, n)
    res = nlsolve(
        (res, xs) -> residual!(res, xs, s, τs, τc, σ),
        coefficients(initial_f),
        iterations = 10)
    return maximum(abs.(res.zero)) / res.residual_norm
end

# example
res = optimize(0.5, 0.9, GoldenSection()) do period
    return -prominence(period, τs, τc, σ, 20)
end
period = Optim.minimizer(res)
display(res)

residuals = Array{Float64}(0)
amplitudes = similar(residuals)
n = 20
periods = linspace(1.3, 1.35, 1_000)
for period in periods
    s = Fourier(0..period)
    initial_f = Fun(x -> 0.25sin(2pi * x / period), s, n)
    res = nlsolve(
        (res, xs) -> residual!(res, xs, s, τs, τc, σ),
        coefficients(initial_f),
        iterations = 10)
    push!(residuals, res.residual_norm)
    push!(amplitudes, maximum(abs.(res.zero)))
end

plot(periods, residuals, label="residual", ".")
plot(periods, amplitudes, label="amplitude", ".")
yscale("log")
legend()
grid()
xlabel("periodicity")
periods[indmin(residuals)] # approximate location of limit cycle

# example
res = optimize(1.2, 1.4, GoldenSection()) do period
    return -prominence(period, τs, τc, σ, 20)
end
period = Optim.minimizer(res)
display(res)

n = 20
#period = 1.3176176176176175
s = Fourier(0..period)
initial_f = Fun(x -> sin(2pi * x / period), s, n)
res = nlsolve(
    (res, xs) -> residual!(res, xs, s, τs, τc, σ),
    coefficients(initial_f),
    ftol = 1e-6)
u1 = Fun(s, res.zero)
display(res)

t2 = linspace(0, 5, 1000)
plot(t2, u1.(t2))

residuals = Array{Float64}(0)
amplitudes = similar(residuals)
n = 20
periods = linspace(0.6875, 0.6885, 1_000)
for period in periods
    s = Fourier(0..period)
    initial_f = Fun(x -> 0.25sin(2pi * x / period), s, n)
    res = nlsolve(
        (res, xs) -> residual!(res, xs, s, τs, τc, σ),
        coefficients(initial_f),
        iterations = 10)
    push!(residuals, res.residual_norm)
    push!(amplitudes, maximum(abs.(res.zero)))
end

plot(periods, residuals, label="residual", ".")
plot(periods, amplitudes, label="amplitude", ".")
yscale("log")
legend()
grid()
xlabel("periodicity")
periods[indmin(residuals)]

# example
res = optimize(0.6, 0.8, GoldenSection()) do period
    return -prominence(period, τs, τc, σ, 20)
end
period = Optim.minimizer(res)
display(res)

n = 20
#period = 0.6880455455455455
s = Fourier(0..period)
initial_f = Fun(x -> 0.25sin(2pi * x / period), s, n)
res = nlsolve(
    (res, xs) -> residual!(res, xs, s, τs, τc, σ),
    coefficients(initial_f))
display(res)
t2 = linspace(0, 5, 1000)
u2 = Fun(s, res.zero)

plot(t2, u1.(t2))
plot(t2, u2.(t2), "--")
#ylim(-0.22, 0.22)
xlabel(L"time $t$")
ylabel(L"process $A(t)$")
tight_layout()

using ApproxFun
using NLsolve
using Optim
using OrdinaryDiffEq
using DelayDiffEq

function residual!(res, xs, basis, τs, τc, σ)
    f = Fun(basis, xs)
    df = f'
    resid = t -> df(t) - rhs(f(t), f, t, τs, τc, σ)
    res .= coefficients(Fun(resid, basis, length(xs)))
    return nothing
end

τs = 0.2
τc = 1.0 # compare 1.0
κ = 5.0
σ = x -> tanh(κ * x)

# example
res = optimize(0.6, 0.9, GoldenSection()) do period
    return -prominence(period, τs, τc, σ, 20)
end
period = Optim.minimizer(res)
display(res)

#period = 0.6880456129134734
n = 20
s = Fourier(0..period)
initial_f = Fun(x -> 0.25sin(2pi * x / period), s, n)
res = nlsolve(
    (res, xs) -> residual!(res, xs, s, τs, τc, σ),
    coefficients(initial_f))
u2 = Fun(s, res.zero)
display(res)

prob = DDEProblem(
    (u, h, p, t) -> rhs(u, t -> h(p, t), t, τs, τc, σ),
    u2(0.), (p, t) -> u2(t), (0., 100.),
    constant_lags = [τs, τs + τc])
sol = solve(prob, MethodOfSteps(Tsit5()));

figure(figsize = (12, 3))
t2 = linspace(0., 100., 10_000)
plot(t2, sol.(t2))
t1 = linspace(-10., 0., 1_000)
plot(t1, u2.(t1))
xlabel(L"time $t$")
ylabel(L"process $A(t)$")
xlim(-10, 100)
axvline(0, color="black", ls="--")

# example
res = optimize(1.0, 1.7, GoldenSection()) do period
    return -prominence(period, τs, τc, σ, 20)
end
period = Optim.minimizer(res)
display(res)

n = 20
#period = 1.317599
s = Fourier(0..period)
initial_f = Fun(x -> 0.25sin(2pi * x / period), s, n)
res = nlsolve(
    (res, xs) -> residual!(res, xs, s, τs, τc, σ),
    coefficients(initial_f),
    iterations = 10,
    ftol = 1e-5)
u1 = Fun(s, res.zero)
display(res)

prob = DDEProblem(
    (u, h, p, t) -> rhs(u, t -> h(p, t), t, τs, τc, σ),
    u1(0.), (p, t) -> u1(t), (0., 100.),
    constant_lags = [τs, τs + τc])
sol = solve(prob, MethodOfSteps(Tsit5()));

figure(figsize = (12, 3))
t2 = linspace(0., 100., 10_000)
plot(t2, sol.(t2))
t1 = linspace(-10., 0., 1_000)
plot(t1, u1.(t1))
xlabel(L"time $t$")
ylabel(L"process $A(t)$")
xlim(-10, 100)
axvline(0, color="black", ls="--")