text:
using Printf
using Logging

using ArgParse
using LoggingExtras
using DataDeps
using Flux
using JLD2
using OrdinaryDiffEq

using Oceananigans
using OceanParameterizations
using FreeConvection

ENV["DATADEPS_ALWAYS_ACCEPT"] = "true"
ENV["GKSwstype"] = "100"

function parse_command_line_arguments()
    settings = ArgParseSettings()

    @add_arg_table! settings begin
        "--grid-points"
            help = "Number of vertical grid points in the trained neural differential equation (LES data will be coarse-grained to this resolution)."
            default = 32
            arg_type = Int

        "--base-parameterization"
            help = "Base parameterization to use for the NDE. Options: nothing, convective_adjustment"
            default = "nothing"
            arg_type = String

        "--time-stepper"
            help = "DifferentialEquations.jl time stepping algorithm to use."
            default = "ROCK4"
            arg_type = String

        "--name"
            help = "Experiment name (also determines name of output directory)."
            default = "layers3_depth4_relu_ROCK4"
            arg_type = String

        "--training-simulations"
            help = "Simulation IDs (list of integers separated by spaces) to train the neural differential equation on." *
                   "All other simulations will be used for testing/validation."
            action = :append_arg
            nargs = '+'
            arg_type = Int
            range_tester = (id -> id in FreeConvection.SIMULATION_IDS)
    end

    return parse_args(settings)
end

@info "Parsing command line arguments..."

args = parse_command_line_arguments()

nde_type = Dict(
    "nothing" => FreeConvectionNDE,
    "convective_adjustment" => ConvectiveAdjustmentNDE
)

Nz = args["grid-points"]
experiment_name = args["name"]
NDEType = nde_type[args["base-parameterization"]]
algorithm = Meta.parse(args["time-stepper"] * "()") |> eval

ids_train = args["training-simulations"][1]
ids_test = setdiff(FreeConvection.SIMULATION_IDS, ids_train)
validate_simulation_ids(ids_train, ids_test)

output_dir = joinpath(@__DIR__, experiment_name)
mkpath(output_dir)

@info "Planting loggers..."

log_filepath = joinpath(output_dir, "$(experiment_name)_testing.log")
TeeLogger(
    OceananigansLogger(),
    MinLevelLogger(FileLogger(log_filepath), Logging.Info)
) |> global_logger

@info "Loading training data..."

data = load_data(ids_train, ids_test, Nz)
coarse_datasets = data.coarse_datasets

@info "Reading neural network from disk..."

nn_history_filepath = joinpath(output_dir, "neural_network_history.jld2")
final_nn_filepath = joinpath(output_dir, "free_convection_final_neural_network.jld2")

final_nn = jldopen(final_nn_filepath, "r")
NN = final_nn["neural_network"]
T_scaling = final_nn["T_scaling"]
wT_scaling = final_nn["wT_scaling"]
close(final_nn)

function free_convection_neural_network(input)
    wT_interior = NN(input.temperature)
    wT = cat(input.bottom_flux, wT_interior, input.top_flux, dims=1)
    return wT
end

@info "Gathering and computing solutions..."

true_solutions = Dict(id => (T=interior(ds["T"])[1, 1, :, :], wT=interior(ds["wT"])[1, 1, :, :]) for (id, ds) in coarse_datasets)
nde_solutions = Dict(id => solve_nde(ds, NN, NDEType, algorithm, T_scaling, wT_scaling) for (id, ds) in coarse_datasets)
kpp_solutions = Dict(id => free_convection_kpp(ds) for (id, ds) in coarse_datasets)
tke_solutions = Dict(id => free_convection_tke_mass_flux(ds) for (id, ds) in coarse_datasets)

convective_adjustment_solutions = Dict()
oceananigans_solutions = Dict()
for (id, ds) in coarse_datasets
    ca_sol, nn_sol = oceananigans_convective_adjustment_nn(ds, output_dir=output_dir, nn_filepath=final_nn_filepath)
    convective_adjustment_solutions[id] = ca_sol
    oceananigans_solutions[id] = nn_sol
end

@info "Plotting loss matrix..."

plot_loss_matrix(coarse_datasets, ids_train, nde_solutions, kpp_solutions, tke_solutions,
                 convective_adjustment_solutions, oceananigans_solutions, T_scaling,
                 filepath = joinpath(output_dir, "loss_matrix_plots.png"))

for (id, ds) in coarse_datasets
    filepath = joinpath(output_dir, "free_convection_comparisons_$id")
    plot_comparisons(ds, id, ids_train, nde_solutions[id], kpp_solutions[id], tke_solutions[id],
                     convective_adjustment_solutions[id], oceananigans_solutions[id], T_scaling,
                     filepath = filepath, frameskip = 5)
end

@info "Animating what the neural network has learned..."

for (id, ds) in coarse_datasets
    filepath = joinpath(output_dir, "learned_free_convection_$id")
    animate_learned_free_convection(ds, NN, free_convection_neural_network, NDEType, algorithm, T_scaling, wT_scaling,
                                    filepath=filepath, frameskip=5)
end

@info "Computing NDE solution history..."

nde_solution_history = compute_nde_solution_history(coarse_datasets, NDEType, algorithm, final_nn_filepath, nn_history_filepath)

@info "Plotting loss(epoch)..."

plot_epoch_loss(ids_train, ids_test, nde_solution_history, true_solutions, T_scaling,
                title = "Free convection loss history",
                filepath = joinpath(output_dir, "free_convection_nde_loss_history.png"))

@info "Plotting loss(time; epoch)..."

animate_nde_loss(coarse_datasets, ids_train, ids_test, nde_solution_history, true_solutions, T_scaling,
                 title = "Free convection loss history",
                 filepath = joinpath(output_dir, "free_convection_nde_loss_evolution"))

#=
@info "Comparing advective fluxes ⟨w'T'⟩ with LES diffusive flux ⟨κₑ∂zT⟩..."

# t = coarse_datasets[1]["T"].times ./ 86400
# p = plot(xlabel="Time (days)", ylabel="|κₑ∂zT| / ( |w'T'| + |κₑ∂zT| )", xlims=extrema(t), grid=false, framestyle=:box,
#          legend=:outertopright, foreground_color_legend=nothing, background_color_legend=nothing, dpi=200)

# for (id, ds) in coarse_datasets
#     advective_heat_flux = sum(ds[:wT].data .|> abs, dims=1)[:]
#     diffusive_heat_flux = sum(ds[:κₑ_∂z_T].data .|> abs, dims=1)[:]
#     total_heat_flux = advective_heat_flux .+ diffusive_heat_flux
#     label = @sprintf("%d W/m²", -ds.metadata[:heat_flux_Wm⁻²])
#     plot!(p, t, diffusive_heat_flux ./ total_heat_flux, linewidth=2, label=label)
# end

# savefig(joinpath(output_dir, "les_flux_contribution.png"))
=#


text:
"""
    haversine(lon1, lat1, lon2, lat2)

Calculate the haversine distance and bearing. Distance is in nm.
"""
@inline @fastmath function haversine(lon1, lat1, lon2, lat2)
    R = 6372.8  # Earth radius in kilometers

    dLat = deg2rad(lat2 - lat1)
    dLon = deg2rad(lon2 - lon1)
    lat1 = deg2rad(lat1)
    lat2 = deg2rad(lat2)
    lon1 = deg2rad(lon1)
    lon2 = deg2rad(lon2)
    a = sin(dLat/2)^2 + cos(lat1)*cos(lat2)*sin(dLon/2)^2
    c = 2*asin(sqrt(a))
    theta = atan(sin(dLon)*cos(lat2),
                 cos(lat1)*sin(lat2)-sin(lat1)*cos(lat2)*cos(dLon))
    theta = (rad2deg(theta) + 360) % 360
    return R*c*0.5399565, theta
end

function rotate_point(ox, oy, px, py, angle,ret_x)
    """
    Rotate a point counterclockwise by a given angle around a given origin.

    The angle should be given in radians.
    """
    qx = ox + cos(angle) * (px - ox) - sin(angle) * (py - oy)
    qy = oy + sin(angle) * (px - ox) + cos(angle) * (py - oy)
    if ret_x == true
        return qx
    else
        return qy
    end
end

function generate_grid(start_lon, start_lat, finish_lon, finish_lat, nodes)
    dist = haversine(start_lon, start_lat, finish_lon, finish_lat)
    spacing = dist[1]/(nodes+1)
    alpha = dist[2]
    x_dist = spacing
    y_dist = spacing
    grid_x = reshape(start_lon.+[i*x_dist for i in range(1, length=nodes) for j in range(0, length=nodes).-(nodes-1)/2], (nodes, nodes))
    grid_y = reshape(start_lat.+[j*y_dist for i in range(1, length=nodes) for j in range(0, length=nodes).-(nodes-1)/2], (nodes, nodes))
    rot_grid_x = [rotate_point(start_lon, start_lat, x, y, alpha, true) for (x, y) in zip(grid_x,grid_y)]
    rot_grid_y = [rotate_point(start_lon, start_lat, x, y, alpha, false) for (x, y) in zip(grid_x,grid_y)]
    return transpose(rot_grid_x), rot_grid_y
end

"""
    euclidean(x1, y1, x2, y2)

Calculate the distance and argument between the vector and north between points 1 and 2.
euclidean(0.0, 0.0, 10.0, 10.0) = (14.142135623730951, 45.0)
euclidean(0.0, 0.0, 10.0, 10.0) = (14.142135623730951, 45.0)
euclidean(0.0, 0.0, 10.0, 10.0) = (14.142135623730951, 45.0)
euclidean(0.0, 0.0, 10.0, 10.0) = (14.142135623730951, 45.0)
"""
@inline @fastmath function euclidean(x1::Float64, y1::Float64, x2::Float64, y2::Float64)
    dx = x2 - x1
    dy = y2 - y1
    dist = (dx^2 + dy^2)^(0.5)
    theta = 0.0
    if dy > 0.0 && dx > 0.0
        theta = 90.0 - rad2deg(atan(dy, dx))
    elseif dy == 0 && dx > 0.0
        theta = 90.0
    elseif dx == 0 && dy < 0.0
        theta = 180.0
    elseif dy < 0.0 && dx > 0.0
        theta = 90.0 + rad2deg(atan(abs(dy), dx))
    elseif dy < 0.0 && dx < 0.0
        theta = 180 + rad2deg(atan(abs(dy), abs(dx)))
    elseif dy > 0.0 && dx < 0.0
        theta = 360 - rad2deg(atan(dy, abs(dx)))
    end
    return dist, theta
end

"""Calculate the number of nodes for a specific distance in nm."""
function calc_nodes(lon1, lon2, lat1, lat2, req_d)
    d = haversine(lon1, lon2, lat1, lat2)[1]
    return round(Int, d/req_d)
end

ext: jl


text:
# I don't want to use map anymore, it's so ambigious, especially to newcomers.
# TODO should this become it's own function?
"""
    lift(f, o1::Observables.AbstractObservable, rest...) -> w

Create a new `w::Observable` by applying `f` to the _values_ of all observables
in `o1` and `rest...`. The initial value of `w` is determined by the first function
evaluation. The observable `w` is updated by calling `f` each time any of the
observables `o1, rest...` are updated.

## Examples
```julia
julia> x = Observable(2)
Observable{Int64} with 0 listeners. Value:
2

julia> y = lift(a -> a^2, x)
Observable{Int64} with 0 listeners. Value:
4

julia> z = lift((a,b) -> a+b, x, y)
Observable{Int64} with 0 listeners. Value:
6
```
"""
function lift(f, o1::Observables.AbstractObservable, rest...; kw...)
    if !isempty(kw)
        error("lift(f, obs...; init=f.(obs...), typ=typeof(init)) is deprecated. Use lift(typ, f, obs...), or map!(f, Observable(), obs...) for different init.")
    end
    init = f(to_value(o1), to_value.(rest)...)
    typ = typeof(init)
    result = Observable{typ}(init)
    map!(f, result, o1, rest...)
    return result
end

function lift(
        f, ::Type{T}, o1::Observables.AbstractObservable, rest...
    ) where {T}
    init = f(to_value(o1), to_value.(rest)...)
    result = Observable{T}(init)
    map!(f, result, o1, rest...)
    return result
end

Base.close(obs::Observable) = empty!(obs.listeners)

"""
Observables.off but without throwing an error
"""
function safe_off(o::Observables.AbstractObservable, f)
    l = listeners(o)
    for i in 1:length(l)
        if f === l[i]
            deleteat!(l, i)
            for g in Observables.removehandler_callbacks
                g(o, f)
            end
            return
        end
    end
end

