text:
# # Instructions for Assignment03
#
# ## Introduction
#
# Note: this file is runnable in its current state,
# but is incomplete.
# You can run the file from the command line,
# or use the VS Code julia extension to run individual lines.

## To view additional instructions for this assignment, #jl
## see https://wellesley-bisc195.github.io/BISC195.jl/stable/Assignments/assignment03 #jl

## To view a rendered version of this document, #jl
## see https://wellesley-bisc195.github.io/BISC195.jl/stable/Assignments/assignment03_code #jl

# ## Writing real Code
#
# In assignments 1 and 2, variable and function names
# were often things like `question1` and `question2`.
# From now on, we'll use more informative function and variable names
# so that our code is "self documenting."
#
# We'll also continue to use doc strings to help understand
# the specifications required for our functions.

# ## Understanding DNA sequences
#
# We've already done a fair amount of work in assignment 2
# and in lesson 3 to make some functions for understanding DNA sequences.
# Below, I've put a couple of function signatures with docstrings,
# but you can (and should!) copy the necessary functionality
# out of the functions you've already defined, if applicable.
#
# !!! tip
#     If you defined those functions in the julia REPL,
#     you can go find them from the command line or the julia REPL!
#     In VS Code, open ~/.julia/logs/repl_history.jl.

# ## Question 1 - a `complement` function

"""
    compliment(base)

Get the DNA compliment of the provided base:

    A <-> T
    G <-> C

Accepts `String` or `Char`, but always returns `Char`.
If a valid base is not provided, the function throws an error.

Examples
≡≡≡≡≡≡≡≡≡≡

    julia> compliment('A')
    'T'

    julia> compliment("G")
    'C'

    julia> compliment("T")
    'A'

    julia> compliment('C')
    'G'
"""
function compliment(base)
    ## See Lesson 3 for more info
end

# ## Question 2 - some boolean functions

"""
    ispurine(base)

A boolean function that returns `true` if the base is a purine (A or G)
and `false` if it is not.
The function only supports bases A, C, G, and T (throws an error for other values).
Accepts `String` or `Char`.

Examples
=========

    julia> ispurine('A')
    true

    julia> ispurine("C")
    false

    julia> ispurine('B')
    Error: "Base B not supported")
"""
function ispurine(base)
    ## We haven't made this before, but you should have all the pieces
end

"""
    ispyrimidine(base)

A boolean function that returns `true` if the base is a pyrimidine (C or T)
and `false` if it is not.
The function only supports bases A, C, G, and T (throws an error for other values).
Accepts `String` or `Char`.

Examples
=========

    julia> ispyrimidine('G')
    false

    julia> ispyrimidine("T")
    true

    julia> ispyrimidine('X')
    Error: "Base X not supported"
"""
function ispyrimidine(base)
    ## This is the strict opposite of `ispurine`.
    ## In principle, you can write this in one line - remember `!` means `NOT`.
    ## Eg `isuppercase(x)` means the same thing as `!islowercase(x)`
end

# ## Question 3 - Using boolean functions for composition
#
# For the following function, you should not need to re-write
# the logic checking what kind of base this is.
# You've already written it, and it's in a convenient function!
#
# A big part of programming is re-use;
# if you find yourself writing the same code multiple times,
# you should probably put it in a function and call that instead!

"""
    base_type(base)

Determines whether a base is a purine (A or G) or pyrimidine (T or C),
and returns a `String`.

Examples
≡≡≡≡≡≡≡≡≡≡

    julia> base_type("G")
    "purine"

    julia> base_type('C')
    "pyrimidine"

    julia> base_type('Z')
    Error: "Base Z not supported"

    julia> x = base_type('A'); println(x)
    purine
"""
function base_type(base)
    ## Note: this is different than the `base_type()` we defined in the lesson.
    ## Here, we want a fruitful function that returns the value rather than `print`ing it.
    ## Also, there's no need to re-write the logic. If your `ispurine` / `ispyrimidine` functions work,
    ## you can use them here.
end

# ## Question 4 - Modifying arguments instead of adding a bunch of logic
#
# One thing that none of our functions can do so far
# is to accept lowercase sequences.
# Most of the time, DNA sequences are written with uppercase letters,
# but we may not be able to count on that.
#
# If we want to be able to accept lowercase strings,
# one possibility would be to add additional logic, Eg
#
# `if base == 'G' || base == "G" || base = 'g' || base == "g"`
#
# But that's a lot of typing -
# especially considering we'd have to do this for every base!
# In programming, it's OK to be lazy (in fact, it's often better)!
# Instead, we can modify the parameters to be formatted the way we expect.
# For example, the `uppercase()` function takes a `String` or `Char`
# and returns the uppercase representation of it.
#
# ```julia
# julia> uppercase('a')
# 'A': ASCII/Unicode U+0041 (category Lu: Letter, uppercase)
#
# julia> uppercase("attc")
# "ATTC"
# ```
#
# Note that, if you simply use this function on a variable or parameter,
# it will remain unchanged (the function doesn't "mutate" its argument),
# but you can re-assign the variable or parameter.
#
# ```julia
# julia> seq = "attcgc"
# "attcgc"

# julia> uppercase(seq)
# "ATTCGC"

# julia> seq
# "attcgc"

# julia> seq = uppercase(seq)
# "ATTCGC"

# julia> seq
# "ATTCGC"
# ```
#
# Some julia functions *can* mutate their aruments -
# we'll encounter some of those soon.

"""
    gc_content(sequence)

Calculates the GC ratio of a DNA sequence.
The GC ratio is the total number of G and C bases divided by the total length of the sequence.
For more info about GC content, [see here](https://en.wikipedia.org/wiki/GC-content).

Example
≡≡≡≡≡≡≡≡≡≡

    julia> gc_content("AATG")
    0.25

    julia> gc_content("cccggg") * 100
    100.0

    julia> gc_content("ATta")
    0.0
"""
function gc_content(sequence)
    ## Start with the same code as `question3()` from assignment 2.
    ## only a small modification is necessary to make this work.
end

# ## Question 5 - Incremental development
#
# Now that you've learned how to do this,
# go back and modify the functions defined in questions 1-3
# so that they are able to take lowercase arguments.
# You should not copy the code below this line to modify it,
# you should modify it in-place.


text:
abstract type ContextDist end
abstract type RewardDist end
abstract type BanditPolicy end

mutable struct Trajectory
    t::Int64
    X::Matrix{Float64}
    a::Vector{Int64}
    r::Vector{Float64}
    total_regret::Vector{Float64}
end

function get_X(traj::Trajectory, a::Int64, view::Bool=false)
    if view
        return traj.X[traj.a .== a, :]
    else
        X = @view traj.X[traj.a .== a, :]
        return X
    end
end

function get_r(traj::Trajectory, a::Int64, view::Bool=false)
    if view
        return traj.r[traj.a .== a]
    else
        r = @view traj.r[traj.a .== a]
        return r
    end
end

function Trajectory(d::Int64)
    X = Matrix{Float64}(undef, 0, d)
    a = Vector{Float64}(undef, 0)
    r = Vector{Float64}(undef, 0)
    total_regret = Vector{Float64}(undef, 0)
    return Trajectory(0, X, a, r, total_regret)
end

function add_memory!(traj::Trajectory, T::Int64)
    X1 = Matrix{Float64}(undef, T, size(traj.X, 2))
    traj.X = vcat(traj.X, X1)
    push!(traj.a, zeros(Int64, T)...)
    push!(traj.r, zeros(Float64, T)...)
    push!(traj.total_regret, zeros(Float64, T)...)
end

function extend!(traj::Trajectory, x::Vector{Float64}, a::Int64, r::Float64, total_regret::Float64)
    traj.t += 1
    traj.X[traj.t, :] = x
    traj.a[traj.t] = a
    traj.r[traj.t] = r
    traj.total_regret[traj.t] = total_regret
end

# function extend!(
#     traj::Trajectory,
#     X::Matrix{Float64},
#     a::Vector{Int64},
#     r::Vector{Float64},
#     total_regret::Vector{Float64}
# )
#     n = length(r)
#     traj.X[(traj.t + 1):(traj.t + n), :] = X
#     traj.a[(traj.t + 1):(traj.t + n)] = a
#     traj.r[(traj.t + 1):(traj.t + n)] = r
#     traj.total_regret[(traj.t + 1):(traj.t + n)] = total_regret
#     traj.t += n
# end

function simulate!(
    T::Int64,
    traj::Trajectory,
    cdist::ContextDist,
    rdist::RewardDist,
    policy::BanditPolicy;
    verbose::Bool = true
)
    add_memory!(traj, T)
    start = traj.t + 1
    finish = traj.t + T
    X = sample(T, cdist)

    for t in start:finish
        if verbose print("\rTime step $t/$finish") end
        a = choose(X[t - start + 1, :], policy)
        r = sample(X[t - start + 1, :], a, rdist)
        regret = expected_regret(X[t - start + 1, :], a, rdist)
        total_regret = t == 1 ? regret : regret + traj.total_regret[t - 1]
        extend!(traj, X[t - start + 1, :], a, r, total_regret)
        update!(policy, X[t - start + 1, :], a, r)
    end
end

ext: jl


text:
for T in (Float32, Float64)
    @testset "$(T): Signed zero" begin
        @test isequal(add_up(zero(T), zero(T)), zero(T))
        @test isequal(add_up(zero(T), -zero(T)), zero(T))
        @test isequal(add_up(-zero(T), zero(T)), zero(T))
        @test isequal(add_up(-zero(T), -zero(T)), -zero(T))

        @test isequal(add_down(zero(T), zero(T)), zero(T))
        @test isequal(add_down(zero(T), -zero(T)), -zero(T))
        @test isequal(add_down(-zero(T), zero(T)), -zero(T))
        @test isequal(add_down(-zero(T), -zero(T)), -zero(T))

        @test isequal(sub_up(zero(T), zero(T)), zero(T))
        @test isequal(sub_up(zero(T), -zero(T)), zero(T))
        @test isequal(sub_up(-zero(T), zero(T)), -zero(T))
        @test isequal(sub_up(-zero(T), -zero(T)), zero(T))

        @test isequal(sub_down(zero(T), zero(T)), -zero(T))
        @test isequal(sub_down(zero(T), -zero(T)), zero(T))
        @test isequal(sub_down(-zero(T), zero(T)), -zero(T))
        @test isequal(sub_down(-zero(T), -zero(T)), -zero(T))

        @test isequal(mul_up(zero(T), zero(T)), zero(T))
        @test isequal(mul_up(zero(T), -zero(T)), -zero(T))
        @test isequal(mul_up(-zero(T), zero(T)), -zero(T))
        @test isequal(mul_up(-zero(T), -zero(T)), zero(T))

        @test isequal(mul_down(zero(T), zero(T)), zero(T))
        @test isequal(mul_down(zero(T), -zero(T)), -zero(T))
        @test isequal(mul_down(-zero(T), zero(T)), -zero(T))
        @test isequal(mul_down(-zero(T), -zero(T)), zero(T))

        @test isequal(div_up(zero(T), zero(T)), T(NaN))
        @test isequal(div_up(zero(T), -zero(T)), T(NaN))
        @test isequal(div_up(-zero(T), zero(T)), T(NaN))
        @test isequal(div_up(-zero(T), -zero(T)), T(NaN))

        @test isequal(div_down(zero(T), zero(T)), T(NaN))
        @test isequal(div_down(zero(T), -zero(T)), T(NaN))
        @test isequal(div_down(-zero(T), zero(T)), T(NaN))
        @test isequal(div_down(-zero(T), -zero(T)), T(NaN))

        @test isequal(sqrt_up(zero(T)), zero(T))
        @test isequal(sqrt_down(-zero(T)), -zero(T))
    end
end

@testset "Corner cases" begin
    # TODO
    # Add tests for Float32

    @testset "twosum intermediate overflow" begin
        # http://verifiedby.me/adiary/09
        a = 3.5630624444874539e+307
        b = -floatmax(Float64)
        x = a + b
        @test isfinite(x)
        tmp = x - a
        @test isinf(tmp)

        @test isequal(add_up(a, b), -1.4413868904135702e308)
        @test isequal(add_down(a, b), -1.4413868904135704e308)
    end

    @testset "twoprod intermediate overflow" begin
        # http://verifiedby.me/adiary/09
        function split(a)
            tmp = a * (2.0^27 + 1.0)
            x = tmp - (tmp - a)
            y = a - x
            x, y
        end
        a = 6.929001713869936e+236
        b = 2.5944475251952003e+71
        x = a * b
        @test isfinite(x)
        a1, _ = split(a)
        b1, _ = split(a)
        tmp = a1 * b1
        @test isinf(tmp)

        @test isequal(mul_up(a, b), floatmax(Float64))
        @test isequal(mul_down(a, b), prevfloat(floatmax(Float64)))
    end
end


text:
function test_convergence_p3_dg()
  @testset "---- P3 Conservative DG Convergence Tests -----" begin
    start_dir = pwd()

    cd("./m1")
    fname = "input_vals_vortex3.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)
    fname = "input_vals_vortex4.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)

    cd("../m2")
    fname = "input_vals_vortex3.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)
    fname = "input_vals_vortex4.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)

    cd("..")

    slope = calc_line()
    println("slope = ", slope)

    data = readdlm("err_data.dat")
    err_vals = data[:, 2]
    #println("err_vals = ", err_vals)

    slope_val = 3.84
    slope_margin = 0.1

    @test  slope  > slope_val - slope_margin
    @test  slope  < slope_val + slope_margin

    err_val = 6.87e-5
    slope_fac = 1.25
    println("err_vals[1] = ", err_vals[1])
    @test  err_vals[1]  > err_val/slope_fac
    @test  err_vals[1]  < err_val*slope_fac

  end

  return nothing
end


text:
@info "Loading DataFrames support for GBIF.jl"

format_gbif_entity(::Missing) = missing
format_gbif_entity(t::Pair{String,Int64}) = t.first

import .DataFrames: DataFrame

function DataFrames.DataFrame(records::GBIFRecords)
  output = DataFrame(
    key = Int64[],
    name = Union{Missing,AbstractString}[],
    dataset = Union{Missing,AbstractString}[],
    published_in = Union{Missing, AbstractString}[],
    observed_in = Union{Missing, AbstractString}[],
    latitude = Union{Missing, AbstractFloat}[],
    longitude = Union{Missing, AbstractFloat}[],
    date = Union{Missing, DateTime}[],
    rank = Union{Missing, AbstractString}[],
    observer = Union{Missing, AbstractString}[],
    license = Union{Missing, AbstractString}[],
    kingdom = Union{Missing,String}[],
    phylum = Union{Missing,String}[],
    class = Union{Missing,String}[],
    order = Union{Missing,String}[],
    family = Union{Missing,String}[],
    genus = Union{Missing,String}[],
    species = Union{Missing,String}[]
  )
  for occ in records
    push!(
      output,
      (
        occ.key, occ.taxon.name, occ.dataset, occ.publishingCountry,
        occ.countryCode, occ.latitude, occ.longitude,
        occ.date, occ.rank, occ.observer, occ.license,
        format_gbif_entity(occ.taxon.kingdom),
        format_gbif_entity(occ.taxon.phylum),
        format_gbif_entity(occ.taxon.class),
        format_gbif_entity(occ.taxon.order),
        format_gbif_entity(occ.taxon.family),
        format_gbif_entity(occ.taxon.genus),
        format_gbif_entity(occ.taxon.species)
      )
    )
  end
  return output
