text:
# This code is modified from part of my master thesis. I tried to use the semiparametric method to
# develop a new Feasible-GLS estimating procedure. In the first step, a standard FGLS estimating
# procedure is implemented with a non-parametrically identified skedactic function Σ. Then, I resemple
# the test statistic using the wild bootstrap method, and then compare the FGLS t-statistic to the
# artificially-generated distribution. This method is proved to have great performance in finite samples
# in terms of lower size distortion and better testing power. The following code demonstrates the method
# using k-nearest neighbors estimator in the first step. Also, a simple leave-one-out cross validation
# is applied to choose the optimal bandwidth.

using LinearAlgebra, Distributions, StatsBase, BenchmarkTools

function KNN(h, X, y, e)
    #-----------------------------------------------------#
    # h : bandwidth                                       #
    # X : Nx5 matrix of constant and explanatory variable #
    # y : Nx1 matrix of explained variable                #
    # e : Nx1 matrix of residuals                         #
    #-----------------------------------------------------#
    N = length(y)
    w = zeros(N, N)
    for pp in 1:N
        x_normsort = sortperm(view.(Ref(w), 1:size(w,1), :))
        v = x_normsort .<= h
        v[pp] = 0
        w[:, pp] = v ./ sum(v)
    end
    s2hat = transpose(w) * e.^2
    small = mean(e.^2) / 10
    s2hat[s2hat .< small] .= small
    return(s2hat)
end

function CrossValidation(H, X, y, e)
    #-----------------------------------------------------#
    # H : vector of bandwidth candidates                  #
    #-----------------------------------------------------#
    CV = []
    for r in 1:length(H)
        s2hat_cv = KNN(H[r], X, y, e)
        error = sum((e.^2 .- s2hat_cv).^2)
        CV = append!(CV, error)
    end
    h = H[argmin(CV)]
    s2hat = KNN(h, X, y, e)
    Σ = Diagonal(s2hat)
    β̂ = (inv(transpose(X) * inv(Σ) * X) * transpose(X) * inv(Σ) * y)
    return(tuple(β̂, Σ))
end

function WildBootstrap(iter, X, y, e, β̂, Σ)
    N = length(y)
    t_memory = []
    for b in 1:iter
        u = float(rand(Bernoulli(0.5), 100))
        u[u.==0] .= -1
        ỹ = X * β̂ + e.*u
        β̃5 = (inv(transpose(X) * inv(Σ) * X) * transpose(X) * inv(Σ) * ỹ)[5, :]
        tstat = β̃5 / sqrt(inv(transpose(X) * inv(Σ) * X)[5, 5])
        t_memory = append!(t_memory, tstat)
    end
    return(sort(t_memory))
end

# DGP
x1 = ones(100)
x25 = rand(Normal(3, 1), 100, 4)
X = hcat(x1, x25)
β = [1, 1, 1, 1, 0]
ε = rand(100)
y = X * β + ε
β̂ = inv(transpose(X) * X) * transpose(X) * y
e = y - X * β̂

# FGLS estimation
H = range(1, stop = 8, length = 10) * (2.50773668907 + 0.021528881383 * (100^(4 / 5)))
FGLS = CrossValidation(H, X, y, e)
β̂5 = FGLS[1][5, :]
Σ = FGLS[2]
tstat = β̂5 / sqrt(inv(transpose(X) * inv(Σ) * X)[5, 5])

# Resempling t-statistic using wild bootstrap
T = WildBootstrap(1000, X, y, e, β̂, Σ)
uppercrit = percentile(T, 75)
lowercrit = percentile(T, 25)
if tstat[1] >= uppercrit || tstat[1] <= lowercrit
    println("Reject the null hypothesis")
else
    println("Do not reject the null hypothesis")
end

# Here are the executing time of the functions
@btime CrossValidation(H, X, y, e)
@btime WildBootstrap(1000, X, y, e, β̂, Σ)
# The CrossValidation function is slow since it contains multi-layers of loops.

# The following functions use multi-threading to improve the speed.
function KNN_multi(h, X, y, e)
    #-----------------------------------------------------#
    # h : bandwidth                                       #
    # X : Nx5 matrix of constant and explanatory variable #
    # y : Nx1 matrix of explained variable                #
    # e : Nx1 matrix of residuals                         #
    #-----------------------------------------------------#
    N = length(y)
    w = zeros(N, N)
    Threads.@threads for pp in 1:N
        x_normsort = sortperm(view.(Ref(w), 1:size(w,1), :))
        v = x_normsort .<= h
        v[pp] = 0
        w[:, pp] = v ./ sum(v)
    end
    s2hat = transpose(w) * e.^2
    small = mean(e.^2) / 10
    s2hat[s2hat .< small] .= small
    return(s2hat)
end

function CrossValidation_multi(H, X, y, e)
    #-----------------------------------------------------#
    # H : vector of bandwidth candidates                  #
    #-----------------------------------------------------#
    CV = []
    Threads.@threads for r in 1:length(H)
        s2hat_cv = KNN_multi(H[r], X, y, e)
        error = sum((e.^2 .- s2hat_cv).^2)
        CV = append!(CV, error)
    end
    h = H[argmin(CV)]
    s2hat = KNN_multi(h, X, y, e)
    Σ = Diagonal(s2hat)
    β̂ = (inv(transpose(X) * inv(Σ) * X) * transpose(X) * inv(Σ) * y)
    return(tuple(β̂, Σ))
end

# It is around 2 times faster when using 4 threads on my computer
@btime CrossValidation_multi(H, X, y, e)

max_stars_repo_path: Problem Sets/knn.jl

max_issues_repo_path: Problem Sets/knn.jl

max_forks_repo_path: Problem Sets/knn.jl


text:
export compensated_sum, recursive_sum, blocked_sum, fab_sum

# Kahan sum
function compensated_sum(A)
    T = eltype(A)
    sum = c = zero(T)
    for a in A
        y = a - c
        t = sum + y
        c = (t - sum) - y
        sum = t
    end
    return sum
end

function recursive_sum(A)
    T = eltype(A)
    sum = zero(T)
    for a in A
        sum += a
    end
    sum
end

function blocked_sum(A, b::Integer=128)
    n = length(A)
    m = (n+b-1) ÷ b
    S = [recursive_sum(view(A,i*b+1:min(i*b+b,n))) for i = 0:m-1]
    return sum(S)
end

function fab_sum(A, b::Integer=128)
    n = length(A)
    m = (n+b-1) ÷ b
    S = [recursive_sum(view(A,i*b+1:min(i*b+b,n))) for i = 0:m-1]
    return compensated_sum(S)
end

ext: jl


text:
# Executable to map MESH descriptors to UMLS concepts. Descriptors are
# assumed to be stored in a SQLite database containing table "mesh_descriptor"
# Date: May 6, 2016
# Authors: Isabel Restrepo
# BCBI - Brown University
# Version: Julia 0.4.5

using ArgParse
using PubMedMiner
using SQLite
using MySQL
import BioMedQuery.UMLS:Credentials

function main(args)

    # initialize the settings (the description is for the help screen)
    s = ArgParseSettings()
    s.description = "This is pubMedMiner.jl - main script"
    s.version = "1.0"

    @add_arg_table s begin
        "--append_results"
            help = "Flag to indicate not to clear the results TABLE:umls2mesh"
            action = :store_true
        "--umls_username"
            help = "UMLS-NLM username"
            arg_type = ASCIIString
            required = true
        "--umls_password"
            help = "UMLS-NLM password"
            arg_type = ASCIIString
            required = true
        "mysql"
            help = "Use MySQL backend"
            action = :command
        "sqlite"
             help = "Use SQLite backend"
             action = :command
    end

    @add_arg_table s["mysql"] begin
     "--host"
         help = "Host where you database lives"
         arg_type = ASCIIString
         default = "localhost"
     "--dbname"
         help = "Database name"
         arg_type = ASCIIString
         required = true
     "--username"
         help = "MySQL username"
         arg_type = ASCIIString
         default = "root"
     "--password"
         help = "MySQL password"
         arg_type = ASCIIString
         default = ""
    end

     @add_arg_table s["sqlite"] begin
     "--db_path"
          help = "Path to SQLite database file to store results"
          arg_type = ASCIIString
          required = true
     end

    parsed_args = parse_args(s)
    println("-------------------------------------------------------------")
    println(s.description)
    println("Parsed args:")
    for (key,val) in parsed_args
        println("  $key   =   $(repr(val))")
    end
    println("-------------------------------------------------------------")

    user = parsed_args["umls_username"]
    psswd = parsed_args["umls_password"]
    credentials = Credentials(user, psswd)
    append = parsed_args["append_results"]
    db = nothing

    if haskey(parsed_args, "sqlite")
        db_path  = parsed_args["sqlite"]["db_path"]
        db = SQLite.DB(db_path)
    elseif haskey(parsed_args, "mysql")
        host = parsed_args["mysql"]["host"]
        dbname = parsed_args["mysql"]["dbname"]
        username = parsed_args["mysql"]["username"]
        pswd = parsed_args["mysql"]["password"]
        db = MySQL.connect(host, username, pswd, dbname)
    else
        error("Unsupported database backend")
    end

    @time begin
        PubMedMiner.map_mesh_to_umls!(db, credentials; append_results=append)
    end
    println("-------------------------------------------------------------")
    println("Done Mapping Mesh to UMLS")
    println("-------------------------------------------------------------")

end

main(ARGS)


text:
using Test

include("zerorange.jl")
include("urange.jl")
include("zeroto.jl")  # deprecated

@test intersect(ZeroRange(4), URange(-1,2)) === 0:2
@test ZeroRange(6)[URange(2,4)] === 1:3

nothing


text:
struct IntVariate{T<:Integer} <: AbstractIntVariate{T}
    name::String
    data::Vector{T}
end

Base.length(intvariate::IntVariate{T}) where {T<:Signed} = length(intvariate.data)

function IntVariate{T}(name::String, length::Integer, datpath::String) where {T<:Signed}
    data = Vector{T}(undef, length)
    open(datpath) do f
        read!(f, data)
    end
    IntVariate{T}(name, data)
end