"""
    map_once(closure, inputs::Observable....)::Observable

Like Reactive.foreach, in the sense that it will be preserved even if no reference is kept.
The difference is, that you can call map once multiple times with the same closure and it will
close the old result Observable and register a new one instead.

```
function test(s1::Observable)
    s3 = map_once(x-> (println("1 ", x); x), s1)
    s3 = map_once(x-> (println("2 ", x); x), s1)

end
test(Observable(1), Observable(2))
>

"""
function map_once(
        f, input::Observable, inputrest::Observable...
    )
    for arg in (input, inputrest...)
        safe_off(arg, f)
    end
    lift(f, input, inputrest...)
end


text:
using DelimitedFiles

struct CharChar
    a::Char
    b::Char
end

struct CharString
    a::Char
    b::String
end

mutable struct Cart
    id::Int8
    position::CartesianIndex
    symbol::Char
    mod_turns::Int8
    alive::Bool
end

function islesscart(c1::Cart,c2::Cart) :: Bool
    return c1.position[1]<c2.position[1] || (c1.position[1]==c2.position[1] && c1.position[2]<c2.position[2])
end

up = '^'
down = 'v'
left = '<'
right = '>'

vertical = '|'
horizontal = '-'
curve_backslash = '\\'
curve_slash = '/'
intersection = '+'
track_symbols = [vertical,horizontal,curve_slash,curve_backslash,intersection]
carts_symbols = [up,down,left,right]
directions = ["left","straight","right"]

straightMoves = Dict{Char,CartesianIndex}(
    up=>CartesianIndex(-1,0),
    down=>CartesianIndex(1,0),
    left=>CartesianIndex(0,-1),
    right=>CartesianIndex(0,1))
curves = Dict{CharChar,Char}(
    CharChar(up,curve_slash)=>right,
    CharChar(up,curve_backslash)=>left,
    CharChar(down,curve_slash)=>left,
    CharChar(down,curve_backslash)=>right,
    CharChar(left,curve_slash)=>down,
    CharChar(left,curve_backslash)=>up,
    CharChar(right,curve_slash)=>up,
    CharChar(right,curve_backslash)=>down)

turns = Dict{CharString,Char}(
    CharString(up,"left")=>left,
    CharString(up,"right")=>right,
    CharString(down,"left")=>right,
    CharString(down,"right")=>left,
    CharString(left,"left")=>down,
    CharString(left,"right")=>up,
    CharString(right,"left")=>up,
    CharString(right,"right")=>down
)

function solve(input)
    raw = readdlm(input,'\n',String)
    n_col = length(raw[1])
    n_row = length(raw)
    grid = Array{Char,2}(undef,n_row,n_col)
    for y in 1:n_row
        for x in 1:n_col
            grid[y,x] = raw[y][x]
        end
    end
    carts=Array{Cart,1}()
    cart_id=1
    for cart_position in findall(c->c in carts_symbols,grid)
        symbol = grid[cart_position]
        push!(carts,Cart(cart_id,cart_position,symbol,0,true))
        if symbol in carts_symbols[1:2]
            grid[cart_position] = vertical
        else
            grid[cart_position] = horizontal
        end
        cart_id+=1
    end
    #@show carts
    #display(grid,carts)
    tick = 0
    while length(carts)>1
        carts=sort(carts,lt=islesscart)
        #@show carts
        #@show tick
        for cart in carts
            #@show (cart.id,cart.position)
            move(cart,grid)
            carts = crash(carts)
        end
        tick += 1
    end
    if length(carts)==1
        pos = carts[1].position
        x=pos[2]-1
        y=pos[1]-1
        println("Last alive cart position: (X,Y)=($x,$y) id $(carts[1].id)")
    end
end

function display(grid,carts)
    displayed = copy(grid)
    for cart in carts
        if displayed[cart.position] in carts_symbols
            displayed[cart.position]='X'
        else
            displayed[cart.position]=cart.symbol
        end
    end
    for row in 1:size(displayed)[1]
        println(join(displayed[row,:]))
    end
    println("")
end

function move(cart::Cart, grid)
    if cart.alive
        cart.position += straightMoves[cart.symbol]
        next_position = grid[cart.position]
        if next_position in track_symbols[3:4]
            cart.symbol = curves[CharChar(cart.symbol,next_position)]
        elseif next_position == intersection
            cart.mod_turns = (cart.mod_turns % length(directions)) + 1
            next_direction = directions[cart.mod_turns]
            if next_direction != "straight"
                cart.symbol = turns[CharString(cart.symbol,next_direction)]
            end
        end
    end
end

function crash(carts)
    num_carts = length(carts)
    for i in 1:num_carts-1
        cart1 = carts[i]
        for j in i+1:num_carts
            cart2 = carts[j]
            if cart1.position == cart2.position
                cart1.alive = false
                cart2.alive = false
                x=cart1.position[2]-1
                y=cart1.position[1]-1
                println("new crash ! (X,Y)=($x,$y) cart ids $(cart1.id) $(cart2.id)")
                return filter(c->c.alive,carts)
            end
        end
    end
    return carts
end

#solve("input_test.txt")
#solve("input_test.2.txt")
solve("input.txt")


text:
using FeedbackParticleFilters, Distributions, Random, PDMats, FillArrays

println("Testing DiffusionStateModel.jl:")

@testset "DiffusionStateModel.jl" begin

    f(x::AbstractVector) = [x[1]*x[2], x[1]-x[3], x[1]+x[2]*x[3]]
    g(x::AbstractVector) = [x[1] x[2]; x[2] x[3]; 0. 1.]
    init                 = MvNormal(3, 1.)

    print("  inner constructor")
    mod = DiffusionStateModel(f, g, init)
    print(".")
    @test_throws BoundsError DiffusionStateModel(f, g, [1., 2.]) # wrong length of initial condition
    print(".")
    @test_throws ErrorException DiffusionStateModel(f, g, [1., 2., 3., 4.]) # wrong length of initial condition
    print(".")
    g2(x::AbstractVector) = [x[1] x[2]; x[2] x[3]]
    @test_throws ErrorException DiffusionStateModel(f, g2, init) # wrong size of diffusion term
    print(".")
    @test_throws ErrorException DiffusionStateModel(f, g, [0, 0, 0]) # wrong type of initial condition
    println("DONE.")

    print("  method drift")
    print(".")
    @test drift(mod)([1., 2., 3.]) == [2., -2., 7.]
    print(".")
    @test hasmethod(drift(mod), Tuple{AbstractMatrix})
    println("DONE.")

    print("  method diffusion")
    print(".")
    @test diffusion(mod)([1., 2., 3.]) == [1. 2.; 2. 3.; 0. 1.]
    print(".")
    @test hasmethod(diffusion(mod), Tuple{AbstractMatrix})
    println("DONE.")

    print("  method initial_condition")
    print(".")
    @test initial_condition(mod) == mod.init
    println("DONE.")

    print("  method state_dim")
    print(".")
    @test state_dim(mod) == 3
    println("DONE.")

    print("  method noise_dim")
    print(".")
    @test noise_dim(mod) == 2
    println("DONE.")

    print("  method initialize")
    print(".")
    Random.seed!(0)
    @test initialize(mod) == [0.6791074260357777, 0.8284134829000359, -0.3530074003005963]
    print(".")
    mod2 = DiffusionStateModel(f, g, [1., 2., 3.])
    @test initialize(mod2) == [1., 2., 3.]
    println("DONE.")

    print("  method drift_function")
    print(".")
    @test drift_function(mod)([1., 2., 3.]) == drift(mod)([1., 2., 3.])
    println("DONE.")

    print("  method diffusion_function")
    print(".")
    @test diffusion_function(mod)([1., 2., 3.]) == diffusion(mod)([1., 2., 3.])
    println("DONE.")

    print("  calling instance of DiffusionStateModel")
    print(".")
    Random.seed!(0)
    x0 = ones(3)
    @test mod(x0, 0.01) == [1.1607520908935813, 1.1507520908935813, 1.1028413482900037]
    print(".")
    Random.seed!(0)
    xx = hcat(ones(3), zeros(3))
    @test mod(xx, 0.01) == [2.0744253554105256 0. ; 2.1415030557533146 0. ; 2.0697335850849417 -0.07628038164104582]
    println("DONE.")

    print("  method propagate!")
    print(".")
    Random.seed!(0)
    x0 = ones(3)
    propagate!(x0, mod, 0.01)
    @test x0 == [1.0744253554105254, 1.1415030557533146, 1.0697335850849417]
    print(".")
    Random.seed!(0)
    xx = hcat(ones(3), zeros(3))
    propagate!(xx, mod, 0.01)
    @test xx == [1.0744253554105254 0. ; 1.1415030557533146 0. ; 1.0697335850849417 -0.07628038164104582]
    println("DONE.")

end; #DiffusionStateModel.jl


text:
using Test, JuMP, Cbc
using JuMP.Containers: DenseAxisArray

include(joinpath(dirname(@__DIR__), "src", "model.jl"))

# --- Test Utility ---

@test data(DenseAxisArray([1, 2], [:a, :b])) == [1, 2]
@test data([1, 2]) == [1, 2]

@test round_int(1.1) == 1
@test convert_int(1.001, Int) == 1
@test convert_int([1.001, 2.003], Array{Int, 1}) == [1, 2]
@test convert_int(1.001, Float64) == 1.001

# --- Test Model ---

m, n, k = 3, 1, 1
a = [-8]
b = [-1]
c = [14, -33, 20]
A = reshape([-1, -4, 2], m, n)
B = reshape([-2, -1, 1], m, k)

specs = Specs(relax_integer=false)
indices = Indices(m, n, k)
params = Params(a, b, c, A, B)
model = ExampleModel(specs, indices, params)

set_optimizer(model, Cbc.Optimizer)
optimize!(model)

@test termination_status(model) != MOI.INFEASIBLE

variables = Variables(model)
objectives = Objectives(model)

@test true

max_stars_repo_name: nnhjy/JuMP-project-template

max_issues_repo_name: nnhjy/JuMP-project-template

max_forks_repo_name: nnhjy/JuMP-project-template


text:
function freeprecess1()

    t = 100
    spin = Spin(1, 1000, 100, 1.25)
    (A1, B1) = freeprecess(spin, t)
    A2 = FreePrecessionMatrix()
    B2 = Magnetization()
    freeprecess!(A2, B2, spin, t)
    @test A1 == A2
    return B1 == B2

end

function freeprecess2()

    s = Spin(Magnetization(1, 0, 0), 1, 1000, 100, 3.75)
    freeprecess!(s, 100)
    M_correct = Magnetization(-0.2601300475114444, -0.2601300475114445, 0.09516258196404048)
    return s.M ≈ M_correct

end

function freeprecess3()

    t = 100
    grad = Gradient(0, 0, 1)
    spin = Spin(1, 1000, 100, 1.25, Position(0, 0, 1))
    (A1, B1) = freeprecess(spin, t, grad)
    A2 = FreePrecessionMatrix()
    B2 = Magnetization(0.0, 0.0, 0.0)
    freeprecess!(A2, B2, spin, t, grad)
    @test A1 == A2
    return B1 == B2

end

function freeprecess4()

    s = Spin(Magnetization(1, 0, 0), 1, 1000, 100, 0, Position(0, 0, 3.75))
    freeprecess!(s, 100, Gradient(0, 0, 1/GAMBAR))
    M_correct = Magnetization(-0.2601300475114444, -0.2601300475114445, 0.09516258196404048)
    return s.M ≈ M_correct

end

function freeprecess5()

    s = Spin(1, 1000, 100, 3.75)
    A = ExcitationMatrix()
    BtoM = Magnetization()
    excite!(A, s, InstantaneousRF(π/2))
    applydynamics!(s, BtoM, A)
    A = FreePrecessionMatrix()
    B = Magnetization()
    freeprecess!(A, B, s, 100)
    applydynamics!(s, BtoM, A, B)
    M_correct = Magnetization(-0.2601300475114444, -0.2601300475114445, 0.09516258196404054)
    return s.M ≈ M_correct

end

function freeprecess6()

    s = Spin(1, 1000, 100, 3.75)
    (A1, B1) = freeprecess(s, 100)
    (A2, B2) = freeprecess(100, s.M0, s.T1, s.T2, s.Δf)

    @test A1 == A2
    return B1 == B2

end

function freeprecessmc1()

    t = 20
    spin = SpinMC(1, (0.7, 0.2, 0.1), (1000, 400, 1000), (100, 20, 0.01), (0, 15, 0), (100, 100, 25, Inf, Inf, Inf))
    (A1, B1) = freeprecess(spin, t)
    A2 = BlochMcConnellMatrix(3)
    B2 = MagnetizationMC(3)
    freeprecess!(A2, B2, spin, t)
    @test A1 == A2
    return B1 == B2

end