end


text:
using Plots, QuasiGeostrophy
# 1D Tests
Ω = S¹(0,2π)
Nx = 2^8;
grid = FourierGrid(Nx, Ω)
ϕ, ϕ2 = create_fields(names = ("ϕ", "ϕ2"), grid = grid)
∂x = create_operators(grid)
#  initialize fields with nontrivial data
#  Heuns Method
x = grid.grid[1]

scale = round(Int, 1/(x[2]-x[1]))
ν = 10^4
κ = 1e-4
iv = @. tanh( ν * (x - π/2) ) * tanh( ν * (x - 3π/2) )
ϕ(iv)
ϕ2(iv)
plot(ϕ)
Δt = 0.01 / scale * 2π
t = randn(1) .* 0

for i in 1:100*scale
    Δ = ∂x^2
    f = ∂x(ϕ) + κ * Δ(ϕ)
    nϕ = ϕ + Δt * f
    g  = ∂x(nϕ)  + κ * Δ(nϕ)
    nϕ = ϕ + Δt/2 * ( f +  g)
    ϕ.data .= nϕ.data
    if i%(scale) == 0
        p1 = plot(ϕ)
        p2 = spectrum(ϕ)
        display(plot(p1,p2))
        sleep(0.01)
    end
    t[1] += Δt
end
norm(ϕ-ϕ2)/norm(ϕ2)

##
# 2D Test

Ω = S¹(0,2π) × S¹(0,2π)
Nx = Ny =  2^6;
grid = FourierGrid((Nx, Ny), Ω)
fieldnames = ("ϕ", "ϕ2", "u", "v", "ψ")
ϕ, ϕ2, u, v, ψ = create_fields(names = fieldnames, grid = grid)
∂x, ∂y =  create_operators(grid)
# initialize fields with nontrivial data
#  Heuns Method
x, y = grid.grid

scale = round(Int, 1/(x[2]-x[1]))
ν = 10000
κ = 0e-2
iv = @. tanh( ν*(x - π/2)) * tanh(ν*(x - 3π/2) ) * tanh( ν*(y - π/2)) * tanh(ν*(y - 3π/2) )
# iv = Reshape(iv, (Nx,Ny))
ϕ(iv)
ϕ2(iv)
plot(ϕ)
Δt = 0.01 / scale * 2π
t = randn(1) .* 0
u( (x .+ y) .* 0 .+1 )
v( (x .+ y) .* 0 .+1 )
ψ(@. sin(x/2) * sin(y/2))
# u = ∂y(ψ)
# v = - ∂x(ψ)
norm(∂x(u) + ∂y(v))
for i in 1:100*scale
    Δ = (∂x^2 + ∂y^2)^1
    f = u * ∂x(ϕ) + v * ∂y(ϕ) + κ * Δ(ϕ)
    nϕ = ϕ + Δt * f
    nϕ = ϕ + Δt/2 * ( f +  u * ∂x(nϕ) + v * ∂y(nϕ)  + κ * Δ(nϕ))
    ϕ.data .= nϕ.data
    if i%(scale) == 0
        display(plot(ϕ))
        sleep(0.01)
    end
    t[1] += Δt
end
for i in 1:100*scale
    Δ = (∂x^2 + ∂y^2)^1
    f =  -u * ∂x(ϕ) + -v * ∂y(ϕ) + κ * Δ(ϕ)
    nϕ = ϕ + Δt * f
    nϕ = ϕ + Δt/2 * ( f +  -u * ∂x(nϕ) + -v * ∂y(nϕ)  + κ * Δ(nϕ))
    ϕ.data .= nϕ.data
    if i%(scale) == 0
        display(plot(ϕ))
        sleep(0.01)
    end
    t[1] += Δt
end
norm(ϕ-ϕ2) / norm(ϕ2)

max_stars_repo_name: upiterbarg/QuasiGeostrophy

max_issues_repo_name: sandreza/QuasiGeostrophy

max_forks_repo_name: upiterbarg/QuasiGeostrophy


text:
using DelimitedFiles
using DifferentialEquations
using ProgressMeter
include("PSS_vern9.jl")
 using PyCall