function slice(intvariate::IntVariate{T}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {T<:Signed}
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    slice(intvariate.data, fromobs, toobs, slicelength)
end

function Base.map(intvariate::IntVariate, dataframe::AbstractDataFrame; permute::Bool = false)
    dataframe[getname(intvariate)]
end

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
module LTI

include("_1_lti_primitives.jl")
include("_2_lti_systems.jl")
include("_3_lti_kernel.jl")
include("_4_lti_freq_resp.jl")
include("_5_lti_pprint.jl")
include("_6_lti_utils.jl")

end # module


text:
function Base.tryparse(::Type{Infinite}, str::AbstractString)
  if str == "∞"
    ∞
  elseif str == "-∞"
    -∞
  else
    nothing
  end
end


text:
using GeneralizedMonteCarlo
using Test, Distributions, Symbolics
@testset "inputPassthrough" begin
    @testset "PMF input" begin
        @test issubset(GeneralizedMonteCarlo.inputPassthrough([.01 .02 .03 .04 .05 .08; .1 .2 .3 .25 .10 .05]), [.01, .02, .03, .04, .05, .08])
        @test isa(GeneralizedMonteCarlo.inputPassthrough([.01 .02 .03 .04 .05 .08; .1 .2 .3 .25 .10 .05]), Float64)
        @test_throws ArgumentError GeneralizedMonteCarlo.inputPassthrough([.01 .02 .03 .04 .05 .08; .1 .2 .3 .25 .10 .05; 1 2 3 4 5 6])
        @test_throws ArgumentError GeneralizedMonteCarlo.inputPassthrough([.01 .02 .03 .04 .05 .08; .1 .2 .3 .25 .10 .15])
    end
    @testset "Sampleable input" begin
        @test 0 ≤ GeneralizedMonteCarlo.inputPassthrough(Distributions.Uniform(0, 1)) ≤ 1
        @test isa(GeneralizedMonteCarlo.inputPassthrough(Distributions.Uniform(0, 1)), Float64)
    end
    @testset "Real number input" begin
        @test GeneralizedMonteCarlo.inputPassthrough(0) == 0.0
        @test isa(GeneralizedMonteCarlo.inputPassthrough(0), Float64)
    end
end

@testset "GMC" begin
    @test begin @variables iws ms p; GeneralizedMonteCarlo.GMC(10_000, [iws, ms, 190 + ms + p], [iws, ms, p], [[10_000_000 10_500_000 11_000_000; 0.2 0.6 0.2], [0.01 0.02 0.03 0.04 0.05 0.08; 0.10 0.20 0.30 0.25 0.10 0.05], [-3.0 3.0; 0.5 0.5]])[1] isa Array end
end


text:
using DifferentialEquations, Plots, LaTeXStrings, LinearAlgebra
gr()
default(titlefont = 8, legendfontsize = 4, guidefont = 5, tickfont = 5, markersize =2.5)

methods = [Rodas4(),Rodas4P(),Rodas4P2(),Rodas5(),Rodas5P()];
s_methods = ["Rodas4","Rodas4P","Rodas4P2","Rodas5","Rodas5P"]
n_method = length(methods)

#---- Problem definitions --------------------------------------------------------------------
title_1 = "Polynomial of order 6"
y0_1 = [0.0; 0.0]; tspan_1 = [0.0, 1.0];  M_1 = zeros(2,2); M_1[1,1] = 1.0; p_1 = []; n_1 = [1,2,4,8,16,32] #- number of time steps
function ode_1(dy, y, p, t)
  dy[1] = 6*t^5
  dy[2] = y[1] - y[2]
  nothing
end
function sol_1(t,p)
  [t.^6, t.^6]
end

title_2 = "Prothero Robinson"
M_2 = zeros(1,1); M_2[1,1] = 1.0; y0_2 = [0.0]; tspan_2 = [0.0, 2.0]; p_2 = []; n_2 = [2,4,8,16,32,64,128,256];
function ode_2(dy, y, p, t)
  lam = 1.0e5;
  g=10.0-(10+t)*exp(-t);  dg = exp(-t)*(9+t);
  dy[1] = -lam*(y[1]-g)+dg;
  nothing
end
function sol_2(t,p)
  [10.0-(10+t)*exp(-t)]
end

title_3 = "Index 2"
M_3 = zeros(2,2); M_3[1,1] = 1.0; y0_3 = [-1.0; 1.0]; tspan_3 = [1.0, 2.0]; p_3 = []; n_3 = [32,64,128,256,512,1024];
function ode_3(du,u,p,t)
  du[1] = u[2]
  du[2] =u[1]^2-1/t^2
end
function sol_3(t,p)
  [-1 / t, 1 / t.^2]
end

title_4 = "Parabol"
nx = 1000; dx = 1.0/nx; n = nx-1; x = zeros(n)
M_4 = Matrix(1.0I, n, n);
for i=1:n
  x[i] = i*dx
end
b1  = zeros(n); b1[1]=3/(4*dx^2); b1[n]=3/(4*dx^2);
J = zeros(n,n)
for i=1:n
  J[i,i] = -2/(dx^2);
  if i>1 J[i,i-1] = 1/(dx^2); end
  if i<n J[i,i+1] = 1/(dx^2); end
end
p_4 = x, b1, J
function ode_4(du,u,p,t)
  x, b1, J = p
  g = -(x .+ 1/2).*(3/2 .- x)./(t+1).^2 .+ 2/(1+t);
  rb = b1/(1+t);
  du[:] = J*u .+ rb .+ g;
end
function sol_4(t,p)
  x, b1, J = p
  u = (x .+ 1/2).*(3/2 .- x)/(1+t);
end
y0_4 = sol_4(0.0,p_4); tspan_4 = [0.0, 0.1]; n_4 = [1,2,4,8];

title_5 = "Hyperbol"
nx = 1000; dx = 1.0/nx; n = nx; x = zeros(n)
M_5 = Matrix(1.0I, n, n);
for i=1:n
  x[i] = i*dx
end
b1  = zeros(n); b1[1]=-1/dx; J = zeros(n,n)
for i=1:n
  J[i,i] = 1/dx;
  if i>1 J[i,i-1] = -1/dx; end
end
p_5 = x, b1, J
function ode_5(du,u,p,t)
  x, b1, J = p
  g = (t .- x)./(t .+ 1).^2;
  rb = b1/(1+t);
  du[:] = -J*u .- rb .+ g;
end
function sol_5(t,p)
  x, b1, J = p
  u = (1 .+ x)./(1 .+ t);
end
y0_5 = sol_5(0.0,p_5); tspan_5 = [0.0, 1]; n_5 = [4,8,16,32,64,128];

#--- problems with constant stepsize -----
title_s = [title_1,title_2,title_3,title_4,title_5]
ode_s = [ode_1,ode_2,ode_3,ode_4,ode_5]
sol_s = [sol_1,sol_2,sol_3,sol_4,sol_5]
y0_s = [y0_1,y0_2,y0_3,y0_4,y0_5]
tspan_s = [tspan_1,tspan_2,tspan_3,tspan_4,tspan_5]
M_s = [M_1,M_2,M_3,M_4,M_5]
param_s = [p_1,p_2,p_3,p_4,p_5]
n_s = [n_1,n_2,n_3,n_4,n_5]

n_probs = length(ode_s); Pl = []

#--- Benchmarks for Rodas-type methods - --
for ip = 1:n_probs
  P_1 = plot(); P_2 = plot()
  f = ODEFunction(ode_s[ip], mass_matrix = M_s[ip]); tspan = tspan_s[ip]; prob = ODEProblem(f, y0_s[ip], tspan, param_s[ip])
  n = n_s[ip]; n_max = length(n); h = (tspan[2] - tspan[1])./n;
  err_end = zeros(n_max); err_dense = zeros(n_max)
  t_dense = range(tspan[1],tspan[2], length = 100)
  for m=1:n_method
    for i=1:n_max
       global sol = solve(prob, methods[m], dt=h[i], adaptive=false)
#       println(sol.destats);
       err_end[i] = max(1.0e-16,maximum(abs.(sol(tspan[2]) .- sol_s[ip](tspan[2],param_s[ip]))))
#       err_end[i] = 1.0e-16
#       for j = 1:length(sol.t)
#         err_end[i] = max(err_end[i], maximum(abs.(sol.u[j] .- sol_s[ip](sol.t[j],param_s[ip]))))
#       end
       err_dense[i] = 1.0e-16
       for j=1:length(t_dense)
         err_dense[i] = max(err_dense[i], maximum(abs.(sol(t_dense[j]) .- sol_s[ip](t_dense[j],param_s[ip]))))
       end
    end
    P_1 = plot!(P_1,h,err_end, xaxis=:log, yaxis=:log, xlabel = "stepsize", ylabel = L"$L_{\infty}$ error at $t_{\textrm{end}}$", label=s_methods[m], marker=([:circle :d]), legend=:bottomright)
    P_2 = plot!(P_2,h,err_dense, xaxis=:log, yaxis=:log, xlabel = "stepsize", ylabel = L"dense $L_{\infty}$ error", label=s_methods[m], marker=([:circle :d]), legend=:bottomright)
  end
  push!(Pl,plot(P_1,P_2,layout=(1,2),title=title_s[ip]))
end
p = plot(Pl[1],Pl[2],layout=(2,1))
savefig(p, "rodas5p_bench_1.pdf")
p = plot(Pl[3],Pl[4],layout=(2,1))
savefig(p, "rodas5p_bench_2.pdf")
p = plot(Pl[5],layout=(1,1))
savefig(p, "rodas5p_bench_3.pdf")


text:
"""
    StreamKmeans(m_target::Int)

Online clustering algorithm [1] to select inducing points in a streaming setting.
Reference :
[1] Liberty, E., Sriharsha, R. & Sviridenko, M. An Algorithm for Online K-Means Clustering. arXiv:1412.5721 [cs] (2015).
"""
mutable struct StreamKmeans{T} <: OnIPSA
    m_target::Int
    m_efficient::Int
    f::T
    q::Int
end

StreamKmeans(m_target::Int) = StreamKmeans(m_target, 0, 0.0, 0)

function Base.show(io::IO, alg::StreamKmeans)
    return print(io, "Streaming Kmeans (m_target=$(alg.m_target))")
end

function initZ(
    rng::AbstractRNG,
    alg::StreamKmeans,
    X::AbstractVector;
    kernel=nothing,
    arraytype=Vector{Float64},
    kwargs...,
)
    length(X) > 10 ||
        throw(ArgumentError("The first batch of data should be bigger than 10 samples"))
    m_efficient = max(1, ceil(Int, (alg.m_target - 15) / 5))
    alg.m_efficient = m_efficient + 10 > length(X) ? 0 : m_efficient
    samp = sample(rng, 1:length(X), alg.m_efficient + 10; replace=false)
    Z = to_vec_of_vecs(X[samp], arraytype)
    m = length(Z)
    w = zeros(m)
    for i in 1:m
        w[i] = 0.5 * find_nearest_center(Z[i], Z[1:end .!= i], kernel)[2]
    end
    alg.f = sum(sort(w)[1:10]) # Take the 10 smallest values
    alg.q = 0
    return Z
end

function add_point!(
    rng::AbstractRNG,
    Z::AbstractVector{T},
    alg::StreamKmeans,
    X::AbstractVector;
    kernel=nothing,
    kwargs...,
) where {T}
    b = length(X) # Size of the input data
    for i in 1:b
        val = find_nearest_center(X[i], Z, kernel)[2]
        if val > (alg.f * rand(rng))
            # new_centers = vcat(new_centers,X[i,:]')
            push!(Z, X[i])
            alg.q += 1
        end
        if alg.q >= alg.m_efficient
            alg.q = 0
            alg.f *= 10
        end
    end
    return Z
end

function add_point(
    rng::AbstractRNG,
    Z::AbstractVector{T},
    alg::StreamKmeans,
    X::AbstractVector;
    kernel=nothing,
    kwargs...,
) where {T}
    b = length(X) # Size of the input data
    for i in 1:b
        val = find_nearest_center(X[i], Z, kernel)[2]
        if val > (alg.f * rand(rng))
            # new_centers = vcat(new_centers,X[i,:]')
            Z = T isa Real ? vcat(Z, T(X[i])) : vcat(Z, T.(X[i:i]))
            alg.q += 1
        end
        if alg.q >= alg.m_efficient
            alg.q = 0
            alg.f *= 10
        end
    end
    return Z
end


text:
macro something(args...)
    if length(args) == 1 && (exprs = statements_in_braces(args[1])) !== nothing
        return something_expr(__module__, __source__, exprs)
    elseif length(args) == 1 && isexpr(args[1], :bracescat)
        return something_expr(__module__, __source__, args[1].args)
    end
    return something_expr(__module__, __source__, args)
end

function something_expr(__module__, __source__, @nospecialize(args))
    args = collect(Any, args)
    lnns = LineNumberNode[]
    ln = __source__
    for x in args
        ln = something(first_line_number_node(is_advancing(ln), x), ln)
        push!(lnns, ln)
    end

    foldr(
        collect(zip(args, lnns));  # `collect` for Julia 1.0
        init = :(throw(ArgumentError("all evaluated as `nothing`"))),
    ) do (x, ln), ex
        block = esc(Expr(:block, ln, x))
        @gensym v
        return Expr(
            :block,
            ln,
            :(local $v = $block),
            :($v !== nothing ? something($v) : $ex),
        )
    end
end

"""
    statements_in_braces(ex)::Union{Vector{Any},Nothing}

This is used only for Julia >= 1.5.

# Examples
```julia
julia> using Maybe.Implementations: statements_in_braces

julia> statements_in_braces(:non_expression)

julia> statements_in_braces(:(@f {a; b; c}).args[end])
3-element Array{Any,1}:
 :a
 :b
 :c
```

# Notes
`@something{a; b; c}` and `@something{a; b; c;}` are parsed as

```lisp
(:macrocall, Symbol("@something"),
             :(#= ... =#),
             (:braces, (:parameters, (:parameters, :c), :b), :a))

(:macrocall, Symbol("@something"),
             :(#= ... =#),
             (:braces, (:parameters, (:parameters, (:parameters,), :c), :b), :a))
```

`parameters_as_statements` called from `statements_in_braces` converts
the `(:braces, ...)` part to `[:a, :b, :c]`.
"""
function statements_in_braces(ex)::Union{Vector{Any},Nothing}
    @nospecialize
    isexpr(ex, :braces) || return nothing
    isempty(ex.args) && return []
    if isexpr(ex.args[1], :parameters)
        return parameters_as_statements(Expr(:parameters, ex.args...))
    elseif length(ex.args) == 1
        return ex.args
    else
        return nothing
    end
end

function parameters_as_statements(ex)::Vector{Any}
    @nospecialize
    isexpr(ex, :parameters) || return Any[ex]
    if length(ex.args) == 0
        return []
    elseif length(ex.args) == 1
        return parameters_as_statements(ex.args[1])
    elseif length(ex.args) == 2 && isexpr(ex.args[1], :parameters)
        return append!(Any[ex.args[2]], parameters_as_statements(ex.args[1]))
    else
        error(
            "Semicolons and commas should not be mixed. ",
            "Use `{a; b; c}` instead of `{a; b, c}`.\n",
            "Error while parsing:\n",
            ex,
        )
    end
end


text:
@testset "right" begin
    _test_binary_op(TimeDag.right, (x, y) -> y)
end

@testset "left" begin
    _test_binary_op(TimeDag.left, (x, y) -> x)
end

@testset "align" begin
    @test _eval(TimeDag.align(n4, n1)) == b1
    @test _eval(TimeDag.align(n1, n4)) == Block([
        DateTime(2000, 1, 1) => 1,
        DateTime(2000, 1, 2) => 2,
        DateTime(2000, 1, 3) => 3,
        DateTime(2000, 1, 4) => 4,
        DateTime(2000, 1, 5) => 4,
        DateTime(2000, 1, 6) => 4,
        DateTime(2000, 1, 7) => 4,
    ])
end

@testset "coalign" begin
    @testset "unary" begin
        @test coalign(n4) === n4
        @test coalign(n4; alignment=INTERSECT) === n4
        @test coalign(n4; alignment=LEFT) === n4
        @test coalign(n4; alignment=UNION) === n4
    end

    @testset "binary" begin
        for alignment in (INTERSECT, LEFT, UNION)
            na, nb = coalign(n1, n2; alignment)
            @test _eval(na) == _eval(left(n1, n2, alignment))
            @test _eval(nb) == _eval(right(n1, n2, alignment))
        end
    end

    @testset "ternary" begin
        for alignment in (INTERSECT, LEFT, UNION)
            na, nb, nc = coalign(n1, n2, n3; alignment)

            l(x, y) = left(x, y, alignment)
            r(x, y) = right(x, y, alignment)

            @test _eval(na) == _eval(l(n1, l(n2, n3)))
            @test _eval(nb) == _eval(r(n1, l(n2, n3)))
            @test _eval(nc) == _eval(r(n1, r(n2, n3)))
        end
    end

    @testset "ordering" begin
        # We don't want to generate new different nodes if we coalign the same nodes in a
        # different order, when alignment is order-independent.
        nodes = [n1, n2, n3]
        for alignment in (INTERSECT, UNION)
            for indices in permutations([1, 2, 3])
                new_nodes = coalign(nodes[indices]...; alignment)
                @test new_nodes === coalign(nodes...; alignment)[indices]
            end
        end

        # For the case of LEFT, the first node must stay in place, but we should be able to
        # permute the order of the others and not change the answer.
        for indices in permutations([1, 2, 3])
            new_nodes = coalign(n4, nodes[indices]...; alignment=LEFT)
            @test new_nodes === coalign(n4, nodes...; alignment=LEFT)[[1; 1 .+ indices]]
        end
    end
end

@testset "first_knot" begin
    @test _eval(first_knot(n4)) == b4[1:1]
    @test _eval(first_knot(n_boolean)) == b_boolean[1:1]
    @test _eval(empty_node(Float64)) == Block{Float64}()
end

@testset "active_count" begin
    @test _eval(active_count(n1)) == Block([b1.times[1]], [1])
    @test _eval(active_count(n1, n2)) == Block([b1.times[1], b2.times[1]], [1, 2])
    @test _eval(active_count(n1, n2, n3)) == Block([b1.times[1], b2.times[1]], [2, 3])
    @test _eval(active_count(n1, n2, n3, n4)) == Block([b1.times[1], b2.times[1]], [3, 4])

    # Testing for optimisations.
    @test active_count(n1, n2) === active_count(n2, n1)
    @test active_count(n1, n2, n3) === active_count(n2, n1, n3)
    @test active_count(n1, n2, n3) === active_count(n3, n1, n2)
    @test active_count(n1, n2, n3) === active_count(n3, n2, n1)
end


text:
# module for representing agent
# decide control order to robot
# log pose and observation into text file

include(joinpath(split(@__FILE__, "src")[1], "src/common/state_transition/state_transition.jl"))

mutable struct LoggerAgent
  speed
  yaw_rate
  time_interval
  estimator
  prev_spd
  prev_yr
  pose
  step
  log
  is_test

  # init
  function LoggerAgent(speed::Float64, yaw_rate::Float64;
                       time_interval::Float64=0.1,
                       estimator=nothing,
                       init_pose=[0.0, 0.0, 0.0],
                       is_test=false)
      self = new()
      self.speed = speed
      self.yaw_rate = yaw_rate
      self.time_interval = time_interval
      self.estimator = estimator
      self.prev_spd = 0.0
      self.prev_yr = 0.0
      self.pose = init_pose
      self.step = 0
      self.is_test = is_test
      if self.is_test == false
        log_path = "src/slam/graph_based_slam/traj_obsrv_edge_input_log.txt"
        self.log = open(joinpath(split(@__FILE__, "src")[1], log_path), "w")
        # record delta time
        write(self.log, "delta $(time_interval)\n")
      end
      return self
  end
end

function close_log_file(self::LoggerAgent)
  if self.is_test == false
    close(self.log)
  end
end

function draw_decision!(self::LoggerAgent, observation)
  if self.is_test == false
    # logging trajectory and observation
    write(self.log, "u $(self.step) $(self.speed) $(self.yaw_rate)\n")
    write(self.log, "x $(self.step) $(self.pose[1]) $(self.pose[2]) $(self.pose[3])\n") # step x y theta
    for obs in observation
      write(self.log, "z $(self.step) $(obs[2]) $(obs[1][1]) $(obs[1][2]) $(obs[1][3])\n") # step id distance direction orientation
    end
    self.step += 1
  end

  # draw estimated pose
  if self.estimator !== nothing
    motion_update(self.estimator, self.prev_spd, self.prev_yr, self.time_interval)
    self.prev_spd, self.prev_yr = self.speed, self.yaw_rate
    observation_update(self.estimator, observation)
    draw!(self.estimator)
  end

  self.pose = state_transition(self.speed, self.yaw_rate, self.time_interval, self.pose)
  return self.speed, self.yaw_rate
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
using Combinatorics

# Checks if val is sum of n values in arr
function is_sum(val, arr; n=2)
    s = with_replacement_combinations(arr, n)
    sums = []
    for i in s
        push!(sums, sum(i))
    end

    if val in sums
        return true
    else
        return false
    end
end

# Check all vals in a window of size preamble
function check_vals(file; preamble = 25)
    i = 1
    while i <= (length(file) - preamble)
        j = i + preamble - 1
        arr = file[i:j]
        val = file[j+1]

        if is_sum(val, arr)
            i += 1
        else
            return val
        end
    end
end

# Function to find the weakness
# Defined as the sum of the min and max
# values of a continuous section of code
# that sum to the invalid value
function find_weakness(file, val)
    i = 1
    di = 1

    # Do forever!
    while true
        j = i + di
        arr = file[i:j]
        s = sum(arr)

        if s == val
            return minimum(arr) + maximum(arr)
        elseif s > val
            i += 1
            di = 1
        else
            di += 1
        end
    end

end

# Test values
test_vals = [35, 20, 15, 25, 47, 40, 62, 55, 65, 95, 102, 117, 150, 182, 127, 219, 299, 277, 309, 576,]

data =  parse.(Int, readlines("./XMAS_data.txt"))
println("Test: ", check_vals(test_vals, preamble=5))
inval = check_vals(data)
println("Real: ", inval)
println("Encryption weakness: ", find_weakness(data, inval))


text:
### UTILS
FloatOrNothing = Union{Float64, Nothing}
IntOrNothing = Union{Int64, Nothing}

### IDA
struct IDA
    name::String # Index Name
    date::Date    # Date
    index::Float64 # Index
    return_1d::FloatOrNothing # 1 Day (%)
    return_month_to_date::FloatOrNothing # Month-to-Date (%)
    return_year_to_date::FloatOrNothing # Year-to-Date (%)
    return_12_months::FloatOrNothing # 12-Months (%)
    return_24_months::FloatOrNothing # 24-Months (%)
    duration::IntOrNothing # Duration (business day)
    weight::FloatOrNothing # Weight (%)
    market_value::FloatOrNothing # Portifolio at Market Value (R$ thousand)
end

struct IDAResult
    elements::Dict{String, IDA}
end

function IDAResult(v::Vector{IDA})
    d = Dict([x.name => x for x in v])
    return IDAResult(d)
end

Base.getindex(x::IDAResult, i::String) = getindex(x.elements, i)
Base.keys(x::IDAResult) = keys(x.elements)

### IMA
struct IMA
    name::String # Index Name
    date::Date    # Date
    index::Float64 # Index
    return_1d::FloatOrNothing # 1 Day (%)
    return_month_to_date::FloatOrNothing # Month-to-Date (%)
    return_year_to_date::FloatOrNothing # Year-to-Date (%)
    return_12_months::FloatOrNothing # 12-Months (%)
    return_24_months::FloatOrNothing # 24-Months (%)
    weight::FloatOrNothing # Weight (%)
    duration::FloatOrNothing # Duration (business day)
    market_value::FloatOrNothing # Portifolio at Market Value (R$ thousand)
    number_of_trades::IntOrNothing # Number of Trades *
    amount_traded_1000_bonds::FloatOrNothing # Amount Traded (1,000 bonds)
    amount_traded_thousand_BRL::FloatOrNothing # Amount Traded (R$ thousand)
    pmr::FloatOrNothing # PMR
    convexity::FloatOrNothing # Convexidade
    yield::FloatOrNothing # Yield
    redemption_yield::FloatOrNothing # Redemption Yield
end

struct IMAResult
    elements::Dict{String, IMA}
end

function IMAResult(v::Vector{IMA})
    d = Dict([x.name => x for x in v])
    return IMAResult(d)
end

Base.getindex(x::IMAResult, i::String) = getindex(x.elements, i)
Base.keys(x::IMAResult) = keys(x.elements)

struct SvenssonCurve
    β1::Float64
    β2::Float64
    β3::Float64
    β4::Float64
    λ1::Float64
    λ2::Float64
end

"Zero rates for a SvenssonCurve. `t` in business days. Returns 0.1 for a 10% p.a. rate."
function zerorate(cv::SvenssonCurve, t::Int64)::Float64
    β1 = cv.β1
    β2 = cv.β2
    β3 = cv.β3
    β4 = cv.β4
    λ1 = cv.λ1
    λ2 = cv.λ2
    τ = t/252
    return (
        β1 +
        β2 *  (1-exp(-λ1*τ)) / (λ1*τ) +
        β3 * ((1-exp(-λ1*τ)) / (λ1*τ) - exp(-λ1*τ)) +
        β4 * ((1-exp(-λ2*τ)) / (λ2*τ) - exp(-λ2*τ))
    )
end

struct NelsonSiegelCurve
    β1::Float64
    β2::Float64
    β3::Float64
    λ::Float64
end

function zerorate(cv::NelsonSiegelCurve, t::Int64)::Float64
    β1 = cv.β1
    β2 = cv.β2
    β3 = cv.β3
    λ = cv.λ
    τ = t/252
    return (
        β1 +
        β2 *  (1-exp(-λ*τ)) / (λ*τ) +
        β3 * ((1-exp(-λ*τ)) / (λ*τ) - exp(-λ*τ))
    )
end


text:
using GeoStats

@info "All tests of the GeoStats.jl framework are available in submodules"


text:
# just visual tests!

using Mellan, Luxor

function alltests()
    mellanize(dirname(@__FILE__) * "/mona.png",
        800,
        lineweight      = 5,
        output          = "mona-mellan-800.pdf",
        foregroundcolor = "gray30",
        backgroundcolor = "antiquewhite2",
        startradius     = 5,
        margin          = 10,
        tightness       = 1,
        chord           = 4)

    mellanize(dirname(@__FILE__) * "/mona.png",
        500,
        lineweight      = 5,
        output          = "mona-mellan-500.pdf",
        foregroundcolor = "gray20",
        backgroundcolor = "antiquewhite2",
        startradius     = 3,
        margin          = 10,
        tightness       = 0.8,
        chord           = 3)

    mellanize(dirname(@__FILE__) * "/mona.png",
        500,
        lineweight      = 3,
        output          = "mona-mellan.svg",
        minlineweight   = 0.1,
        foregroundcolor =  "grey5",
        backgroundcolor =  "gold",
        startradius     = 5,
        margin          = 0,
        tightness       = 0.5,
        chord           = 3,
        annotation      = true)
end

if get(ENV, "MELLAN_KEEP_TEST_RESULTS", false) == "true"
        # they changed mktempdir in v1.3
        if VERSION <= v"1.2"
            cd(mktempdir())
        else
            cd(mktempdir(cleanup=false))
        end
        @info("...Keeping the Mellan test output in: $(pwd())")
        alltests()
        @info("Mellan test images were saved in: $(pwd())")
else
    mktempdir() do tmpdir
        cd(tmpdir) do
            @info("running Mellan tests in: $(pwd())")
            @info("but not keeping the results")
            @info("because you didn't do: ENV[\"MELLAN_KEEP_TEST_RESULTS\"] = \"true\"")
            alltests()
            @info("Test images weren't saved. To see the test images, next time do this before running:")
            @info(" ENV[\"MELLAN_KEEP_TEST_RESULTS\"] = \"true\"")
        end
    end
end

max_stars_repo_name: cormullion/Mellan

max_issues_repo_name: cormullion/Mellan

max_forks_repo_name: cormullion/Mellan


text:
# matplotlib.use("macosx")

@testset "eval_f_on_gauss_pts" begin
    f = (x,y)->x^2+y^2
    F = eval_f_on_gauss_pts(f, m, n, h)
    pts = AdFem.get_gauss_points(m, n, h)
    close("all")
    scatter3D(pts[:,1], pts[:,2], F, marker=".")
    mesh(X, Y, (@. f(X, Y)), color="orange", alpha=0.5)
    savefig("eval_f_on_gauss_pts.png")
end

@testset "compute_fvm_mechanics_term" begin
    u = zeros(2*(m+1)*(n+1))
    k = 0
    for j = 1:n+1
        for i = 1:m+1
            x = (i-1)*h
            y = (j-1)*h
            k += 1
            u[k] = x^2+y^2
            u[k + (m+1)*(n+1)] = x^2-y^2
        end
    end
    eps_v = compute_fvm_mechanics_term(u, m, n, h)
    close("all")
    scatter3D(Xv, Yv, eps_v, marker=".")
    scatter3D(Xv, Yv, 2h^2 * (@. Xv-Yv), marker="+")
    savefig("compute_fvm_mechanics_term.png")
end

@testset "compute_interaction_matrix" begin
    A = compute_interaction_matrix(m, n, h)
    u = zeros(2*(m+1)*(n+1))
    k = 0
    for j = 1:n+1
        for i = 1:m+1
            x = (i-1)*h
            y = (j-1)*h
            k += 1
            u[k] = x^2+y^2
            u[k + (m+1)*(n+1)] = x^2-y^2
        end
    end
    p = zeros(m*n)
    k = 0
    for j = 1:n
        for i = 1:m
            k += 1
            x = (i-1/2)*h
            y = (j-1/2)*h
            p[k] =  x+y
        end
    end
    val = p'*A*u
    @test val ≈ 64
end

@testset "compute_fvm_source_term" begin
    f = (x,y)->x^2+y^2
    F = eval_f_on_gauss_pts(f, m, n, h)
    S = compute_fvm_source_term(F, m, n, h)
    Fs = f.(Xv, Yv)*h^2
    close("all")
    scatter3D(Xv, Yv, S, marker = ".")
    scatter3D(Xv, Yv, Fs, marker = "+")
    savefig("compute_fvm_source_term.png")
end

@testset "compute_fvm_tpfa_matrix" begin
    M = compute_fvm_tpfa_matrix(m, n, h)
    p = Xv.^2 + Yv.^2
    V = M * p
    close("all")
    scatter3D(Xv, Yv, V, marker = ".")
    scatter3D(Xv, Yv, 4h^2 * ones(m*n), marker = "+")
    savefig("compute_fvm_tpfa_matrix.png")
end

@testset "static linear elasticity" begin
    bdedge = []
    for i = 1:m
        push!(bdedge, [i i+1])
    end
    bdedge = vcat(bdedge...)

    bdnode = Int64[]
    for j = 1:n+1
        push!(bdnode, (j-1)*(m+1)+1)
        push!(bdnode, (j-1)*(m+1)+m+1)
    end
    for i = 2:m
        push!(bdnode, n*(m+1)+i)
    end

    F1 = eval_f_on_gauss_pts((x,y)->3.0, m, n, h)
    F2 = eval_f_on_gauss_pts((x,y)->-1.0, m, n, h)
    F = compute_fem_source_term(F1, F2, m, n, h)

    t1 = eval_f_on_boundary_edge((x,y)->-x-y, bdedge, m, n, h)
    t2 = eval_f_on_boundary_edge((x,y)->2y, bdedge, m, n, h)
    T = compute_fem_traction_term([t1 t2], bdedge, m, n, h)

    D = diagm(0=>[1,1,0.5])
    K = compute_fem_stiffness_matrix(D, m, n, h)
    rhs = T - F
    bdval = [eval_f_on_boundary_node((x,y)->x^2+y^2, bdnode, m, n, h);
            eval_f_on_boundary_node((x,y)->x^2-y^2, bdnode, m, n, h)]
    rhs[[bdnode;bdnode .+ (m+1)*(n+1)]] = bdval
    K, Kbd = fem_impose_Dirichlet_boundary_condition(K, bdnode, m, n, h)
    u = K\(rhs-Kbd*bdval)
    X, Y, U, V = plot_u(u, m, n, h)

    figure(figsize=[10,4])
    subplot(121)
    pcolormesh(X, Y, (@. X^2+Y^2-U), alpha=0.6); xlabel("x"); ylabel("y"); title("Error for u")
    colorbar()
    subplot(122)
    pcolormesh(X, Y, (@. X^2-Y^2-V), alpha=0.6); xlabel("x"); ylabel("y"); title("Error for v")
    colorbar()

end

@testset "Poisson" begin
    m = 40
    n = 20
    h = 0.1
    # bdedge = []
    # for i = 1:m
    #     push!(bdedge, [i i+1])
    # end
    # bdedge = vcat(bdedge...)

    bdnode = Int64[]
    for j = 1:n+1
        push!(bdnode, (j-1)*(m+1)+1)
        push!(bdnode, (j-1)*(m+1)+m+1)
    end
    for i = 2:m
        push!(bdnode, n*(m+1)+i)
        push!(bdnode, i)
    end

    K_ = [2.0 1.0
        1.0 2.0]
    K = compute_fem_stiffness_matrix1(K_, m, n, h)

    K, Kbd = fem_impose_Dirichlet_boundary_condition1(K, bdnode, m, n, h)
    F = eval_f_on_gauss_pts((x,y)->-8, m, n, h)
    rhs = compute_fem_source_term1(F, m, n, h)
    bdval = eval_f_on_boundary_node( (x,y)->(x^2+y^2), bdnode, m, n, h)
    rhs[bdnode] = bdval
    U = K\(rhs-Kbd*bdval)

    Uexact = zeros(n+1,m+1)
    for j = 1:n+1
        for i = 1:m+1
            x = (i-1)*h; y = (j-1)*h
            Uexact[j, i] = (x^2+y^2)
        end
    end
    pcolormesh(reshape(U, m+1, n+1)'-Uexact); colorbar()
end

@testset "heat equation" begin
    m = 40
    n = 20
    h = 0.01
    bdedge = []
    for i = 1:m
        push!(bdedge, [i i+1])
    end
    bdedge = vcat(bdedge...)

    bdnode = Int64[]
    for j = 1:n+1
        push!(bdnode, (j-1)*(m+1)+1)
        push!(bdnode, (j-1)*(m+1)+m+1)
    end
    for i = 2:m
        push!(bdnode, n*(m+1)+i)
    end

    ρ = eval_f_on_gauss_pts((x,y)->1+x^2+y^2, m, n, h)
    M = compute_fem_mass_matrix1(ρ, m, n, h)

    K_ = [2.0 1.0
        1.0 2.0]
    K = compute_fem_stiffness_matrix1(K_, m, n, h)

    NT = 200
    Δt = 1/NT
    A = M/Δt+K
    A, Abd = fem_impose_Dirichlet_boundary_condition1(A, bdnode, m, n, h)
    U = zeros((m+1)*(n+1), NT+1)
    for i = 1:m+1
        for j = 1:n+1
            x = (i-1)*h; y = (j-1)*h
            U[(j-1)*(m+1)+i, 1] = x^2+y^2
        end
    end
    for i = 1:NT
        F = eval_f_on_gauss_pts((x,y)->(-(1+x^2+y^2)*(x^2+y^2)-8)*exp(-i*Δt), m, n, h)
        F = compute_fem_source_term1(F, m, n, h)

        T = eval_f_on_boundary_edge((x,y)->-(2*x+4*y)*exp(-i*Δt), bdedge, m, n, h)
        T = compute_fem_flux_term1(T, bdedge, m, n, h)

        rhs = F  + M*U[:,i]/Δt + T
        bdval = eval_f_on_boundary_node( (x,y)->(x^2+y^2)*exp(-i*Δt), bdnode, m, n, h)
        rhs[bdnode] = bdval
        U[:,i+1] = A\(
            rhs - Abd*bdval
        )
    end

    Uexact = zeros(n+1,m+1)
    for j = 1:n+1
        for i = 1:m+1
            x = (i-1)*h; y = (j-1)*h
            Uexact[j, i] = (x^2+y^2)*exp(-1)
        end
    end
    pcolormesh(reshape(U[:,end], m+1, n+1)'-Uexact); colorbar()
end

@testset "elastodynamics" begin
    # alpha method
    # table 9.1.1, page 493
    β = 1/4; γ = 1/2
    a = b = 0.1
    m = 40
    n = 20
    h = 0.01
    NT = 200
    Δt = 1/NT
    bdedge = []
    for j = 1:n
        push!(bdedge, [(j-1)*(m+1)+m+1 j*(m+1)+m+1])
    end
    bdedge = vcat(bdedge...)

    bdnode = Int64[]
    for j = 1:n+1
        push!(bdnode, (j-1)*(m+1)+1)
    end

    M = compute_fem_mass_matrix1(m, n, h)
    S = spzeros((m+1)*(n+1), (m+1)*(n+1))
    M = [M S;S M]

    D = diagm(0=>[1,1,0.5])
    K = compute_fem_stiffness_matrix(D, m, n, h)
    C = a*M + b*K # damping matrix

    L = M + γ*Δt*C + β*Δt^2*K
    L, Lbd = fem_impose_Dirichlet_boundary_condition(L, bdnode, m, n, h)

    a = zeros(2(m+1)*(n+1))
    v = zeros(2(m+1)*(n+1))
    d = zeros(2(m+1)*(n+1))
    U = zeros(2(m+1)*(n+1),NT+1)
    for i = 1:NT
        T = eval_f_on_boundary_edge((x,y)->0.01, bdedge, m, n, h)
        # T = eval_f_on_boundary_edge((x,y)->0.0, bdedge, m, n, h)
        T = [zeros(length(T)) -T]
        T = compute_fem_traction_term(T, bdedge, m, n, h)
        f1 = eval_f_on_gauss_pts((x,y)->0., m, n, h)
        f2 = eval_f_on_gauss_pts((x,y)->0., m, n, h)
        # f2 = eval_f_on_gauss_pts((x,y)->0.1, m, n, h)
        F = compute_fem_source_term(f1, f2, m, n, h)

        rhs = F+T

        td = d + Δt*v + Δt^2/2*(1-2β)*a
        tv = v + (1-γ)*Δt*a
        rhs = rhs - C*tv - K*td
        rhs[[bdnode; bdnode.+(m+1)*(n+1)]] .= 0.0

        a = L\rhs
        d = td + β*Δt^2*a
        v = tv + γ*Δt*a
        U[:,i+1] = d
    end

    x = []
    y = []
    for j= 1:n+1
        for i = 1:m+1
            push!(x, (i-1)*h)
            push!(y, (j-1)*h)
        end
    end
    for i = 1:10:NT+1
        close("all")
        scatter(x+U[1:(m+1)*(n+1), i], y+U[(m+1)*(n+1)+1:end, i])
        xlabel("x")
        ylabel("y")
        k = string(i)
        k = repeat("0", 3-length(k))*k
        title("t = $k")
        ylim(-0.05,0.25)
        gca().invert_yaxis()
        savefig("u$k.png")
    end
end

@testset "tpfa" begin

    n = 50
    m = 2n
    h = 1/n
    X = zeros((m+1)*(n+1))
    Y = zeros((m+1)*(n+1))
    for i = 1:m+1
        for j = 1:n+1
            k = i+(j-1)*(m+1)
            X[k] = (i-1)*h
            Y[k] = (j-1)*h
        end
    end
    bd = bcedge("all", m, n, h)
    bdval = Float64[]
    for i = 1:size(bd,1)
        x = X[bd[i,1]]; y = Y[bd[i,1]]
        p1 = x^2+y^2
        x = X[bd[i,2]]; y = Y[bd[i,2]]
        p2 = x^2+y^2
        push!(bdval, (p1+p2)/2)
    end

    A, rhs = compute_fvm_tpfa_matrix(ones(m*n), bd, bdval, m, n, h)
    f = compute_fvm_source_term(4*ones(4*m*n), m, n, h)
    f -= rhs
    sol = A\f
    figure(figsize=(4,10))
    subplot(211)
    pcolormesh(LinRange(h,2,m)|>collect, LinRange(h,1,n)|>collect, reshape(sol, m, n)',
        vmin = -1.0, vmax = 3.5)
    gca().invert_yaxis()
    axis("equal")
    colorbar()
    subplot(212)
    x0 = LinRange(h,2,m)|>collect
    y0 = LinRange(h,1,n)|>collect
    X, Y = np.meshgrid(x0, y0)
    Z = @. X^2+Y^2
    pcolormesh(x0, y0, Z,
        vmin = -1.0, vmax = 3.5)
    gca().invert_yaxis()
    axis("equal")
    colorbar()

    figure()
    pcolormesh(log.(abs.(reshape(sol, m, n)'-Z))/log(10))
    colorbar()
end


text:
struct PropertyArray{T, N, name, A <: AbstractArray{<: Any, N}} <: AbstractArray{T, N}
    parent::A
end

PropertyArray{T, N, name}(parent::A) where {T, N, name, A} = PropertyArray{T, N, name, A}(parent)

Base.size(x::PropertyArray) = size(x.parent)
Base.axes(x::PropertyArray) = axes(x.parent)
Base.parent(x::PropertyArray) = x.parent
rootnode(x::PropertyArray) = rootnode(parent(x))
Base.IndexStyle(::Type{<: PropertyArray{<: Any, <: Any, <: Any, A}}) where {A} = IndexStyle(A)

@inline function Base.getindex(x::PropertyArray{<: Any, <: Any, name}, i::Int...) where {name}
    @boundscheck checkbounds(x, i...)
    @inbounds getproperty(x.parent[i...], name)
end

@inline function unsafe_getindex(x::PropertyArray{<: Any, <: Any, name}, i::Int...) where {name}
    @boundscheck checkbounds(x, i...)
    @inbounds getproperty(unsafe_getindex(x.parent, i...), name)
end

@inline function Base.setindex!(x::PropertyArray{<: Any, <: Any, name}, v, i::Int...) where {name}
    @boundscheck checkbounds(x, i...)
    @inbounds allocated = allocate!(x.parent, i...)
    @inbounds set!(allocated, name, v)
    x
end

@inline function isactive(x::PropertyArray, i::Int...)
    @boundscheck checkbounds(x, i...)
    @inbounds isactive(x.parent, i...)
end


text:
using Test, EqualizerFilters, DSP

@testset "sampling_rate setting" begin
    @test EqualizerFilters.sampling_rate() == 48e3 #test default
    @test EqualizerFilters.sampling_rate(1.) == 1.
    @test EqualizerFilters.sampling_rate() == 1.
    @test EqualizerFilters.sampling_rate(48e3) == 48e3
    @test EqualizerFilters.sampling_rate() == 48e3
end

fTest = 1e3
Qtest = .7
dbGainTest = -10

@testset "individual filter functions" begin
    #missing validation measurements from real EQ APO
    @test_nowarn LP(fTest, Qtest)
    @test_nowarn HP(fTest, Qtest)
    @test_nowarn BP(fTest, Qtest)
    @test_nowarn NO(fTest, Qtest)
    @test_nowarn AP(fTest, Qtest)
    @test_nowarn PK(fTest, dbGainTest, Qtest)
    @test_nowarn LS(fTest, dbGainTest, Qtest)
    @test_nowarn HS(fTest, dbGainTest, Qtest)
end

@testset "Biquad(Symbol, ...)" begin
    @test Biquad(:LP, fTest, Qtest) == LP(fTest, Qtest)
    @test Biquad(:HP, fTest, Qtest) == HP(fTest, Qtest)
    @test Biquad(:BP, fTest, Qtest) == BP(fTest, Qtest)
    @test Biquad(:NO, fTest, Qtest) == NO(fTest, Qtest)
    @test Biquad(:AP, fTest, Qtest) == AP(fTest, Qtest)
    @test Biquad(:PK, fTest, dbGainTest, Qtest) == PK(fTest, dbGainTest, Qtest)
    @test Biquad(:LS, fTest, dbGainTest, Qtest) == LS(fTest, dbGainTest, Qtest)
    @test Biquad(:HS, fTest, dbGainTest, Qtest) == HS(fTest, dbGainTest, Qtest)
end

@testset "Biquad(Tuple-Format)" begin
    @test Biquad((:LP, fTest, Qtest)) == LP(fTest, Qtest)
    @test Biquad((:HP, fTest, Qtest)) == HP(fTest, Qtest)
    @test Biquad((:BP, fTest, Qtest)) == BP(fTest, Qtest)
    @test Biquad((:NO, fTest, Qtest)) == NO(fTest, Qtest)
    @test Biquad((:AP, fTest, Qtest)) == AP(fTest, Qtest)
    @test Biquad((:PK, fTest, dbGainTest, Qtest)) == PK(fTest, dbGainTest, Qtest)
    @test Biquad((:LS, fTest, dbGainTest, Qtest)) == LS(fTest, dbGainTest, Qtest)
    @test Biquad((:HS, fTest, dbGainTest, Qtest)) == HS(fTest, dbGainTest, Qtest)
end

function sos_approx(x::SecondOrderSections,y::SecondOrderSections)
    if !(x.g ≈ y.g)
        return false
    elseif !(length(x.biquads) == length(y.biquads))
        return false
    else
        for i in eachindex(x.biquads)
            if biquad_approx(x.biquads[i],y.biquads[i]) == false
                return false
            end
        end
    end
    return true
end

function biquad_approx(x::Biquad, y::Biquad)
    x.b0 ≈ y.b0 &&
    x.b1 ≈ y.b1 &&
    x.b2 ≈ y.b2 &&
    x.a1 ≈ y.a1 &&
    x.a2 ≈ y.a2
end

@testset "SecondOrderSections(Array{Tuple-Format})" begin
    @test sos_approx(
    SecondOrderSections([(:LP, fTest, Qtest)]),
    LP(fTest, Qtest)|>SecondOrderSections
    )

    testArray = [(:LP, fTest, Qtest), (:HP, fTest, Qtest), (:BP, fTest, Qtest),
    (:NO, fTest, Qtest), (:AP, fTest, Qtest), (:PK, fTest, dbGainTest, Qtest),
    (:LS, fTest, dbGainTest, Qtest), (:HS, fTest, dbGainTest, Qtest)]
    @test sos_approx(
    SecondOrderSections(testArray),
    LP(fTest, Qtest)*HP(fTest, Qtest)*BP(fTest, Qtest)*
    NO(fTest, Qtest)*AP(fTest, Qtest)*PK(fTest, dbGainTest, Qtest)*
    LS(fTest, dbGainTest, Qtest)*HS(fTest, dbGainTest, Qtest)
    )
end

@testset "eqAPOstring(Symbol, ...)" begin
    @test_nowarn eqAPOstring(:LP, fTest, Qtest)
    @test_nowarn eqAPOstring(:HP, fTest, Qtest)
    @test_nowarn eqAPOstring(:BP, fTest, Qtest)
    @test_nowarn eqAPOstring(:NO, fTest, Qtest)
    @test_nowarn eqAPOstring(:AP, fTest, Qtest)
    @test_nowarn eqAPOstring(:PK, fTest, dbGainTest, Qtest)
    @test_nowarn eqAPOstring(:LS, fTest, dbGainTest, Qtest)
    @test_nowarn eqAPOstring(:HS, fTest, dbGainTest, Qtest)
end

@testset "eqAPOstring(Tuple-Format)" begin
    @test_nowarn eqAPOstring((:LP, fTest, Qtest))
    @test_nowarn eqAPOstring((:HP, fTest, Qtest))
    @test_nowarn eqAPOstring((:BP, fTest, Qtest))
    @test_nowarn eqAPOstring((:NO, fTest, Qtest))
    @test_nowarn eqAPOstring((:AP, fTest, Qtest))
    @test_nowarn eqAPOstring((:PK, fTest, dbGainTest, Qtest))
    @test_nowarn eqAPOstring((:LS, fTest, dbGainTest, Qtest))
    @test_nowarn eqAPOstring((:HS, fTest, dbGainTest, Qtest))
end

@testset "eqAPOstring(Array{Tuple-Format})" begin
    @test_nowarn eqAPOstring([(:LP, fTest, Qtest)])

    testArray = [(:LP, fTest, Qtest), (:HP, fTest, Qtest), (:BP, fTest, Qtest),
    (:NO, fTest, Qtest), (:AP, fTest, Qtest), (:PK, fTest, dbGainTest, Qtest),
    (:LS, fTest, dbGainTest, Qtest), (:HS, fTest, dbGainTest, Qtest)]
    @test_nowarn eqAPOstring(testArray)
end

@testset "eqAPOstring(x::FilterCoefficients)" begin
    @test_nowarn eqAPOstring(Biquad(1,2,3,4,5))
    @test_nowarn eqAPOstring(LP(fTest, Qtest))
    @test_nowarn eqAPOstring(PK(fTest, dbGainTest, Qtest)*AP(fTest, Qtest))
end


text:
module HashFunctions

export FNV1, FNV1A, calc

include("fnv.jl")

end


text:
using FluxArchitectures
using Documenter

DocMeta.setdocmeta!(FluxArchitectures, :DocTestSetup, :(using FluxArchitectures); recursive=true)

makedocs(;
    modules=[FluxArchitectures],
    authors="Sören Dobberschütz and contributors",
    repo="https://github.com/sdobber/FluxArchitectures.jl/blob/{commit}{path}#{line}",
    sitename="FluxArchitectures.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://sdobber.github.io/FluxArchitectures.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Examples" => "examples/examples.md",
        "Datasets" => "datasets/datasets.md",
        "Exported Functions" => "functions.md",
        "Models" =>
                    ["DARNN" => "models/darnn.md",
                     "DSANet" => "models/dsanet.md",
                     "LSTnet" => "models/lstnet.md",
                     "TPALSTM" => "models/tpalstm.md"],
        "Benchmarks" => "benchmark.md",
        "Reference" => "reference.md",
    ],
)

deploydocs(;
    repo="github.com/sdobber/FluxArchitectures.jl",
)

max_stars_repo_name: sdobber/FluxArchitectures

max_issues_repo_name: sdobber/FluxArchitectures

max_forks_repo_name: sdobber/FluxArchitectures


text:
"""
Provides unofficial implementations (in the form of `ODESystem` and
`ODEFunction` instances) of a polynomial approximation for longitudinal
aircraft dynamics. Specifically, these models approximate NASA's
Generic Transport Model – a radio-controlled, sub-scale model aircraft which
is used flight control research. These publicly available equations were published by
[Chakraborty et al](https://www.sciencedirect.com/science/article/abs/pii/S0967066110002595).

# Extended help

## License

$(LICENSE)

## Exports

$(EXPORTS)

## Imports

$(IMPORTS)
"""
module PolynomialGTM

using Memoize
using Symbolics
using LinearAlgebra
using ModelingToolkit
using DocStringExtensions

@template (FUNCTIONS, METHODS, MACROS) =
    """
    $(SIGNATURES)

    $(DOCSTRING)
    """

@template (TYPES, CONSTANTS) =
    """
    $(TYPEDEF)

    $(DOCSTRING)
    """

export GTM, GTMFunction

"""
NASA's Generic Transport Model can be approximated
(near select flight conditions) as a polynomial model.
Here, `GTM` is a `ModelingToolkit.ODESystem` which
provides publicly available polynomial approximations
for GTM longitudinal flight dynamics.

# Extended Help

## Initial Conditions

The default initial conditions are one trim condition.
Two trim conditions for these polynomial-approximated
dynamics are shown below.

```julia
trim₁ = [[29.6, deg2rad(9), 0.0, deg2rad(9)],   [deg2rad(0.68), 12.7]]
trim₂ = [[25.0, deg2rad(18), 0.0, deg2rad(18)], [deg2rad(-7.2), 59]]
```

## Usage

```julia
model = GTM()
```
## References:

- [Chakraborty et al](https://www.sciencedirect.com/science/article/abs/pii/S0967066110002595)
- [Joe Carpinelli](https://github.com/cadojo/Replicated-ROA-Analysis)
"""
@memoize function GTM(; stm=false, structural_simplify=true, name=:GTM)

    # First, let's define the necessary parameters
    # and variables for modeling the longitudinal
    # dynamics for NASA's Generic Transport Model

    @parameters t δₑ δₜ
    @variables V(t) α(t) q(t) θ(t)
    δ = Differential(t)
    x = [V, α, q, θ]
    p = [δₑ, δₜ]

    # Great! Now, we'll need to hard-code in the polynomial
    # approximations the longitudinal dynamics.

    eqs = [

        # Damn Unicode and its variable character font width!!!

        δ(V) ~  1.233e-8*V^4*q^2 + 4.853e-9*α^3*δₜ^3 	  	+
                3.705e-5*V^3*α*q - 2.184e-6*V^3*q^2 		+
                2.203e-2*V^2*α^3 - 2.836e-6*α^3*δₜ^2 	  	+
                3.885e-7*α^2*δₜ^3 - 1.069e-6*V^3*q 	  		-
                4.517e-2*V^2*α^2 - 2.140e-3*V^2*α*δₑ 		-
                3.282e-3*V^2*α*q - 8.901e-4*V^2*δₑ^2 		+
                9.677e-5*V^2*q^2 - 2.037e-4*α^3*δₜ 	  	 	 -
                2.270e-4*α^2*δₜ^2 - 2.912e-8*α*δₜ^3 	     +
                1.591e-3*V^2*α - 4.077e-4*V^2*δₑ 		  	+
                9.475e-5*V^2*q - 1.637*α^3 				  	-
                1.631e-2*α^2*δₜ + 4.903*α^2*θ 		  		-
                4.903*α*θ^2 + 1.702e-5*α*δₜ^2 		  		-
                7.771e-7*δₜ^3 + 1.634*θ^3 - 4.319e-4*V^2    -
                2.142e-1*α^2 + 1.222e-3*α*δₜ 			  	+
                4.541e-4*δₜ^2 + 9.823*α + 3.261e-2*δₜ 	    -
                9.807*θ + 4.282e-1,
        δ(α) ~ -3.709e-11*V^5*q^2 + 6.869e-11*V*α^3*δₜ^3 	+
                7.957e-10*V^4*α*q + 9.860e-9*V^4*q^2 		+
                1.694e-5*V^3*α^3 - 4.015e-8*V*α^3*δₜ^2 		 -
                7.722e-12*V*α^2*δₜ^3 - 6.086e-9*α^3*δₜ^3 	 -
                2.013e-8*V^4*q - 5.180e-5*V^3*α^2 			-
                2.720e-6*V^3*α*δₑ - 1.410e-7*V^3*α*q 		+
                7.352e-7*V^3*δₑ^2 - 8.736e-7*V^3*q^2 		-
                1.501e-3*V^2*α^3 - 2.883e-6*V*α^3*δₜ 		+
                4.513e-9*V*α^2*δₜ^2 - 4.121e-10*V*α*δₜ^3 	+
                3.557e-6*α^3*δₜ^2 + 6.841e-10*α^2*δₜ^3 		+
                4.151e-5*V^3*α + 3.648e-6*V^3*δₑ 		    +
                3.566e-6*V^3*q + 6.246e-6*V^2*α*q 			+
                4.589e-3*V^2*α^2 + 2.410e-74*V^2*α*δₑ 		-
                6.514e-5*V^2*δₑ^2 + 2.580e-5*V^2*q^2 		-
                3.787e-5*V*α^3 + 3.241e-7*V*α^2*δₜ 			+
                2.409e-7*V*α*δₜ^2 + 1.544e-11*V*δₜ^3 		+
                2.554e-4*α^3*δₜ - 3.998e-7*α^2*δₜ^2 		+
                3.651e-8*α*δₜ^3 + 4.716e-7*V^3 				-
                3.677e-3*V^2*α - 3.231e-4*V^2*δₑ 			-
                1.579e-4*V^2*q + 2.605e-3*V*α^2 			+
                1.730e-5*V*α*δₜ - 5.201e-3*V*α*θ 			-
                9.026e-9*V*δₜ^2 + 2.601e-3*V*θ^2 			+
                3.355e-3*α^3 - 2.872e-5*α^2*δₜ 				-
                2.134e-5*α*δₜ^2 - 1.368e-9*δₜ^3 		    -
                4.178e-5*V^2 + 2.272e-4*V*α 				-
                6.483e-7*V*δₜ - 2.308e-1*α^2 				-
                1.532e-3*α*δₜ + 4.608e-1*α*θ - 2.304e-1*θ^2 +
                7.997e-7*δₜ^2 - 5.210e-3*V - 2.013e-2*α 	+
                5.744e-5*δₜ + q + 4.616e-1,
        δ(q) ~ -6.573e-9*V^5*q^3 + 1.747e-6*V^4*q^3 	    -
                1.548e-4*V^3*q^3 - 3.569e-3*V^2*α^3 	    +
                4.571e-3*V^2*q^3 + 4.953e-5*V^3*q 		    +
                9.596e-3*V^2*α^2 + 2.049e-2*V^2*α*δₑ 	    -
                2.431e-2*V^2*α - 3.063e-2*V^2*δₑ 		    -
                4.388e-3*V^2*q - 2.594e-7*δₜ^3 			    +
                2.461e-3*V^2 + 1.516e-4*δₜ^2 + 1.089e-2*δₜ 	 +
                1.430e-1,
        δ(θ) ~ q
    ]

    # Let's set some default values, so we can plug this system directly
    # into `ODEProblem` and perform quick analysis. Note that what we're
    # doing here is setting one equilibrium position as a default
    # initial condition if no arguments are provided to `ODEProblem`!
    defaults = Dict(
        V  => 29.6,
        α  => deg2rad(9),
        q  => 0.0,
        θ  => deg2rad(0),
        δₑ => deg2rad(0.68),
        δₜ => 12.7
    )

    # If state transition matrix dynamics are enabled, append
    # the dynamics to our equations of motion, append the
    # state variables to our state vector, and append default
    # values (the identity matrix) to the defaults field!
    if stm
        @variables Φ[1:4,1:4](t)
        Φ = Symbolics.scalarize(Φ)
        A = Symbolics.jacobian(map(el -> el.rhs, eqs), x)

        LHS = map(δ, Φ)
        RHS = map(simplify, A * Φ)

        eqs = vcat(eqs, [LHS[i] ~ RHS[i] for i in 1:length(LHS)])

        for ϕ in vec(Φ)
            push!(x, ϕ)
        end

        for (ϕ, i) in zip(vec(Φ), vec(Matrix(I(4))))
            defaults[ϕ] = i
        end
    end

    # Model name
    if string(name) == "GTM" && stm
        modelname = Symbol("GTMWithSTM")
    else
        modelname = name
    end

    # Make the model and return
    model = ODESystem(eqs, t, x, p; defaults=defaults, name=Symbol(modelname))
    return structural_simplify ? ModelingToolkit.structural_simplify(model) : model
end

"""
Returns a `DifferentialEquations`-compatible `ODEFunction` for GTM dynamics.
The `stm`, `structural_simplify`, and `name` keyword arguments are passed
to `GTM`. All other keyword arguments are passed directly to `ODEFunction`.

# Extended Help

## Usage

Note that this `ODEFunction` output has several methods, including an in-place
method! Function signatures follow `ModelingToolkit` and `DifferentialEquations`
conventions.

```julia
f = GTMFunction()
let u = randn(4), p = randn(2), t = rand()
    f(u,p,t)
end
```
"""
@memoize function GTMFunction(; stm=false, structural_simplify=true, name=:GTM, kwargs...)
    defaults = (; jac = true)
    options  = merge(defaults, kwargs)

    model = GTM(; stm=stm, structural_simplify=structural_simplify, name=name)
    return ODEFunction(model; options...)
end

end # module

max_stars_repo_name: cadojo/PolynomialGTMModel

max_issues_repo_name: cadojo/PolynomialGTMModel

max_forks_repo_name: cadojo/PolynomialGTMModel


text:
tests = ["lreg",
  		 "whiten",
         "pca",
         "cca",
         "cmds",
         "lda",
         "mclda",
         "ica"]

println("Running tests:")

for t in tests
    fp = string(t, ".jl")
    println(" * $(fp)")
    include(fp)
end


text:
function read_hdrgen_file(path::String)
# description: read hdrgen file and parse the lines as (image, shutter_time)
#              pair.
# input: path to hdrgen file
# return value: vector of the (image, shutter_time) pairs.
    hdrgen_file = open(path)
    hdrgen_lines = readlines(hdrgen_file)
    base_dir = dirname(path)

    #Enhancement: take general image types instead of only RGB{N0f8}
    hdr_pairs = HDRSequence()

    for line in hdrgen_lines
        image_name, exposure = split(line)
        image_path = joinpath(base_dir, image_name)
        println("reading image ", image_path)
        image = load(image_path)

        push!(hdr_pairs, (image, 1.0 / float(exposure)))
    end

    return hdr_pairs
end


text:
@testset "BoolVar.jl" begin
    @testset "isbound()" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.BoolVar("x", trailer)

        @test !SeaPearl.isbound(x)

        SeaPearl.assign!(x, false)

        @test SeaPearl.isbound(x)
    end

    @testset "assignedValue()" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.BoolVar("x", trailer)

        @test_throws AssertionError SeaPearl.assignedValue(x)

        SeaPearl.assign!(x, false)

        @test !SeaPearl.assignedValue(x)
    end

    @testset "BoolVar()" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.BoolVar("x", trailer)

        @test length(x.domain) == 2
        @test true in x.domain && false in x.domain
    end
end


text:
module Preconditioners

using JuLIP: AbstractAtoms, Preconditioner, JVecs, JVecsF, Dofs, maxdist,
            constraint, pairs, cutoff, positions, defm, JVecF, forces, mat, vecs,
            set_positions!, julipwarn, chemical_symbols, rnn, JVec, JMat ,
            AbstractCalculator, calculator

using JuLIP.Potentials: @analytic, evaluate, evaluate_d, PairPotential, HS,
         SitePotential, sites, C0Shift, _precon_or_hessian_pos

using JuLIP.Constraints: project!, FixedCell, _pos_to_alldof

using AMG

import JuLIP: update!
import JuLIP.Potentials: precon, cutoff
import Base: A_ldiv_B!, A_mul_B!, dot, *, \, size

export Exp, FF

"""
`IPPrecon`: the standard preconditioner type for JuLIP

`IPPrecon` stores a field `p <: AbstractCalculator` which is
should implement `precon(p, r, R)` which is a block matrix analogue of
`hess(p, r, R)` and is used to assemble the preconditioner matrix via
```julia
precon_matrix(P, at::AbstractAtoms)
```
If `p` should be updated in response to an update of `at` then
the type can overload `JuLIP.Preconditioners.update_inner!`.

## Constructor:
```
IPPrecon(p::AbstractCalculator, at::AbstractAtoms; kwargs...)
```

## Keyword arguments:

* `updatedist`: determines after how much movement of the atoms, the
preconditioner is updated
* `tol`: solver tolerance if AMG is used
* `updatefreq`: after how many iterations do we force an update {10}
* `stab`: stabilisation constant {0.01}, add `stab * I` to `P`
* `solve`: which solver to used, `:amg` or `:chol`, default is `:amg`
"""
type IPPrecon{TV, TS, T <: AbstractFloat, TI <: Integer} <: Preconditioner
   p::TV
   solver::TS
   A::SparseMatrixCSC{T, TI}
   oldX::Vector{JVec{T}}
   updatedist::T
   tol::T
   updatefreq::TI
   skippedupdates::TI
   stab::T
end

function IPPrecon(p::AbstractCalculator, at::AbstractAtoms;
         updatedist=0.2 * rnn(at), tol=1e-7, updatefreq=10, stab=0.01,
         solver = :chol)
   # make sure we don't use this in a context it is not intended for!
   @assert isa(constraint(at), FixedCell)
   A = AMG.poisson(12)
   if solver == :amg
      solver = ruge_stuben(A)
   elseif solver in [:chol, :direct]
      solver = cholfact(A)
   else
      error("`IPPrecon` : unknown solver $(solver)")
   end
   P = IPPrecon(p, solver, A, positions(at), updatedist, tol, updatefreq, 0, stab)
   # the force_update! makes the first assembly pass.
   return force_update!(P, at)
end

# some standard Base functionality lifted to IPPrecon
A_ldiv_B!(out::Dofs, P::IPPrecon{TV, TS}, f::Dofs) where TV where TS <: AMG.MultiLevel =
   copy!(out, solve(P.solver, f, 200, AMG.V(), P.tol))
A_ldiv_B!(out::Dofs, P::IPPrecon{TV, TS}, f::Dofs) where TV where TS <: Base.SparseArrays.CHOLMOD.Factor =
   copy!(out, P.solver \ f)
A_mul_B!(out::Dofs, P::IPPrecon, x::Dofs) = A_mul_B!(out, P.A, x)
dot(x, P::IPPrecon, y) = dot(x, P * y)
*(P::IPPrecon, x::AbstractVector) = P.A * x
\(P::IPPrecon, x::AbstractVector) = A_ldiv_B!(zeros(length(x)), P, x)
Base.size(P::IPPrecon) = size(P.A)

need_update(P::IPPrecon, at::AbstractAtoms) =
   (P.skippedupdates > P.updatefreq) ||
   (maxdist(positions(at), P.oldX) >= P.updatedist)

update!(P::IPPrecon, at::AbstractAtoms) =
   need_update(P, at) ? force_update!(P, at) : (P.skippedupdates += 1; P)

update_solver!(P::IPPrecon{TV, TS}) where TV where TS <: Base.SparseArrays.CHOLMOD.Factor =
   cholfact(Symmetric(P.A))

update_solver!(P::IPPrecon{TV, TS}) where TV where TS <: AMG.MultiLevel =
   ruge_stuben(P.A)

function force_update!(P::IPPrecon, at::AbstractAtoms)
   # perform updates of the potential p (if needed; usually not)
   P.p = update_inner!(P.p, at)
   # construct the preconditioner matrix ...
   Pmat = precon_matrix(P.p, at)
   Pmat = Pmat + P.stab * speye(size(Pmat, 1))
   Pmat = project!(constraint(at), Pmat)
   # and the AMG solver
   P.A = Pmat
   P.solver = update_solver!(P)
   # remember the atom positions
   copy!(P.oldX, positions(at))
   # and remember that we just did a full update
   P.skippedupdates = 0
   return P
end

# ============== Implementation of some concrete types ======================

# default implementation of update_inner!;
# most of the time no inner update is required

"update `IPPrecon.p`, must return updated `p`"
update_inner!(p::Any, at::AbstractAtoms) = p

"return the three linear indices associated with an atom index"
atind2lininds(i::Integer) = (i-1) * 3 + [1;2;3]

"""
build the preconditioner matrix associated with the potential V
"""
precon_matrix(V::AbstractCalculator, at::AbstractAtoms) =
   _pos_to_alldof(_precon_or_hessian_pos(V, at, precon), at)

"""
A variant of the `Exp` preconditioner; see

### Constructor: `Exp(at::AbstractAtoms)`

Keyword arguments:

* `A`: stiffness of potential {default 3.0}
* `r0`: nn distance estimate {default is an automicatic estimate}
* `cutoff_mult`: cut-off is `r0 * cutoff_mult` {default 2.2}
* `kwargs...`: `IPPrecon` parameters

### Reference

      D. Packwood, J. Kermode, L. Mones, N. Bernstein, J. Woolley, N. I. M. Gould,
      C. Ortner, and G. Csanyi. A universal preconditioner for simulating condensed
      phase materials. J. Chem. Phys., 144, 2016.
"""
mutable struct Exp{T, TV} <: PairPotential
   Vexp::TV
   energyscale::T
end

cutoff(P::Exp) = cutoff(P.Vexp)

precon{T}(P::Exp{T}, r, R) = (P.energyscale * P.Vexp(r)) * one(JMat{T})

function Exp(at::AbstractAtoms;
             A=3.0, r0=rnn(at), cutoff_mult=2.2, energyscale = 1.0,
             kwargs...)
   e0 = energyscale == :auto ? 1.0 : energyscale
   rcut = r0 * cutoff_mult
   Vexp = let A=A, r0=r0
      @analytic r -> exp( - A * (r/r0 - 1.0))
   end * C0Shift(rcut)
   P = IPPrecon(Exp(Vexp, e0), at; kwargs...)
   if energyscale == :auto
      P.p.energyscale = estimate_energyscale(at, P)
      force_update!(P, at)
   end
   return P
end

# want μ * <P v, v> ~ <∇E(x+hv) - ∇E(x), v> / h
function estimate_energyscale(at, P)
   # get the P-matrix at current configuration
   A = precon_matrix(P.p, at)
   # determine direction in which the cell is maximal
   F = Matrix(defm(at))
   X0 = positions(at)
   _, i = findmax( norm(F[:, j]) for j = 1:3 )
   Fi = JVecF(F[:, i])
   # an associated perturbation
   V = [sin(2*pi * dot(Fi, x)/dot(Fi,Fi)) * Fi for x in X0]
   # compute gradient at current positions
   f0 = forces(at) |> mat
   # compute gradient at perturbed positions
   h = 1e-3
   set_positions!(at, X0 + h * V)
   f1 = forces(at) |> mat
   # return the estimated value for the energyscale
   V = mat(V)
   μ = - vecdot((f1 - f0)/h, V) / dot(A * V[:], V[:])
   if μ < 0.1 || μ > 100.0
      warn("""
      e0 = $(μ) in `estimate_energyscale`; this is likely due to a poor
      starting guess in an optimisation, probably best to set the
      energyscale manually, or use CG instead of LBFGS.
      """)
   end
   return μ
end

"""
`FF`: defines a preconditioner based on a force-field;

TODO: thorough documentation and reference once the paper is finished
"""
FF(at::AbstractAtoms, V::AbstractCalculator; kwargs...) =
   IPPrecon(V, at; kwargs...)

FF(at::AbstractAtoms) = FF(at, calculator(at))

end # end module Preconditioners


text:
module ThingSpeak
using Requests
immutable TSChannel
  writekey::ASCIIString
end
const updateapiurl = "https://api.thingspeak.com/update"
const updateheaders = {"Content-type"=> "application/x-www-form-urlencoded","Accept"=> "text/plain"}
defaultchannel = TSChannel("")
setdefaultchannel(s::ASCIIString) = (global defaultchannel = TSChannel(s))

function updatedata(channel::TSChannel, status::ASCIIString, fields::Dict)
  data = "key=$(channel.writekey)"
  for (i,f) in fields
    data *= "&field$i="*@sprintf("%f",f)
  end
  if status != ""
    length(status)>140 && error("status is limited to 140 characters in ThingSpeak API")
    data *= "&status=$status"
  end
  data
end
function tsupdate(fields::Dict;status::ASCIIString="", channel::TSChannel=defaultchannel)
  response = post(updateapiurl, data=updatedata(channel, status, fields), headers=updateheaders)
end
function tsupdate(fields...;status::ASCIIString="", channel::TSChannel=defaultchannel)
  response = post(updateapiurl, data=updatedata(channel, status, {k=>v for (k,v) in enumerate(fields)}), headers=updateheaders)
end

tssuccess(r::Requests.Response) = r.data != "0"

export tsupdate, tssuccess
end # module


text:
#!/usr/bin/env julia
using ArgParse

s = ArgParseSettings()

@add_arg_table s begin
    "--geometry"
    "--N"
    "--d"
    "--gamma"
        default = "1.0"
    "--edipole"
    "--phi"
    "--theta"
    "--o"
    "--T"
end

parameters = parse_args(s)

# Output
const odir = parameters["o"]
@assert isdir(odir)

# Integration
const T = float(eval(Meta.parse(parameters["T"])))
@assert T[1]<T[end]

# System parameters
const γ = float(parameters["gamma"])
@assert 0<γ

# Initial state (Bloch state)
const phi = float(parameters["phi"])
const theta = float(parameters["theta"])

# System geometry
using QuantumOptics, CollectiveSpins
const cs = CollectiveSpins

const edipole = float(eval(Meta.parse(parameters["edipole"])))
const geomN = int(parameters["N"])
const d = float(parameters["d"])
const geomstring = parameters["geometry"]
if geomstring=="chain"
    const systemgeometry = cs.geometry.chain(d, geomN)
elseif geomstring=="square"
    const systemgeometry = cs.geometry.square(d; Nx=geomN, Ny=geomN)
elseif geomstring=="cube"
    const systemgeometry = cs.geometry.cube(d)
else
    error("Unknown geometry: $geomstring")
end

const system = SpinCollection(systemgeometry, edipole, γ)
const N = length(system.spins)

const td_ind = Float64[]
const td_mf = Float64[]
const td_mpc = Float64[]

# Independent
state0 = cs.independent.blochstate(phi,theta,N)
tout, state_ind_t = cs.independent.timeevolution(T, system, state0)

# Meanfield
state0 = cs.meanfield.blochstate(phi,theta,N)
tout, state_mf_t = cs.meanfield.timeevolution(T, system, state0)

# Meanfield + Correlations
state0 = cs.mpc.blochstate(phi,theta,N)
tout, state_mpc_t = cs.mpc.timeevolution(T, system, state0)

# Quantum: master equation

function fout(t, rho::Operator)
    i = findfirst(T, t)
    rho_ind = cs.independent.densityoperator(state_ind_t[i])
    rho_mf = cs.meanfield.densityoperator(state_mf_t[i])
    rho_mpc = cs.mpc.densityoperator(state_mpc_t[i])
    push!(td_ind, tracedistance(rho, rho_ind))
    push!(td_mf, tracedistance(rho, rho_mf))
    push!(td_mpc, tracedistance(rho, rho_mpc))
end

Ψ₀ = cs.quantum.blochstate(phi,theta,N)
ρ₀ = Ψ₀⊗dagger(Ψ₀)
cs.quantum.timeevolution(T, system, ρ₀, fout=fout)

keyparameters = Dict(
    "theta"=>parameters["theta"],
    "edipole"=>parameters["edipole"],
    "d"=>parameters["d"]
    )

name = QuantumOptics.io.dict2filename(keyparameters)
f = open(joinpath(odir, name), "w")
cs.io.write_head(f, system, parameters)
cs.io.write_state(f, "timevector", T; time=false)
cs.io.write_state(f, "td_ind", td_ind; time=false)
cs.io.write_state(f, "td_mf", td_mf; time=false)
cs.io.write_state(f, "td_mpc", td_mpc; time=false)
close(f)


text:
using Plots

function bracket_sign_change(f′, a, b; k=2)
    if a > b; a,b = b,a; end # ensure a < b
    center, half_width = (b+a)/2, (b-a)/2
    while f′(a)*f′(b) > 0
        half_width *= k
        a = center - half_width
        b = center + half_width
    end
    return (a,b)
end

function bisection(f′, a, b, ϵ)
    if a > b; a,b = b,a; end # ensure a < b

    ya, yb = f′(a), f′(b)
    if ya == 0; b = a; end
    if yb == 0; a = b; end

    while b - a > ϵ
        x = (a+b)/2
        y = f′(x)
        if y == 0
            a, b = x, x
        elseif sign(y) == sign(ya)
            a = x
        else
            b = x
        end
    end
    return (a,b)
end

#definicao da funcao
f = x-> -(0.4/((x^2+1))^(1/2))+((x^2+1)^(1/2))*(1-0.4/(x^2+1))-x    #multiplicada por -1 para encontrar o minimo
g = x-> -(-(0.4/((x^2+1))^(1/2))+((x^2+1)^(1/2))*(1-0.4/(x^2+1))-x) #funcao original

#definicao da primeira derivada
df(x) = 1.2*x/((x^2+1)^(3/2))+(1-0.4/(x^2+1))*x/((x^2+1)^(1/2))-1

#dois pontos iniciais
x0 = 3.4
x1 = 3.8

#convergencia
c_tol = 1.0e-8

as,bs = bracket_sign_change(df,x0,x1)
as,df(as),bs,df(bs)

af,bf = bisection(df,as,bs,c_tol)

println("x* = ", af)
println("df = ", -df(af))
println("f = ", g(af))

plot(g,as-1.,bs+1.,lw=2,draw_arrow="true",label="Function")
scatter!([as,bs],[g(as),g(bs)],label="interval")
scatter!([af,bf],[g(af),g(bf)],label="x*")