function freeprecessmc2()

    t = 20
    M0 = 1
    (fa, fb) = (0.5, 0.5)
    (T1a, T2a, Δfa) = (1000, 100, 0)
    (T1b, T2b, Δfb) = (400, 20, 15)
    spina = Spin(fa * M0, T1a, T2a, Δfa)
    spinb = Spin(fb * M0, T1b, T2b, Δfb)
    spinmc = SpinMC(M0, (fa, fb), (T1a, T1b), (T2a, T2b), (Δfa, Δfb), (Inf, Inf))
    Aa = FreePrecessionMatrix()
    Ba = Magnetization()
    Ab = FreePrecessionMatrix()
    Bb = Magnetization()
    Amc = BlochMcConnellMatrix(2)
    Bmc = MagnetizationMC(2)
    freeprecess!(Aa, Ba, spina, t)
    freeprecess!(Ab, Bb, spinb, t)
    freeprecess!(Amc, Bmc, spinmc, t)
    @test Matrix(Amc) ≈ [Matrix(Aa) zeros(3, 3); zeros(3, 3) Matrix(Ab)]
    return Bmc ≈ MagnetizationMC(Ba, Bb)

end

function freeprecessmc3()

    t = 20
    grad = Gradient(0, 0, 1)
    spin = SpinMC(1, (0.8, 0.2), (1000, 400), (100, 20), (0, 15), (100, 25), Position(0, 0, 1))
    (A1, B1) = freeprecess(spin, t, grad)
    A2 = BlochMcConnellMatrix(2)
    B2 = MagnetizationMC(2)
    freeprecess!(A2, B2, spin, t, grad)
    @test A1 ≈ A2
    return B1 ≈ B2

end

function freeprecessmc4()

    s = SpinMC(1, (0.8, 0.2), (1000, 400), (100, 20), (0, 15), (Inf, Inf))
    t = 100
    (A1, B1) = freeprecess(s, t)
    (A2, B2) = freeprecess(s, t, nothing)

    @test A1 ≈ A2
    return B1 ≈ B2

end

function combine1()

    A = BlochMatrix()
    B = Magnetization()
    A1 = BlochMatrix(1, 2, 3, 4, 5, 6, 7, 8, 9)
    B1 = nothing
    A2 = BlochMatrix(1, 2, 3, 4, 5, 6, 7, 8, 9)
    B2 = Magnetization(1, 1, 1)
    combine!(A, B, A1, B1, A2, B2)
    Ac = A2 * A1
    Bc = B2

    @test A == Ac
    return B == Bc

end

function combine2()

    A = BlochMatrix()
    B = Magnetization()
    A1 = BlochMatrix(1, 2, 3, 4, 5, 6, 7, 8, 9)
    B1 = Magnetization(1, 1, 1)
    A2 = nothing
    B2 = nothing
    combine!(A, B, A1, B1, A2, B2)
    Ac = A1
    Bc = B1

    @test A == Ac
    return B == Bc

end

function combine3()

    A = BlochMatrix()
    B = Magnetization()
    A1 = nothing
    B1 = nothing
    A2 = BlochMatrix(1, 2, 3, 4, 5, 6, 7, 8, 9)
    B2 = Magnetization(1, 1, 1)
    combine!(A, B, A1, B1, A2, B2)
    Ac = A2
    Bc = B2

    @test A == Ac
    return B == Bc

end

function combine4()

    A = BlochMatrix()
    A1 = BlochMatrix(1, 1, 1, 1, 1, 1, 1, 1, 1)
    A2 = BlochMatrix(1, 2, 3, 4, 5, 6, 7, 8, 9)
    combine!(A, A1, A2)
    Ac = A2 * A1

    return A == Ac

end

@testset "Free Precession" begin

    @test freeprecess1()
    @test freeprecess2()
    @test freeprecess3()
    @test freeprecess4()
    @test freeprecess5()
    @test freeprecess6()
    @test freeprecessmc1()
    @test freeprecessmc2()
    @test freeprecessmc3()
    @test freeprecessmc4()

end

@testset "combine!" begin

    @test combine1()
    @test combine2()
    @test combine3()
    @test combine4()

end


text:
using JinjaTemplates
""" JinjaTemplates.renders a given element tree by passing elements thru Jinja Templates """
function jinja(env, md::Markdown.MD)
    env.scratch[:settings_cache] = flatten(env.settings)
    stream = IOBuffer()
    for element in md.content
        println(stream, jinja(env, element))
    end
    text = readstring(seek(stream, 0))
    return text
end

function jinja(env, md::String)
    # println("ESCAPING ", md, " to ", Markdown.latexesc(md))
    latexesc(md)
end
jinja(env, md::Number) = string(md)
jinja(env, md::Array) = join([jinja(env, i) for i in md], "")

function jinja(env, md)
    fields = Dict{Symbol, Any}()
    for i in fieldnames(md)
        fields[i] = jinja(env, getfield(md, i))
    end
    return JinjaTemplates.render(env.templates, "elements/latex/$(get_template_name(md)).tex"; settings=env.scratch[:settings_cache], fields...)
end

# Elements which need special attention
jinja(env, md::Marble.Document) = ""

function jinja(env, md::Base.Markdown.List)
    listitems = map(md.items) do item
        return jinja(env, item)
    end
    return JinjaTemplates.render(env.templates, "elements/latex/$(get_template_name(md)).tex";
        settings=env.scratch[:settings_cache],
        ordered=md.ordered,
        items=listitems)
end

jinja(env, md::Base.Markdown.Code) = JinjaTemplates.render(
    env.templates,
    "elements/latex/$(get_template_name(md)).tex";
    settings=env.scratch[:settings_cache],
    language=md.language,
    code=md.code
    )

function jinja(env, md::Base.Markdown.Table)
    rows = map(md.rows) do row
        return [jinja(env, col) for col in row]
    end
    headers = rows[1]
    content = rows[2:end]
    alignment = map(md.align) do a
        if a == :l
            return "left"
        elseif a == :r
            return "right"
        else
            return "center"
        end
    end
    return JinjaTemplates.render(env.templates, "elements/latex/$(get_template_name(md)).tex";
        settings=env.scratch[:settings_cache],
        headers = headers,
        content = content,
        alignment = alignment)
end

function jinja(env, md::Marble.Table)
    # TODO: Decide if we want latex or MD inlines within tabels
    # content = map(md.content) do row
    #     return [jinja(env, col) for col in row]
    # end
    # headers = map(md.headers) do row
    #     return [jinja(env, col) for col in row]
    # end
    return JinjaTemplates.render(env.templates, "elements/latex/$(get_template_name(md)).tex";
        settings=env.scratch[:settings_cache],
        path=md.path,
        name=md.name,
        content=md.content,
        caption=jinja(env, md.caption),
        headers=md.headers,
        alignment=md.alignment)
end

function jinja(env, md::Marble.Figure)
    return JinjaTemplates.render(env.templates, "elements/latex/$(get_template_name(md)).tex";
        settings=env.scratch[:settings_cache],
        name=md.name,
        path=isempty(md.path) ? "" : abspath(md.path),
        caption=jinja(env, md.caption),
        tex=md.tex,
        width=md.width)
end

function jinja(env, md::Marble.Equation)
    return JinjaTemplates.render(env.templates, "elements/latex/$(get_template_name(md)).tex";
        settings=env.scratch[:settings_cache],
        name=md.name,
        content=md.content)
end

function jinja(env, md::Marble.Tex)
    return JinjaTemplates.render(env.templates, "elements/latex/$(get_template_name(md)).tex";
        settings=env.scratch[:settings_cache],
        content=md.content)
end

function jinja(env, md::Marble.InlineCite)
    return JinjaTemplates.render(env.templates, "elements/latex/$(get_template_name(md)).tex";
        settings=env.scratch[:settings_cache],
        sources=md.sources)
end

function jinja(env, md::Marble.InlineUnit)
    return JinjaTemplates.render(env.templates, "elements/latex/$(get_template_name(md)).tex";
        settings=env.scratch[:settings_cache],
        units=md.units)
end

function jinja(env, md::Marble.InlineData)
    text = format(md.format, env.scratch[:analysis][md.text])
    return JinjaTemplates.render(env.templates, "elements/latex/$(get_template_name(md)).tex";
        settings=env.scratch[:settings_cache],
        text=text,
        format=md.format,
        original=md.text)
end

function jinja(env, md::Marble.InlineMath)
    return JinjaTemplates.render(env.templates, "elements/latex/$(get_template_name(md)).tex";
        settings=env.scratch[:settings_cache],
        text=md.text)
end

function jinja(env, md::Marble.InlineCE)
    return JinjaTemplates.render(env.templates, "elements/latex/$(get_template_name(md)).tex";
        settings=env.scratch[:settings_cache],
        text=md.text)
end

function jinja(env, md::Marble.InlineTex)
    return JinjaTemplates.render(env.templates, "elements/latex/$(get_template_name(md)).tex";
        settings=env.scratch[:settings_cache],
        text=md.text)
end

##################### SUPPORT FUNCTIONS #######################
const _latexescape_chars = Dict{Char, AbstractString}(
   '~'=>"{\\textasciitilde}", '^'=>"{\\textasciicircum}", '\\'=>"{\\textbackslash}")
for ch in "&%\$#_{}"
    _latexescape_chars[ch] = "\\$ch"
end

function latexesc(io, s::String)
    for ch in s
        print(io, get(_latexescape_chars, ch, ch))
    end
end
latexesc(s::String) = sprint(latexesc, s)


text:
function hunt(;nx::Int=3, ny::Int=3, Re::Float64 = 10.0, Ha::Float64 = 10.0,
    f_u = nothing, resultsfile = nothing)

  reset_timer!()

  @timeit "model" begin
  L = 1.0
  N = Ha^2/Re
  K = Ha / (1-0.95598*Ha^(-1/2)-Ha^(-1))
  ∂p∂z = -Re * K / L^3

  if f_u === nothing
    _f_u = VectorValue(0.0,0.0, -∂p∂z) * L/Re^2 # Assumed ν=L=1 by default
  else
    _f_u = f_u
  end
  g_u = VectorValue(0.0,0.0,0.0)
  g_j = VectorValue(0.0,0.0,0.0)
  g_φ = 0.0
  B = VectorValue(0.0,1.0,0.0)

  # Discretization
  order = 2
  domain = (-1.0,1.0,-1.0,1.0,0.0,0.1)
  map(x) = VectorValue(sign(x[1])*(abs(x[1])*0.5)^0.5,
                       sign(x[2])*(abs(x[2])*0.5)^0.5, x[3]*sqrt(2)/2)*2/sqrt(2)

  dirichlet_tags_u = append!(collect(1:20),[23,24,25,26])
  dirichlet_tags_j = append!(collect(1:20),[25,26])

  partition=(nx,ny,3)
  model = CartesianDiscreteModel(domain,partition;
    isperiodic=(false,false,true), map=map)

  labels = get_face_labeling(model)
  add_tag_from_tags!(labels,"dirichlet_u",dirichlet_tags_u)
  add_tag_from_tags!(labels,"dirichlet_j",dirichlet_tags_j)
  end
  @timeit "FE spaces" begin

  Vu = FESpace(model, ReferenceFE(lagrangian,VectorValue{3,Float64},order);
      conformity=:H1, dirichlet_tags="dirichlet_u")

  Vp = FESpace(model, ReferenceFE(lagrangian,Float64,order-1,space=:P);
      conformity=:L2)

  Vj = FESpace(model, ReferenceFE(raviart_thomas,Float64,order-1);
      conformity=:Hdiv, dirichlet_tags="dirichlet_j")

  Vφ = FESpace(model, ReferenceFE(lagrangian,Float64,order-1,space=:Q);
      conformity=:L2)

  U = TrialFESpace(Vu,g_u)
  P = TrialFESpace(Vp)
  J = TrialFESpace(Vj,g_j)
  Φ = TrialFESpace(Vφ)

  Y = MultiFieldFESpace([Vu, Vp, Vj, Vφ])
  X = MultiFieldFESpace([U, P, J, Φ])
  end

  @timeit "Integration" begin
  # Integration
  trian = Triangulation(model)
  degree = 2*(order)
  dΩ = Measure(trian,degree)

  res(x,y) = InductionlessMHD.dimensionless_residual(x, y, Re, N, B, _f_u, dΩ)
  jac(x,dx,y) = InductionlessMHD.dimensionless_jacobian(x, dx, y, Re, N, B, dΩ)

  op  = FEOperator(res,jac,X,Y)
  end

  @timeit "Setup solver" begin
  # Solver
  nls = NLSolver(GmresSolver(preconditioner=ilu,τ=1e-6);
    show_trace=true, method=:newton, linesearch=BackTracking())
  solver = FESolver(nls)
  end

  @timeit "solve" xh = solve(solver,op)

  print_timer()
  println()

  if resultsfile != nothing
    uh, ph, jh, φh = xh
    writevtk(trian, resultsfile,
      cellfields=["uh"=>uh, "ph"=>ph, "jh"=>jh, "φh"=>φh])
  end
  (xh, trian, dΩ)
end