pygui(:qt)
using PyPlot
pygui(true)
max_hit=2
Hamil(XX,YY,QQ,PP)=( (QQ-XX)^2+(PP-YY)^2 )*( (QQ+XX)^2+(PP+YY)^2 )/((PP^4+2*PP^2*(XX^2-1)+(1+XX^2)^2 )*(QQ^4+2*QQ^2*(YY^2+1)+(YY^2-1)^2 ))
ODE2(z,w)=conj(  im * z.*( 1 ./(w.^2-z.^2)+1 ./(1+z.^2) ))
ODE1(z,w)=conj(  im * w.*( 1 ./(z.^2-w.^2)+1 ./(1+w.^2) ))
function Eq_of_M(du,u,p,t)
    du[1]=real(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
    du[2]=imag(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
    du[3]=real(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
    du[4]=imag(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
    du[5]=0
return
end
condition(u,t,integrator)= u[5]>max_hit
function condition2(u,t,integrator) # Event when event_f(u,t) == 0
u[1]
end
affect!(integrator) = terminate!(integrator)
function affect2!(integrator)
integrator.u[5]=integrator.u[5]+1
end
cb2 = ContinuousCallback(condition2,affect2!,nothing)
cb1 = DiscreteCallback(condition,affect!)

cb=CallbackSet(cb2,cb1)
H_test(u)=( (u[3]-u[1])^2+(u[2]-u[4])^2 )*( (u[3]+u[1])^2+(u[2]+u[4])^2 )/ ((u[2]^4+2*u[2]^2*(u[1]^2-1)+(1+u[1]^2)^2 )*(u[3]^4+2*u[3]^2*(u[4]^2+1)+(u[4]^2-1)^2 ))

Energy=.25
t_end=100.
n_iter_P=500
n_iter_Q=501

# ArrP=range(1.1,stop=2,length=n_iter_P)
# ArrQ=range(1e-10,stop=1,length=n_iter_Q)
# ArrP=range(0,stop=.1,length=n_iter_P)
# ArrQ=range(1e-10,stop=.1,length=n_iter_Q)
ArrP=zeros(n_iter_P)
ArrQ=zeros(n_iter_Q)
Q_bound=zeros(0)
P_bound=zeros(0)

Q_0=zeros(0)
P_0=zeros(0)

Q_1=zeros(0)
P_1=zeros(0)

Q_2=zeros(0)
P_2=zeros(0)

Q_3=zeros(0)
P_3=zeros(0)

Q_4=zeros(0)
P_4=zeros(0)

Q_5=zeros(0)
P_5=zeros(0)

@showprogress 1 "Computing..." for j=1:n_iter_P
    for k=1:n_iter_Q
        # ArrQ[k]=.1*rand()
        # ArrP[j]=.01*rand()
        K=PSS_vern9(ArrQ[k],ArrP[j], Energy, t_end)
        if K==0
            push!(Q_0,ArrQ[k])
            push!(P_0,ArrP[j])
        elseif K==1
            push!(Q_1,ArrQ[k])
            push!(P_1,ArrP[j])
        elseif K==2
            push!(Q_2,ArrQ[k])
            push!(P_2,ArrP[j])
        elseif K==3
            push!(Q_3,ArrQ[k])
            push!(P_3,ArrP[j])
        elseif K==4
            push!(Q_4,ArrQ[k])
            push!(P_4,ArrP[j])

        elseif K==-1
            push!(Q_bound,ArrQ[k])
            push!(P_bound,ArrP[j])
        end
    end
end

figure()
# axis([ -2.5, 2.5,-1.1, 1.1])
Qbound=vcat(Q_bound,-Q_bound,Q_bound,-Q_bound)
Pbound=vcat(P_bound,P_bound,-P_bound,-P_bound)
Tbound=hcat(Qbound,Pbound)
plot(Pbound,Qbound ,color="k",".",markersize=1, markeredgewidth=.1)

Q0=vcat(Q_0,-Q_0,Q_0,-Q_0)
P0=vcat(P_0,P_0,-P_0,-P_0)
T0=hcat(Q0,P0)
plot(P0,Q0 ,color="r",".",markersize=1, markeredgewidth=.1)

Q1=vcat(Q_1,-Q_1,Q_1,-Q_1)
P1=vcat(P_1,P_1,-P_1,-P_1)
T1=hcat(Q1,P1)
plot(P1,Q1 ,color="r",".",markersize=1, markeredgewidth=.1)

Q2=vcat(Q_2,-Q_2,Q_2,-Q_2)
P2=vcat(P_2,P_2,-P_2,-P_2)
T2=hcat(Q2,P2)
plot(P2,Q2 ,color="g",".",markersize=1, markeredgewidth=.1)

Q3=vcat(Q_3,-Q_3,Q_3,-Q_3)
P3=vcat(P_3,P_3,-P_3,-P_3)
T3=hcat(Q3,P3)
plot(P3,Q3 ,color="m",".",markersize=1, markeredgewidth=.1)

Q4=vcat(Q_4,-Q_4,Q_4,-Q_4)
P4=vcat(P_4,P_4,-P_4,-P_4)
T4=hcat(Q4,P4)
plot(P4,Q4 ,color="c",".",markersize=1, markeredgewidth=.1)

Q5=vcat(Q_5,-Q_5,Q_5,-Q_5)
P5=vcat(P_5,P_5,-P_5,-P_5)
T5=hcat(Q5,P5)
plot(P5,Q5 ,color="y",".",markersize=1, markeredgewidth=.1)

outfile_0 = "Vern9_outfile0_"*string(Energy)*"_nP_"*string(n_iter_P)*"_nQ_"*string(n_iter_Q)*".dat"
outfile_1 = "Vern9_outfile1_"*string(Energy)*"_nP_"*string(n_iter_P)*"_nQ_"*string(n_iter_Q)*".dat"
outfile_2 = "Vern9_outfile2_"*string(Energy)*"_nP_"*string(n_iter_P)*"_nQ_"*string(n_iter_Q)*".dat"
outfile_3 = "Vern9_outfile3_"*string(Energy)*"_nP_"*string(n_iter_P)*"_nQ_"*string(n_iter_Q)*".dat"
outfile_4 = "Vern9_outfile4_"*string(Energy)*"_nP_"*string(n_iter_P)*"_nQ_"*string(n_iter_Q)*".dat"
outfile_5 = "Vern9_outfile5_"*string(Energy)*"_nP_"*string(n_iter_P)*"_nQ_"*string(n_iter_Q)*".dat"
outfile_bound = "Vern9_outfile1_bound"*string(Energy)*"_nP_"*string(n_iter_P)*"_nQ_"*string(n_iter_Q)*".dat"

a = open(outfile_0, "w")
b = open(outfile_1, "w")
c = open(outfile_2, "w")
d = open(outfile_3, "w")
e = open(outfile_4, "w")
f = open(outfile_5, "w")
g = open(outfile_bound, "w")

writedlm(outfile_0, T0)
writedlm(outfile_1, T1)
writedlm(outfile_2, T2)
writedlm(outfile_3, T3)
writedlm(outfile_4, T4)
writedlm(outfile_5, T5)

writedlm(outfile_bound, Tbound)

close(a)
close(b)
close(c)
close(d)
close(e)
close(f)
close(g)


text:
module vfld

using Base
import Dates
import DataFrames
import SQLite

include("arg_handler.jl")
include("vfld_to_sqlite.jl")
include("date_naming.jl")

import .arg_handler
import .vfld_to_sqlite
import .date_naming

function __init__()

end

function main(args)
    cmd_message = arg_handler.main_args(args)
    cmd_message[1]==="vfld_to_sqlite" ? vfld_to_sqlite.make_sqlite(cmd_message) : nothing
    cmd_message[1]==="date_naming" ? date_naming.rename(cmd_message) : nothing
end

end


text:
using Base.Test
using Klara

fnames = (
  :value,
  :loglikelihood,
  :logprior,
  :logtarget
)

println("    Testing BasicDiscUnvParameterState constructors and methods...")

v = Int64(-2)
s = BasicDiscUnvParameterState(v)

@test eltype(s) == (Int64, Float64)
@test s.value == v
for i in 2:4
  @test isnan(getfield(s, fnames[i]))
end
@test s.diagnosticvalues == []
@test s.diagnostickeys == Symbol[]

@test diagnostics(s) == Dict{Symbol, Any}()

v = Int32(4)
dv = [true]
s = BasicDiscUnvParameterState(v, [:accept], Float32, dv)

@test eltype(s) == (Int32, Float32)
@test s.value == v
for i in 2:4
  @test isnan(getfield(s, fnames[i]))
end
@test s.diagnosticvalues == dv
@test s.diagnostickeys == [:accept]

@test diagnostics(s) == Dict(:accept => dv[1])

z = BasicDiscUnvParameterState(Int64(-3), [:accept], Float64, [false])
w = deepcopy(z)
@test isequal(z, w)


text:
module CirculantCov

using Dierckx: Spline1D
using ApproxFun: Fun, Jacobi
using FastTransforms: jac2cheb
using FFTW: plan_fft
using LinearAlgebra: mul!

# J
# ==================================================

"""
Jperm(ℓ::Int, n::Int) return the column number in the J matrix U^2
where U is unitary FFT. The J matrix looks like this:

|1   0|
|  / 1|
| / / |
|0 1  |

"""
function Jperm end

function Jperm(ℓ::Int, n::Int)
    @assert 1 <= ℓ <= n
    ℓ==1 ? 1 : n - ℓ + 2
end

# geometric methods with angles and geodesics and periodize
# ==================================================

function periodize(f::Vector{T}, freq_mult::Int) where {T}
    n = length(f)
    nfm = n÷freq_mult
    @assert nfm == n//freq_mult
    f′ = sum( circshift(f, k*nfm) for k=0:freq_mult-1)
    f′[1:nfm]
end

in_0_2π(φ) = rem2pi(φ, RoundDown)

counterclock_Δφ(φstart, φstop) = in_0_2π(φstop - φstart)

function fullcircle(φ::AbstractVector)
    Δφpix  = counterclock_Δφ(φ[1], φ[2])
    Δφspan = counterclock_Δφ(φ[1], φ[end]) + Δφpix
    # The extra Δφpix makes Δφspan measure the (angular) distance between the
    # left boundary of the starting pixel and the right boundary of the ending pixel

    @assert div(2π, Δφspan, RoundNearest) ≈ 2π / Δφspan
    freq_mult = Int(div(2π, Δφspan, RoundNearest))

    nφ2π = length(φ)*freq_mult
    φ2π  = @. in_0_2π(φ[1] + 2π * (0:nφ2π-1) / nφ2π)

    return φ2π, freq_mult
end

"""
# slated for removal ...
`fraccircle(∂φstart, ∂φstop, nφ)` specifies a uniform grid on a
contiguous interval of azumimuth. `∂φstart` begins the interval. Moving counter
clockwise (looking down from the north pole) to `∂φstop`.
Only integer fractions are allowed and both `∂φstart`, `∂φstop` must be `≥ 0`.

Note: `(∂φstart, ∂φstop) = (5.3, 1.0) ≡ (5.3, 1.0 + 2π)`
"""
# function fraccircle(∂φstart, ∂φstop, nφ)
#     ∂φstart′, ∂φstop′ = in_0_2π(∂φstart), in_0_2π(∂φstop)
#     Δφspan = counterclock_Δφ(∂φstart′, ∂φstop′)
#     φ  = @. in_0_2π(∂φstart′ + Δφspan * (0:nφ-1) / nφ)
#     return φ
# end

"""
`geoβ(θ₁, θ₂, φ₁, φ₂)` -> Geodesic between two spherical points at
(θ₁, φ₁) and (θ₂, φ₂).
"""
function geoβ(θ₁, θ₂, φ₁, φ₂)
    sΔθ½, sΔφ½ = sin((θ₁ - θ₂)/2), sin((φ₁ - φ₂)/2)
    2asin(√(sΔθ½^2 + sin(θ₁)*sin(θ₂) * sΔφ½^2))
end

"""
`cosgeoβ(θ₁, θ₂, φ₁, φ₂)` -> Cosine of the geodesic between two spherical points
at (θ₁, φ₁) and (θ₂, φ₂).
"""
function cosgeoβ(θ₁, θ₂, φ₁, φ₂)
    cos(θ₁-θ₂) - sin(θ₁)*sin(θ₂)*(1-cos(φ₁-φ₂))/2
end

"""
`βsingrid(ngrid, period)` creats a grid of points on [0,period) ⊂ [0,π)
with extra density near the endpoints
"""
βsingrid(ngrid, period)   = shift_scale_sin.(range(0,π,length=ngrid), period)

"""
`shift_scale_sin(β,period)` internal function used in `βsingrid`
"""
shift_scale_sin(β,period) = period * (sin(π*β/period - π/2) + 1) / 2

# Methods for generating some useful polar grids
# ================================================

function θ_healpix_j_Nside(j_Nside)
    0 < j_Nside < 1  ? acos(1-abs2(j_Nside)/3)      :
    1 ≤ j_Nside ≤ 3  ? acos(2*(2-j_Nside)/3)        :
    3 < j_Nside < 4  ? acos(-(1-abs2(4-j_Nside)/3)) :
    error("argument ∉ (0,4)")
end

θ_healpix(Nside) = θ_healpix_j_Nside.((1:4Nside-1)/Nside)

function θ_grid(;θspan::Tuple{T,T}, N::Int, type=:equiθ) where T<:Real
    @assert N > 0
    @assert 0 < θspan[1] < θspan[2] < π

    if type==:equiθ
        θgrid″ = range(θspan[1], θspan[2], length=N+2)
    elseif type==:equicosθ
        θgrid″ = acos.(range(cos(θspan[2]), cos(θspan[1]), length=N+2)[end:-1:1])
    elseif type==:healpix
        @warn """
            When `type` argument is set to `:healpix` the parameter `N` corresponds
            to Healpix `Nside`, _not_ the number of θ grid points within the interval
            specified by `θspan` as it does when `type ∈ {:equiθ, :equicosθ}`.
            """
        θgrid′ = θ_healpix(N)
        θgrid″ = θgrid′[θspan[1] .≤ θgrid′ .≤ θspan[2]]
    else
        error("`type` argument variable is not a valid option. Choose from `type ∈ {:equiθ, :equicosθ, :healpix}`")
    end

    # θgrid″ subsets θgrid′ to be within θspan
    # δ½south″ and δ½north″ are the arclength midpoints to the adjacent pixel
    δ½south″ = (circshift(θgrid″,-1)  .- θgrid″) ./ 2
    δ½north″ = (θgrid″ .- circshift(θgrid″,1)) ./ 2

    # now restrict to the interior of the range of θgrid″
    θ       = θgrid″[2:end-1]
    δ½south = δ½south″[2:end-1]
    δ½north = δ½north″[2:end-1]

    # These are the pixel boundaries along polar
    # so length(θ∂) == length(θ)+1
    θ∂ = vcat(θ[1] .- δ½north[1], θ .+ δ½south)

    return θ, θ∂
end

"""
`φ_grid(;φspan::Tuple{T,T}, N::Int)` specifies a uniform grid on a
contiguous interval of azumimuth. `∂φstart` begins the interval. Moving counter
clockwise (looking down from the north pole) to `∂φstop`.
Only integer fractions are allowed and both `∂φstart`, `∂φstop` must be `≥ 0`.

Note: `(∂φstart, ∂φstop) = (5.3, 1.0) ≡ (5.3, 1.0 + 2π)`
"""
function φ_grid(;φspan::Tuple{T,T}, N::Int) where T<:Real
    ∂φstart′, ∂φstop′ = in_0_2π(φspan[1]), in_0_2π(φspan[2])
    Δφspan = counterclock_Δφ(∂φstart′, ∂φstop′)
    φ∂  = @. in_0_2π(∂φstart′ + Δφspan * (0:N) / N)
    Δφ  = Δφspan / N
    φ   = φ∂[1:end-1] .+ Δφ / 2
    return φ, φ∂
end

# these are the generic versions ...
# you just need to define types for Γθ₁θ₂φ₁φ⃗, Cθ₁θ₂φ₁φ⃗
# for dispatch
# ==================================================

# Note: the reason we restrict to φ::AbstractVector is that Spline1D's are optimized for
# Union{Vector, Number} so in general it is better to broadcast
# via via whole columns

# overload this
function γθ₁θ₂φ⃗_ξθ₁θ₂φ⃗(θ₁::Real, θ₂::Real, φ::AbstractVector, Γθ₁θ₂φ₁φ⃗, Cθ₁θ₂φ₁φ⃗)

    φ2π, freq_mult = fullcircle(φ)
    covPP̄  = Γθ₁θ₂φ₁φ⃗(θ₁, θ₂, φ2π[1], φ2π)
    covPP  = Cθ₁θ₂φ₁φ⃗(θ₁, θ₂, φ2π[1], φ2π)
    covPP̄′ = periodize(covPP̄, freq_mult)
    covPP′ = periodize(covPP, freq_mult)

    return covPP̄′, covPP′
end

# overload this
function γθ₁θ₂φ⃗(θ₁::Real, θ₂::Real, φ::AbstractVector, Γθ₁θ₂φ₁φ⃗)

    φ2π, freq_mult = fullcircle(φ)
    covPP̄  = Γθ₁θ₂φ₁φ⃗(θ₁, θ₂, φ2π[1], φ2π)
    covPP̄′ = periodize(covPP̄, freq_mult)

    return covPP̄′
end

# behavior comes directly from γθ₁θ₂φ⃗_ξθ₁θ₂φ⃗
function γθ₁θ₂ℓ⃗_ξθ₁θ₂ℓ⃗(
    θ₁::Real, θ₂::Real, φ::AbstractVector, Γθ₁θ₂φ₁φ⃗, Cθ₁θ₂φ₁φ⃗,
    planFFT = FFTW.plan_fft(Vector{ComplexF64}(undef,length(φ))),
    γstorage = Vector{ComplexF64}(undef,length(φ)),
    ξstorage = Vector{ComplexF64}(undef,length(φ)),
    )

    γ₁₂φ⃗, ξ₁₂φ⃗ =  γθ₁θ₂φ⃗_ξθ₁θ₂φ⃗(θ₁, θ₂, φ, Γθ₁θ₂φ₁φ⃗, Cθ₁θ₂φ₁φ⃗)

    mul!(γstorage, planFFT, γ₁₂φ⃗)
    mul!(ξstorage, planFFT, ξ₁₂φ⃗)

    return γstorage, ξstorage
end

# behavior comes directly from γθ₁θ₂φ⃗_ξθ₁θ₂φ⃗
function γθ₁θ₂ℓ⃗(
    θ₁::Real, θ₂::Real, φ::AbstractVector, Γθ₁θ₂φ₁φ⃗,
    planFFT  = FFTW.plan_fft(Vector{ComplexF64}(undef,length(φ))),
    γstorage = Vector{ComplexF64}(undef,length(φ)),
    )

    γ₁₂φ⃗ = γθ₁θ₂φ⃗(θ₁, θ₂, φ, Γθ₁θ₂φ₁φ⃗)
    mul!(γstorage, planFFT, γ₁₂φ⃗)

    return γstorage
end

# structs Γθ₁θ₂φ₁φ⃗_CMBpol and Cθ₁θ₂φ₁φ⃗_CMBpol
# =====================================================

"""
IAU uses rotation around outward normal to the sphere => Q + iU is spin (+2)
"""
struct Γθ₁θ₂φ₁φ⃗_CMBpol
    IAU::Bool
    premult_spln::Spline1D
end

struct Cθ₁θ₂φ₁φ⃗_CMBpol
    IAU::Bool
    premult_spln::Spline1D
end

# Constructor for both Γ and C
function ΓCθ₁θ₂φ₁φ⃗_CMBpol(
        ℓ, eeℓ, bbℓ;
        IAU = false,
        ngrid::Int = 100_000,
        βgrid = βsingrid(ngrid, π),
    )
    @assert ℓ[1] == 0
    @assert ℓ[2] == 1
    @assert IAU == false # TODO remove this an impliment the spin(+2) version
    nℓ = @. (2ℓ+1)/(4π)
    ## ↓ starts at 2 since the Jacobi expansion goes like J^(a,b)_{ℓ-2}
    j2⁺2ℓ = (@. (eeℓ + bbℓ) * nℓ)[2:end]
    j2⁻2ℓ = (@. (eeℓ - bbℓ) * nℓ)[2:end]
    ## ↓  TODO: check the a,b swap
    f2⁺2  = ((a,b,jℓ)=(0,4,j2⁺2ℓ);  Fun(Jacobi(b,a),jℓ))
    f2⁻2  = ((a,b,jℓ)=(4,0,j2⁻2ℓ);  Fun(Jacobi(b,a),jℓ))
    # !! pre-canceled out cos β½ and sin β½ in the denom
    covPP̄ = x-> f2⁺2(cos(x))
    covPP = x-> f2⁻2(cos(x))
    β2covPP̄ = Spline1D(βgrid, covPP̄.(βgrid), k=3)
    β2covPP = Spline1D(βgrid, covPP.(βgrid), k=3)
    Γθ₁θ₂φ₁φ⃗_CMBpol(IAU, β2covPP̄), Cθ₁θ₂φ₁φ⃗_CMBpol(IAU, β2covPP)
end

# Hook into method γθ₁θ₂φ⃗_ξθ₁θ₂φ⃗
function γθ₁θ₂φ⃗_ξθ₁θ₂φ⃗(
    θ₁::Real, θ₂::Real, φ::AbstractVector,
    Γθ₁θ₂φ₁φ⃗::Γθ₁θ₂φ₁φ⃗_CMBpol,
    Cθ₁θ₂φ₁φ⃗::Cθ₁θ₂φ₁φ⃗_CMBpol,
    )

    φ2π, freq_mult = fullcircle(φ)
    β      = geoβ.(θ₁, θ₂, φ2π[1], φ2π)
    covPP̄  = Γθ₁θ₂φ₁φ⃗.premult_spln(β) .* multPP̄.(θ₁, θ₂, φ2π[1], φ2π)
    covPP  = Cθ₁θ₂φ₁φ⃗.premult_spln(β) .* multPP.(θ₁, θ₂, φ2π[1], φ2π)
    covPP̄′ = periodize(covPP̄, freq_mult)
    covPP′ = periodize(covPP, freq_mult)

    return covPP̄′, covPP′
end

# for isotropic spin 0
# =======

struct Γθ₁θ₂φ₁φ⃗_Iso
    spln::Spline1D
end

# constructor
function Γθ₁θ₂φ₁φ⃗_Iso(
        ℓ, ttℓ;
        ngrid::Int = 100_000,
        βgrid = βsingrid(ngrid,π),
    )
    @assert ℓ[1] == 0
    @assert ℓ[2] == 1
    nℓ = @. (2ℓ+1)/(4π)
    j0⁺0tℓ = @. ttℓ * nℓ
    f0⁺0t = ((a,b,jℓ)=(0,0,j0⁺0tℓ); Fun(Jacobi(b,a),jℓ))
    covtt = x-> f0⁺0t(cos(x))
    β2covtt = Spline1D(βgrid, covtt.(βgrid), k=3)
    Γθ₁θ₂φ₁φ⃗_Iso(β2covtt)
end

# Hook into method γθ₁θ₂φ⃗
function γθ₁θ₂φ⃗(
    θ₁::Real, θ₂::Real, φ::AbstractVector,
    Γθ₁θ₂φ₁φ⃗::Γθ₁θ₂φ₁φ⃗_Iso,
    )

    φ2π, freq_mult = fullcircle(φ)
    β      = geoβ.(θ₁, θ₂, φ2π[1], φ2π)
    covPP̄  = Γθ₁θ₂φ₁φ⃗.spln(β)
    covPP̄′ = periodize(covPP̄, freq_mult)

    return complex(covPP̄′)
end

# Multipliers needed to convert the isotropic parts to full polarization cov
# =====================================================

function sincosΔθpθΔφ(θ₁, θ₂, φ₁, φ₂)
    𝓅θ½ = (θ₁ + θ₂)/2
    Δθ½ = (θ₁ - θ₂)/2
    Δφ½ = (φ₁ - φ₂)/2
    s𝓅θ½, c𝓅θ½ = sincos(𝓅θ½)
    sΔθ½, cΔθ½ = sincos(Δθ½)
    sΔφ½, cΔφ½ = sincos(Δφ½)
    return sΔθ½, sΔφ½, cΔθ½, cΔφ½, s𝓅θ½, c𝓅θ½
end

function multPP̄(θ₁, θ₂, φ₁, φ₂)
    sΔθ½, sΔφ½, cΔθ½, cΔφ½, s𝓅θ½, c𝓅θ½ = sincosΔθpθΔφ(θ₁, θ₂, φ₁, φ₂)
    return complex(sΔφ½ * c𝓅θ½,   cΔφ½ * cΔθ½)^4
end

function multPP(θ₁, θ₂, φ₁, φ₂)
    sΔθ½, sΔφ½, cΔθ½, cΔφ½, s𝓅θ½, c𝓅θ½ = sincosΔθpθΔφ(θ₁, θ₂, φ₁, φ₂)
    return complex(sΔφ½ * s𝓅θ½, - cΔφ½ * sΔθ½)^4
end

## multII(θ₁, θ₂, φ₁, φ₂) = 1

Q1Q2(covPP̄, covPP) = ( real(covPP̄) + real(covPP) ) / 2

U1U2(covPP̄, covPP) = ( real(covPP̄) - real(covPP) ) / 2

Q1U2(covPP̄, covPP) = ( imag(covPP̄) + imag(covPP) ) / 2

U1Q2(covPP̄, covPP) = (- imag(covPP̄) + imag(covPP) ) / 2

#
# ==================================================

"""
`cheb2spherecov(β, θs)` compute auto-covariance function at
angular distances `θs` from chebyshev expansion.
"""
function cheb2spherecov(β, θs)
    cθ = zero(θs)
    n  = length(β)
    for i ∈ eachindex(cθ)
        for (k,kp1) ∈ zip(0:n-1, 1:n)
            @inbounds cθ[i] += cos(k * θs[i]) * β[kp1]
        end
    end
    return cθ
end

"""
`spec2spherecov(cl, θs)` compute auto-covariance function at
angular distances `θs` from spherical spectral density `cl`.
Note: `cl` is assumed to be a vector whose values correspond to
multipoles `l = [0,1,2 ... ,lmax]`
"""
function spec2spherecov(cl, θs)
    lmax = length(cl)-1
    l    = 0:lmax
    j00l = @. cl * (2l + 1) / (4π)
    β    = jac2cheb(j00l, 0, 0)
    cheb2spherecov(β, θs)
end

# These are slated for removal ...
# ==================================================

struct βcovSpin2
    covPP̄_premult_spln::Spline1D
    covPP_premult_spln::Spline1D
end

struct βcovSpin0
    covII_premult_spln::Spline1D
end

function βcovSpin2(
        ℓ, eeℓ, bbℓ;
        ngrid::Int = 100_000,
        βgrid = βsingrid(ngrid, π),
    )
    @assert ℓ[1] == 0
    @assert ℓ[2] == 1
    nℓ = @. (2ℓ+1)/(4π)
    ## ↓ starts at 2 since the Jacobi expansion goes like J^(a,b)_{ℓ-2}
    j2⁺2ℓ = (@. (eeℓ + bbℓ) * nℓ)[2:end]
    j2⁻2ℓ = (@. (eeℓ - bbℓ) * nℓ)[2:end]
    ## ↓  TODO: check the a,b swap
    f2⁺2  = ((a,b,jℓ)=(0,4,j2⁺2ℓ);  Fun(Jacobi(b,a),jℓ))
    f2⁻2  = ((a,b,jℓ)=(4,0,j2⁻2ℓ);  Fun(Jacobi(b,a),jℓ))
    # !! pre-canceled out cos β½ and sin β½ in the denom
    covPP̄ = x-> f2⁺2(cos(x))
    covPP = x-> f2⁻2(cos(x))
    β2covPP̄ = Spline1D(βgrid, covPP̄.(βgrid), k=3)
    β2covPP = Spline1D(βgrid, covPP.(βgrid), k=3)
    βcovSpin2(β2covPP̄, β2covPP)
end

function βcovSpin0(
        ℓ, ttℓ;
        ngrid::Int = 100_000,
        βgrid = βsingrid(ngrid,π),
    )
    @assert ℓ[1] == 0
    @assert ℓ[2] == 1
    nℓ = @. (2ℓ+1)/(4π)
    j0⁺0tℓ = @. ttℓ * nℓ
    f0⁺0t = ((a,b,jℓ)=(0,0,j0⁺0tℓ); Fun(Jacobi(b,a),jℓ))
    covtt = x-> f0⁺0t(cos(x))
    β2covtt = Spline1D(βgrid, covtt.(βgrid), k=3)
    βcovSpin0(β2covtt)
end

function (covP::βcovSpin2)(β::Matrix)
    rtnPP̄ = similar(β)
    rtnPP = similar(β)
    for (col, cβ) ∈ enumerate(eachcol(β))
        rtnPP̄[:,col] = covP.covPP̄_premult_spln(cβ)
        rtnPP[:,col] = covP.covPP_premult_spln(cβ)
    end
    return complex(rtnPP̄), complex(rtnPP)
end
function (covP::βcovSpin2)(β::Union{Vector, Number})
    rtnPP̄ = covP.covPP̄_premult_spln(β)
    rtnPP = covP.covPP_premult_spln(β)
    return complex(rtnPP̄), complex(rtnPP)
end

# Note: the reason we have different methods for Matrix vrs
# Union{Vector, Number} is that Spline1D's are optimized for
# Union{Vector, Number} so in general it is better to broadcast
# via via whole columns

# Also note: the only reason we make the return argument complex
# is that is the eltype the planned FFT will expect.

function (covP::βcovSpin0)(β::Matrix)
    rtn = similar(β)
    for (col, cβ) ∈ enumerate(eachcol(β))
        rtn[:,col] = covP.covII_premult_spln(cβ)
    end
    return complex(rtn)
end
function (covP::βcovSpin0)(β::Union{Vector, Number})
    return complex(covP.covII_premult_spln(β))
end

end


text:
# for stationary problems
function add_stationarity_constraints!(model::Model, MP::MarkovProcess, v::Int, ::Partition, domain::AbstractSemialgebraicSet, w, rhs)
    @constraint(model, inf_generator(MP, w[v]) >= rhs, domain = domain)
end

function add_stationarity_constraints!(model::Model, MP::MarkovProcess, v::Int, ::Partition, domain::Vector{<:AbstractSemialgebraicSet}, w, rhs)
    cons = []
    for X in domain
        push!(cons, @constraint(model, inf_generator(MP, w[v]) >= rhs, domain = X))
    end
    return cons
end

function add_stationarity_constraints!(model::Model, MP::MarkovProcess, v::Int, P::Partition, ::Singleton, w, rhs::APL)
    @constraint(model, inf_generator(MP, w, v, P) >= rhs(MP.x => props(P.graph, v)[:cell].x))
end

function add_stationarity_constraints!(model::Model, MP::MarkovProcess, v::Int, P::Partition, ::Singleton, w, rhs)
	@constraint(model, inf_generator(MP, w, v, P) >= rhs)
end

function add_coupling_constraints!(model::Model, MP::MarkovProcess, e::Edge, P::Partition, w)
    if props(P.graph, e.dst)[:cell] isa Singleton
         @constraint(model, w[e.dst] - w[e.src](MP.x => props(P.graph, e.dst)[:cell].x) == 0)
    else
        Xs = props(P.graph, e)[:interface]
        if Xs isa AbstractVector
            for X in Xs # TODO FIX THIS
                @constraint(model, w[e.dst] - w[e.src] == 0, domain = X)
            end
        else
            @constraint(model, w[e.dst] - w[e.src] == 0, domain = Xs)
        end
    end
end

# for dynamic problems
function add_dynamics_constraints!(model::Model, MP::MarkovProcess, v::Int, P::Partition,
                                   space_domain::Singleton,
                                   time_domain::AbstractSemialgebraicSet,
                                   Δt::AbstractVector, w, rhs; ρ::Real = 0)
    nT = length(Δt)
    t = MP.iv
    model.obj_dict[Symbol("dynamics_$v")] = @constraint(model, [k in 1:nT], extended_inf_generator(MP, w, (k, v), P; scale = Δt[k]) - ρ*w[k, v] >= Δt[k]*rhs, domain = time_domain, base_name = "dynamics_$v")
    model.obj_dict[Symbol("temporal_coupling_$v")] = @constraint(model, [k in 2:nT], w[k, v](t => 0) - w[k-1, v](t => 1) >= 0, base_name = "temporal_coupling_$v")
end

function add_dynamics_constraints!(model::Model, MP::MarkovProcess, v::Int, P::Partition,
                                   space_domain::Singleton,
                                   time_domain::AbstractSemialgebraicSet,
                                   Δt::AbstractVector, w, rhs::APL; ρ::Real = 0)
    nT = length(Δt)
    t = MP.iv
    model.obj_dict[Symbol("dynamics_$v")] = @constraint(model, [k in 1:nT], extended_inf_generator(MP, w, (k, v), P; scale = Δt[k])  - ρ*w[k, v] >= Δt[k]*subs(rhs, MP.x => space_domain.x), domain = time_domain, base_name = "dynamics_$v")
    model.obj_dict[Symbol("temporal_coupling_$v")] = @constraint(model, [k in 2:nT], w[k, v](t => 0) - w[k-1,v](t => 1) >= 0, base_name = "temporal_coupling_$v")
end

function add_dynamics_constraints!(model::Model, MP::MarkovProcess, v::Int, P::Partition,
                                   space_domain::AbstractSemialgebraicSet,
                                   time_domain::AbstractSemialgebraicSet,
                                   Δt::AbstractVector, w, rhs; ρ::Real = 0)
    nT = length(Δt)
    XT = intersect(space_domain,time_domain)
    t = MP.iv
    model.obj_dict[Symbol("dynamics_$v")] = @constraint(model, [k in 1:nT], extended_inf_generator(MP, w[k, v]; scale = Δt[k]) - ρ*w[k, v] >= Δt[k]*rhs, domain = XT, base_name = "dynamics_$v")
    model.obj_dict[Symbol("temporal_coupling_$v")] = @constraint(model, [k in 2:nT], subs(w[k, v], t => 0) - subs(w[k-1, v], t => 1) >= 0, domain = space_domain, base_name = "temporal_coupling_$v")
end

function add_dynamics_constraints!(model::Model, MP::MarkovProcess, v::Int, P::Partition,
                                   space_domain::Vector{<:AbstractSemialgebraicSet},
                                   time_domain::AbstractSemialgebraicSet,
                                   Δt::AbstractVector, w, rhs; ρ::Real = 0)
    nT = length(Δt)
    t = MP.iv
    model.obj_dict[Symbol("dynamics_$v")] = @constraint(model, [k in 1:nT, X in space_domain], extended_inf_generator(MP, w[k, v]; scale = Δt[k]) - ρ*w[k, v] >= Δt[k]*rhs, domain = XT = intersect(X,time_domain), base_name = "dynamics_$v")
    model.obj_dict[Symbol("temporal_coupling_$v")] = @constraint(model, [k in 2:nT, X in space_domain], subs(w[k, v], t => 0) - subs(w[k-1,v], t => 1) >= 0, domain = X, base_name = "temporal_coupling_$v")
end

function add_transversality_constraints!(model::Model, MP::MarkovProcess, space_domain::Singleton, w, rhs::Real, v::Int)
    model.obj_dict[Symbol("transversality_$v")] = @constraint(model, subs(w, MP.iv => 1) <= rhs, base_name = "transversality_$v")
end

function add_transversality_constraints!(model::Model, MP::MarkovProcess, space_domain::Singleton, w, rhs::APL, v::Int)
    model.obj_dict[Symbol("transversality_$v")] = @constraint(model, subs(w, MP.iv => 1) <= subs(rhs, MP.x => space_domain.x), base_name = "transversality_$v")
end

function add_transversality_constraints!(model::Model, MP::MarkovProcess, space_domain::AbstractSemialgebraicSet, w, rhs, v::Int)
    model.obj_dict[Symbol("transversality_$v")] = @constraint(model, subs(w, MP.iv => 1) <= rhs, domain = space_domain, base_name = "transversality_$v")
end

function add_transversality_constraints!(model::Model, MP::MarkovProcess, space_domain::Singleton, w, rhs, v::Int)
    model.obj_dict[Symbol("transversality_$v")] = @constraint(model, subs(w, MP.iv => 1) <= rhs, base_name = "transversality_$v")
end

function add_transversality_constraints!(model::Model, MP::MarkovProcess, space_domain::Vector{<:AbstractSemialgebraicSet}, w, rhs, v::Int)
    model.obj_dict[Symbol("transversality_$v")] = @constraint(model, [X in space_domain], subs(w, MP.iv => 1) <= rhs, domain = X, base_name = "transversality_$v")
end

function add_coupling_constraints!(model::Model, MP::MarkovProcess, e::Edge, P::Partition,
                                   time_domain::AbstractSemialgebraicSet, Δt::AbstractVector, w)
    nT = length(Δt)
    if props(P.graph, e.dst)[:cell] isa Singleton
        @constraint(model, [k in 1:nT], w[k, e.dst] - subs(w[k, e.src], MP.x => props(P.graph, e.dst)[:cell].x) == 0, domain = time_domain)
    else
        Xs = props(P.graph, e)[:interface]
        for X in Xs
            XT = intersect(X,time_domain)
            @constraint(model, [k in 1:nT], w[k, e.dst] - w[k, e.src] == 0, domain = XT)
        end
    end
end

function add_path_chance_constraint!(model::Model, MP::MarkovProcess, PC::ChanceConstraint, time_domain, Δt::AbstractVector, w_pc, s_pc)
	∂X = ∂(PC.X)
	t = MP.iv
	nT = length(Δt)
	@constraint(model, [i in 1:nT], extended_inf_generator(MP, w_pc[i]; scale = Δt[i]) >= 0, domain = intersect(PC.X, time_domain))
	@constraint(model, [i in 2:nT], subs(w_pc[i], t => 0) - subs(w_pc[i-1], t => 1) >= 0, domain = PC.X)
	@constraint(model, [i in 1:nT, k in 1:length(∂X)], w_pc[i] >= 0, domain = ∂X[k])
	@constraint(model, -s_pc - subs(w_pc[nT], t => 1) >= 0, domain = PC.X)
	@constraint(model, s_pc >= 0)
end

function add_boundary_constraints!(model::Model, MP::MarkovProcess, nT::Int, v::Int, P::Partition, space_domain::AbstractSemialgebraicSet,
								   time_domain::AbstractSemialgebraicSet, w, X::AbstractSemialgebraicSet, ∂X::AbstractVector{<:AbstractSemialgebraicSet})
	@constraint(model, [i in 1:nT, k in 1:length(∂X)], w[i, v] >= 0, domain = intersect(∂X[k], space_domain, T))
	@constraint(model, - 1 - subs(w[nT, v], MP.iv => 1) >= 0, domain = intersect(space_domain, X))
end

function add_terminal_chance_constraint!(model::Model, MP::MarkovProcess, TC::ChanceConstraint, v::Int, P::Partition, space_domain::AbstractSemialgebraicSet, w, rhs)
	@constraint(model, subs(w[nT, v], MP.iv => 1) <= rhs, domain = intersect(TC.X,space_domain))
end


text:
module Cryptowatch

    using DataFrames
    using Dates

    import HTTP: get
    import JSON2: read

    export assets,
           pairs,
           exchanges,
           markets,
           prices,
           summary,
           orderbook,
           ohlc

    include("public.jl")

    const URL = "https://api.cryptowat.ch"

end


text:
module SymmetricSparseMatrices

import Base: size, getindex, show, *
import SparseArrays: nnz, sparse

export SymmetricSparseMatrix
export nnz, getindex, show, size, symsparse, *

using SparseArrays

"""
    SymmetricSparseMatrix{T<:Real, I<:Integer}
Symmetric Sparse matrix implementation.
Extension from two-dimensional arrays (or array-like types) with
elements of type `T` and `I`. Alias for [`AbstractArray{T,2}`](@ref).
"""
    struct SymmetricSparseMatrix{T<:Real, I<:Integer} <: AbstractMatrix{T}
        m      ::Int                # Number of rows/columns
        rowptr ::Vector{I}          # Row i is in rowptr[i]:(rowptr[i+1]-1)
        colval ::Vector{I}          # Col indices of stored values
        nzval  ::Vector{T}          # Stored values, typically nonzeros

        function SymmetricSparseMatrix{T,I}(m::Integer, rowptr::Vector{I}, colval::Vector{I}, nzval::Vector{T}) where {T<:Real,I<:Integer}
            rowptr[end]-1 == length(nzval) == length(colval) || throw(ArgumentError(string("Length of colval and nzval arrays does not match with rowptr[end]-1 ", "(",rowptr[end]-1,")")))
            rowptr[end]-1 == length(nzval) == length(colval) || throw(ArgumentError(string("Length of colval and nzval arrays does not match with rowptr[end]-1 ", "(",rowptr[end]-1,")")))
            new(Int(m), rowptr, colval, nzval)
        end
    end # SymmetricSparseMatrix

"""
    SymmetricSparseMatrix(m::Integer, n::Integer, rowptr::Vector, colval::Vector, nzval::Vector)
SymetricSparseMatrix outter constructor method with full signature.
"""
    function SymmetricSparseMatrix(m::Integer, rowptr::Vector, colval::Vector, nzval::Vector)
        Tv = eltype(nzval)
        Ti = promote_type(eltype(rowptr), eltype(colval))
        SymmetricSparseMatrix{Tv,Ti}(m, rowptr, colval, nzval)
    end

"""
    SymmetricSparseMatrix()
SymetricSparseMatrix outter constructor method to create empty .
"""
    SymmetricSparseMatrix() = SymmetricSparseMatrix(0,[1],Vector{Int}(),Vector{Int}())

"""
    convertCOOrowstoCSR(m, I)
Convert ordered COO rows (I) to CSR rowptr given its number of rows.
"""
    function convertCOOrowstoCSR(m, nnz, I)
        rowptr=zeros(eltype(I), m+1)
        rowptr[1]=1

        if length(I) > 0
            i = 1
            j = 1
            irw = I[j]

            while i<=m
                while i<irw; rowptr[i+1] = rowptr[i]; i+=1; end

                j+=1
                j>=nnz && break
                if I[j] != irw
                    rowptr[i+1] = j
                    irw = I[j]
                    i+=1
                end
            end
            while i <= m; rowptr[i+1] = j; i+=1; end
        end
        rowptr[m+1]=nnz+1
        return rowptr
    end

"""
    symsparse(S::SparseMatrixCSC{T,Ti}) where {T, Ti}
Convert SparseMatrixCSC{T,Ti} to SymmetricSparseMatrix.
"""
    function symsparse(S::SparseMatrixCSC{T,Ti}) where {T, Ti}
        S.m == S.n || throw(ArgumentError(string("Only square input matrix are supported")))
        m  = S.m                     # Number of rows
        nz = nnz(S)                  # Number of columns

        for i in 1:m setindex!(S, getindex(S, i, i), i, i) end

        # Sort by rows
        I,J,K = findnz(S)
        permutation = sortperm(I)

        # Return SymmetricSparseMatrix CSR
        return SymmetricSparseMatrix(m, convertCOOrowstoCSR(m, nz, I[permutation]), J[permutation], K[permutation])
    end

"""
    symsparse(A::AbstractMatrix{Tv}) where {Tv}
Convert AbstractMatrix{Tv} to SymmetricSparseMatrix.
"""
    symsparse(A::AbstractMatrix{Tv}) where {Tv} = symsparse(convert(SparseMatrixCSC{Tv,Int}, A))

"""
    show(io::IO, A::SymmetricSparseMatrix)
SymetriSparseMatrix show method implementation.
"""
    show(io::IO, A::SymmetricSparseMatrix) = dump(A::SymmetricSparseMatrix)
#show(io::IO, m::MIME"text/plain", A::SymmetricSparseMatrix) = show(io::IO, A::SymmetricSparseMatrix)

"""
    size(A::SymmetricSparseMatrix)
Returns the size (#rows, #columns) of the SymmetricSparseMatrix.
"""
    size(A::SymmetricSparseMatrix) = (A.m, A.m)

"""
    nnz(A::SymmetricSparseMatrix)
Returns the number of stored (filled) elements in a SymmetricSparseMatrix.
"""
    nnz(A::SymmetricSparseMatrix) = A.rowptr[end]-1

"""
    getindex(A::SymmetricSparseMatrix{T}, x, y) where {T<:Real}
Return the (x, y) value of the SymmetricSparseMatrix.
"""
    function getindex(A::SymmetricSparseMatrix{T,I}, x::Integer, y::Integer) where {T, I}
        if !(x in 1:A.m) || !(y in 1:A.m) throw(BoundsError(string(size(A), " SymmetricSparseMatrix at index ", "[",x,",", y,"]"))) end
        x_= min(x,y)
        y_= max(x,y)
        lbound=A.rowptr[x_]
        ubound=A.rowptr[x_+1]-1
        if ubound - lbound  < 0
            return zero(T)
        else
            index = findfirst(isequal(y_), A.colval[lbound:ubound])
            return index == nothing ? zero(T) : A.nzval[index+lbound-1]
        end
    end # getindex

"""
    *(A::SymmetricSparseMatrix, v::Vector{T}) where {T<:Real}
Return the A*v product where `A` is SymmetricSparseMatrix{T,I} and `v` a Vector{T}.
"""
    function *(A::SymmetricSparseMatrix, v::Vector{T}) where {T<:Real}
        b = length(v) == A.m ? zero(v) : throw(ArgumentError(string("v::Vector{T} size must match with SymmetricSparseMatrix number of rows", "(",A.m,")")))
        for x in 1:A.m
            lbound=A.rowptr[x]
            ubound=A.rowptr[x+1]-1
            for i in lbound:ubound
                y = A.colval[i]
                if x==y
                    b[x] += A.nzval[i] * v[x]
                elseif y>x
                    b[x] += A.nzval[i] * v[x]
                    b[y] += A.nzval[i] * v[y]
                end
            end
        end
        return b
    end # *

end # module


text:
using Documenter, GridapDistributed

pages = [
  "Home" => "index.md",
  "Getting Started" => "getting-started.md",
  "GridapDistributed" => "GridapDistributed.md",
  "Algebra" => "Algebra.md",
  "CellData" => "CellData.md",
  "DivConformingFESpaces" => "DivConformingFESpaces.md",
  "FESpaces" => "FESpaces.md",
  "Geometry" => "Geometry.md",
  "MultiField" => "MultiField.md",
  "Visualization" => "Visualization.md",
 ]

makedocs(;
    modules=[GridapDistributed],
    format=Documenter.HTML(),
    pages=pages,
    repo="https://github.com/gridap/GridapDistributed.jl/blob/{commit}{path}#L{line}",
    sitename="GridapDistributed.jl",
    authors="S. Badia <santiago.badia@monash.edu>, A. F. Martin <alberto.martin@monash.edu>, F. Verdugo <fverdugo@cimne.upc.edu>",
)

deploydocs(;
    repo="github.com/gridap/GridapDistributed.jl",
)


text:
# Replacement for deprecated functions:

function linearspace(a,b,n)
if VERSION >= v"0.7"
  return range(a,stop=b,length=n)
else
  return linspace(a,b,n)
end
return
end

function logarithmspace(a,b,n)
if VERSION >= v"0.7"
  return 10.0 .^range(a,stop=b,length=n)
else
  return logspace(a,b,n)
end
return
end


text:
#=
    based on original work by John Myles White and Toivo Henningsson
    (see end of file for source code refs)
=#

"""
    macros for doing delegation

    import Base: length, last, (abs), (+)

    Given these types

       type MyInt          type MyInts                   type MyNums{T}
          i::Int              elems::Vector{Int}           elems::T
       end                 end                           end

    These macro calls

       @delegateTyped MyInt.i      [ abs, ]
       @delegateTyped2 MyInt.i     [ (+), ]

       @delegate MyInts.elems [ length,  last ]
       @delegate MyNums.elems [ length,  last ]

    produces these blocks of expressions

      abs(a::MyInt)           = MyInt( abs( getfield(a, :i) ) )
      (+)(a::MyInt, b::MyInt) = MyInt( (+)( getfield(a, :i), getfield(b, :i) ) )

      last( a::MyInts)  = last( getfield(a, :elems) )
      length(a::MyInts) = length( getfield(a, :elems) )

      last( a::MyNums)  = last( getfield(a, :elems) )
      length(a::MyNums) = length( getfield(a, :elems) )

    and allows

      myFirstInt = MyInt(-1)
      mySecondInt = MyInt(2)

      abs(myFirstInt)          # MyInt(1)
      myFirstInt + mySecondInt # MyInt(1)

      myInts = MyInts([5, 4, 3, 2, 1])
      myNums = MyNums([1.0, 2.0, 3.0])

      length(myInts) # 5
      length(myNums) # 3

      last(myInts)  # 1
      last(myNums)  # 3.0
"""
macro delegate(source, targets)
  typename = esc(source.args[1])
  fieldname = esc(Expr(:quote, source.args[2].args[1]))
  funcnames = targets.args
  n = length(funcnames)
  fdefs = Array(Any, n)
  for i in 1:n
    funcname = esc(funcnames[i])
    fdefs[i] = quote
                 ($funcname)(a::($typename), args...) = ($funcname)(getfield(a,($fieldname)), args...)
               end
    end
  return Expr(:block, fdefs...)
end

# for methods that take two equi-typed source arguments

macro delegate2(sourceExemplar, targets)
  typesname = esc(sourceExemplar.args[1])
  fieldname = esc(Expr(:quote, sourceExemplar.args[2].args[1]))
  funcnames = targets.args
  n = length(funcnames)
  fdefs = Array(Any, n)
  for i in 1:n
    funcname = esc(funcnames[i])
    fdefs[i] = quote
                 ($funcname)(a::($typesname), b::($typesname), args...) =
                   ($funcname)(getfield(a,($fieldname)), getfield(b,($fieldname)), args...)
               end
    end
  return Expr(:block, fdefs...)
end

# for methods that take one typed argument and return an iso-typed result

macro delegateTyped(source, targets)
  typename = esc(source.args[1])
  fieldname = esc(Expr(:quote, source.args[2].args[1]))
  funcnames = targets.args
  n = length(funcnames)
  fdefs = Array(Any, n)
  for i in 1:n
    funcname = esc(funcnames[i])
    fdefs[i] = quote
                 ($funcname)(a::($typename), args...) =
                   ($typename)( ($funcname)(getfield(a,($fieldname)), args...) )
               end
    end
  return Expr(:block, fdefs...)
end

# for methods that take two equi-typed source arguments) and return an iso-typed result

macro delegateTyped2(sourceExemplar, targets)
  typesname = esc(sourceExemplar.args[1])
  fieldname = esc(Expr(:quote, sourceExemplar.args[2].args[1]))
  funcnames = targets.args
  n = length(funcnames)
  fdefs = Array(Any, n)
  for i in 1:n
    funcname = esc(funcnames[i])
    fdefs[i] = quote
                 ($funcname)(a::($typesname), b::($typesname), args...) =
                   ($typesname)( ($funcname)(getfield(a,($fieldname)), getfield(b,($fieldname)), args...) )
               end
    end
  return Expr(:block, fdefs...)
end

#=
    initial implementation
    (description and logic from https://gist.github.com/johnmyleswhite/5225361)

    additional macro text from
      https://github.com/JuliaLang/DataStructures.jl/blob/master/src/delegate.jl

    and from Toivo for delegation with nary ops
    (https://groups.google.com/forum/#!msg/julia-dev/MV7lYRgAcB0/-tS50TreaPoJ)

    julia> type T
               x
           end

    julia> import Base.sin, Base.cos

    julia> for f in (:+, :- )    # delegate binary + and - to T.x
               @eval $f(a::T, b::T) = $f(a.x, b.x)
           end

    julia> for f in (:sin, :cos) # delegate sin and cos
               @eval $f(a::T) = $f(a.x)
           end
=#


text:
using Sobol, StatsFuns

function simulate!(X, mcProcess::HestonProcess, rfCurve::ZeroRate, mcBaseData::MonteCarloConfiguration{type1, type2, type3, SerialMode, type5}, T::numb) where {numb <: Number, type1 <: Integer, type2 <: Integer, type3 <: SobolMode, type5 <: Random.AbstractRNG}
    r = rfCurve.r
    S0 = mcProcess.underlying.S0
    d = dividend(mcProcess)
    Nsim = mcBaseData.Nsim
    Nstep = mcBaseData.Nstep
    σ = mcProcess.σ
    σ_zero = mcProcess.σ_zero
    λ1 = mcProcess.λ
    κ = mcProcess.κ
    ρ = mcProcess.ρ
    θ = mcProcess.θ
    @assert T > 0.0

    ####Simulation
    ## Simulate
    κ_s = κ + λ1
    θ_s = κ * θ / κ_s
    seq = SobolSeq(2 * Nstep)
    skip(seq, Nstep * Nsim)
    dt = T / Nstep
    isDualZero = T * r * σ_zero * θ_s * κ_s * σ * ρ * 0.0 * S0
    view(X, :, 1) .= isDualZero
    for i = 1:Nsim
        v_m = σ_zero^2
        vec = norminvcdf.(next!(seq))
        for j = 1:Nstep
            @views e1 = vec[j]
            @views e2 = e1 * ρ + vec[Nstep+j] * sqrt(1 - ρ * ρ)
            @views X[i, j+1] = X[i, j] + ((r - d) - 0.5 * v_m) * dt + sqrt(v_m) * sqrt(dt) * e1
            v_m += κ_s * (θ_s - v_m) * dt + σ * sqrt(v_m) * sqrt(dt) * e2
            #when v_m = 0.0, the derivative becomes NaN
            v_m = max(v_m, isDualZero)
        end
    end
    ## Conclude
    X .= S0 .* exp.(X)
    return
end

function simulate!(X, mcProcess::HestonProcess, rfCurve::ZeroRateCurve, mcBaseData::MonteCarloConfiguration{type1, type2, type3, SerialMode, type5}, T::numb) where {numb <: Number, type1 <: Integer, type2 <: Integer, type3 <: SobolMode, type5 <: Random.AbstractRNG}
    r = rfCurve.r
    S0 = mcProcess.underlying.S0
    d = dividend(mcProcess)
    Nsim = mcBaseData.Nsim
    Nstep = mcBaseData.Nstep
    σ = mcProcess.σ
    σ_zero = mcProcess.σ_zero
    λ1 = mcProcess.λ
    κ = mcProcess.κ
    ρ = mcProcess.ρ
    θ = mcProcess.θ
    @assert T > 0.0

    ####Simulation
    ## Simulate
    κ_s = κ + λ1
    θ_s = κ * θ / (κ + λ1)
    r_d = r - d
    dt = T / Nstep
    seq = SobolSeq(2 * Nstep)
    skip(seq, Nstep * Nsim)
    zero_drift = r_d(dt * 0.0, dt)
    isDualZero = S0 * T * σ_zero * κ * θ * λ1 * σ * ρ * 0.0 * zero_drift
    #X=Matrix{typeof(isDualZero)}(undef,Nsim,Nstep+1);
    X[:, 1] .= isDualZero
    tmp_ = [r_d((j - 1) * dt, dt) for j = 1:Nstep]
    for i = 1:Nsim
        v_m = σ_zero^2
        vec = norminvcdf.(next!(seq))

        for j = 1:Nstep
            @views e1 = vec[j]
            @views e2 = e1 * ρ + vec[Nstep+j] * sqrt(1 - ρ * ρ)
            @views X[i, j+1] = X[i, j] + (tmp_[j] - 0.5 * v_m) * dt + sqrt(v_m) * sqrt(dt) * e1
            v_m += κ_s * (θ_s - v_m) * dt + σ * sqrt(v_m) * sqrt(dt) * e2
            #when v_m = 0.0, the derivative becomes NaN
            v_m = max(v_m, isDualZero)
        end
    end
    ## Conclude
    X .= S0 .* exp.(X)
    return
end


text:
using Parameters, UUIDs

abstract type Component end

@with_kw struct InputComponent{T} <: Component
    hash::UUID = uuid4()
    component::String
    default::T
    parameters::Union{Dict{String, Any}, Nothing}
end

@with_kw struct StaticComponent <: Component
    hash::UUID = uuid4()
    component::String
    parameters::Union{Dict{String, Any}, Nothing}
    children::Union{Vector{Component}, Nothing}
end

@with_kw struct Bind
    hash::UUID = uuid4()
    variable::String
    component::InputComponent
end

@with_kw struct BindSet
    hash::UUID = uuid4()
    binds::Array{Bind}
end

@with_kw struct OutputCallback
    hash::UUID = uuid4()
    bind_set::BindSet
    fn::Function
end

@with_kw struct UpdateScope
    hash::UUID = uuid4()
    binds::BindSet
end

@with_kw struct OutputComponent{T} <: Component
    hash::UUID = uuid4()
    component::String = "Text"
    callback::OutputCallback
    parameters::Union{Dict{String, Any}, Nothing}
end

macro bind_set(binds)
    if typeof(binds) == Symbol
        local bind_component = eval(binds)
        local component_type = typeof(bind_component)

        if !(component_type <: InputComponent)
            throw(ArgumentError(bind_variable * " is expected to be `InputComponent`, but " * string(component_type) * " given"))
        end

        local bind_set = [Bind(
            variable = string(binds),
            component = bind_component
        )]

        return BindSet(binds = bind_set)
    end

    if binds.head == :...
        local bind_set = eval(binds.args[1])
        local bind_set_type = typeof(bind_set)

        if !(bind_set_type <: BindSet)
            throw(ArgumentError("`binsSet` expected, but " * string(returned_bind) * " given"))
        end

        return bind_set
    end

    if binds.head !== :tuple
        throw(ArgumentError("binds should be tuple, " * string(binds.head) * " provided"))
    end

    local parsed_binds::Vector{Bind} = []

    for i in 1:length(binds.args)
        if typeof(binds.args[i]) == Expr
            if binds.args[i].head !== :...
                throw(ArgumentError("unknown syntax"))
            end

            local bind_set = eval(binds.args[i].args[1])

            if typeof(bind_set) != BindSet
                throw(ArgumentError("`BindSet` expected, but " * string(typeof(bind_set)) * " given"))
            end

            for bind in bind_set.binds
                append!(parsed_binds, [Bind(
                    variable = bind.variable,
                    component = bind.component
                )])
            end
        elseif typeof(binds.args[i]) == Symbol
            local bind_variable = string(binds.args[i])
            local bind_component = eval(binds.args[i])
            local component_type = typeof(bind_component)

            if !(component_type <: InputComponent)
                throw(ArgumentError(bind_variable * " is expected to be `InputComponent`, but " * string(component_type) * " given"))
            end

            append!(parsed_binds, [Bind(
                variable = bind_variable,
                component = bind_component
            )])
        else
            throw(ArgumentError("element of " * string(binds.args[i]) * " should be `Expr` or `Symbol`, " * string(typeof(binds.args[i])) * " provided"))
        end
    end

    return BindSet(binds = parsed_binds)
end

macro output_fn(definition, block)
    if definition.head !== :(::)
        throw(ArgumentError("bind definition must have two parts, concate with `::`"))
    end

    local binds = definition.args[1]
    local output_type = definition.args[2]

    local bind_set = eval(:(@bind_set $binds))

    local call_paras::Vector{Expr} = [Expr(:(::), :__update_hash, UUID)]

    for i in 1:length(bind_set.binds)
        local para_name = Symbol(bind_set.binds[i].variable)
        local para_type = typeof(bind_set.binds[i].component).parameters[1]
        local var_def = Expr(:(::), para_name, para_type)

        append!(call_paras, [var_def])
    end

    local fn_paras = Expr(
        :(::), Expr(
            :call, :__internal_output_callback, Expr(
                :parameters, call_paras...
            )
        ), Symbol(output_type))

        local fn_expr = Expr(:function, fn_paras, block)
        local fn = eval(fn_expr)

        OutputCallback(
            bind_set = bind_set,
            fn = fn
        )
end

function Toggle(
    label::String = "",
    default::Bool = false;
    inline_label::Bool = false,
    on_text::String = "On",
    off_text::String = "Off")
    InputComponent{Bool}(
        component = "Toggle",
        default = default,
        parameters = Dict(
            "label"       => label,
            "inlineLabel" => inline_label,
            "onText"      => on_text,
            "offText"     => off_text
        )
    )
end

@with_kw struct ChoiceOption
    key::String = ""
    option::String = ""
end

function Choice(
    label::String = "",
    default::String = "";
    options::Dict{String, String} = Dict("" => ""))
    @assert haskey(options, default) "default option is not in `options`"

    options_para::Vector{ChoiceOption} = []

    for (k, v) in options
        append!(options_para, [ChoiceOption(k, v)])
    end

    InputComponent{String}(
        component = "ChoiceGroup",
        default = default,
        parameters = Dict(
            "label"              => label,
            "options"            => options_para,
            "defaultSelectedKey" => default
        )
    )
end

function Slider(
    label::String = "",
    default::Float64 = 0.0;
    min::Float64 = 0.0,
    max::Float64 = 10.0,
    step::Float64 = 1.0,
    show_value::Bool = true,
    origin_from_zero::Bool = false,
    vertical::Bool = false)
    InputComponent{Float64}(
        component = "Slider",
        default = default,
        parameters = Dict(
            "min"            => min,
            "max"            => max,
            "step"           => step,
            "showValue"      => show_value,
            "originFromZero" => origin_from_zero,
            "vertical"       => vertical
        )
    )
end

function TextOutput(
    callback::OutputCallback;
    monospace::Bool = false)
    OutputComponent{String}(
        callback = callback,
        component = "Text",
        parameters = Dict(
            "monospace" => monospace
        )
    )
end

function SplitView(
    children::Vector{T} where T <: Component;
    widths::Union{Vector{Int64}, Nothing} = nothing)
    StaticComponent(
        component = "SplitView",
        parameters = Dict(
            "widths" => widths
        ),
        children = children
    )
end

function StackView(
    children::Vector{T} where T <: Component)
    StaticComponent(
        component = "StackView",
        parameters = nothing,
        children = children
    )
end


text:
println("Hello! This is a test on the basic Assoc Array construction.\n")
println("Associative array takes on entries of triplets, and it will parse an array of substrings that is divided by char divider.\n")
println("--Please note that this divider is indicated as the last char in the string.\n")

row = "a,a,a,a,a,a,a,aa,aaa,b,bb,bbb,a,aa,aaa,b,bb,bbb,"
column = "a,aa,aaa,b,bb,bbb,a,a,a,a,a,a,a,aa,aaa,b,bb,bbb,"
vals = "a-a,a-aa,a-aaa,a-b,a-bb,a-bbb,a-a,aa-a,aaa-a,b-a,bb-a,bbb-a,a-a,aa-aa,aaa-aaa,b-b,bb-bb,bbb-bbb,"

# Create assoc array and list triples.
A = Assoc(row,column,vals)

# Show data structure.
println("This is the data structure of the Associative Array Class")
dump(A)

# Display in tabular form.
println("Print allows it to be printed in a tabular form.")
printFull(A)

# Write to a CSV file.
println("When written into CSV form, the data is stored in the tabular form")
WriteCSV(A,joinpath(Base.source_dir(),"A.csv"))


text:
# Automatically generated using Clang.jl

# Skipping MacroDefinition: DEFINE_C_API_STRUCT ( name , storage ) struct name { storage * ptr ; } ; typedef struct name name

struct MlirAffineExpr
    ptr::Ptr{Cvoid}
end

struct MlirAffineMap
    ptr::Ptr{Cvoid}
end

struct MlirDiagnostic
    ptr::Ptr{Cvoid}
end

@cenum MlirDiagnosticSeverity::UInt32 begin
    MlirDiagnosticError = 0
    MlirDiagnosticWarning = 1
    MlirDiagnosticNote = 2
    MlirDiagnosticRemark = 3
end

const MlirDiagnosticHandlerID = UInt64
const MlirDiagnosticHandler = Ptr{Cvoid}

struct MlirContext
    ptr::Ptr{Cvoid}
end

struct MlirDialect
    ptr::Ptr{Cvoid}
end

struct MlirOperation
    ptr::Ptr{Cvoid}
end

struct MlirOpPrintingFlags
    ptr::Ptr{Cvoid}
end

struct MlirBlock
    ptr::Ptr{Cvoid}
end

struct MlirRegion
    ptr::Ptr{Cvoid}
end

struct MlirAttribute
    ptr::Ptr{Cvoid}
end

struct MlirIdentifier
    ptr::Ptr{Cvoid}
end

struct MlirLocation
    ptr::Ptr{Cvoid}
end

struct MlirModule
    ptr::Ptr{Cvoid}
end

struct MlirType
    ptr::Ptr{Cvoid}
end

struct MlirValue
    ptr::Ptr{Cvoid}
end

struct MlirNamedAttribute
    name::MlirIdentifier
    attribute::MlirAttribute
end

struct MlirStringRef
    data::Cstring
    length::Csize_t
end

mutable struct MlirOperationState
    name::MlirStringRef
    location::MlirLocation
    nResults::intptr_t
    results::Ptr{MlirType}
    nOperands::intptr_t
    operands::Ptr{MlirValue}
    nRegions::intptr_t
    regions::Ptr{MlirRegion}
    nSuccessors::intptr_t
    successors::Ptr{MlirBlock}
    nAttributes::intptr_t
    attributes::Ptr{MlirNamedAttribute}
end

struct MlirPass
    ptr::Ptr{Cvoid}
end

struct MlirPassManager
    ptr::Ptr{Cvoid}
end

struct MlirOpPassManager
    ptr::Ptr{Cvoid}
end

# Skipping MacroDefinition: MLIR_CAPI_EXPORTED __attribute__ ( ( visibility ( "default" ) ) )

const MlirStringCallback = Ptr{Cvoid}

struct MlirLogicalResult
    value::Int8
end


text:
using DRMSMIP
using Test

@testset "DRMSMIP.jl" begin
    # Write your tests here.
end


text:
typealias NoData Union(Nothing,Number)

const FILE_TYPE_NPY                 = 1
const FILE_TYPE_AAGRID              = 2
const FILE_TYPE_TXTLIST             = 3
const FILE_TYPE_INCL_PAIRS_AAGRID   = 4
const FILE_TYPE_INCL_PAIRS          = 5

const STR_FILE_TYPES = ["numpy", "aagrid", "text list", "include grid", "include pairs"]

const FILE_HDR_GZIP                 = [0x1f, 0x8b, 0x08]
const FILE_HDR_NPY                  = [0x93, "NUMPY".data]
const FILE_HDR_AAGRID               = ["ncols".data]
const FILE_HDR_INCL_PAIRS_AAGRID    = ["min".data]
const FILE_HDR_INCL_PAIRS           = ["mode".data]

const iologger = getlogger("circuitscape.io")

MSG_RESAMPLE(raster_name, attrib) = "$raster_name raster has different $attrib than habitat raster. Circuitscape will try to crudely resample the raster. We recommend using the \"Export to Circuitscape\" ArcGIS tool to create ASCII grids with compatible cell size and extent."

MSG_NO_RESAMPLE(raster_name, attrib) = "$raster_name raster must have same $attrib as habitat raster"

#######################################
# UTILITY METHODS
#######################################
_check_file_exists(filename::String) = ((stat(filename).ctime == 0) && error("File $filename does not exist"); true)

function _open_auto_uncompress(filename::String)
    f = open(filename, "r")
    try
        hdr = read!(f, Array(Uint8, 3))
        seek(f, 0)
        if hdr == FILE_HDR_GZIP
            close(f)
            f = GZip.open(filename, "r")
        end
    catch
        close(f)
        f = nothing
    end
    f
end

function _guess_file_type(filename::String)
    io = _open_auto_uncompress(filename)
    try
        return _guess_file_type(io)
    catch ex
        close(io)
        rethrow(ex)
    end
end

function _guess_file_type(io::IO)
    hdr = read!(io, Array(Uint8, 10))
    seek(io, 0)

    beginswith(hdr, FILE_HDR_NPY)               && return FILE_TYPE_NPY
    beginswith(hdr, FILE_HDR_AAGRID)            && return FILE_TYPE_AAGRID
    beginswith(hdr, FILE_HDR_INCL_PAIRS_AAGRID) && return FILE_TYPE_INCL_PAIRS_AAGRID
    beginswith(hdr, FILE_HDR_INCL_PAIRS)        && return FILE_TYPE_INCL_PAIRS
    return FILE_TYPE_TXTLIST
end

#######################################
# AAGRID
#######################################
type FileAAGridHeader
    ncols::Int
    nrows::Int
    xllcorner::Float64
    yllcorner::Float64
    cellsize::Float64
    nodata::NoData

    file_type::Int
    file_name::String

    FileAAGridHeader() = new(0, 0, 0.0, 0.0, 0.0, nothing, [], 0, "")

    function FileAAGridHeader(filename::String="")
        ncols, nrows, xllcorner, yllcorner, cellsize, nodata, file_type = _ascii_grid_read_header(filename)
        new(ncols, nrows, xllcorner, yllcorner, cellsize, nodata, file_type, filename)
    end

    function _ascii_grid_read_header(filename::String)
        _check_file_exists(filename)

        io = _open_auto_uncompress(filename)

        try
            ncols = nrows = 0
            xllcorner = yllcorner = cellsize = 0.0
            nodata = nothing

            file_type = _guess_file_type(io)
            if file_type == FILE_TYPE_NPY
                file_base, file_extension = splitext(filename)
                filename = file_base + ".hdr"
                ncols, nrows, xllcorner, yllcorner, cellsize, nodata, _file_type = _ascii_grid_read_header(filename)
            else
                try
                    ncols = int(split(readline(io))[2])
                    nrows = int(split(readline(io))[2])
                    xllcorner = float(split(readline(io))[2])
                    yllcorner = float(split(readline(io))[2])
                    cellsize = float(split(readline(io))[2])
                catch e
                    err(iologger, "Error reading ASCII grid. $(e)")
                    error("Unable to read ASCII grid: \"$(filename)\".")
                end

                next_line = split(readline(io))
                if length(next_line) == 2
                    try
                        nodata = int(next_line[2])
                    catch
                        try
                            nodata = float(next_line[2])
                        end
                    end
                end
            end

            return ncols, nrows, xllcorner, yllcorner, cellsize, nodata, file_type

        finally
            close(io)
        end
    end
end

function show(io::IO, hdr::FileAAGridHeader)
    println("$(STR_FILE_TYPES[hdr.file_type]) $(hdr.nrows)x$(hdr.ncols) ($(hdr.file_name))")
    println("llcorner: ($(hdr.xllcorner),$(hdr.yllcorner)). cellsize: $(hdr.cellsize). nodata: $(hdr.nodata)")
end

type FileAAGrid
    hdr::FileAAGridHeader
    data::Array
    data_type::String

    FileAAGrid() = new(FileAAGridHeader(), [])

    function FileAAGrid(filename::String, data_type="float")
        hdr = FileAAGridHeader(filename)

        if hdr.file_type == FILE_TYPE_NPY
            data = np.load(filename, mmap_mode=nothing)
            data = convert(Array{Float64, 2}, data)
        else
            if hdr.nodata == nothing
                data = np.loadtxt(filename, skiprows=5, dtype=data_type)
            else
                data = np.loadtxt(filename, skiprows=6, dtype=data_type)
                data[find(x->(x==hdr.nodata), data)] = -9999
            end
        end

        (hdr.nrows == 1) && (data = reshape(data, 1, length(data)))
        (hdr.ncols == 1) && (data = reshape(data, length(data), 1))
        new(hdr, data, data_type)
    end
end

function show(io::IO, aag::FileAAGrid)
    print(io, "$(aag.data_type) ")
    show(io, aag.hdr)
    show(io, aag.data)
end

function write(filename::String, aagrid::FileAAGrid; file_type=aagrid.hdr.file_type, compress=false)
    if file_type == FILE_TYPE_NPY
        np.save(filename, aagrid.data)
    else
        io = compress ? GZip.open(filename, "w") : open(filename, "w")
        try
            write(io, aagrid)
        finally
            close(io)
        end
    end
end

function write(io::IO, aagrid::FileAAGrid)
    hdr = aagrid.hdr
    println(io, "ncols         $(hdr.ncols)")
    println(io, "nrows         $(hdr.nrows)")
    println(io, "xllcorner     $(hdr.xllcorner)")
    println(io, "yllcorner     $(hdr.yllcorner)")
    println(io, "cellsize      $(hdr.cellsize)")
    (hdr.nodata != nothing) && println(io, "NODATA_value  $(hdr.nodata)")

    writedlm(io, aagrid.data, ' ')
end

#######################################
# TXT LIST
#######################################

type FileTxtList
#    data::Array
#    data_type::String
#
#    file_name::String
#
#    function FileTxtList(filename::String, habitat_size::FileAAGridHeader, data_type::String="float")
#        typ = (data_type == "float") ? Float64 : (data_type == "int") ? Int : error("invalid data type $data_type")
#        data = readdlm(filename)
#
#        data2 = ceil(habitat_size.nrows - (data[:,3] - habitat_size.yllcorner) / habitat_size.cellsize) - 1
#        data3 = ceil((data[:,2] - habitat_size.xllcorner) / habitat_size.cellsize) - 1
#        data[:,2] = data2
#        data[:,3] = data3
#        data = convert(Array{typ, 2}, data)
#
#        new(data, data_type, filename)
#    end
end


text:
# using Pkg
# Pkg.activate(".")
#
#
#
#
using Pkg
Pkg.activate(".")

Pkg.test()

using AbstractLogic

Pkg.test()

# Pkg.test()

# Main.include("AbstractLogic.jl")
#
# import AbstractLogic: abstractlogic,
#        checkfeasible,
#        help,
#        LogicalCombo,
#        logicalparse,
#        search,
#        showfeasible,
#        discover,
#        dashboard!,
#        dashboard,
#        expand,
#        help,
#        nfeasible,
#        replset,
#        returnreplset,
#        showlogichistory,
#        showcommandhistory,
#        showuserinput,
#        showsetlocation,
#        showcmdlocation,
#        showcommandlist
#
# import AbstractLogic
#
# using MyPackage: myfunction1, myfunction2, myfunction3
#
# # Make the documents
# include("C:\\Users\\francis.smart.ctr\\GitDir\\AbstractLogic.jl\\docs\\make.jl")
#
# using SparseArrays, Serialization

# v7 = logicalparse("a,b,c,d,e,f,g ∈ 1:7 || {{i}} != {{>i}}")
# sv7 = sparsevec(v7[:])
# vector(sv7)
#
# serialize("src/sv/sv7.ser", sv7)
#
# v8 = logicalparse("a,b,c,d,e,f,g,h ∈ 1:8 || {{i}} != {{>i}}")
# sv8 = sparsevec(v8[:])
#
# serialize("src/sv/sv8.ser", sv8)
#
# v9 = logicalparse("a,b,c,d,e,f,g,h,i ∈ 1:9 || {{i}} != {{>i}}")
# sv9 = sparsevec(v9[:])

# function f1()
#     j = 0
#     for i in (1:length(v9[:]))[v9[:]]
#         j += 1
#     end
# end
#
# function f2()
#     j = 0
#     for i in (sv9)
#         j += 1
#     end
# end
#
#
# serialize("src/sv/sv9.ser", sv9)
#

# n = 6
# factorial(n)
# 20
#
# function geti(x, n, k)
#     fset = 1:n
#     y = x
#     z = 0
#     setout = Integer[]
#     for j in 1:k
#         fset, y, z, setout
#         K = factorial(n-j)
#         z = Integer(ceil(y/ K))
#         y = y - (z-1) * K
#         push!(setout, fset[z])
#         fset = fset[[(1:(z-1))..., ((z+1):length(fset))...]]
#     end
#     setout[k]
# end
#
# X = [geti(i,6,j) for i in 1:720, j in 1:n]
#
# unique(X, dims=1)
#
# X = LogicalCombo([:a,:b,:c,:d], 1:4, fill(true, factorial(4)), permutationuniquelookup)
# size(X)
# X[:,:]


text:
using CountDownLatches
using Test
using ThreadPools

@testset "create a latch with a negative count" begin
    @test_throws ArgumentError CountDownLatch(-1)
end

@testset "count down from 1" begin
    latch = CountDownLatch(1)
    count_down(latch)
    await(latch)
    @test get_count(latch) == 0
end

@testset "count down from 1 two times" begin
    latch = CountDownLatch(1)
    count_down(latch)
    count_down(latch)
    await(latch)
    @test get_count(latch) == -1
end

@testset "wait for another thread to count down from 1" begin
    latch = CountDownLatch(1)

    t = Threads.@spawn begin
        count_down(latch)
    end

    await(latch)

    @test get_count(latch) == 0
    @test istaskdone(t)
end

@testset "mutual count down" begin
    latch1 = CountDownLatch(1)
    latch2 = CountDownLatch(1)
    latch3 = CountDownLatch(1)

    t = Threads.@spawn begin
        count_down(latch2)
        await(latch1)
        count_down(latch3)
    end

    await(latch2)
    count_down(latch1)
    await(latch3)

    @test get_count(latch1) == 0
    @test get_count(latch2) == 0
    @test get_count(latch3) == 0
    @test istaskdone(t)
end

@testset "wait for another thread to count down from 2" begin
    latch = CountDownLatch(2)

    t = Threads.@spawn begin
        count_down(latch)
        count_down(latch)
    end

    await(latch)

    @test get_count(latch) == 0
    @test istaskdone(t)
end

@testset "wait for this thread to count down from 2 using ThreadPools" begin
    latch = CountDownLatch(2)

    this_thread_id = Threads.threadid()
    @tspawnat this_thread_id begin
        # I would like to assert that await() was called but this is the best I can come up with
        @test get_count(latch) == 1
        count_down(latch)
    end
    count_down(latch)
    await(latch)

    @test get_count(latch) == 0
end


text:
n = 10 # rand(1:10)
A = matrixdepot("pascal", n)

@test isposdef(A)
@test issymmetric(A)

P = diagm(0 => (-1).^[0:n-1;])
P[:,1] = ones(n, 1)

for j = 2: n-1
    for i = j+ 1:n
        P[i,j] = P[i-1, j] - P[i-1, j-1]
    end
end

@test A ≈ P*P'
println("'pascal' passed test...")


text:
module SCAIterators

export StateIterator, ActionIterator

using SCAs

type StateIterator

    xs::Vector{Float64}
    ys::Vector{Float64}
    bearings::Vector{Float64}
    speedsOwnship::Vector{Float64}
    speedsIntruder::Vector{Float64}

    xCounter::Int64
    yCounter::Int64
    bearingCounter::Int64
    speedOwnshipCounter::Int64
    speedIntruderCounter::Int64

    nStates::Int64

    function StateIterator(
            xs::Vector{Float64},
            ys::Vector{Float64},
            bearings::Vector{Float64},
            speedsOwnship::Vector{Float64},
            speedsIntruder::Vector{Float64})

        nStates =
            length(xs) *
            length(ys) *
            length(bearings) *
            length(speedsOwnship) *
            length(speedsIntruder) +
            1  # +1 for clear of conflict state

        return new(
            xs,  # xs
            ys,  # ys
            bearings,  # bearings
            speedsOwnship,  # speedsOwnship
            speedsIntruder,  # speedsIntruder
            1,  # xCounter
            1,  # yCounter
            1,  # bearingCounter
            1,  # speedOwnshipCounter
            1,  # speedIntruderCounter
            nStates)  # nStates

    end # function StateIterator

end # type StateIterator

Base.start(states::StateIterator) = 1
Base.done(states::StateIterator, stateCounter::Int64) = stateCounter - 1 == states.nStates

function Base.next(states::StateIterator, stateCounter::Int64)

    nextState = State(
        states.xs[states.xCounter],  # x
        states.ys[states.yCounter],  # y
        states.bearings[states.bearingCounter],  # bearing
        states.speedsOwnship[states.speedOwnshipCounter],  # speedOwnship
        states.speedsIntruder[states.speedIntruderCounter],  # speedIntruder
        false)  # clearOfConflict

    updateCounters!(states)

    if stateCounter - 1 == states.nStates - 1
        nextState.clearOfConflict = true
    end # if

    return nextState, stateCounter + 1

end # function Base.next

function updateCounters!(states::StateIterator)

    if (states.xCounter) % length(states.xs) == 0
        states.xCounter = 1
        if (states.yCounter) % length(states.ys) == 0
            states.yCounter = 1
            if (states.bearingCounter) % length(states.bearings) == 0
                states.bearingCounter = 1
                if (states.speedOwnshipCounter) % length(states.speedsOwnship) == 0
                    states.speedOwnshipCounter = 1
                    if (states.speedIntruderCounter) % length(states.speedsIntruder) == 0
                        states.speedIntruderCounter = 1  # we are done
                    else
                        states.speedIntruderCounter += 1
                    end # if: speedIntruder
                else
                    states.speedOwnshipCounter += 1
                end # if: speedOwnship
            else
                states.bearingCounter += 1
            end # if: bearing
        else
            states.yCounter += 1
        end # if: y
    else
        states.xCounter += 1
    end # if: x

end # function updateCounters!

type ActionIterator

    actions::Vector{Symbol}

    ownshipCounter::Int64
    intruderCounter::Int64

    nActions::Int64

    function ActionIterator(actions::Vector{Symbol})

        return new(
            actions,  # actions
            1,  # ownshipCounter
            1,  # intruderCounter
            length(actions)^2)  # nActions

    end # function ActionIterator

end # type ActionIterator

Base.start(actions::ActionIterator) = 1
Base.done(actions::ActionIterator, actionCounter::Int64) = actions.nActions == actionCounter - 1

function Base.next(actions::ActionIterator, actionCounter::Int64)

    nextAction = Action(
        actions.actions[actions.ownshipCounter],  # actionOwnship
        actions.actions[actions.intruderCounter])  # actionIntruder

    updateCounters!(actions)

    return nextAction, actionCounter + 1

end # function Base.next

function updateCounters!(actions::ActionIterator)

    if (actions.ownshipCounter) % length(actions.actions) == 0
        actions.ownshipCounter = 1
        if (actions.intruderCounter) % length(actions.actions) == 0
            actions.intruderCounter = 1
        else
            actions.intruderCounter += 1
        end # if: intruder
    else
        actions.ownshipCounter += 1
    end # if: ownship

end # function updateCounters!

end # module


text:
using TeukEvolution
using Test

include("TestRadial.jl")
include("TestSphere.jl")

import .TestRadial, .TestSphere

#============================================================================
## Testing Radial functions
============================================================================#
# Test 4th order convergence of derivatives

n11, n12 = TestRadial.norms_diff(128,48)
n21, n22 = TestRadial.norms_diff(256,60)

@test n11/n21 >= 15.5
@test n12/n22 >= 15.5

# Test filter is TVD

n1, n2 = TestRadial.norms_filter(128,48)
n3, n4 = TestRadial.norms_filter(180,90)

@test n2 < n1
@test n4 < n3

#============================================================================
## Testing Spherical functions
============================================================================#

println("Testing spherical spin-weighted spherical harmonic inner product")
TestSphere.test_swal_inner_product(48,-2,-3)
TestSphere.test_swal_inner_product(48, 0,-3)
TestSphere.test_swal_inner_product(48, 1,-3)

TestSphere.test_swal_inner_product(48,-1,-2)
TestSphere.test_swal_inner_product(48, 0,-2)
TestSphere.test_swal_inner_product(48, 2,-2)

TestSphere.test_swal_inner_product(48,-2,-1)
TestSphere.test_swal_inner_product(48, 0,-1)
TestSphere.test_swal_inner_product(48, 1,-1)

TestSphere.test_swal_inner_product(48,-1,0)
TestSphere.test_swal_inner_product(48, 0,0)
TestSphere.test_swal_inner_product(48, 2,0)

TestSphere.test_swal_inner_product(48,-2,1)
TestSphere.test_swal_inner_product(48, 0,1)
TestSphere.test_swal_inner_product(48, 1,1)

TestSphere.test_swal_inner_product(48,-1,2)
TestSphere.test_swal_inner_product(48, 0,2)
TestSphere.test_swal_inner_product(48, 2,2)

TestSphere.test_swal_inner_product(48,-2,3)
TestSphere.test_swal_inner_product(48, 0,3)
TestSphere.test_swal_inner_product(48, 1,3)

println("Testing spherical spin-weighted laplacian operator")
TestSphere.test_norm_swal_lap(32, 0,0,2)

TestSphere.test_norm_swal_lap(32, 0,1,2)
TestSphere.test_norm_swal_lap(32, 1,1,2)

TestSphere.test_norm_swal_lap(32,-1,-2,2)
TestSphere.test_norm_swal_lap(32, 0,-2,2)
TestSphere.test_norm_swal_lap(32, 2,-2,2)

TestSphere.test_norm_swal_lap(32,-1,2,2)
TestSphere.test_norm_swal_lap(32, 0,2,2)
TestSphere.test_norm_swal_lap(32, 2,2,2)

TestSphere.test_norm_swal_lap(48,-1,5,23)
TestSphere.test_norm_swal_lap(48, 0,5,23)
TestSphere.test_norm_swal_lap(48, 2,5,23)

println("Testing spherical spin-weighted raising operator")
TestSphere.test_norm_swal_raising(32, 0,0,2)

TestSphere.test_norm_swal_raising(32, 0,1,2)
TestSphere.test_norm_swal_raising(32, 1,1,2)

TestSphere.test_norm_swal_raising(32,-1,-2,2)
TestSphere.test_norm_swal_raising(32, 0,-2,2)
TestSphere.test_norm_swal_raising(32, 2,-2,2)

TestSphere.test_norm_swal_raising(32,-1,2,2)
TestSphere.test_norm_swal_raising(32, 0,2,2)
TestSphere.test_norm_swal_raising(32, 2,2,2)

TestSphere.test_norm_swal_raising(48,-1,5,23)
TestSphere.test_norm_swal_raising(48, 0,5,23)
TestSphere.test_norm_swal_raising(48, 2,5,23)

println("Testing spherical spin-weighted lowering operator")
TestSphere.test_norm_swal_lowering(32, 0,0,2)

TestSphere.test_norm_swal_lowering(32, 0,1,2)
TestSphere.test_norm_swal_lowering(32, 1,1,2)

TestSphere.test_norm_swal_lowering(32,-1,-2,2)
TestSphere.test_norm_swal_lowering(32, 0,-2,2)
TestSphere.test_norm_swal_lowering(32, 2,-2,2)

TestSphere.test_norm_swal_lowering(32,-1,2,2)
TestSphere.test_norm_swal_lowering(32, 0,2,2)
TestSphere.test_norm_swal_lowering(32, 2,2,2)

TestSphere.test_norm_swal_lowering(48,-1,5,23)
TestSphere.test_norm_swal_lowering(48, 0,5,23)
TestSphere.test_norm_swal_lowering(48, 2,5,23)

println("Testing spherical spin-weighted filter operator")
TestSphere.test_norm_filter(32,-2,2,2)
TestSphere.test_norm_filter(32,-1,3,4)
TestSphere.test_norm_filter(32, 0,1,5)
TestSphere.test_norm_filter(32, 1,0,2)
TestSphere.test_norm_filter(32, 2,-1,2)


text:
const colnames_style = Dict(:color=>:yellow,)
const rownames_style = Dict(:color=>:cyan,)

# Base.show
function Base.show(io::IO, colored::ColoredString)
    if Base.have_color == true && colored.style !== nothing
         printstyled(io, colored.data, color=colored.style[:color])
    else
         print(io, colored.data)
    end
end

function Base.show(io::IO, v::Vertical)
    print(io, v.data)
end

function Base.show(io::IO, dash::Dash)
    print(io, repeat(dash.data, dash.repeat))
end

function Base.show(io::IO, connector::Connector)
    print(io, connector.data)
end

function print_plate(io::IO, horizon::Horizontal, tablemode::TableMode, islast::Bool)
    @inbounds for j=1:length(horizon)
        hr = horizon[j]
        print(io, hr)
    end
    !islast && println(io)
end

function print_plate(io::IO, linear::Linear, tablemode::TableMode, islast::Bool)
    isempty(linear) && return
    firstcell = first(linear)
    height = firstcell.height
    @inbounds for i=1:height
        for j=1:length(linear)
            cv = linear[j]
            if isa(cv, Vertical)
                print(io, cv)
            else
                cel = cv
                rows,cols = arraysize(cel.data)
                if rows >= i
                    for n=1:cols
                        print(io, cel.data[i,n])
                    end
                else
                    print(io, repeat(" ", cel.width))
                end
            end
        end
        if tablemode.style === :markdown && i==height
            !islast && println(io)
        else
            println(io)
        end
    end
end

function print_plates(io::IO, plates::PlateVector, tablemode::TableMode)
    len = length(plates)
    @inbounds for i=1:len
        plate = plates[i]
        print_plate(io, plate, tablemode, i==len)
    end
end

function Base.show(io::IO, mill::Mill)
    prepend_newline = display_style[:prepend_newline]
    prepend_newline && println(io)
    tablemode = display_style[:table_mode]
    plates = decking(mill, tablemode)
    print_plates(io, plates, tablemode)
end

# Base.isless
Base.isless(::String, ::Symbol) = false

# arraysize
arraysize(a::AbstractArray) = size(a[:,:])

# precell
precell(::Nothing, style=nothing) = DataCell([], 0, 0, style)

function precell(el::String, style=nothing)
    if occursin("\n", el)
        a = map(split(el, "\n")) do x
            width = textwidth(x)
            string(x, repeat(" ", width-textwidth(x)))
        end
        m,n = arraysize(a)
        DataCell(a, maximum(map(length, a)), m, style)
    else
        width = textwidth(el)
        DataCell([el], width, 1, style)
    end
end

function precell(el::AbstractArray, style=nothing)
    rows,cols = arraysize(el)
    if 0==rows
        DataCell([""], 0, 0, style)
    else
        prelines = Array{String}(undef, size(el))
        widths = zeros(Int, rows, cols)
        @inbounds for i=1:rows
            for j=1:cols
                s = replace(string(el[i,j]), "\n" => " ")
                widths[i,j] = length(s)
                prelines[i,j] = s
            end
        end
        maxwidths = zeros(Int, cols)
        @inbounds for j=1:cols
            maxwidths[j] = maximum(widths[:,j])
        end

        lines = AbstractArray{String}([])
        @inbounds for i=1:rows
            line = AbstractArray{String}([])
            for (j,x) in enumerate(prelines[i,:])
                push!(line, lpad(x, maxwidths[j]))
            end
            push!(lines, join(line, " "))
        end
        DataCell(lines, sum(maxwidths) + (cols-1), rows, style)
    end
end

function precell(el::Any, style=nothing)
    str = sprint(show, el, context = :color => false)
    l = split(str, '\n')
    m = maximum(length.(l))
    n = length(l)
    str_color = sprint(show, el, context = :color => true)
    l_color = split(str_color, '\n')
    DataCell(rpad.(l_color, m), m, n, style)
end

# postcell
function postcell(cell::DataCell, width::Int, height::Int, margin::Margin)
    data = cell.data
    rows,cols = arraysize(data)
    A = Array{ColoredString}(undef, rows, cols)
    @inbounds for i=1:rows
        for j=1:cols
            x = data[i,j]
            str = replace(x, r"\e[[0-9;]*[mGKF]" => "")
            padding = margin.leftside + width - textwidth(str) - 1
            prep = repeat(' ', padding > 0 ? padding : 0)
            A[i,j] = ColoredString(string(' ', x, prep, repeat(" ", margin.rightside)), cell.style)
        end
    end
    Cell(A, width+margin.leftside+margin.rightside, height)
end

function horizontal(maxwidths::Vector{Int}, cols::Int, margin::Margin, tablemode::TableMode; dash="-")
    corner_corner = string(tablemode.corner_corner)
    horizon = Horizontal{Union{Dash,Connector}}([])
    push!(horizon, Connector(corner_corner))
    @inbounds for j=1:cols
        width = maxwidths[j]
        push!(horizon, Dash(dash, width + margin.leftside + margin.rightside))
        push!(horizon, Connector(corner_corner))
    end
    horizon
end

function marginal(option::Dict{Symbol,Any})
    Margin(1, 1) # fixed
end

function decking(mill::Mill, tablemode::TableMode)
    board = mill.board
    option = mill.option
    header_fillchar = string(tablemode.header_fillchar)

    rows,cols = arraysize(board)
    plates = PlateVector{Union{Linear,Horizontal}}([])
    if 0==rows
        return plates
    end

    widths = zeros(Int, rows+1, cols+1)
    heights = zeros(Int, rows+1, cols+1)
    # prepare ells
    preplates = Vector{Vector{DataCell}}([])
    headlinear = Vector{Union{DataCell,Vertical}}([])
    zerocolname = nothing
    if haskey(option, :zerocolname)
        zerocolname = option[:zerocolname]
    end
    cel = precell(zerocolname, colnames_style)
    widths[1,1] = cel.width
    heights[1,1] = cel.height
    push!(headlinear, cel)
    ols = Vector{DataCell}(undef, cols)
    @inbounds for j=1:cols
        ols[j] = precell(j, colnames_style)
    end
    if haskey(option, :colnames)
        @inbounds for (j,name) in enumerate(option[:colnames])
            if cols >= j
                ols[j] = precell(name, colnames_style)
            end
        end
    end
    @inbounds for j=1:cols
        cel = ols[j]
        widths[1,j+1] = cel.width
        heights[1,j+1] = cel.height
        push!(headlinear, cel)
    end
    push!(preplates, headlinear)

    prerows = Vector{DataCell}(undef, rows)
    @inbounds for i=1:rows
        prerows[i] = precell(i, rownames_style)
    end
    if haskey(option, :rownames)
        @inbounds for (i,name) in enumerate(option[:rownames])
            if rows >= i
                prerows[i] = precell(name, rownames_style)
            end
        end
    end
    @inbounds for i=1:rows
        linear = Vector{Union{DataCell,Vertical}}([])
        rownamecell = prerows[i]
        widths[1+i,1] = rownamecell.width
        heights[1+i,1] = rownamecell.height
        push!(linear, rownamecell)
        for j=1:cols
            cel = precell(board[i,j])
            widths[1+i,1+j] = cel.width
            heights[1+i,1+j] = cel.height
            push!(linear, cel)
        end
        push!(preplates, linear)
    end

    maxwidths = zeros(Int, cols+1)
    @inbounds for j=1:cols+1
        maxwidths[j] = maximum(widths[:,j])
    end
    maxheights = zeros(Int, rows+1)
    @inbounds for i=1:rows+1
        maxheights[i] = maximum(heights[i,:])
    end
    margin = marginal(option)
    if 0==rows
        cols = 1
    else
        rows += 1
        cols += 1
    end

    isgridmode = tablemode.style === :grid
    isgridmode && push!(plates, horizontal(maxwidths, cols, margin, tablemode, dash="-"))
    @inbounds for i=1:rows
        linear = Linear{Union{Cell,Vertical}}([])
        preplate = preplates[i]
        push!(linear, Vertical(maxheights[i]))
        maxheight = maxheights[i]
        for j=1:cols
            el = preplate[j]
            cel = postcell(el, maxwidths[j], maxheight, margin)
            if 1==j && cel.height > 1
                prep = [repeat(" ", cel.width) for x in 1:(cel.height-length(cel.data))]
                cel.data = vcat(prep, cel.data)
            end
            push!(linear, cel)
            push!(linear, Vertical(maxheight))
        end
        push!(plates, linear)
        (1==i || isgridmode) && push!(plates, horizontal(maxwidths, cols, margin, tablemode, dash = 1==i ? header_fillchar : "-"))
    end
    plates
end

# table
function table(board::Type{T}; options...) where T
    zerocolname = (String ∘ nameof)(T)
    options = merge(Dict(options), Dict(:zerocolname => zerocolname, :rownames => string.(fieldnames(T)), :colnames => [""]))
    Mill(hcat(collect(fieldtypes(T))); options...)
end

function table(board::T; options...) where T
    names = fieldnames(T)
    if !isempty(names)
        zerocolname = (String ∘ nameof)(T)
        options = merge(Dict(options), Dict(:zerocolname => zerocolname, :rownames => [""], :colnames => string.(names)))
        a = collect(getfield.(Ref(board), names))
        board = reshape(a, (1,length(a)))
        Mill(board; options...)
    else
        Mill([board]; options...)
    end
end
function table(board::AbstractArray; options...)
    Mill(board; options...)
end
function table(board::Vector; options...)
    isempty(board) && return Mill(board; options...)
    first_eltype = (typeof ∘ first)(board)
    if all(x -> x isa first_eltype, board)
        if first_eltype <: NamedTuple
            options = merge(Dict(options), Dict(:colnames => string.(keys(first(board)))))
            board = hcat(vcat.(values.(board)...)...)
        else
            names = fieldnames(first_eltype)
            if !(first_eltype <: Tuple) && !isempty(names)
                zerocolname = string(:Vector, '{', (String ∘ nameof)(first_eltype), '}')
                options = merge(Dict(options), Dict(:zerocolname => zerocolname, :colnames => string.(names)))
                board = hcat(vcat.(map(x -> getfield.(Ref(x), names), board)...)...)
            end
        end
    end
    Mill(board; options...)
end
table(board::Tuple; options...) = Mill(collect(board); options...)

function table(board::Dict; options...)
    dict = Dict{Symbol, Any}(options)
    if !haskey(dict, :colnames)
        dict[:zerocolname] = "KEY"
        dict[:colnames] = ["VALUE"]
    end
    rownames = sort(collect(keys(board)))
    dict[:rownames] = rownames
    Mill(getindex.(Ref(board), rownames); tuple(dict...)...)
end

table(board::T; options...) where {T<:Base.Enums.Enum} = table(typeof(board), options...)
function table(board::Type{T}; options...) where {T<:Base.Enums.Enum}
    dict = Dict{Symbol, Any}(options)
    enums = instances(board)
    dict[:rownames] = (String ∘ Symbol).(enums)
    if !haskey(dict, :colnames)
        dict[:zerocolname] = (String ∘ nameof)(board)
        dict[:colnames] = ["VALUE"]
    end
    Mill(collect(Int.(enums)); tuple(dict...)...)
end