function analytical_hunt_u(a::Float64,       # semi-length of side walls
                           b::Float64,       # semi-length of Hartmann walls
                           μ::Float64,       # fluid viscosity
                           grad_pz::Float64, # presure gradient
                           Ha::Float64,      # Hartmann number
                           n::Int,           # number of sumands included in Fourier series
                           x)                # evaluation point
  l = b/a
  ξ = x[1]/a
  η = x[2]/a

  V = 0.0; V0=0.0;
  for k in 0:n
    α_k = (k + 0.5)*π/l
    N = (Ha^2 + 4*α_k^2)^(0.5)
    r1_k = 0.5*( Ha + N)
    r2_k = 0.5*(-Ha + N)

    num = exp(-r1_k*(1-η))+exp(-r1_k*(1+η))
    den = 1+exp(-2*r1_k)
    V2 = (r2_k/N)*(num/den)

    num = exp(-r2_k*(1-η))+exp(-r2_k*(1+η))
    den = 1+exp(-2*r2_k)
    V3 = (r1_k/N)*(num/den)

    V += 2*(-1)^k*cos(α_k * ξ)/(l*α_k^3) * (1-V2-V3)
  end
  u_z = V/μ * (-grad_pz) * a^2

  return VectorValue(0.0*u_z,0.0*u_z,u_z)
end

function analytical_hunt_j(a::Float64,       # semi-length of side walls
                           b::Float64,       # semi-length of Hartmann walls
                           σ::Float64,       # fluid conductivity
                           μ::Float64,       # fluid viscosity
                           grad_pz::Float64, # presure gradient
                           Ha::Float64,      # Hartmann number
                           n::Int,           # number of sumands included in Fourier series
                           x)                # evaluation point
  l = b/a
  ξ = x[1]/a
  η = x[2]/a

  H_dx = 0.0; H_dy = 0.0
  for k in 0:n
    α_k = (k + 0.5)*π/l
    N = sqrt(Ha^2 + 4*α_k^2)
    r1_k = 0.5*( Ha + N)
    r2_k = 0.5*(-Ha + N)

    num = exp(-r1_k*(1-η)) - exp(-r1_k*(1+η))
    num_dy = exp(-r1_k*(1-η))*(r1_k/a) + exp(-r1_k*(1+η))*(r1_k/a)
    den = 1+exp(-2*r1_k)
    H2 = (r2_k/N)*(num/den)
    H2_dy = (r2_k/N)*(num_dy/den)

    num = exp(-r2_k*(1-η)) - exp(-r2_k*(1+η))
    num_dy = exp(-r2_k*(1-η))*(r2_k/a) + exp(-r2_k*(1+η))*(r2_k/a)
    den = 1+exp(-2*r2_k)
    H3 = (r1_k/N)*(num/den)
    H3_dy = (r1_k/N)*(num_dy/den)

    H_dx += -2*(-1)^k * sin(α_k * ξ)/(a*l*α_k^2) * (H2 - H3)
    H_dy += 2*(-1)^k * cos(α_k * ξ)/(l*α_k^3) * (H2_dy - H3_dy)
  end
  j_x = a^2*σ^0.5 / μ^0.5 * (-grad_pz) * H_dy
  j_y = a^2*σ^0.5 / μ^0.5 * (-grad_pz) * (-H_dx)

  return VectorValue(j_x,j_y,0.0)
end


text:
# En Pfam 30.0 PF00400 has 268378 sequences
@benchgroup "Residue conversions" ["IO", "MSA"] begin

    Random.seed!(1)
    chars = rand(['.','-','a':'z'...,'A':'Z'...], 268378*2)
    residues = Residue[ char for char in chars ]
    ints = Int[ res for res in residues ]

    @bench "char2res" begin
                        for char in $chars
                            Residue(char)
                        end
                      end
    @bench "res2char" begin
                        for res in $residues
                            Char(res)
                        end
                      end
    @bench "int2res"  begin
                        for i in $ints
                            Residue(i)
                        end
                      end
    @bench "res2int"  begin
                        for res in $residues
                            Int(res)
                        end
                      end
end


text:
using Base.Test
using MPFR

# constructors
x = MPFRFloat{53}()
x = MPFRFloat(12)
y = MPFRFloat(x)
@test x == y
y = MPFRFloat(0xc)
@test x == y
y = MPFRFloat(12.)
@test x == y
y = MPFRFloat(BigInt(12))
@test x == y
y = MPFRFloat(BigFloat(12))
@test x == y
y = MPFRFloat("12")
@test x == y
y = MPFRFloat(float32(12.))
@test x == y
y = MPFRFloat(12//1)
@test x == y

# +
x = MPFRFloat(12)
y = MPFRFloat(30)
@test x + y == MPFRFloat(42)

# -
x = MPFRFloat(12)
y = MPFRFloat(-30)
@test x - y == MPFRFloat(42)

# *
x = MPFRFloat(6)
y = MPFRFloat(9)
@test x * y != MPFRFloat(42)
@test x * y == MPFRFloat(54)

# /
x = MPFRFloat(9)
y = MPFRFloat(6)
@test x / y == MPFRFloat(9/6)

# < / > / <= / >=
x = MPFRFloat(12)
y = MPFRFloat(42)
z = MPFRFloat(30)
@test y > x
@test y >= x
@test y > z
@test y >= z
@test x < y
@test x <= y
@test z < y
@test z <= y
@test y - x >= z
@test y - x <= z
@test !(x >= z)
@test !(y <= z)

# ^
x = MPFRFloat(12)
y = MPFRFloat(4)
@test x^y == MPFRFloat(20736)

# ceil
x = MPFRFloat(12.042)
@test MPFRFloat(13) == ceil(x)

# copysign
x = MPFRFloat(1)
y = MPFRFloat(-1)
@test copysign(x, y) == y
@test copysign(y, x) == x

# isfinite / isinf
x = MPFRFloat(Inf)
y = MPFRFloat(1)
@test isinf(x) == true
@test isinf(y) == false
@test isfinite(x) == false
@test isinf(x) == true

# isnan
x = MPFRFloat(NaN)
y = MPFRFloat(1)
@test isnan(x) == true
@test isnan(y) == false

# convert to
@test convert(MPFRFloat{53}, 1//2) == MPFRFloat("0.5")
@test convert(MPFRFloat{53}, 0.5) == MPFRFloat("0.5")
@test convert(MPFRFloat{53}, 40) == MPFRFloat("40")
@test convert(MPFRFloat{53}, float32(0.5)) == MPFRFloat("0.5")

# convert from
@test convert(Float64, MPFRFloat(0.5)) == 0.5
@test convert(Float32, MPFRFloat(0.5)) == float32(0.5)

# exponent
x = MPFRFloat(0)
@test_fails exponent(x)
x = MPFRFloat(Inf)
@test_fails exponent(x)
x = MPFRFloat(15.674)
@test exponent(x) == exponent(15.674)

# sqrt DomainError
@test_fails sqrt(MPFRFloat(-1))

# precision
old_precision = get_default_precision()
x = MPFRFloat(0)
@test prec(x) == old_precision
set_default_precision(256)
x = MPFRFloat(0)
@test prec(x) == 256
set_default_precision(old_precision)
z = with_precision(240) do
    z = x + 20
    return z
end
@test float(z) == 20.
@test prec(z) == 240
x = MPFRFloat(12)
@test prec(x) == old_precision
@test_fails set_default_precision(1)

# integer_valued
@test integer_valued(MPFRFloat(12))
@test !integer_valued(MPFRFloat(12.12))

# nextfloat / prevfloat
with_precision(53) do
    x = MPFRFloat(12.12)
    @test MPFRFloat(nextfloat(12.12)) == nextfloat(x)
    @test MPFRFloat(prevfloat(12.12)) == prevfloat(x)
end
@test isnan(nextfloat(MPFRFloat(NaN)))
@test isnan(prevfloat(MPFRFloat(NaN)))

# comparisons
x = MPFRFloat(1)
y = MPFRFloat(-1)
z = MPFRFloat(NaN)
ipl = MPFRFloat(Inf)
imi = MPFRFloat(-Inf)
@test x > y
@test x >= y
@test x >= x
@test y < x
@test y <= x
@test y <= y
@test x < ipl
@test x <= ipl
@test x > imi
@test x >= imi
@test imi == imi
@test ipl == ipl
@test imi < ipl
@test z != z
@test !(z == z)
@test !(z <= z)
@test !(z < z)
@test !(z >= z)
@test !(z > z)

# modf
x = MPFRFloat(12)
y = MPFRFloat(0.5)
@test modf(x+y) == (y, x)
x = MPFRFloat(NaN)
@test map(isnan, modf(x)) == (true, true)
x = MPFRFloat(Inf)
y = modf(x)
@test (isnan(y[1]), isinf(y[2])) == (true, true)

# rem
with_precision(53) do
    x = MPFRFloat(2)
    y = MPFRFloat(1.67)
    @test rem(x,y) == rem(2, 1.67)
    y = MPFRFloat(NaN)
    @test isnan(rem(x,y))
    @test isnan(rem(y,x))
    y = MPFRFloat(Inf)
    @test rem(x,y) == x
    @test isnan(rem(y,x))
end

# min/max
x = MPFRFloat(4)
y = MPFRFloat(2)
@test max(x,y) == x
@test min(x,y) == y
y = MPFRFloat(NaN)
@test max(x,y) == x
@test min(x,y) == x
@test isnan(max(y,y))
@test isnan(min(y,y))

# sum
x = MPFRFloat(1)
y = MPFRFloat(2)
z = MPFRFloat(3)
w = MPFRFloat(4)
@test sum([x,y,z,w]) == MPFRFloat(10)
big_array = ones(MPFRFloat, 100)
@test sum(big_array) == MPFRFloat(100)

# promotion
# the array converts everyone to the DEFAULT_PRECISION!
x = MPFRFloat(12)
y = with_precision(60) do
    MPFRFloat(42)
end
@test [x,y] == [MPFRFloat(12), MPFRFloat(42)]

# log / log2 / log10
x = MPFRFloat(42)
@test log(x) == log(42)
@test isinf(log(MPFRFloat(0)))
@test_fails log(MPFRFloat(-1))
@test log2(x) == log2(42)
@test isinf(log2(MPFRFloat(0)))
@test_fails log2(MPFRFloat(-1))
@test log10(x) == log10(42)
@test isinf(log10(MPFRFloat(0)))
@test_fails log10(MPFRFloat(-1))

# exp / exp2 / exp10
x = MPFRFloat(10)
@test exp(x) == exp(10)
@test exp2(x) == 1024
@test exp10(x) == 10000000000

# convert to integer types
x = MPFRFloat(12.1)
y = MPFRFloat(42)
@test_fails convert(Int32, x)
@test_fails convert(Int64, x)
@test_fails convert(BigInt, x)
@test_fails convert(Uint32, x)
@test_fails convert(Uint32, x)
@test convert(Int32, y) == 42
@test convert(Int64, y) == 42
@test convert(BigInt, y) == 42
@test convert(Uint32, y) == 42
@test convert(Uint32, y) == 42

# iround
x = MPFRFloat(42.42)
y = with_precision(256) do
    MPFRFloat("9223372036854775809.2324")
end
z = BigInt("9223372036854775809")
@test iround(x) == 42
@test iround(y) == z
@test typeof(iround(Uint8, x)) == Uint8 && iround(Uint8, x) == 0x2a
@test typeof(iround(Uint16, x)) == Uint16 && iround(Uint16, x) == 0x2a
@test typeof(iround(Uint32, x)) == Uint32 && iround(Uint32, x) == 0x2a
@test typeof(iround(Uint64, x)) == Uint64 && iround(Uint64, x) == 0x2a
@test typeof(iround(Int64, x)) == Int64 && iround(Int64, x) == 42
@test typeof(iround(Int, x)) == Int && iround(Int, x) == 42
@test typeof(iround(Uint, x)) == Uint && iround(Uint, x) == 0x2a

# factorial
with_precision(256) do
    x = MPFRFloat(42)
    @test factorial(x) == factorial(BigInt(42))
    x = MPFRFloat(10)
    @test factorial(x) == factorial(10)
    @test_fails factorial(MPFRFloat(-1))
    @test_fails factorial(MPFRFloat(331.3))
end

# bessel functions
@test_approx_eq besselj(4, MPFRFloat(2)) besselj(4, 2.)
@test_approx_eq besselj0(MPFRFloat(2))  besselj0(2.)
@test_approx_eq besselj1(MPFRFloat(2))  besselj1(2.)
@test_approx_eq bessely(4, MPFRFloat(2))  bessely(4, 2.)
@test_approx_eq bessely0(MPFRFloat(2))  bessely0(2.)
@test_approx_eq bessely1(MPFRFloat(2))  bessely1(2.)

# trigonometric functions
for f in (:sin,:cos,:tan,:sec,:csc,:cot,:acos,:asin,:atan,
        :cosh,:sinh,:tanh,:sech,:csch,:coth,:asinh),
    j in (-1., -0.5, -0.25, .25, .5, 1.)
    @eval begin
        @test_approx_eq ($f)(MPFRFloat($j)) ($f)($j)
    end
end
for f in (:acos,:asin,:acosh,:atanh),
    j in (-2, -1.5)
    @eval begin
        @test_fails ($f)(MPFRFloat($j))
    end
end
for f in (:sin,:cos,:tan,:sec,:csc,:cot,:cosh,:sinh,:tanh,
        :sech,:csch,:coth,:acosh,:asinh),
    j in (1., 1.5, 1.9)
    @eval begin
        @test_approx_eq ($f)(MPFRFloat($j)) ($f)($j)
    end
end
for j in (.25, .5)
    @test_approx_eq atanh(MPFRFloat(j)) atanh(j)
end


text:
using Flux, GraphNeuralNetworks, Graphs, BenchmarkTools, CUDA
using DataFrames, Statistics, JLD2, SparseArrays
CUDA.device!(2)
CUDA.allowscalar(false)

BenchmarkTools.ratio(::Missing, x) = Inf
BenchmarkTools.ratio(x, ::Missing) = 0.0
BenchmarkTools.ratio(::Missing, ::Missing) = missing

function run_single_benchmark(N, c, D, CONV; gtype=:lg)
    data = erdos_renyi(N, c / (N-1), seed=17)
    X = randn(Float32, D, N)

    g = GNNGraph(data; ndata=X, graph_type=gtype)
    g_gpu = g |> gpu

    m = CONV(D => D)
    ps = Flux.params(m)

    m_gpu = m |> gpu
    ps_gpu = Flux.params(m_gpu)

    res = Dict()

    res["CPU_FWD"] = @benchmark $m($g)
    res["CPU_GRAD"] = @benchmark gradient(() -> sum($m($g).ndata.x), $ps)

    try
        res["GPU_FWD"] = @benchmark CUDA.@sync($m_gpu($g_gpu)) teardown=(GC.gc(); CUDA.reclaim())
    catch
        res["GPU_FWD"] = missing
    end

    try
        res["GPU_GRAD"] = @benchmark CUDA.@sync(gradient(() -> sum($m_gpu($g_gpu).ndata.x), $ps_gpu)) teardown=(GC.gc(); CUDA.reclaim())
    catch
        res["GPU_GRAD"] = missing
    end

    return res
end

"""
    run_benchmarks(;
        Ns = [10, 100, 1000, 10000],
        c = 6,
        D = 100,
        layers = [GCNConv, GraphConv, GATConv]
        )

Benchmark GNN layers on Erdos-Renyi ranomd graphs
with average degree `c`. Benchmarks are perfomed for each graph size in the list `Ns`.
`D` is the number of node features.
"""
function run_benchmarks(;
        Ns = [10, 100, 1000, 10000],
        c = 6,
        D = 100,
        layers = [GCNConv, GATConv],
        gtypes = [:coo, :sparse, :dense],
        )

    df = DataFrame(N=Int[], c=Float64[], layer=String[], gtype=Symbol[],
                   time_cpu=Any[], time_gpu=Any[]) |> allowmissing

    for gtype in gtypes
        for N in Ns
            println("## GRAPH_TYPE = $gtype  N = $N")
            for CONV in layers
                res = run_single_benchmark(N, c, D, CONV; gtype)
                row = (;layer = "$CONV",
                        N = N,
                        c = c,
                        gtype = gtype,
                        time_cpu = ismissing(res["CPU"]) ? missing : median(res["CPU"]),
                        time_gpu = ismissing(res["GPU"]) ? missing : median(res["GPU"]),
                    )
                push!(df, row)
            end
        end
    end

    df.gpu_to_cpu = ratio.(df.time_gpu, df.time_cpu)
    sort!(df, [:layer, :N, :c, :gtype])
    return df
end

# df = run_benchmarks()
# for g in groupby(df, :layer); println(g, "\n"); end

# @save "perf/perf_master_20210803_carlo.jld2" dfmaster=df
## or
# @save "perf/perf_pr.jld2" dfpr=df

function compare(dfpr, dfmaster; on=[:N, :c, :gtype, :layer])
    df = outerjoin(dfpr, dfmaster; on=on, makeunique=true, renamecols = :_pr => :_master)
    df.pr_to_master_cpu = ratio.(df.time_cpu_pr, df.time_cpu_master)
    df.pr_to_master_gpu = ratio.(df.time_gpu_pr, df.time_gpu_master)
    return df[:,[:N, :c, :gtype, :layer, :pr_to_master_cpu, :pr_to_master_gpu]]
end

# @load "perf/perf_pr.jld2" dfpr
# @load "perf/perf_master.jld2" dfmaster
# compare(dfpr, dfmaster)


text:
import Pkg; Pkg.activate(@__DIR__); Pkg.instantiate();
using RobotDynamics, Rotations
using TrajectoryOptimization
using StaticArrays, LinearAlgebra
#using TrajOptPlots
#using MeshCat
#using Plots
#using FileIO, MeshIO
using Altro
#using PyPlot
#using Distributions

function Rotations.getproperty(q::QuatRotation, f::Symbol)
    if f == :w
        getfield(q,:q).s
    elseif f == :x
        getfield(q,:q).v1
    elseif f == :y
        getfield(q,:q).v2
    elseif f == :z
        getfield(q,:q).v3
    else
        getfield(q,f)
    end
end

struct Quadrotor{R} <: RigidBody{R}
    n::Int
    m::Int
    mass::Float64
    J::Diagonal{Float64,SVector{3,Float64}}
    Jinv::Diagonal{Float64,SVector{3,Float64}}
    gravity::SVector{3,Float64}
    motor_dist::Float64
    kf::Float64
    km::Float64
    bodyframe::Bool  # velocity in body frame?
    info::Dict{Symbol,Any}
end

function Quadrotor{R}(;
        mass=0.5,
        J=Diagonal(@SVector [0.0023, 0.0023, 0.004]),
        gravity=SVector(0,0,-9.81),
        motor_dist=0.1750,
        kf=1.0,
        km=0.0245,
        bodyframe=false,
        info=Dict{Symbol,Any}()) where R
    Quadrotor{R}(13,4,mass,J,inv(J),gravity,motor_dist,kf,km,bodyframe,info)
end

(::Type{Quadrotor})(;kwargs...) = Quadrotor{UnitQuaternion{Float64}}(;kwargs...)

RobotDynamics.control_dim(::Quadrotor) = 4

function RobotDynamics.forces(model::Quadrotor, x, u)
    q = orientation(model, x)
    kf = model.kf
    g = model.gravity
    m = model.mass

    # Extract motor speeds
    w1 = u[1]
    w2 = u[2]
    w3 = u[3]
    w4 = u[4]

    # Calculate motor forces
    F1 = max(0,kf*w1);
    F2 = max(0,kf*w2);
    F3 = max(0,kf*w3);
    F4 = max(0,kf*w4);
    F = @SVector [0., 0., F1+F2+F3+F4] #total rotor force in body frame

    m*g + q*F # forces in world frame
end

function RobotDynamics.moments(model::Quadrotor, x, u)

    kf, km = model.kf, model.km
    L = model.motor_dist

    # Extract motor speeds
    w1 = u[1]
    w2 = u[2]
    w3 = u[3]
    w4 = u[4]

    # Calculate motor forces
    F1 = max(0,kf*w1);
    F2 = max(0,kf*w2);
    F3 = max(0,kf*w3);
    F4 = max(0,kf*w4);

    # Calculate motor torques
    M1 = km*w1;
    M2 = km*w2;
    M3 = km*w3;
    M4 = km*w4;
    tau = @SVector [L*(F2-F4), L*(F3-F1), (M1-M2+M3-M4)] #total rotor torque in body frame
end

RobotDynamics.inertia(model::Quadrotor) = model.J
RobotDynamics.inertia_inv(model::Quadrotor) = model.Jinv
RobotDynamics.mass(model::Quadrotor) = model.mass

model = Quadrotor();
n,m = size(model)
N = 51             # number of knot points
tf = 10.0               # total time (sec)
dt = tf/(N-1)          # time step (sec)

x0_pos = SA[0, -10, 1.]
xf_pos = SA[0, +10, 1.]
x0 = RobotDynamics.build_state(model, x0_pos, UnitQuaternion(I), zeros(3), zeros(3))
xf = RobotDynamics.build_state(model, xf_pos, UnitQuaternion(I), zeros(3), zeros(3));

x0 = RBState([0,0,1.5], UnitQuaternion(I), zeros(3), zeros(3))
xf = RBState([1,3,2.5], UnitQuaternion(I), zeros(3), zeros(3))

# objective
Q = Diagonal(@SVector fill(0.1, n))
R = Diagonal(@SVector fill(0.01, m))
Qf = Diagonal(@SVector fill(100.0, n))
obj = LQRObjective(Q,R,Qf,xf,N)

# constraints
cons = ConstraintList(n,m,N)
add_constraint!(cons, BoundConstraint(n,m, u_min=zeros(4), u_max=fill(10.0,4)), 1:N-1)
#add_constraint!(cons, CircleConstraint(n, SA_F64[1,2], SA_F64[1,2], SA[0.1,0.1]), 1:N-1)
add_constraint!(cons, GoalConstraint(xf, SA[1,2,3]), N)

# problem
prob = Problem(model, obj, xf, tf, x0=x0, constraints=cons);

u0 = @SVector fill(0.5*model.mass/m, m)
U_hover = [copy(u0) for k = 1:N-1]; # initial hovering control trajectory

initial_controls!(prob, U_hover)
rollout!(prob);

opts = SolverOptions(
    penalty_scaling=100.,
    penalty_initial=0.1,
)

solver = ALTROSolver(prob, opts);
solve!(solver)
println("Cost: ", cost(solver))
println("Constraint violation: ", max_violation(solver))
println("Iterations: ", iterations(solver))

X = states(solver);
U = controls(solver);
X = hcat(Vector.(X)...);
X[1:3,:]
#x = X[1,:]
#y = X[2,:]
#z = X[3,:];

max_stars_repo_name: labicon/crazyswarm-labicon

max_issues_repo_name: labicon/crazyswarm-labicon

max_forks_repo_name: labicon/crazyswarm-labicon


text:
# compare the solvers for all cases
using Distributed;

addprocs(20);
@everywhere using PowerModels, Gurobi, Ipopt, Random, Distributions;
@everywhere using JLD,HDF5;
@everywhere include("solver_run.jl");
@everywhere const GUROBI_ENV = Gurobi.Env();

list_file = readdir("./data");
result_Gurobi = Dict();
result_Ipopt = Dict();

itema = "case5.m";
pmap(i -> run_diff_solver("./data/",itema), 1:20);

# nominal condition run
list_run = [];
list_error = [];
objDiff_dict = Dict();
pgDiff_dict = Dict();
timeDiff_dict = Dict();
cond_dict = Dict();

list_run = [itema for itema in list_file if occursin(".m",itema)];
solver_result = pmap(fn -> run_diff_solver("./data/",fn), list_run);
for j in 1:length(solver_result)
    if length(solver_result[j]) == 2
        objDiff_dict[solver_result[j][1]] = solver_result[j][2][1];
        pgDiff_dict[solver_result[j][1]] = solver_result[j][2][2];
        timeDiff_dict[solver_result[j][1]] = solver_result[j][2][3];
        cond_dict[solver_result[j][1]] = solver_result[j][2][4];
    end
end
save("nominal.jld","results", [objDiff_dict, pgDiff_dict, timeDiff_dict, cond_dict]);


text:
const _js_loaded = [false]

function ijulia()

# if we're in IJulia call setup the notebook js interop
if isdefined(Main, :IJulia) && Main.IJulia.inited
    # borrowed from https://github.com/plotly/plotly.py/blob/2594076e29584ede2d09f2aa40a8a195b3f3fc66/plotly/offline/offline.py#L64-L71
    # and https://github.com/JuliaLang/Interact.jl/blob/cc5f4cfd34687000bc6bc70f0513eaded1a7c950/src/IJulia/setup.jl#L15
    if !_js_loaded[1]
        session = Session()
        d = Dict("session_id" => session.id,"port" => port)
        template = Mustache.template_from_file(joinpath(dirname(@__FILE__),"..","res","PagesJL.js"))
        const _PagesJL_js = render(template,d)

        # println(_PagesJL_js)

        display("text/html", """

PagesJL.js loaded.

         """)
        _js_loaded[1] = true
    end
end

end


text:
# Autogenerated wrapper script for GAP_pkg_orb_jll for aarch64-linux-gnu
export orb

using GAP_jll
using GAP_lib_jll
JLLWrappers.@generate_wrapper_header("GAP_pkg_orb")
JLLWrappers.@declare_file_product(orb)
function __init__()
    JLLWrappers.@generate_init_header(GAP_jll, GAP_lib_jll)
    JLLWrappers.@init_file_product(
        orb,
        "lib/gap/orb.so",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using VerilogFIR
using Test

@testset "Testing Auxiliary Functions" begin
    # maxrep
    for (i,j) in zip((2,3,4,5,6),(1,3,7,15,31))
        @test VerilogFIR.maxrep(i) == j
    end
    # minrep
    for (i,j) in zip((2,3,4,5,6),(2,4,8,16,32))
        @test VerilogFIR.minrep(i) == -j
    end
    # fminabs
    let
        t1 = (1.0, 10, -400, -0.1)
        @test VerilogFIR.fminabs(t1...) == 0
        # Test case when a division by 0 occurs, make sure result is still correct.
        t2 = (Inf, 1.1)
        @test VerilogFIR.fminabs(t2...) == 1
        t3 = (-Inf, Inf, -2.3)
        @test VerilogFIR.fminabs(t3...) == 2
    end
    # clog2
    let
        tests       = (1,2,3,2.5,100.23)
        expected    = (0,1,2,2,7)
        for (t,e) in zip(tests, expected)
            @test VerilogFIR.clog2(t) == e
        end
    end
    # flog2
    let
        tests       = (1,2,3,2.5,100.23)
        expected    = (0,1,1,1,6)
        for (t,e) in zip(tests, expected)
            @test VerilogFIR.flog2(t) == e
        end
    end
    # nbits_u
    let
        tests    = (6,7,8,9)
        expected = (3,3,4,4)
        for (t,e) in zip(tests, expected)
            @test VerilogFIR.nbits_u(t) == e
        end
    end
    # nbits_s
    let
        tests    = (7,8,9,-7,-8,-9)
        expected = (4,5,5,4,4,5)
        for (t,e) in zip(tests, expected)
            @test VerilogFIR.nbits_s(t) == e
        end
    end
    # tab
    let
        @test VerilogFIR.tab(1) == "    "
        @test VerilogFIR.tab(2) == "        "
    end
    # makedict
    let
        a = 1
        b = [10,20,30]
        c = "hello"
        d = 5
        f = 2.21
        dict = VerilogFIR.@makedict a b c d f
        expected = Dict(
            :a => 1,
            :b => [10,20,30],
            :c => "hello",
            :d => 5,
            :f => 2.21,
        )
        @test dict == expected
    end
end

#=
Tests for functions that scale coefficients.
=#
@testset "Testing Analysis Functions" begin
    # Test coefficient scaling
    let
        # Test with all positive numbers

        ## TEST 1
        coeffs = [2,1,1]
        w_coeff = 3
        @test VerilogFIR.scale_coefficients(coeffs, w_coeff) == ([2,1,1],0)

        ## TEST 2
        w_coeff = 4
        @test VerilogFIR.scale_coefficients(coeffs, w_coeff) == ([4,2,2],1)

        ## TEST 3
        # Test with negative numbers
        # With the same values as TEST 1, should scale since negative numbers
        # have one extra number to work with on the low end.
        coeffs = [-2,-1,-1]
        w_coeff = 3
        @test VerilogFIR.scale_coefficients(coeffs, w_coeff) == ([-4,-2,-2],1)

        ## TEST 4
        # Test with positive and negative numbers
        coeffs = [2,0,-2]
        w_coeff = 3
        @test VerilogFIR.scale_coefficients(coeffs, w_coeff) == ([2,0,-2],0)

        ## TEST 5
        w_coeff = 4
        @test VerilogFIR.scale_coefficients(coeffs, w_coeff) == ([4,0,-4],1)

        ## TEST 6
        coeffs = [1,0,-2]
        w_coeff = 3
        @test VerilogFIR.scale_coefficients(coeffs, w_coeff) == ([2,0,-4],1)
    end
    # Test accumulator width
    let
        coeffs = [1,-2,0,-1,3]
        w_input = 4
        # Get the maximum positive and negative values that can be represented
        # with 4 bits
        big     = 7
        little  = -8
        # Compute maximum and minimum values for the filter
        maxval = 7*(1+3) + -8*(-1-2)
        minval = -8*(1+3) + 7*(-1-2)
        # Get the number of bits for each
        bits_max = VerilogFIR.nbits_s(maxval)
        bits_min = VerilogFIR.nbits_s(minval)
        nbits = min(bits_max, bits_min)
        @test VerilogFIR.acc_width(coeffs, w_input) == nbits
    end
end

@testset "Testing complete run" begin
    generate_fir(stdout, ones(4)/4, w_input = 16, w_output = 16, w_coeff = 12)
    generate_fir(stdout, ones(4)/4, w_input = 8,  w_output = 8, w_coeff = 21)
    generate_fir(stdout, ones(4)/4, w_input = 16, w_output = 8, w_coeff = 9)
    # Test generating saturating logic
    generate_fir(ones(8))

    @test_throws InexactError generate_fir(stdout, ones(4)/4, w_coeff = 1)
    # Test with file
    file = "test15713209.v"
    generate_fir(file, ones(4)/4)
    rm(file)
end


text:
import Base: show, ==, hash

export Qubit, Edge, Marker

import .config.get_qubit_params
import .config.get_marker_params
import .config.get_edge_params
import .config.get_instrument_params

@compat abstract type Channel end
show(io::IO, c::Channel) = print(io, c.label)

# TODO: is there a benefit to having the concrete channel types immutable?
"""
Channel representing single qubit drive.
"""
immutable Qubit <: Channel
	label::String
	awg_channel::String
	gate_channel::String
	shape_params::Dict{Symbol,Any}
	frequency::Float64
end

function Qubit(label)
	channel_params = get_qubit_params()

	if label in keys(channel_params)
		q_params = channel_params[label]["control"]
		phys_chan = get(q_params, "AWG", "")
		gate_chan = get(q_params, "gate", "")

		# pull out shape_params and convert keys to symbols for splatting into shape function
		shape_params = Dict{Symbol, Any}(Symbol(k) => v for (k,v) in q_params["pulse_params"])

		# translate pulse function in shape params from a string to a function handle and snakeify key
		shape_params[:shape_function] = getfield(QGL.PulseShapes, Symbol(pop!(shape_params, :shape_fun)))

		# snakeify some other parameters
		if :dragScaling in keys(shape_params)
			shape_params[:drag_scaling] = pop!(shape_params, :dragScaling)
		end

		Qubit(label, phys_chan, gate_chan, shape_params, q_params["frequency"])
	else
		warn("Unable to find Qubit label $label in channel json file. Creating default Qubit channel")
		Qubit(label,"", "", Dict{String,Any}(), 0.0)
	end

end

"""
Channel representing a digital output marker line.
"""
immutable Marker <: Channel
	label::String
	awg_channel::String
	shape_params::Dict{Any, Any}
end

function Marker(label)
	phys_chan = get_marker_params()[label]
	m_params = get_instrument_params()[split(phys_chan)[1]]["markers"][split(phys_chan)[2]]
	# translate pulse function in shape params from a string to a function handle and snakeify key
	shape_params =  Dict{Symbol, Any}(Symbol(k) => v for (k,v) in m_params)
	shape_params[:shape_function] = getfield(QGL.PulseShapes,  Symbol(pop!(shape_params, :shape_fun)))
	Marker(label, phys_chan, shape_params)
end

# NOTE is this mutable on purpose??
type QuadratureAWGChannel
	awg::String
	delay::Real
	mixer_correction::Matrix{Real}
end

==(a::Channel, b::Channel) = a.label == b.label
hash(c::Channel, h::UInt) = hash(c.label, h)

"""
Channel representing qubit measurement drive
"""
immutable Measurement <: Channel
	label::String
	awg_channel::String
	gate_channel::String
	trigger_channel::String
	shape_params::Dict{Any,Any}
	frequency::Real
end

"""
	measurement_channel(l:String)

Looks us the measurement channel associated with qubit.
"""
function measurement_channel(l::String)
	m_params = get_qubit_params()[l]["measure"]
	phys_chan = get(m_params, "AWG", "")
	gate_chan = get(m_params, "gate", "")
	trig_chan = get(m_params, "trigger", "")
	# pull out shape_params and convert keys to symbols for splatting into shape function
	shape_params = Dict{Symbol, Any}(Symbol(k) => v for (k,v) in m_params["pulse_params"])
	# translate pulse function in shape params from a string to a function handle and snakeify key
	shape_params[:shape_function] = getfield(QGL.PulseShapes,  Symbol(pop!(shape_params, :shape_fun)))
	# inject autodyne frequence into the `shape_params` where it should be
	shape_params[:autodyne_freq] = m_params["autodyne_freq"]
	Measurement(string("M-", l), phys_chan, gate_chan, trig_chan, shape_params, 0)
end

"""
Channel representing an interaction drive.
"""
immutable Edge <: Channel
	label::String
	source::Qubit
	target::Qubit
	awg_channel::String
	gate_channel::String
	shape_params::Dict{Symbol,Any}
	frequency::Float64
end

"""
	Edge(source::Qubit, target::Qubit)

Create the edge representing interaction drive from `source` to `target`.
"""
function Edge(source::Qubit, target::Qubit)
	# look up whether we have an edge connecting source -> target
	channel_params = get_edge_params()

	edges = filter(
		(k,v) -> v["source"] == source.label && v["target"] == target.label,
		channel_params)

	@assert length(edges) == 1 "Found $(length(edges)) matching edges for $source → $target"

	e_params = collect(values(edges))[1]

	phys_chan = get(e_params, "AWG", "")
	gate_chan = get(e_params, "gate", "")

	# pull out shape_params and convert keys to symbols for splatting into shape function
	shape_params = e_params["pulse_params"]
	shape_params = Dict{Symbol, Any}(Symbol(k) => v for (k,v) in shape_params)

	# translate pulse function in shape params from a string to a function handle and snakeify key
	shape_params[:shape_function] = getfield(QGL.PulseShapes, Symbol(pop!(shape_params, :shape_fun)))

	Edge(collect(keys(edges))[1], source, target, phys_chan, gate_chan, shape_params, e_params["frequency"])
end


text:
"""
nk_test()

CI for nsoli
Testing Eisenstat-Walker and functions witout precomputed data
"""
function nk_test()
    passsimple = nksimple()
    passsimple || println("nksimple fails")
    jvpass = jacvec2d()
    jvpass || println("jacvec2d fails")
    nkpass = passsimple && jvpass
    return nkpass
end
"""
nksimple()

Test nsoli with the simple 2D problem and line search failure and success.
"""
function nksimple()
    x0 = [2.0; 0.5]
    FPS = zeros(2, 3)
    FPJ = zeros(2, 2)
    FS = copy(x0)
    #
    # For the easy problem we will do analytic Jacobians for
    # Newton and forward difference directional derivatives for Newton-GMRES
    #
    dout = nsol(simple!, x0, FS, FPJ, jsimple!; sham = 1, keepsolhist = true)
    kout = nsoli(
        simple!,
        x0,
        FS,
        FPS;
        eta = 1.e-10,
        keepsolhist = true,
        fixedeta = false,
    )
    dsolhist = norm(kout.solhist - dout.solhist, Inf)
    shpass = (dsolhist < 1.e-7)
    shpass || println("solhist compare fails in easy nksimple", dsolhist)
    vconverge = krstest(dout, kout)
    #
    # For the stagnating problem we will do analytic Jacobians for
    # Newton and analytic Jacobian-vector products for Newton-GMRES
    #
    x0 = [3.0; 5.0]
    dout = nsol(simple!, x0, FS, FPJ, jsimple!; sham = 1)
    kout = nsoli(simple!, x0, FS, FPS, JVsimple; fixedeta = true,
                eta = 1.e-10)
    vdiverge = krstest(dout, kout)
    vdiverge || println("failure hard nksimple problem")
    return vconverge && vdiverge && shpass
end

function krstest(dout, kout)
    hdiff = norm(kout.history - dout.history, Inf)
    hpass = (hdiff < 5.e-7)
    hpass || println("history compare fails in nksimple")
    #
    adiff = kout.stats.iarm - dout.stats.iarm
    apass = (sum(adiff) == 0)
    apass || println("line search compare fails in nksimple")
    #
    fdiff = kout.stats.ifun - dout.stats.ifun
    fpass = (sum(fdiff) == 0)
    fpass || println("function value compare fails in nksimple")
    #
    soldiff = kout.solution - dout.solution
    solpass = (norm(soldiff, Inf) < 1.e-9)
    solpass || println("solution compare fails in nksimple  ", norm(soldiff, Inf))
    krpass = (fpass && apass && hpass && solpass)
end

"""
jacvec2d()

Analytic Jacobian-vector product. Compare Eisenstat-Walker to
fixed eta. Test precomputed data support.
"""
function jacvec2d()
    x0 = ones(2)
    fv = zeros(2)
    jv = zeros(2, 2)
    jvs = zeros(2, 3)
    pdata = zeros(2)
    nout = nsol(f!, x0, fv, jv; sham = 1, pdata = pdata)
    kout = nsoli(f!, x0, fv, jvs, JVec; fixedeta = false, eta = 0.9, lmaxit = 2, pdata = pdata)
    kout2 = nsoli(fv2!, x0, fv, jvs, JVecv2; fixedeta = true, eta = 0.1,
                 lmaxit = 2, Pvec=PVecv2)
    histdiff = norm(nout.history - kout2.history)
    histpass = (histdiff < 1.e-5)
    histpass || println("hist test fails in jacvec2d")
    ncost = funcost(nout)
    nplot = acost(nout)
    kcost = funcost(kout)
    kplot = acost(kout)
    kcost2 = funcost(kout2)
    kplot2 = acost(kout2)
    costpass = (ncost == 10) && (kcost == 15) && (kcost2 == 14)
    costpass || println("cost compare fails in jacvec2d")
    costpass || println(ncost,"  ", kcost, "  ", kcost2)
    soldiff = (
        norm(kout.solution - nout.solution, Inf) +
        norm(kout2.solution - nout.solution, Inf)
    )
    solpass = (soldiff < 1.e-7)
    solpass || println("solution compare fails in jacvec2d")
    jvpass = histpass && costpass && solpass
    return jvpass
end

function f!(fv, x, pdata)
    fv[1] = x[1] + sin(x[2])
    fv[2] = cos(x[1] + x[2])
    return fv
end

"""
fv2!(fv, x)

Function evaluation witout precomputed data for testing.
"""
function fv2!(fv, x)
    fv[1] = x[1] + sin(x[2])
    fv[2] = cos(x[1] + x[2])
    return fv
end

"""
PVecv2(v, x)

Here's a preconditioner that does not need procomputed data and
does not do anything.
"""
function PVecv2(v, x)
return v
end

"""
JVecv2(v, fv, x)

Precondition without precomputed/stored data
"""
function JVecv2(v, fv, x)
    jvec = zeros(2)
    p = -sin(x[1] + x[2])
    jvec[1] = v[1] + cos(x[2]) * v[2]
    jvec[2] = p * (v[1] + v[2])
    return jvec
end

"""
JVec(v, fv, x, pdata)

Precondition with precomputed/stored data
"""
function JVec(v, fv, x, pdata)
    jvec = zeros(2)
    p = -sin(x[1] + x[2])
    pdata[1] = v[1] + cos(x[2]) * v[2]
    pdata[2] = p * (v[1] + v[2])
    return pdata
end

function funcost(itout)
    netcost = itout.stats.ifun + itout.stats.ijac
    cost = sum(netcost)
end

function acost(itout)
    netcost = itout.stats.ifun + itout.stats.ijac
    cost = cumsum(netcost)
end


text:
# Affinity propagation
#
#   Reference:
#       Clustering by Passing Messages Between Data Points.
#       Brendan J. Frey and Delbert Dueck
#       Science, vol 315, pages 972-976, 2007.
#

#### Interface

mutable struct AffinityPropResult <: ClusteringResult
    exemplars::Vector{Int}      # indexes of exemplars (centers)
    assignments::Vector{Int}    # assignments for each point
    counts::Vector{Int}         # number of samples in each cluster
    iterations::Int             # number of iterations executed
    converged::Bool             # converged or not
end

const _afp_default_maxiter = 200
const _afp_default_damp = 0.5
const _afp_default_tol = 1.0e-6
const _afp_default_display = :none

function affinityprop(S::AbstractArray{T};
                      maxiter::Integer=_afp_default_maxiter,
                      tol::Real=_afp_default_tol,
                      damp::Real=_afp_default_damp,
                      display::Symbol=_afp_default_display) where T<:AbstractFloat

    # check arguments
    n = size(S, 1)
    size(S, 2) == n || error("S must be a square matrix.")
    n >= 2 || error("the number of samples must be at least 2.")
    tol > 0 || error("tol must be a positive value.")
    0 <= damp < 1 || error("damp must be a non-negative real value below 1.")

    # invoke core implementation
    _affinityprop(S, round(Int, maxiter), tol, convert(T, damp), display_level(display))
end

#### Implementation

function _affinityprop(S::AbstractArray{T},
                       maxiter::Int,
                       tol::Real,
                       damp::T,
                       displevel::Int) where T<:AbstractFloat
    n = size(S, 1)
    n2 = n * n

    # initialize messages
    R = zeros(T, n, n)  # responsibilities
    A = zeros(T, n, n)  # availabilities

    # prepare storages
    Rt = Matrix{T}(undef, n, n)
    At = Matrix{T}(undef, n, n)

    if displevel >= 2
        @printf "%7s %12s | %8s \n" "Iters" "objv-change" "exemplars"
        println("-----------------------------------------------------")
    end

    t = 0
    converged = false
    while !converged && t < maxiter
        t += 1

        # compute new messages
        _afp_compute_r!(Rt, S, A)
        _afp_dampen_update!(R, Rt, damp)

        _afp_compute_a!(At, R)
        _afp_dampen_update!(A, At, damp)

        # determine convergence
        ch = max(Linfdist(A, At), Linfdist(R, Rt)) / (one(T) - damp)
        converged = (ch < tol)

        if displevel >= 2
            # count the number of exemplars
            ne = _afp_count_exemplars(A, R)
            @printf("%7d %12.4e | %8d\n", t, ch, ne)
        end
    end

    # extract exemplars and assignments
    exemplars = _afp_extract_exemplars(A, R)
    assignments, counts = _afp_get_assignments(S, exemplars)

    if displevel >= 1
        if converged
            println("Affinity propagation converged with $t iterations: $(length(exemplars)) exemplars.")
        else
            println("Affinity propagation terminated without convergence after $t iterations: $(length(exemplars)) exemplars.")
        end
    end

    # produce output struct
    return AffinityPropResult(exemplars, assignments, counts, t, converged)
end

# compute responsibilities
function _afp_compute_r!(R::Matrix{T}, S::AbstractArray{T}, A::Matrix{T}) where T
    n = size(S, 1)

    I1 = Vector{Int}(undef, n)  # I1[i] is the column index of the maximum element in (A+S)[i,:]
    Y1 = Vector{T}(undef, n)    # Y1[i] is the maximum element in (A+S)[i,:]
    Y2 = Vector{T}(undef, n)    # Y2[i] is the second maximum element in (A+S)[i,:]

    # Find the first and second maximum elements along each row
    @inbounds for i = 1:n
        v1 = A[i,1] + S[i,1]
        v2 = A[i,2] + S[i,2]
        if v1 > v2
            I1[i] = 1
            Y1[i] = v1
            Y2[i] = v2
        else
            I1[i] = 2
            Y1[i] = v2
            Y2[i] = v1
        end
    end
    @inbounds for j = 3:n, i = 1:n
        v = A[i,j] + S[i,j]
        if v > Y2[i]
            if v > Y1[i]
                Y2[i] = Y1[i]
                I1[i] = j
                Y1[i] = v
            else
                Y2[i] = v
            end
        end
    end

    # compute R values
    @inbounds for j = 1:n, i = 1:n
        mv = (j == I1[i] ? Y2[i] : Y1[i])
        R[i,j] = S[i,j] - mv
    end

    return R
end

# compute availabilities
function _afp_compute_a!(A::Matrix{T}, R::Matrix{T}) where T
    n = size(R, 1)
    z = zero(T)
    for j = 1:n
        @inbounds rjj = R[j,j]

        # compute s <- sum_{i \ne j} max(0, R(i,j))
        s = z
        for i = 1:n
            if i != j
                @inbounds r = R[i,j]
                if r > z
                    s += r
                end
            end
        end

        for i = 1:n
            if i == j
                @inbounds A[i,j] = s
            else
                @inbounds r = R[i,j]
                u = rjj + s
                if r > z
                    u -= r
                end
                A[i,j] = ifelse(u < z, u, z)
            end
        end
    end
    return A
end

# dampen update
function _afp_dampen_update!(x::Array{T}, xt::Array{T}, damp::T) where T
    ct = one(T) - damp
    for i = 1:length(x)
        @inbounds x[i] = ct * xt[i] + damp * x[i]
    end
    return x
end

# count the number of exemplars
function _afp_count_exemplars(A::Matrix, R::Matrix)
    n = size(A,1)
    c = 0
    for i = 1:n
        @inbounds if A[i,i] + R[i,i] > 0
            c += 1
        end
    end
    return c
end

# extract all exemplars
function _afp_extract_exemplars(A::Matrix, R::Matrix)
    n = size(A,1)
    r = Int[]
    for i = 1:n
        @inbounds if A[i,i] + R[i,i] > 0
            push!(r, i)
        end
    end
    return r
end

# get assignments
function _afp_get_assignments(S::AbstractArray, exemplars::Vector{Int})
    n = size(S, 1)
    k = length(exemplars)
    Se = S[:, exemplars]
    a = Vector{Int}(undef, n)
    cnts = zeros(Int, k)
    for i = 1:n
        p = 1
        v = Se[i,1]
        for j = 2:k
            s = Se[i,j]
            if s > v
                v = s
                p = j
            end
        end
        a[i] = p
    end
    for i = 1:k
        a[exemplars[i]] = i
    end
    for i = 1:n
        @inbounds cnts[a[i]] += 1
    end
    return (a, cnts)
end


text:
module LevenshteinToolkit

# ----------------------------------------
# EXPORTED INTERFACE
# ----------------------------------------

export distance_matrix
export distance_row

export dfa
export nfa

export check
export draw

include("matrix.jl")
include("row.jl")
include("automata.jl")

end

max_stars_repo_name: jlegare/LevenshteinToolkit

max_issues_repo_name: jlegare/LevenshteinToolkit

max_forks_repo_name: jlegare/LevenshteinToolkit


text:
# This file was generated, do not modify it. # hide
r1 = range(xgb, :max_depth, lower=3, upper=10)
r2 = range(xgb, :num_round, lower=1, upper=25);


text:
using Base.Test

@testset "Study" begin
  @test true
  @test_throws DomainError throw(DomainError())
  @test_throws ErrorException error("Error Message")
end

try
  error("Error 3")
catch e
  println("Error type: $(typeof(e)), message: $(e.msg)")
end

#error("Error 2")
function f()
end


text:
# Base.@kwdef struct KBOptions{T}
#     β::T = 0.5
#     solver::F = IPOT
#     tol::Float64 = 1e-6
#     iters::Int = 1000
# end

struct KBResult{T} <: Clustering.ClusteringResult
    assignments::Vector{Int}
    barycenters::T
    cost::Float64
end

Clustering.nclusters(R::KBResult) = length(R.barycenters)
Clustering.counts(R::KBResult) = [count(R.assigments .== i) for i in 1:nclusters(R)]
Clustering.assignments(R::KBResult) = R.assignments

"""
    kbarycenters(X, p, k; seed = :rand, kiters = 10, verbose = false, output = :best, kwargs...)

Clustering using K-barycenters. If you want to cluster spectra, consider the method that accepts models instead.

# Arguments:
- `X`: Support of input measures, Vector{Matrix} where each matrix is n_dims×n_atoms
- `p`: Weights of input measures Vector{Vector} where each matrix is of length n_atoms and should sum to 1.
- `k`: number of clusters
- `seed`: :rand or :eq
- `kiters`: number of iterations
- `verbose`: print stuff?
- `output`: output lowest cost clustering or :latest?
- `kwargs`: are sent to the inner solvers.
"""
function kbarycenters(X,p,k; seed=:rand, kiters=10, verbose=false, output=:best, kwargs...)
    N = length(X)
    @assert length(p) == N
    @assert k < N "number of clusters must be smaller than number of points"
    if seed === :rand
        Q = perturb!.(copy.(X[randperm(N)[1:k]]), Ref(X))
    elseif seed === :eq
        Q = perturb!.(copy.(X[1:N÷k:end]), Ref(X))
    else
        throw(ArgumentError("Unknown symbol for seed: $seed"))
    end
    C = distmat_euclidean(Q[1],Q[1])
    q = ones(size(X[1],2)) |> s1
    λ = ones(N) |> s1
    ass, cost = assignments(C,X,Q,p,q;kwargs...)
    ass_old = copy(ass)
    verbose && @info "kbarycenters: Iter: 0 cost: $(cost)"
    bestcost = cost
    bestass = ass
    for iter = 1:kiters
        # @show iter
        Q = barycenters(C,X,p,q,λ,ass,k;kwargs...)
        ass, cost = assignments(C,X,Q,p,q;kwargs...)
        # @show ass
        if cost < bestcost
            bestcost = cost
            bestass = ass
        end
        if ass == ass_old
            verbose && @info "kbarycenters: Iter: $iter converged"
            break
        end
        verbose && @info "kbarycenters: Iter: $iter cost: $(cost) num changes: $(count(ass .!= ass_old))"
        ass_old = ass
        yield()
    end
    if ass != bestass && output === :best
        verbose && @info "kbarycenters: Best cost: $bestcost"
        cost = bestcost
        Q = barycenters(C,X,p,q,λ,bestass,k;kwargs...)
    end
    KBResult(ass, Q, cost)
end

function barycenters(C,X,p,q,λ,ass,k;kwargs...)
    N = length(X)
    unnull!(ass,k)
    Q = tmap(1:k) do i
        inds = ass .== i
        if count(inds) == 0
            @warn "null cluster"
            inds = randperm(N)[1:2]
        end
        barycenter(X[inds],p[inds], s1(λ[inds]); kwargs...)[1] # TODO: not considering the weights
    end
end

function kwcostfun(C,X,Q,p,q; solver=sinkhorn_log!, kwargs...)
    C = distmat_euclidean!(C, X,Q)
    sum(solver(C, p, q; kwargs...)[1] .* C)
end

function assignments(C,X,Q,p,q;kwargs...)
    k = length(Q)
    c = 0.0
    ass = map(1:length(X)) do i
        dists = map(1:k) do j
            kwcostfun(C,X[i],Q[j],p[i],q;kwargs...)
        end
        ind = argmin(dists)
        c += dists[ind]
        ind
    end
    ass, c
end

function unnull!(ass,k)
    N = length(ass)
    change = true
    while change
        change = false
        for i = 1:k
            if i ∉ ass
                change = true
                ass[randperm(N)[1:1]] .= i
            end
        end
    end
end

"""
    kbarycenters(d::OptimalTransportRootDistance, models::Vector{<:AbstractModel}, k; normalize = true, kwargs...)

This function is only available if `using Clustering`.

# Example:
```julia
clusterresult = kbarycenters(
    dist,
    models,
    k,         # number of clusters
    seed       = :rand,
    solver     = sinkhorn_log!,
    tol        = 2e-6,
    innertol   = 2e-6,
    iters      = 100000,
    inneriters = 100000,
    verbose    = true,
    output     = :best,
    uniform    = true,
    kiters     = 10
)
```
The docs contain [a more detailed example](https://baggepinnen.github.io/SpectralDistances.jl/latest/interpolations/#K-Barycenters-1)

# Arguments:
- `models`: A vector of models
- `k`: number of clusters
- `normalize`: Whether or not to normalize the weight vectors (recommended)
- `kwargs`: are sent to inner solvers, (`solver,tol,iters`)
"""
function kbarycenters(d::OptimalTransportRootDistance, models::Vector{<:AbstractModel}, k; normalize=true, kwargs...)
    d.p == 2 || throw(ArgumentError("p must be 2"))
    X, w, realpoles = barycenter_matrices(d, models, normalize)
    res = kbarycenters(X, w, k; β=d.β, kwargs...)
    KBResult(res.assignments, bc2model.(res.barycenters,realpoles), res.cost)
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.3'
#       jupytext_version: 0.8.6
#   kernel_info:
#     name: julia-1.0
#   kernelspec:
#     display_name: Julia 1.0.0
#     language: julia
#     name: julia-1.0
# ---

# +
using Pkg
Pkg.activate(".")
using SemanticModels
using SemanticModels.Unitful: DomainError, s, d, C, uconvert, NoUnits
using DifferentialEquations
using DataFrames
using Unitful
using Test

using Distributions: Uniform
using GLM
using DataFrames

# -

using Plots

stripunits(x) = uconvert(NoUnits, x)

# +
function flusim(tfinal)
    # annual cycle of temperature control flu infectiousness
    springmodel = SpringModel([u"(1.0/(365*8))d^-2"], # parameters (frequency)
                              (u"0d",tfinal), # time domain
                              [u"25.0C", u"0C/d"]) # initial_conditions T, T'
    function create_sir(m, solns)
        sol = solns[1]
        initialS = u"10000person"
        initialI = u"1person"
        initialpop = [initialS, initialI, u"0.0person"]
        β = u"1.0/18"/u"d*C" * sol(sol.t[end-2])[1] #infectiousness
        @show β
        sirprob = SIRSimulation(initialpop, #initial_conditions S,I,R
                                (u"0.0d", u"20d"), #time domain
                                SIRParams(β, u"40.0person/d")) # parameters β, γ
        return sirprob
    end

    function create_flu(cm, solns)
        sol = solns[1]
        finalI = stripunits(sol(u"8.0d")[2]) # X
        population = stripunits(sol(sol.t[end])[2])
        # population = stripunits(sum(sol.u[end]))
        df = SemanticModels.generate_synthetic_data(population, 0,100)
        f = @formula(vaccines_produced ~ flu_patients)
        model =  lm(f,
            df[2:length(df.year),
            [:year, :flu_patients, :vaccines_produced]])
        println("GLM Model:")
        println(model)

        year_to_predict = 1
        num_flu_patients_from_sim = finalI
        vaccines_produced = missing
        targetDF = DataFrame(year=year_to_predict,
            flu_patients=num_flu_patients_from_sim,
            vaccines_produced=missing)
        @show targetDF

        return RegressionProblem(f, model, targetDF, missing)
    end
    cm = CombinedModel([springmodel], create_sir)
    flumodel = CombinedModel([cm], create_flu)
    return flumodel
end

tfinal = 240π*u"d" #(~2 yrs)
flumodel = flusim(tfinal)

# -

springmodel = flumodel.deps[1].deps[1]
sirmodel = flumodel.deps[1]
sol = solve(springmodel)
plot(sol.t./d, map(x->x[1], sol.u) ./ C)

sirsol = solve(sirmodel)

plot(sirsol.t./d,map(x->stripunits.(x)[2], sirsol.u))

sol = solve(flumodel)

print(typeof(sol.u))


text:
include("../../src/SPICE.jl")
using .SPICE

# ---- State-change matrix
const ν = Array([1 -1 1 -1]')

# ---- Initial state-space
x0 = [250]

# ---- Highest orders of reaction for tau-leaping
const hor = [33]

# ---- Bounds on the initial search space for parameters
# ---- Real parameters are [3e-7,1e-4,1e-3,3.5]
const bounds = [[1e-9,1e-6,1e-5,1e-2],[1e-5, 1e-2, 1e-1, 100.0]]

function F(p)
	H(p)
	p.a .*= p.θ
end

# ---- Hazard functions
function H(p)
    p.a[1] = p.x[1] * (p.x[1]-1) / 2 * 100000
    p.a[2] = p.x[1] * (p.x[1]-1) * (p.x[1] - 2) / 6
    p.a[3] = 200000
    p.a[4] = p.x[1]
end

# how many times to run?
system = System(Model(x0,F,H, ν,bounds,hor=hor,obsname=[:X]), "./data/10",routine=CEM(ssa=:Tau, nElite = 10, nRepeat=10, nSamples=1000, maxIter=250, mSamples=20000, shoot=false, splitting=false, tauOpt=TauOpt(ϵ=0.1)))

estimate(system, 100, "schlogl")

max_stars_repo_name: pzuliani/SPICE

max_issues_repo_name: pzuliani/SPICE

max_forks_repo_name: pzuliani/SPICE


text:
__precompile__()

module GLM
    using Reexport
    @reexport using Distributions
    using Base.LinAlg.LAPACK: potrf!, potrs!
    using Base.LinAlg.BLAS: gemm!, gemv!
    using Base.LinAlg: QRCompactWY, Cholesky, BlasReal
    using StatsBase: StatsBase, CoefTable, StatisticalModel, RegressionModel
    using StatsFuns: logit, logistic
    using Distributions: sqrt2, sqrt2π
    using Compat

    import Base: (\), cholfact, convert, cor, show, size
    import StatsBase: coef, coeftable, confint, deviance, nulldeviance, dof, dof_residual, loglikelihood, nullloglikelihood, nobs, stderr, vcov, residuals, predict, fit, model_response, r2, r², adjr2, adjr², PValue
    import StatsFuns: xlogy
    import SpecialFunctions: erfc, erfcinv
    export coef, coeftable, confint, deviance, nulldeviance, dof, dof_residual, loglikelihood, nobs, stderr, vcov, residuals, predict, fit, fit!, model_response, r2, r², adjr2, adjr²

    export                              # types
        CauchitLink,
        CloglogLink,
        DensePred,
        DensePredQR,
        DensePredChol,
        GeneralizedLinearModel,
        GlmResp,
        IdentityLink,
        InverseLink,
        InverseSquareLink,
        LinearModel,
        Link,
        LinPred,
        LinPredModel,
        LogitLink,
        LogLink,
        LmResp,
        ProbitLink,
        SqrtLink,

                                        # functions
        canonicallink,  # canonical link function for a distribution
        delbeta!,       # evaluate the increment in the coefficient vector
        deviance,       # deviance of fitted and observed responses
        devresid,       # vector of squared deviance residuals
        formula,        # extract the formula from a model
        glm,            # general interface
        glmvar,         # the variance function
        inverselink,    # returns μ, dμ/dη and, when appropriate, μ*(1-μ)
        linkfun,        # link function mapping mu to eta, the linear predictor
        linkinv,        # inverse link mapping eta to mu
        linpred,        # linear predictor
        linpred!,       # update the linear predictor in place
        lm,             # linear model
        logistic,
        logit,
        mueta,          # derivative of inverse link
        mustart,        # derive starting values for the mu vector
        nobs,           # total number of observations
        predict,        # make predictions
        updateμ!,       # update the response type from the linear predictor
        wrkresp,        # working response
        ftest           # compare models with an F test

    const FP = AbstractFloat
    const FPVector{T<:FP} = AbstractArray{T,1}

    abstract type ModResp end                         # model response

    abstract type LinPred end                         # linear predictor in statistical models
    abstract type DensePred <: LinPred end            # linear predictor with dense X
    abstract type LinPredModel <: RegressionModel end # model based on a linear predictor

    include("linpred.jl")
    include("lm.jl")
    include("glmtools.jl")
    include("glmfit.jl")
    include("ftest.jl")

end # module


text:
export Study, getStudy, getStudies, validate, enforceStudy, changeStudy

# The following are base types describing
# the dataset store at a certain level
struct Study{D<:DatasetStore}
  store::D
  name::String
  foldername::String
  date::DateTime
  subject::String
end

function Study(store::DatasetStore, name::String; foldername::String="",
         date::DateTime = now(), subject::String="", createDir=true)
  #remove Milliseconds
  date_ = trunc(date, Dates.Second)
  if foldername == ""
    foldername = getMDFStudyFolderName(name,date_)
  end
  path = joinpath( studydir(store), foldername )
  if createDir
    mkpath(path)
  end
    return Study(store, name, foldername, date_, subject)
end

path(s::Study) = joinpath( studydir(s.store), s.foldername )
@mustimplement path(s::Study{D}, numExp::Integer) where D<:DatasetStore

function getStudies(d::DatasetStore, name::String="")
  s = Study[]

  files = readdir( studydir(d) )
  for file in files
    fullpath = joinpath(studydir(d),file)
    if isdir(fullpath) && !ishidden(fullpath)
      try
      study = getStudy(d, file)
      # return only studies where name is a substring of the studyname
      if study !== nothing && match(Regex("(?i)"*name), study.name) !== nothing
        push!(s, study)
      end
      catch e
        @warn e
      end
    end
  end
  return s
end

function remove(study::Study)
  if isdir(path(study))
    rm(path(study), recursive=true)

    #TODO remove recos!
  end
end

function getNewExperimentNum(s::Study)
  return getNewNumInFolder(path(s))
end

function getNewExperimentPath(s::Study)
  addStudy(s.store, s)
  expNum = getNewExperimentNum(s)
  path_ = joinpath(path(s),string(expNum)*".mdf")
  # touch new mdf file
  touch(path_)
  try_chmod(path_, 0o660)
  return path_
end

function changeParam(s::Study, paramName::AbstractString, paramValue)
  exps = getExperiments(s)
  for e in exps
    changeParam(e, paramName, paramValue)
  end
end

function changeStudy(s::Study, newStudyName::AbstractString; date::DateTime = now())
  changeParam(s, "/study/name", newStudyName)
  changeParam(s, "/study/time", string(date))
  sNew = Study(s.store, newStudyName, date=date, subject=s.subject, createDir=false)
  mv(path(s), path(sNew))
  return
end

function validate(s::Study)
  exps = getExperiments(s)
  valid = true
  for e in exps
    MPIFile(path(e), fastMode=true) do f
      date1 = trunc(studyTime(f), Dates.Second)
      date2 = trunc(s.date, Dates.Second)
      if studyName(f) != s.name || date1 != date2
        valid = false
        @info "file $path(e) is not valid"
        @show studyName(f), s.name, date1, date2
      #else
      #  @info "file $path(e) is valid"
      end
    end
  end
  return valid
end

function enforceStudy(s::Study)
  changeParam(s, "/study/name", s.name)
  changeParam(s, "/study/time", string(s.date))
  return
end