text:
using Plots
using Random
Random.seed!(1)

include_ddp()

# Model
include_model("double_integrator")

function f(model::DoubleIntegratorContinuous, x, u, w)
    [x[2] + w[1]; (1.0 + w[3]) * u[1] + w[2]]
end

function fd(model::DoubleIntegratorContinuous{Midpoint, FixedTime}, x, u, w, h, t)
	if t == 1
		x1 = u[2:3]
    	return x1 + h * f(model, x1 + 0.5 * h * f(model, x1, u, w), u, w)
	else
		return x + h * f(model, x + 0.5 * h * f(model, x, u, w), u, w)
	end
end

model = DoubleIntegratorContinuous{Midpoint, FixedTime}(2, 1, 3)
n = [model.n for t = 1:T]
m = [t == 1 ? model.m + model.n : model.m for t = 1:T]

# Time
T = 101
h = 0.1

z = range(0.0, stop = 3.0 * 2.0 * π, length = T)
p_ref = 1.0 * cos.(1.0 * z)
plot(z, p_ref)

# Initial conditions, controls, disturbances
x1 = [p_ref[1]; 0.0]
x1_alt = [0.0; 0.0]
xT = [[p_ref[t]; 0.0] for t = 1:T]
ū = [rand(m[t]) for t = 1:T-1]
w = [zeros(model.d) for t = 1:T-1]

# Rollout
x̄ = rollout(model, x1, ū, w, h, T)

# Objective
Q = [(t < T ? h : 1.0) * Diagonal([100.0; 0.1]) for t = 1:T]
q = [-2.0 * Q[t] * xT[t] for t = 1:T]
R = [h * Diagonal(0.01 * ones(m[t])) for t = 1:T-1]
r = [zeros(m[t]) for t = 1:T-1]
obj = StageCosts([QuadraticCost(Q[t], q[t],
	t < T ? R[t] : nothing, t < T ? r[t] : nothing) for t = 1:T], T)

function g(obj::StageCosts, x, u, t)
	T = obj.T
    if t < T
		Q = obj.cost[t].Q
		q = obj.cost[t].q
	    R = obj.cost[t].R
		r = obj.cost[t].r
        return x' * Q * x + q' * x + u' * R * u + r' * u
    elseif t == T
		Q = obj.cost[T].Q
		q = obj.cost[T].q
        return x' * Q * x + q' * x
    else
        return 0.0
    end
end

g(obj, x̄[1], ū[1], 1)
g(obj, x̄[2], ū[2], 2)
g(obj, x̄[T], nothing, T)
objective(obj, x̄, ū)

# Constraints
p = [t == 1 ? n[1] : 0 for t = 1:T]
info_t = Dict(:x1 => x1_alt)#:ul => [-5.0], :uu => [5.0], :inequality => (1:2 * m))
info_T = Dict()#:xT => xT)
con_set = [StageConstraint(p[t], t < T ? info_t : info_T) for t = 1:T]

function c!(c, cons::StageConstraints, x, u, t)
	T = cons.T
	if t == 1
		c .= view(u, 2:3) - cons.con[t].info[:x1]
	end
end

prob = problem_data(model, obj, con_set, copy(x̄), copy(ū), w, h, T, n = n, m = m)

# Solve
@time constrained_ddp_solve!(prob,
    max_iter = 1000, max_al_iter = 6,
	ρ_init = 1.0, ρ_scale = 10.0)

# # Solve
# @time ddp_solve!(prob,
#     max_iter = 100, verbose = true)

x, u = current_trajectory(prob)
x̄, ū = nominal_trajectory(prob)

x[1] = u[1][2:3]
x̄[1] = ū[1][2:3]

# Visualize
using Plots
plot(hcat([[p_ref[t]; 0.0] for t = 1:T]...)',
    width = 2.0, color = :black, label = "")
plot!(hcat(x...)[:, 1:T]', color = :magenta, label = "")
# plot(hcat(u..., u[end])', linetype = :steppost)


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

# Compute X <- a X + b Y.
function matfun_axpby!(X,a,b,Y)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:m
        @inbounds for j=1:n
            if (b isa ValueOne)
                X[i,j]+=Y[i,j]
            else
                X[i,j]+=b*Y[i,j]
            end
        end
    end
end

@inline function exp_native_jl_rho0_25(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_native_jl_rho0_25!(A_copy)
end

@inline function exp_native_jl_rho0_25!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=5
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    value_one=ValueOne()
    # Computation order: A2 A4 Ua U V Z X P
    # Computing A2 with operation: mult
    mul!(memslots2,memslots1,memslots1)
    # Computing A4 with operation: mult
    mul!(memslots3,memslots2,memslots2)
    # Computing Ua = x*I+x*A2+x*A4
    coeff1=15120.0
    coeff2=420.0
    coeff3=1.0
    memslots4 .= coeff2.*memslots2 .+ coeff3.*memslots3
    mul!(memslots4,true,I*coeff1,true,true)
    # Computing U with operation: mult
    mul!(memslots5,memslots4,memslots1)
    # Deallocating Ua in slot 4
    # Deallocating A in slot 1
    # Computing V = x*I+x*A2+x*A4
    coeff1=30240.0
    coeff2=3360.0
    coeff3=30.0
    # Smart lincomb recycle A2
    memslots2 .= coeff2.*memslots2 .+ coeff3.*memslots3
    mul!(memslots2,true,I*coeff1,true,true)
    # Deallocating A4 in slot 3
    # Computing Z = x*V+x*U
    coeff1=1.0
    coeff2=-1.0
    memslots1 .= coeff1.*memslots2 .+ coeff2.*memslots5
    # Computing X = x*V+x*U
    coeff1=1.0
    coeff2=1.0
    # Smart lincomb recycle V
    memslots2 .= coeff1.*memslots2 .+ coeff2.*memslots5
    # Deallocating U in slot 5
    # Computing P with operation: ldiv
    memslots3 .=memslots1\memslots2
    # Deallocating Z in slot 1
    # Deallocating X in slot 2
    return memslots3 # Returning P
end

ext: jl

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

include("default_initval_algorithm.jl")
include("default_transform_algorithm.jl")
include("default_sampling_algorithm.jl")
include("default_autocor_len.jl")
include("default_eff_sample_size.jl")
include("default_mode_estimator.jl")
include("default_differentiation_algorithm.jl")
include("default_integration_algorithm.jl")


text:
#!/usr/bin/env julia

using DataFrames

import Celeste.AccuracyBenchmark
import Celeste.ArgumentParse
import Celeste: Config
import Celeste.Infer
import Celeste.ParallelRun
import Celeste.SDSSIO

const OUTPUT_DIRECTORY = joinpath(splitdir(Base.source_path())[1], "output")

parser = Celeste.ArgumentParse.ArgumentParser()
ArgumentParse.add_argument(parser, "--joint", help="Use joint inference", action=:store_true)
ArgumentParse.add_argument(
    parser,
    "--use-full-initialization",
    help="Use all information from initialization catalog (otherwise will just use noisy position",
    action=:store_true,
)
ArgumentParse.add_argument(
    parser,
    "--limit-num-sources",
    help="Target only the given number of sources, for quicker testing",
    arg_type=Int,
)
ArgumentParse.add_argument(
    parser,
    "--min-radius-px",
    help="Override Config.min_radius_px with the given value",
    arg_type=Float64,
)
ArgumentParse.add_argument(
    parser,
    "--image-fits",
    help="FITS file containing synthetic imagery; if not specified, will use SDSS RCF",
)
ArgumentParse.add_argument(
    parser,
    "--run",
    help="SDSS run #",
    arg_type=Int,
    default=AccuracyBenchmark.STRIPE82_RCF.run,
)
ArgumentParse.add_argument(
    parser,
    "--camcol",
    help="SDSS camcol #",
    arg_type=Int,
    default=AccuracyBenchmark.STRIPE82_RCF.camcol,
)
ArgumentParse.add_argument(
    parser,
    "--field",
    help="SDSS field #",
    arg_type=Int,
    default=AccuracyBenchmark.STRIPE82_RCF.field,
)
ArgumentParse.add_argument(
    parser,
    "catalog_csv",
    help="CSV catalog for initialization",
    required=true,
)
parsed_args = ArgumentParse.parse_args(parser, ARGS)

srand(12345)

run_camcol_field = nothing
if haskey(parsed_args, "image-fits")
    extensions = AccuracyBenchmark.read_fits(parsed_args["image-fits"])
    images = AccuracyBenchmark.make_images(extensions)
    catalog_label = splitext(basename(parsed_args["image-fits"]))[1]
else
    rcf = SDSSIO.RunCamcolField(
        parsed_args["run"],
        parsed_args["camcol"],
        parsed_args["field"],
    )
    strategy = SDSSIO.PlainFITSStrategy(AccuracyBenchmark.SDSS_DATA_DIR)
    images = SDSSIO.load_field_images(strategy, [rcf])
    catalog_label = @sprintf("sdss_%s_%s_%s", rcf.run, rcf.camcol, rcf.field)
end
@assert length(images) == 5

catalog_data = AccuracyBenchmark.read_catalog(parsed_args["catalog_csv"])
if parsed_args["use-full-initialization"]
    @printf("Using full initialization from %s\n", parsed_args["catalog_csv"])
end
catalog_entries = AccuracyBenchmark.make_initialization_catalog(
    catalog_data,
    parsed_args["use-full-initialization"],
)
@printf("Loaded %d sources...\n", length(catalog_entries))

if haskey(parsed_args, "limit-num-sources")
    target_sources = collect(1:parsed_args["limit-num-sources"])
else
    target_sources = collect(1:length(catalog_entries))
end
neighbor_map = Infer.find_neighbors(target_sources, catalog_entries, images)

config = Config()
if haskey(parsed_args, "min-radius-px")
    config.min_radius_pix = parsed_args["min-radius-px"]
end
results = AccuracyBenchmark.run_celeste(
    config,
    catalog_entries,
    target_sources,
    images,
    use_joint_inference=parsed_args["joint"],
)
results_df = AccuracyBenchmark.celeste_to_df(results)

# save results in both JLD and CSV format

jld_filename = joinpath(
    OUTPUT_DIRECTORY,
    @sprintf(
        "celeste-%.4f-%.4f-%.4f-%.4f.jld",
        minimum(results_df[:right_ascension_deg]),
        maximum(results_df[:right_ascension_deg]),
        minimum(results_df[:declination_deg]),
        maximum(results_df[:declination_deg]),
    )
)
println("Writing $(jld_filename)...")
JLD.save(jld_filename, "results", results)

csv_filename = joinpath(OUTPUT_DIRECTORY, @sprintf("%s_predictions.csv", catalog_label))
AccuracyBenchmark.write_catalog(csv_filename, results_df)
AccuracyBenchmark.append_hash_to_file(csv_filename)


text:
using LinearAlgebra, ForwardDiff, SparseArrays, Optim, LineSearches
# include("al.jl")
prob = prob
function f(x)
    MOI.eval_objective(prob, x)
end

function g!(G, x)
    MOI.eval_objective_gradient(prob, G, x)
    nothing
end

function c!(c, x)
    MOI.eval_constraint(prob, c, x)
    c .*= -1.0
    nothing
end

spar = sparsity_jacobian(prob)
global jac = zeros(length(spar))
global ii = [s[1] for s in spar]
global jj = [s[2] for s in spar]
function d!(D, x)
    MOI.eval_constraint_jacobian(prob, jac, x)
    D .= sparse(ii, jj, -1.0 .* jac)
end

function f_al(x, al::AugmentedLagrangian)
    # evaluate constraints
    c!(al.c, x)
    bounds!(al, x)
    active_set_update!(al)

    # compute objective
    J = f(x)

    # add augmented Lagrangian terms
    J += al.λ' * al.c + 0.5 * sum(al.as .* al.ρ .* (al.c.^2.0))
    J += al.λl' * al.cl + 0.5 * sum(al.asl .* al.ρl .* (al.cl.^2.0))
    J += al.λu' * al.cu + 0.5 * sum(al.asu .* al.ρu .* (al.cu.^2.0))
end

function g_al!(G, x, al::AugmentedLagrangian)
    # compute objective gradient
    g!(G, x)

    # evaluate constraints
    # ForwardDiff.jacobian!(al.∇c, c!, al.c, x)
    c!(al.c, x)
    d!(al.∇c, x)

    bounds!(al, x)
    active_set_update!(al)

    # add augmented Lagrangian gradient terms
    G .+= al.∇c' * (al.λ + al.as .* al.ρ .* al.c)
    G[al.idx_l] -= (al.λl + al.asl .* al.ρl .* al.cl)
    G[al.idx_u] += (al.λu + al.asu .* al.ρu .* al.cu)
    return nothing
end

function solve(x, al; alg = :LBFGS, max_iter = 5, c_tol = 1.0e-3)
    # reset augmented Lagrangian
    reset!(al)
    println("solving...")
    for i = 1:max_iter
        # x̄, ū = unpack(x, prob)
        #
        # visualize!(vis, model, state_to_configuration(x̄), Δt = h)

        # update augmented Lagrangian methods
        _f(z) = f_al(z, al)
        _g!(G, z) = g_al!(G, z, al)

        # solve
        opt = Optim.Options()
                # x_tol = 1.0e-6,
                # f_tol = 1.0e-6,
                # g_tol = 1.0e-2)
                # iterations = convert(Int,1e4))

        sol = optimize(_f, _g!, x,
                        eval(alg)(),#linesearch = LineSearches.BackTracking()),
                        opt)

        # evaluate constraints
        x = sol.minimizer
        c_max = constraint_violation(al, x)
        println("iter: $i")
        println("c_max: $c_max")

        # check for convergence -> update augmented Lagrangian
        if c_max < c_tol
            return x, sol
        else
            c!(al.c, x)
            bounds!(al, x)
            update!(al)
        end
    end

    return x, sol
end

n = prob.num_var
m = prob.num_con
xl, xu = prob.primal_bounds
cl, cu = prob.constraint_bounds
idx_ineq = (1:m)[cu .> cl]
sum(isfinite.(xl))
sum(isfinite.(xu))
sum(cu .> cl)

al = augmented_lagrangian(n, m,
    xl = xl, xu = xu, ρ0 = 10.0, s = 10.0,
    idx_ineq = idx_ineq)

@time x_sol_al, sol , info = solve(copy(z0), al,
    alg = :MomentumGradientDescent, max_iter = 20, c_tol = 1.0e-2)

# Visualize
# using Plots
# x̄, ū = unpack(x_sol_al, prob)
# visualize!(vis, model, state_to_configuration(x̄), Δt = h)

# plot(hcat(x̄...)', width = 2.0)
# plot(hcat(ū...)[1:model.m, :]', width = 2.0, linetype = :steppost)
# s̄ = [ū[t][model.idx_s] for t = 1:T-1]
# plot(hcat(s̄...)', width = 2.0, linetype = :steppost)


text:
@testset "EEMD" begin
    @testset "EEMD Test 1" begin
        n = 200
        s = EEMDSetting(n, 50, 10, 2, 10, 0.2, 450)

        t = collect(range(0,1, length=n))
        x = zeros(n)

        for i=1:n
            x[i] = cos(22*pi*t[i]^2) + 6*t[i]^2
        end

        output = eemd(x, s)

        @test size(output, 2) == s.emd_setting.m

        recovered = zeros(n)

        for i=1:size(output, 2)
            recovered += output[:, i]
        end

        @test recovered == x || isapprox(recovered, x, rtol=0.05)
    end

    @testset "EEMD Test 2" begin
        N = 1024
        s = EEMDSetting(N, 50, 10, 10, 250, 0.3, 200)

        function input_signal(x::Float64)
            omega = 2*pi/(N-1)
            return sin(17*omega*x)+0.5*(1.0-exp(-0.002*x))*sin(51*omega*x+1)
        end

        input = zeros(N)

        for i=1:N
            input[i] = input_signal(float(i))
        end

        output = eemd(input, s)

        @test size(output, 2) == s.emd_setting.m

        recovered = zeros(N)

        for i=1:size(output, 2)
            recovered += output[:, i]
        end

        @test isapprox(recovered,input, rtol=0.03)
    end
end


text:
#!/usr/bin/env bash
    #=
    exec julia --project="$(realpath $(dirname $0))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \
    "${BASH_SOURCE[0]}" "$@"
    =#

include(joinpath(dirname(dirname(@__FILE__)), "src", "ComputabilityTheory.jl"))

using .ComputabilityTheory
using Plots
using CSV
using DataFrames
using ProgressMeter: @showprogress
using CurveFit
using GLM

function graphing_time()
    stop_m_at = 100
    stop_n_at = 4
    num_of_datapoints = 50
    data = Matrix{Float64}(undef,num_of_datapoints,4)# needs 5 columns: n; m; bonk time; doov time; iagueqnar time
    m_data = []

    @showprogress for i in 1:num_of_datapoints
        n = rand(2:stop_n_at)
        m = rand(1:stop_m_at)
        bonks = @elapsed π(m, n)
        alg = @elapsed π(m, n, algebraic)
        data[i,:] .= [n, m, bonks, alg]
    end

    save_path = joinpath(dirname(@__FILE__), "unpairing", "search_v_algebraic,n=$stop_n_at,m=$stop_m_at,i=$num_of_datapoints.pdf")
    theme(:solarized)

    println(typeof(data[:,1]))
    println(typeof(data[:,2]))

    fit = curve_fit(ExpFit, data[:,1], data[:,2])
    y0b = fit.(data[:,1])

    plot_points = scatter(data[:,1], data[:,3:4], fontfamily = font("Times"), xlabel = "n", ylabel = "Time elapsed during calculation [seconds]", label = ["Bonk's Iterative Search" "Bonk's Algebraic Method"])#, xlims = (0, stop_n_at)) # smooth = true
    # plot = plot!(plot_points, model(1:stop_m_at, params_search))

    # plot = plot!(plot_points, data[:,1], data[:,2]
    plot = plot!(data[:,1], y0b, label = "model")
    savefig(plot, save_path)

    println("Plot saved at $save_path.")
end

function programme_ui()

	programmed_programmes = [
		121, # nothing
		5780, # increments R0
		363183787614755732766753446033240, # R0 + R1
		15064444022569784320075110954333881106157772192293890937642519911533240246282230807658733450927595112277856126952898743268714565377889514636997, # 2x
		70412762139173751174325247391799151776337005517147051043631834182952422253779831226208595923309323737831968573315841128818823056224690660136852, # halts ⟺ R0 > R1
		972292871301644916468488152875266508938968846389326007980307063346008398713128885682044504108288931767348821063618087715644933567266540511345568504718733339523678538338052787779884557674350959673597803113281693069940562881722205193604550737455583875504348606989700013337656597740101535, # relation "x is even"
		223927124717252638681796127568111496742698512180714163023526568831822132049774230195118493600338381851704911289089017787857413597404988541338847992210485499035130781598573348126975813033599979584118468386578559370872281223796019975651453204160222340161200928955988080303016030336930443304245145601448154569733221264381538305168568931449480758528673246332402357557389885850502823987502500615923705666157574953869234044687656342114443578201036399996781757409020565065471713601342262949677684562201881239545438799475865994538241042725394500589505167297761345293434717297622500387774051925726286270828011593028154888482643933182720230237091007716467114407940246910876104242298811097230924566373379040422996695018496814405564315525869126758556165426204979060785448299907433924693776629082547220185023385651497716000521725878783901307004251127081977484444374187001323436144140037580903329365965139055311902469340438083050750301104320244848087790925135435440549816132311216931423038414958594054194345955328449597351606723033103611762253425583323756613052352836016212463541320408142893340287911691519461708609068873926465524275772411714118028019844718904105500330674632838771428681137130136100546871186789563992080600787879878511406660580151025959186841418554827311891320184308989618516409485943704466045877998837588192158884659614365477933514118459878126838662936672543012162679280965785850875755241288227892173373302073145637378088093219495972690064022128760997033867243024673762202229065063901594672092113076799336225728489852391219624477091216808040617586878843148705687256752473769380170895833165256985381066579262059980097322456692517960764426357076378502591620840598428480692392558582556549972917361099374457224379447165159150577431975140314355645616431652416541261437750090494682945620970384916782903782096424825827978169475334019068642656156261072845169384253772103615589882161531092319774765872915674167557806826727743079891261233740545808889798017329652560892724042630851029267797497619411669453984450602449597986084188806823478595595515128790674171852156626442453849830854084257473658507273475076486871236759558604076795418209173234381881185600648269038607854092549797177142124163026021172162187852101251819011578114092027496262838478199601979508935066069585104555431366780222617619635159394639488403252338359673576720653539411491577093278404645966204761032796579698188013739359885785188427988372869573092916550596059755937168888105513627640203861038587598828520438054944864945385514172661237373005642977222691283341092758484095459330166556676864303459850601868153088510594207374543818165331580579808638399043608636650159746931808355843721671440780553546753959760350975906249580764849993052921544867678610217996851589828713473915270585295877597236593571404899502354444405479601753716517182781778676830983428681476532848100787997408528892345104613523387370533475767447648602910646760355660958540505065155007961785878896275786367616723562594625648464444056163755337577182731373365691826878791057479596578563844746371661, # relation "x = 3"
	]

	global programme_counter = -1
	programme_description(description::AbstractString) = (global programme_counter; programme_counter += 1; println("\t$programme_counter) $description"))

	println("The following are the list of programmes to choose from:")

	programme_description("Any ol' programme will do!")
	programme_description("A halting programme")
	programme_description("Increments R0 by one")
	programme_description("Computes R0 + R1")
	programme_description("Computes 2 × R1")
	programme_description("Halts if and only if R0 > R1")
	programme_description("Computes the relation \"x is even\"")
	programme_description("Computes the relation \"x = 3\"")

	println()
	println("Please choose a number as defined above.")

	chosen_programme = parse(Int, readline())

	println("-"^displaysize(stdout)[2], "\n")

	iszero(chosen_programme) ? test_random(rand(1:20)) : show_programme(programmed_programmes[chosen_programme])

end

# helper
function binary_search(a::BigInt, b::BigInt) # useful for bug fixing
    while a != b
        centre = BigInt(round((b + a) / 2))
        println(a, "\n", b, "\n")
        try
            π(centre, algebraic)
            a = centre
        catch
            b = centre
        end
    end
    return a
end

# programme_ui()


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   ModiaMath.Variables (ModiaMath/Variables/_module.jl)
#

# Definitions of Real variables

#=
The dynamic system is described as index-1 DAE:

    wr := wr(der(x),x,t)      # Additional variables, that are explicitly computed but are not reported to the integrator (if used in fc(...), wr is not allowed to depend on der(x)
    0   = fd(der(x),x,t,wr)
    0   = fc(x,t,wr)
    wc := wc(der(x),x,t,wr)   # Additional variables, that need to be only computed at communication points.

such that (assuming wr is inlined in fd/fc)

   [der(fd, der(x));
    der(fc, x)     ] is regular

The vectors are partitioned as:

   x      = [xd_exp;        # XD_EXP: variables appearing differentiated in the model and the model explicitly computes der(xd_exp), so der(xd_exp) is an output of the model (used in model)
             xd_imp;        # XD_IMP: variables appearing differentiated in the model and the model does not compute der(xd_imp), so der(xd_imp) is an input to the model (used in model)
             xa;            # XA    : variables not appearing differentiated in the model, so algebraic variables (used in model)
             lambda_int;    #         integral of index-1 algebraic variables (not used in model; used only by integrator)
             mue_int   ]    #         integral of stabilizing variables der(mue_int) (not used in model; used only by integrator)

   der(x) = [der(xd_exp);   # DER_XD_EXP: derivative of xd_exp; is computed explicitly in the model (used in model).
             der(xd_imp);   # DER_XD_IMP: derivative of xd_imp; is provided by integrator (used in model)
             der(xa);       #             derivative of xa (not used in model; used only in integrator)
             lambda;        # LAMBDA    : index-1 algebraic variable; is provided by integrator (used in model)
             mue        ]   # MUE       : stabilizing variable to reduce higher index to index 1; is provided by integrator (used in model)

   residue = [fd_exp;       #           : residue of explicitly solved fd equations; computed by code below (e.g. if x[1] = phi and x[2] = v (= der(phi)), then residue[1] = derx[1] - v)
              fd_imp;       # FD_IMP    : residue of implicitly solved fd equations (is provided by model)
              fc]           # FC        : residue of fc equations (is provided by model)

If x consists only of xd_exp, then the DAE is an ODE and an ODE integration method can be used.

All variables x and der(x) are subtypes of type Real (typically Float64).
These variables are below defined as subtypes of type ContinuousTimeVariable

Types WR are auxiliary variables that are explicitly computed as functions of t,x,der(x) and are used to compute residues
(so these variables are computed at every call of the model)

Types WC are additional variables that are computed only at communication points.

Type TIME is just used to mark which variable is the independent time variable.
=#

isNothing(T) = typeof(T) == NOTHING
hasValue(T)  = typeof(T) != NOTHING

#------------------------------ Generic Real Variable ----------------------------------------
"""
    v = RealVariable{ValueType,ElementType}(...)

Generate a new variable with `v.value::ValueType`, `v.nominal::ElementType`.
The argument list is described in module [`Variables`](@ref).
"""
mutable struct RealVariable{ValueType,ElementType} <: ModiaMath.AbstractRealVariable
    # Generic Attributes of every AbstractVariable (partially based on FMI 2.0)
    _internal::ComponentInternal
    value::ValueType                 # Actual value of variable (can be a scalar or an array)
    info::AbstractString
    causality::Causality
    variability::Variability
    start::ValueType                 # Initial value of variable
    fixed::Bool                      # = true, start is fixed during initialization; = false, start is a guess value (can be changed during initialization)
    analysis::VariableAnalysisType   # Analysis for which the Variable is used

    # Attributes specific to Real variables
    min::ElementType
    max::ElementType
    nominal::ElementType                                              # nominal value; is used to compute absolute tolerances and might be used for scaling
    flow::Bool
    numericType::NumericType                                          # how the variable is used in the equations
    integral::Union{NOTHING,RealVariable{ValueType,ElementType}}  # if present, integral is the variable that represents the integral of the actual variable (so variable = d(integral)/dt     derivative::Union{NOTHING, RealVariable{ValueType, ElementType}}  # if present, derivative is the variable that represents the derivative of the actual variable (so derivative = d(variable)/dt)
    derivative::Union{NOTHING,RealVariable{ValueType,ElementType}}  # if present, derivative is the variable that represents the derivative of the actual variable (so derivative = d(variable)/dt)
    unit::String                                                      # unit of the variable (temporal solution until package Unitful is supported in Julia v0.7)

    # How the variable is stored in vectors
    ivar::Int        # if value is stored in x/derx/residue vector: x/derx/residue[ ivar ] = first value of variable.value
    iresult::Int     # if value is stored in result vector: result[ iresult ] = first value of variable.value

    function RealVariable{ValueType,ElementType}(name, within, info, start, unit, fixed, min, max, nominal, flow, causality, variability,
      numericType, integral, analysis) where {ValueType, ElementType}

        variable = new(ComponentInternal(name, within), deepcopy(start), info, causality, variability,
                     start, fixed, analysis, min, max, nominal, flow, numericType, integral, nothing, unit, 0, 0)

        if typeof(within) != NOTHING
            setfield!(within, Symbol(name), variable)
        end

        if hasValue(integral)
            integral.derivative = variable
        end

        return variable
    end
end

RealVariable{ValueType,ElementType}(name=NoNameDefined,
    within=nothing;
    info="",
    start=nothing,
    unit=NoUnit,
    fixed=false,
    min=-ElementType(Inf),
    max=ElementType(Inf),
    nominal=ElementType(1.0),
    flow=false,
    causality=Local,
    variability=Continuous,
    integral=nothing,
    numericType=WR,
    analysis=ModiaMath.AllAnalysis) where {ValueType, ElementType} =
        RealVariable{ValueType,ElementType}(name, within, info, start, unit, fixed, min, max, nominal, false, causality, variability, numericType, integral, analysis)

function Base.show(io::IO, r::RealVariable{ValueType,ElementType}) where {ValueType, ElementType}
    print(r.value, " ", r.unit)
    if r.info != ""
      # print(" # ", r.info)
   end
end

#------------------------------ Concrete Real Variable types ----------------------------------------

"""
    v = RealScalar(..)

Generate a variable of type [`RealVariable`](@ref)`{Float64,Float64}` to define scalar, Float64, real variables.
The argument list is described in module [`Variables`](@ref).
"""
const RealScalar = RealVariable{Float64,Float64}

"""
    v = RealSVector{Size}(..)

Generate a variable of type [`RealVariable`](@ref)`{SVector{Size,Float64},Float64}` where the values have
type StaticArrays.SVector{Size,Float64}.
The argument list is described in module [`Variables`](@ref).
"""
const RealSVector{Size} = RealVariable{SVector{Size,Float64},Float64}

"""
    v = RealSVector3(..)

Generate a variable of type [`RealVariable`](@ref)`{SVector{3,Float64},Float64}` where the values have
type StaticArrays.SVector{3,Float64}.
The argument list is described in module [`Variables`](@ref).
"""
const RealSVector3 = RealSVector{3}

RealScalar(name=NoNameDefined,
    within=nothing;
    info="",
    start=0.0,
    unit=NoUnit,
    fixed=false,
    min=-Inf,
    max=Inf,
    nominal=1.0,
    flow=false,
    causality=Local,
    variability=Continuous,
    integral=nothing,
    numericType=WR,
    analysis=ModiaMath.AllAnalysis) =
         RealScalar(name, within, info, start, unit, fixed, min, max, nominal, flow, causality, variability, numericType, integral, analysis)

RealSVector{Size}(name=NoNameDefined,
    within=nothing;
    info="",
    start=zeros(Size),
    unit=NoUnit,
    fixed=false,
    min=-Inf,
    max=Inf,
    nominal=1.0,
    flow=false,
    causality=Local,
    variability=Continuous,
    integral=nothing,
    numericType=WR,
    analysis=ModiaMath.AllAnalysis) where {Size} =
        RealSVector{Size}(name, within, info, start, unit, fixed, min, max, nominal, flow, causality, variability, numericType, integral, analysis)


text:
export iauBp06
"""
Frame bias and precession, IAU 2006.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   date1,date2  double         TT as a 2-part Julian Date (Note 1)

Returned:
   rb           double[3][3]   frame bias matrix (Note 2)
   rp           double[3][3]   precession matrix (Note 3)
   rbp          double[3][3]   bias-precession matrix (Note 4)

Notes:

   1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

            date1         date2

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

   2. The matrix rb transforms vectors from GCRS to mean J2000.0 by
   applying frame bias.

   3. The matrix rp transforms vectors from mean J2000.0 to mean of
   date by applying precession.

   4. The matrix rbp transforms vectors from GCRS to mean of date by
   applying frame bias then precession.  It is the product rp x rb.

   5. It is permissible to re-use the same array in the returned
   arguments.  The arrays are filled in the order given.

Called:
   iauPfw06     bias-precession F-W angles, IAU 2006
   iauFw2m      F-W angles to r-matrix
   iauPmat06    PB matrix, IAU 2006
   iauTr        transpose r-matrix
   iauRxr       product of two r-matrices
   iauCr        copy r-matrix

References:

   Capitaine, N. & Wallace, P.T., 2006, Astron.Astrophys. 450, 855

   Wallace, P.T. & Capitaine, N., 2006, Astron.Astrophys. 459, 981

This revision:  2013 August 21

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauBp06(double date1, double date2,
#              double rb[3][3], double rp[3][3], double rbp[3][3])

function iauBp06(date1::Real, date2::Real)
   # Allocate return value
   rb  = zeros(Float64, 3, 3)
   rp  = zeros(Float64, 3, 3)
   rbp = zeros(Float64, 3, 3)

   ccall((:iauBp06, libsofa_c), Cvoid,
       (Cdouble, Cdouble,
       Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
       convert(Float64, date1), convert(Float64, date2),
       rb, rp, rbp)

   return SMatrix{3,3}(rb'), SMatrix{3,3}(rp'), SMatrix{3,3}(rbp')
end


text:
@testcase "record creation" begin
    @lift struct Record_1
        a
        b
    end

    @lift @as_record Record_1

    @lift struct Record_2{A}
        x::A
        y::Int
    end

    @lift @as_record Record_2

    @lift @as_record struct AsRecordInDeclaration
        a
        b
    end

    @testset "matching" begin
        @test @match Record_1(1, 2) begin
            Record_1(1, 2) => true
            _ => false
        end
    end

    @testset "field extraction" begin

        @test @match Record_1(1, 2) begin
            Record_1(a=1) => true
            _ => false
        end

        @test @match Record_1(1, 2) begin
            Record_1(b=2) => true
            _ => false
        end
    end

    @testset "field punning" begin
        @test @match Record_1(1, 2) begin
            Record_1(;a=1) => true
            _ => false
        end

        @test @match Record_1(1, 2) begin
            Record_1(;b=2) => true
            _ => false
        end

        @test @match Record_1(1, 2) begin
            Record_1(;b) => b == 2
            _ => false
        end
    end

    @testset "parametric record" begin
        @test @match Record_2(1, 2) begin
            Record_2{A}(_) where A => A == typeof(1)
            _ => false
        end
    end

    @testset "macro in declaration" begin
        @test @match AsRecordInDeclaration(1, 2) begin
            AsRecordInDeclaration(1, 2) => true
            _ => false
        end
    end
end

@testcase "expression problem" begin
    # test whether addition of a new record in an old ADT works as expected
    # see https://en.wikipedia.org/wiki/Expression_problem for details

    # data structure compiled in some module M
    @lift @data A{B} begin
        C1(x::Int, y::Int, z::B)
        C2(x::Real, y::Real, z::B)
        C3(x::Complex, y::Complex, z::B)
    end

    f(x::A) = @match x begin
        C1(;z) => z
        C2(;z) => z
        C3(;z) => z
        _ => false
    end

    c₁ = C1(0, 0, 1)
    c₂ = C2(0.0, 0.0, 1.0)
    c₃ = C3(0.0 + im, 0.0 + im, 1.0 + im)

    # programmer tests his code, compiles and releases software:

    # commented for performance
    # @testcase "basic match" begin
    #     @test f(c₁) == 1
    #     @test f(c₂) == 1.0
    #     @test f(c₃) == 1.0 + im
    #     @test f(0) == false
    # end

    # programmer goes drink beer, later finds out he forgot to add UInts!

    @lift struct C0{B} <: A{B}
        x::UInt
        y::UInt
        z::B
    end

    # if programmer forgets this, pattern matching will not work
    @lift @as_record C0

    c₀ = C0(UInt(0x000), UInt(0xFFF), UInt(0x111))

    @testset "new records match" begin
        @test f(c₀) == false
        @test f(c₁) == 1
        @test f(c₂) == 1.0
        @test f(c₃) == 1.0 + im
    end

    # Unfortunately the function f is too generic to treat the C0 case,
    # but no problem:

    f(x::C0) = @match x begin
        C0(;z) => z
        _ => false
    end

    @testset "new implementation matches" begin
        @test f(c₀) == 0x111
        @test f(c₁) == 1
        @test f(c₂) == 1.0
        @test f(c₃) == 1.0 + im
    end

    # programmer is happy!

end


text:
using Documenter
using Codex

makedocs(
    sitename = "Codex.jl",
    pages = [
        "Codex" => "index.md",
        "Backup" => "backup.md",
        "GitLab" => "gitlab.md",
        "TransformExport" => "transformexport.md",
        "Questionnaires" => "questionnaires.md"
    ],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true")
)

repo = "github.com/rikhuijzer/Codex.jl.git"
devbranch = "main"
deploydocs(; repo, devbranch)


text:
"""
    MahalanobisKernel(P::AbstractMatrix)

Mahalanobis distance-based kernel given by
```math
    κ(x,y) =  exp(-r^2), r^2 = maha(x,P,y) = (x-y)'*inv(P)*(x-y)
```
where the matrix P is the metric.

"""
struct MahalanobisKernel{T<:Real, A<:AbstractMatrix{T}} <: SimpleKernel
    P::A
    function MahalanobisKernel(P::AbstractMatrix{T}) where {T<:Real}
        LinearAlgebra.checksquare(P)
        new{T,typeof(P)}(P)
    end
end

kappa(κ::MahalanobisKernel, d::T) where {T<:Real} = exp(-d)

metric(κ::MahalanobisKernel) = SqMahalanobis(κ.P)

Base.show(io::IO, κ::MahalanobisKernel) = print(io, "Mahalanobis Kernel (size(P) = ", size(κ.P), ")")


text:
"""
    acquisition_duration(x) -> F64Sec

Duration (in seconds) of volume acquisition. This field is REQUIRED for
sequences that are described with the volume_timingfield and that do not have the
slice_timing field set to allowed for accurate calculation of "acquisition time".
This field is mutually exclusive with repetition_time.

Corresponds to BIDS "AcquisitionDuration".
"""
acquisition_duration(x) = getter(x, :acquisition_duration, F64Sec, i -> 1.0u"s")
acquisition_duration!(x, val) = setter!(x, :acquisition_duration, F64Sec, val)

"""
    anatomical_landmark_coordinates(x) -> CoordinateList

Key:value pairs of any number of additional anatomical landmarks and their coordinates
in voxel units (where first voxel has index 0,0,0) relative to the associated anatomical
MRI, (e.g. Dict(:AC => (127.0,119.0,149.0), :PC=> (128.0,93.0,141.0),
:IH=> (131.0114.0,206.0)).

Corresponds to BIDS "AnatomicalLandmarkCoordinates".
"""
anatomical_landmark_coordinates(x) = getter(x, :anatomical_landmark_coordinates, CoordinateList, i -> CoordinateList())
anatomical_landmark_coordinates!(x, val) = setter!(x, :anatomical_landmark_coordinates, CoordinateList, val)

"""
    anatomical_landmark_coordinate_system(x)

Defines the coordinate system for the anatomical landmarks. See [`CoordinateSystem`](@ref).

Corresponds to BIDS "AnatomicalLandmarkCoordinateSystem".
"""
anatomical_landmark_coordinate_system(x) = coordinate_system(x)

"""
    anatomical_landmark_coordinate_units(x)

Units of the coordinates of anatomical_landmark_coordinate_system. MUST be m, cm, or mm.

Corresponds to BIDS "AnatomicalLandmarkCoordinateUnits"
"""
anatomical_landmark_coordinate_units(x) = spatial_units(x)

"""
anatomical_landmark_coordinate_systemDescription

Freeform text description or link to document describing the anatomical coordinate
system detail.

Corresponds to BIDS "AnatomicalLandmarkCoordinateSystemDescription"
"""
anatomical_landmark_coordinate_system_description(x) = getter(x, :anatomical_landmark_coordinate_system_description, String, i -> "")
anatomical_landmark_coordinate_system_description!(x, val) = setter!(x, :anatomical_landmark_coordinate_system_description, String, val)

"""
    coil_combination_method(x) -> String

Returns the coil combination method.

Almost all fMRI studies using phased-array coils use root-sum-of-squares (rSOS)
combination, but other methods exist. The image reconstruction is changed by
the coil combination method (as for the matrix coil mode above), so anything
non-standard should be reported.

Corresponds to BIDS "CoilCombinationMethod"
"""
coil_combination_method(x) =  getter(x, :coil_combination_method, String, i -> "rSOS")
coil_combination_method!(x, val) =  setter!(x, :coil_combination_method, String, val)

"""
    CogAtlasID(x) -> String

URL of the corresponding [Cognitive Atlas Task](https://www.cognitiveatlas.org/) term.
"""
CogAtlasID(x) =  getter(x, :CogAtlasID, String, i -> "")
CogAtlasID!(x, val) =  setter!(x, :CogAtlasID, String, val)

"""
    CogPOID(x) -> String

URL of the corresponding [CogPO](http://www.cogpo.org/) term.
"""
CogPOID(x) =  getter(x, :CogPOID, String, i -> "")
CogPOID!(x, val) =  setter!(x, :CogPOID, String, val)

"""
    contrast_bolus_ingredient(x) -> String

Return active ingredient of constrast agent. See [`ContrastIngrediant`](@ref) for
more details.

Corresponds to BIDS "ContrastBolusIngredient".
"""
contrast_bolus_ingredient(x) = getter(x, :contrast_bolus_ingredient, ContrastIngrediant, i -> UnkownContrast)
contrast_bolus_ingredient!(x, val) = setter!(x, :contrast_bolus_ingredient, ContrastIngrediant, val)

"""
    delay_time(x) -> F64Sec

Returns the user specified time (in seconds) to delay the acquisition of data for
the following volume. If the field is not present it is assumed to be set to zero.
Corresponds to Siemens CSA header field ldelay_timeInTR. This field is REQUIRED
for sparse sequences using the repetition_time field that do not have the
slice_timing field set to allowed for accurate calculation of "acquisition time".
This field is mutually exclusive with volume_timing.

Corresponds to BIDS "DelayTime".
"""
delay_time(x) = getter(x, :delay_time, F64Sec, i -> 1.0u"s")
delay_time!(x, val) = setter!(x, :delay_time, F64Sec, val)

"""
    delay_after_trigger(x) -> F64Sec

Returns duration (in seconds) from trigger delivery to scan onset. This delay is
commonly caused by adjustments and loading times. This specification is entirely
independent of number_of_volumes_discarded_by_scanner or number_of_volumes_discarded_by_user,
as the delay precedes the acquisition.

Corresponds to BIDS "DelayAfterTrigger".
"""
delay_after_trigger(x) = getter(x, :delay_after_trigger, F64Sec, i -> 1.0u"s")
delay_after_trigger!(x, val) = setter!(x, :delay_after_trigger, F64Sec, val)

"""
    device_serial_number(x) -> String

The serial number of the equipment that produced the composite instances.
A pseudonym can also be used to prevent the equipment from being identifiable,
so long as each pseudonym is unique within the dataset.

Corresponds to BIDS "DeviceSerialNumber".
"""
device_serial_number(x) = getter(x, :device_serial_number, String, i -> "00000")
device_serial_number!(x, val) = setter!(x, :device_serial_number, String, val)

"""
    dewar_position(x)

Position of the dewar during the MEG scan: upright, supine or degrees of angle
from vertical: for example on CTF systems, upright=15°, supine = 90°.

Corresponds to BIDS "DewarPosition".
"""
dewar_position(x) = getter(x, :dewar_position, IntDeg, i -> OneIntDeg)
dewar_position!(x, val) = setter!(x, :dewar_position, IntDeg, val)

"""
    dwell_time(x) -> F64Sec

Actual dwell time (in seconds) of the receiver per point in the readout
direction, including any oversampling. For Siemens, this corresponds to DICOM
field (0019,1018) (in ns). This value is necessary for the optional readout
distortion correction of anatomicals in the HCP Pipelines. It also usefully
provides a handle on the readout bandwidth, which isn’t captured in the other
metadata tags. Not to be confused with `effective_echo_spacing`, and the frequent
mislabeling of echo spacing (which is spacing in the phase encoding direction)
as "dwell time" (which is spacing in the readout direction).

Corresponds to BIDS "DwellTime".
"""
dwell_time(x) = getter(x, :dwell_time, F64Sec, i -> 1.0u"s")
dwell_time!(x, val) = setter!(x, :dwell_time, F64Sec, val)

"""
    echo_time(x) -> F64Sec

The echo time (TE) for the acquisition. This parameter is REQUIRED if corresponding
fieldmap data is present or the data comes from a multi echo sequence.

Corresponds to BIDS "EchoTime".
"""
echo_time(x) = getter(x, :echo_time, F64Sec, i -> 1.0u"s")
echo_time!(x, val) = setter!(x, :echo_time, F64Sec, val)

"""
    effective_echo_spacing(x) -> F64Sec

Returns the effective echo spacing.

The "effective" sampling interval, specified in seconds, between lines in the
phase-encoding direction, defined based on the size of the reconstructed image
in the phase direction. It is frequently, but incorrectly, referred to as
"dwell time" (see `dwell_time` parameter below for actual dwell time). It is
required for unwarping distortions using field maps. Note that beyond just
in-plane acceleration, a variety of other manipulations to the phase encoding
need to be accounted for properly, including partial fourier, phase
oversampling, phase resolution, phase field-of-view and interpolation.2
This parameter is REQUIRED if corresponding fieldmap data is present.

2Conveniently, for Siemens’ data, this value is easily obtained as
1/[`BWPPPE` * `ReconMatrixPE`], where BWPPPE is the
"`BandwidthPerPixelPhaseEncode` in DICOM tag (0019,1028) and ReconMatrixPE is
the size of the actual reconstructed data in the phase direction (which is NOT
reflected in a single DICOM tag for all possible aforementioned scan
manipulations). See [here](https://lcni.uoregon.edu/kb-articles/kb-0003) and
[here](https://github.com/neurolabusc/dcm_qa/tree/master/In/total_readout_time)

Corresponds to BIDS "EffectiveEchoSpacing".
"""
effective_echo_spacing(x) = getter(x, :effective_echo_spacing, F64Sec, i -> 1.0u"s")
effective_echo_spacing!(x, val) = setter!(x, :effective_echo_spacing, F64Sec, val)

"""
    eeg_electrode_groups(x) -> String

Field to describe the way electrodes are grouped into strips, grids or depth probes
e.g., `Dict(:grid1 => "10x8 grid on left temporal pole", :strip2 => "1x8 electrode strip on xxx")`.

Corresponds to BIDS "EEGElectrodeGroups".
"""
eeg_electrode_groups(x) = getter(x, :eeg_electrode_groups, Dict{Symbol,String}, i -> Dict{Symbol,String}())
eeg_electrode_groups!(x, val) = setter!(x, :eeg_electrode_groups, Dict{Symbol,String}, val)

"""
    eeg_ground(x) -> String

Description of the location of the ground electrode ("placed on right mastoid (M2)").

Corresponds to BIDS "EEGGround".
"""
eeg_ground(x) = getter(x, :eeg_ground, String, i -> "")
eeg_ground!(x, val) = setter!(x, :eeg_ground, String, val)

"""
    eeg_placement_scheme(x)

Freeform description of the placement of the iEEG electrodes.
Left/right/bilateral/depth/surface (e.g., "left frontal grid and bilateral
hippocampal depth" or "surface strip and STN depth" or "clinical indication
bitemporal, bilateral temporal strips and left grid").

Corresponds to BIDS "EEGPlacementScheme".
"""
eeg_placement_scheme(x) = getter(x, :eeg_placement_scheme, String, i -> "")
eeg_placement_scheme!(x, val) = setter!(x, :eeg_placement_scheme, String, val)

"""
    electrical_stimulation(x) -> Bool

Boolean field to specify if electrical stimulation was done during the recording
(options are `true` or `false`). Parameters for event-like stimulation should be
specified in the _events.tsv file.

Corresponds to BIDS "ElectricalStimulation".
"""
electrical_stimulation(x) = getter(x, :electrical_stimulation, Bool, i -> false)
electrical_stimulation!(x, val) = setter!(x, :electrical_stimulation, Bool, val)

"""
    electrical_stimulation_parameters(x) -> String

Free form description of stimulation parameters, such as frequency, shape etc.
Specific onsets can be specified in the _events.tsv file. Specific shapes can be
described here in freeform text.

Corresponds to BIDS "ElectricalStimulationParameters".
"""
electrical_stimulation_parameters(x) = getter(x, :electrical_stimulation_parameters, String, i -> "")
electrical_stimulation_parameters!(x, val) = setter!(x, :electrical_stimulation_parameters, String, val)

"""
    epoch_length(x) -> F64Sec

Duration of individual epochs in seconds (e.g., 1) in case of epoched data. If
recording was continuous or discontinuous, leave out the field.

Corresponds to BIDS "EpochLength".
"""
epoch_length(x) = getter(x, :epoch_length, F64Sec, OneF64Sec)
epoch_length!(x, val) = setter!(x, :epoch_length, F64Sec, val)

"""
    fiducial_description(x) -> String

A freeform text field documenting the anatomical landmarks that were used and how
the head localization coils were placed relative to these. This field can describe,
for instance, whether the true anatomical locations of the left and right pre-auricular
points were used and digitized, or rather whether they were defined as the intersection
between the tragus and the helix (the entry of the ear canal), or any other anatomical
description of selected points in the vicinity of the ears.

Corresponds to BIDS "FiducialDescription".
"""
fiducial_description(x) = getter(x, :fiducial_description, String, i -> "")
fiducial_description!(x, val) = setter!(x, :fiducial_description, String, val)

"""
    flip_angle(x) -> IntDeg

Returns the flip angle for the acquisition in degrees.

Corresponds to BIDS "FlipAngle".
"""
flip_angle(x) = getter(x, :flip_angle, IntDeg, i -> OneIntDeg)
flip_angle!(x, val) = setter!(x, :flip_angle, IntDeg, val)

"""
    gradient_set_type(x) -> String

Returns the gradient set type. It should be possible to infer the gradient coil
from the scanner model. If not, e.g. because of a custom upgrade or use of a
gradient insert set, then the specifications of the actual gradient coil should
be reported independently.

Corresponds to BIDS "GradientSetType".
"""
gradient_set_type(x) = getter(x, :gradient_set_type, String, i -> "")
gradient_set_type!(x, val) = setter!(x, :gradient_set_type, String, val)

"""
    head_coil_coordinates(x) -> CoordinateList

Key:value pairs describing head localization coil labels and their coordinates,
interpreted following the head_coil_coordinate_system, e.g., {NAS: [12.7,21.3,13.9],
LPA: [5.2,11.3,9.6], RPA: [20.2,11.3,9.1]}. Note that coils are not always placed
at locations that have a known anatomical name (e.g. for Elekta, Yokogawa systems);
in that case generic labels can be used (e.g. {coil1: [12.2,21.3,12.3],
coil2: [6.7,12.3,8.6], coil3: [21.9,11.0,8.1]} ).

Corresponds to BIDS "HeadCoilCoordinates".
"""
head_coil_coordinates(x) = getter(x, :head_coil_coordinates, CoordinateList, i -> CoordinateList())
head_coil_coordinates!(x, val) = setter!(x, :head_coil_coordinates, CoordinateList, val)

"""
    head_coil_coordinate_system(x) -> CoordinateSystem

Defines the coordinate system for the coils.

Corresponds to BIDS "HeadCoilCoordinateSystem".
"""
head_coil_coordinate_system(x) = getter(x, :head_coil_coordinate_system, CoordinateSystem, i -> UnknownSpace)

"""
    head_coil_coordinate_system_description(x) -> String

Freeform text description or link to document describing the Head Coil coordinate
system system in detail.

Corresponds to BIDS "HeadCoilCoordinateSystemDescription".
"""
head_coil_coordinate_system_description(x) = getter(x, :head_coil_coordinate_system_description, String, i -> "")
head_coil_coordinate_system_description!(x, val) = setter!(x, :head_coil_coordinate_system_description, String, val)

"""
    head_coil_coordinate_units(x) -> m, cm, or mm

Units of the coordinates of head_coil_coordinate_system.

Corresponds to BIDS "HeadCoilCoordinateUnits".
"""
head_coil_coordinate_units(x) = spatial_units(x)

"""
    institution_address(x) -> String

Return the address of the institution in charge of the equipment that produced
the composite instances.

Corresponds to BIDS "InstitutionAddress".
"""
institution_address(x) = getter(x, :institution_address, String, i -> "")
institution_address!(x, val) = setter!(x, :institution_address, String, val)

"""
    institutional_department_name(x) -> String

Return the department in the institution in charge of the equipment that produced
the composite instances.

Corresponds to BIDS "InstitutionalDepartmentName".
"""
institutional_department_name(x) = getter(x, :institutional_department_name, String, i -> "")
institutional_department_name!(x, val) = setter!(x, :institutional_department_name, String, val)

"""
    institution_name(x) -> String

Returns the name of the institution in charge of the equipment that produced the
composite instances.

Corresponds to BIDS "InstitutionName".
"""
institution_name(x) = getter(x, :institution_name, String, i -> "")
institution_name!(x, val) = setter!(x, :institution_name, String, val)

"""
    instructions(x) -> String

Text of the instructions given to participants before the scan. This is especially
important in context of resting state fMRI and distinguishing between eyes open and
eyes closed paradigms.

Corresponds to BIDS "Instructions".
"""
instructions(x) = getter(x, :instructions, String, i -> "")
instructions!(x, val) = setter!(x, :instructions, String, val)

"""
    intended_for(x) -> Vector{String}

Path or list of path relative to the subject subfolder pointing to the structural
MRI, possibly of different types if a list is specified, to be used with the MEG
recording. The path(s) need(s) to use forward slashes instead of backward slashes
(e.g. ses-/anat/sub-01_T1w.nii.gz).

Corresponds to BIDS "IntendedFor".
"""
intended_for(x) = getter(x, :intended_for, Vector{String}, i -> String[])
intended_for!(x, val) = setter!(x, :intended_for, Vector{String}, val)

"""
    inversion_time(x) -> F64Sec

Returns the inversion time (TI) for the acquisition, specified in seconds.
Inversion time is the time after the middle of inverting RF pulse to middle of
excitation pulse to detect the amount of longitudinal magnetization.

Corresponds to BIDS "InversionTime".
"""
inversion_time(x) = getter(x, :inversion_time, F64Sec, i -> 1.0u"s")
inversion_time!(x, val) = setter!(x, :inversion_time, F64Sec, val)

"""
    magnetic_field_strength(x) -> F64Tesla

The nominal field strength of MR magnet in Tesla.

Corresponds to BIDS "MagneticFieldStrength".
"""
magnetic_field_strength(x) = getter(x, :magnetic_field_strength, F64Tesla, i -> 3.0u"T")
magnetic_field_strength!(x, val) = setter!(x, :magnetic_field_strength, F64Tesla, val)

"""
    manufacturer(x) -> String

The manufacturer of the equipment that produced the composite instances.

Corresponds to BIDS "Manufacturer".
"""
manufacturer(x) = getter(x, :manufacturer, String, i -> "")
manufacturer!(x, val) = setter!(x, :manufacturer, String, val)

"""
    manufacturer_model_name(x) -> String

The manufacturer's model name of the equipment that produced the composite instances.

Corresponds to BIDS "ManufacturerModelName".
"""
manufacturer_model_name(x) = getter(x, :manufacturer_model_name, String, i -> "")
manufacturer_model_name!(x, val) = setter!(x, :manufacturer_model_name, String, val)

"""
    matrix_coil_mode(x) -> String

Returns the matrix coil mode. A method for reducing the number of independent
channels by combining in analog the signals from multiple coil elements. There
are typically different default modes when using un-accelerated or accelerated
(e.g. GRAPPA, SENSE) imaging.

Corresponds to BIDS "MatrixCoilMode".
"""
matrix_coil_mode(x) = getter(x, :matrix_coil_mode, String, i -> "")
matrix_coil_mode!(x, val) = setter!(x, :matrix_coil_mode, String, val)

"""
    mr_transmit_coil_sequence(x) -> String

A sequence that provides information about the transmit coil used. This is a
relevant field if a non-standard transmit coil is used.

Corresponds to BIDS "MRTransmitCoilSequence".
"""
mr_transmit_coil_sequence(x) = getter(x, :mr_transmit_coil_sequence, String, i -> "")
mr_transmit_coil_sequence!(x, val) = setter!(x, :mr_transmit_coil_sequence, String, val)

"""
    multiband_acceleration_factor(x) -> String

Returns the multiband factor, for multiband acquisitions.

Corresponds to BIDS "MultibandAccelerationFactor".
"""
multiband_acceleration_factor(x) = getter(x, :multiband_acceleration_factor, String, i -> "")
multiband_acceleration_factor!(x, val) = setter!(x, :multiband_acceleration_factor, String, val)

"""
    negative_contrast(x) -> Bool

Specifies whether increasing voxel intensity (within sample voxels)
denotes a decreased value with respect to the contrast suffix. This is commonly
the case when Cerebral Blood Volume is estimated via usage of a contrast agent
in conjunction with a T2* weighted acquisition protocol.

Corresponds to BIDS "NegativeContrast".
"""
negative_contrast(x) = getter(x, :negative_contrast, Bool, i -> false)
negative_contrast!(x, val) = setter!(x, :negative_contrast, Bool, val)

"""
    nonlinear_gradient_correction(x) -> Bool

Returns `Bool` stating if the image saved has been corrected for gradient
nonlinearities by the scanner sequence. Default is `false`.

Corresponds to BIDS "".
"""
nonlinear_gradient_correction(x) = getter(x, :nonlinear_gradient_correction, Bool, i -> false)
nonlinear_gradient_correction!(x, val) = setter!(x, :nonlinear_gradient_correction, Bool, val)

"""
    number_of_volumes_discarded_by_scanner(x) -> Int

Returns the number of volumes ("dummy scans") discarded by the scanner (as opposed
to those discarded by the user post hoc) before saving the imaging file. For example,
a sequence that automatically discards the first 4 volumes before saving would
have this field as 4. A sequence that doesn't discard dummy scans would have
this set to 0. Please note that the onsets recorded in the _event.tsv file
should always refer to the beginning of the acquisition of the first volume in
the corresponding imaging file - independent of the value of
number_of_volumes_discarded_by_scanner field.

Corresponds to BIDS "NumberOfVolumesDiscardedByScanner".
"""
number_of_volumes_discarded_by_scanner(x) = getter(x, :number_of_volumes_discarded_by_scanner, Int, i -> 0)
number_of_volumes_discarded_by_scanner!(x, val) = setter!(x, :number_of_volumes_discarded_by_scanner, Int, val)

"""
    number_of_volumes_discarded_by_user(x) -> Int

Returns the  number of volumes ("dummy scans") discarded by the user before
including the file in the dataset. If possible, including all of the volumes is
strongly recommended. Please note that the onsets recorded in the _event.tsv file
should always refer to the beginning of the acquisition of the first volume in
the corresponding imaging file - independent of the value of
number_of_volumes_discarded_by_user field.

Corresponds to BIDS "NumberOfVolumesDiscardedByUser".
"""
number_of_volumes_discarded_by_user(x) = getter(x, :number_of_volumes_discarded_by_user, Int, i -> 0)
number_of_volumes_discarded_by_user!(x, val) = setter!(x, :number_of_volumes_discarded_by_user, Int, val)

"""
    number_shots(x) -> Int

Returns the number of RF excitations needed to reconstruct a slice or volume.
Please mind that this is not the same as Echo Train Length which denotes the
number of lines of k-space collected after an excitation.

Corresponds to BIDS "NumberShots".
"""
number_shots(x) = getter(x, :number_shots, Int, i -> 0)
number_shots!(x, val) = setter!(x,  :number_shots, Int, val)

"""
    parallel_acquisition_technique(x) -> String

Returns the type of parallel imaging used (e.g. GRAPPA, SENSE).

Corresponds to BIDS "ParallelAcquisitionTechnique".
"""
parallel_acquisition_technique(x) = getter(x, :parallel_acquisition_technique, String, i -> "")
parallel_acquisition_technique!(x, val) = setter!(x, :parallel_acquisition_technique, String, val)

"""
    parallel_reduction_factor(x) -> Int

The parallel imaging (e.g, GRAPPA) factor. Use the denominator of the fraction
of k-space encoded for each slice. For example, 2 means half of k-space is
encoded.                                                                                |

Corresponds to BIDS "ParallelReductionFactor".
"""
parallel_reduction_factor(x) = getter(x, :parallel_reduction_factor, Int, i -> 0)
parallel_reduction_factor!(x, val) = setter!(x, :parallel_reduction_factor, Int, val)

"""
    partial_fourier(x) -> Float64

Returns the fraction of partial Fourier information collected.

Corresponds to BIDS "PartialFourier".
"""
partial_fourier(x) = getter(x, :partial_fourier, Float64, i -> 1.0)
partial_fourier!(x, val) = setter!(x, :partial_fourier, Float64, val)

"""
    partial_fourier_direction(x) -> String

Returns the direction where only partial Fourier information was collected.

Corresponds to BIDS "PartialFourierDirection".
"""
partial_fourier_direction(x) = getter(x, :partial_fourier_direction, String, i -> "")
partial_fourier_direction!(x, val) = setter!(x, :partial_fourier_direction, String, val)

"""
    phase_encoding_direction(x) -> EncodingDirection

Returns the phase encoding direction.

`phase_encoding_direction` is defined as the direction along which phase is was
modulated which may result in visible distortions. Note that this is not the
same as the DICOM term `InPlanephase_encoding_directiong` which can have `ROW` or
`COL` values. This parameter is REQUIRED if corresponding fieldmap data is present
or when using multiple runs with different phase encoding directions (which can
be later used for field inhomogeneity correction).

Corresponds to BIDS "PhaseEncodingDirection".
"""
phase_encoding_direction(x) = EncodingDirection(phasedim(x))
phase_encoding_direction!(x, val) = phasedim!(x, val)

"""
    power_line_frequency(x) -> F64Hz

Frequency (in Hz) of the power grid at the geographical location of the EEG
instrument (i.e., 50 or 60).

Corresponds to BIDS "PowerLineFrequency".
"""
power_line_frequency(x) = getter(x, :power_line_frequency, F64Hz, i -> 1.0u"Hz")
power_line_frequency!(x, val) = setter!(x, :power_line_frequency, F64Hz, val)

"""
    pulse_sequence(x) -> String

General description of the pulse sequence used for the scan (i.e. MPRAGE,
Gradient Echo EPI, Spin Echo EPI, Multiband gradient echo EPI).

Corresponds to BIDS "PulseSequence".
"""
pulse_sequence(x) = getter(x, :pulse_sequence, String, i -> "")
pulse_sequence!(x, val) = setter!(x, :pulse_sequence, String, val)

"""
    pulse_sequence_details(x) -> String

Information beyond pulse sequence type that identifies the specific pulse
sequence used (i.e. "Standard Siemens Sequence distributed with the VB17
software," "Siemens WIP ### version #.##," or "Sequence written by X using a
version compiled on MM/DD/YYYY").

Corresponds to BIDS "PulseSequenceDetails".
"""
pulse_sequence_details(x) = getter(x, :pulse_sequence_details, String, i -> "")
pulse_sequence_details!(x, val) = setter!(x, :pulse_sequence_details, String, val)

"""
    pulse_sequence_type(x) -> String

A general description of the pulse sequence used for the scan (i.e. MPRAGE,
Gradient Echo EPI, Spin Echo EPI, Multiband gradient echo EPI).

Corresponds to BIDS "PulseSequenceType".
"""
pulse_sequence_type(x) = getter(x, :pulse_sequence_type, String, i -> "")
pulse_sequence_type!(x, val) = setter!(x, :pulse_sequence_type, String, val)

"""
    receive_coil_name(x) -> String

The information describing the receiver coil. Corresponds to DICOM Tag
Receive Coil Name, although not all vendors populate that DICOM Tag, in which
case this field can be derived from an appropriate private DICOM field.

Corresponds to BIDS "ReceiveCoilName".
"""
receive_coil_name(x) = getter(x, :receive_coil_name, String, i -> "")
receive_coil_name!(x, val) = setter!(x, :receive_coil_name, String, val)

"""
    receive_coil_active_elements(x) -> String

Information describing the active/selected elements of the receiver coil. This
doesn’t correspond to a tag in the DICOM ontology. The vendor-defined
terminology for active coil elements can go in this field. As an example, for
Siemens, coil channels are typically not activated/selected individually, but
rather in pre-defined selectable "groups" of individual channels, and the list
of the "groups" of elements that are active/selected in any given scan
populates the Coil String entry in Siemens’ private DICOM fields (e.g., HEA;HEP
for the Siemens standard 32 ch coil when both the anterior and posterior groups
are activated). This is a flexible field that can be used as most appropriate
for a given vendor and coil to define the "active" coil elements. Since
individual scans can sometimes not have the intended coil elements selected, it
is preferable for this field to be populated directly from the DICOM for each
individual scan, so that it can be used as a mechanism for checking that a
given scan was collected with the intended coil elements selected.

Corresponds to BIDS "ReceiveCoilActiveElements".
"""
receive_coil_active_elements(x) = getter(x, :receive_coil_active_elements, String, i -> "")
receive_coil_active_elements!(x, val) = setter!(x, :receive_coil_active_elements, String, val)

"""
    recording_duration(x) -> F64Sec

Length of the recording in seconds (e.g., 3600)

Corresponds to BIDS "RecordingDuration".
"""
recording_duration(x) = getter(x, :recording_duration, F64Sec, OneF64Sec)
recording_duration!(x, val) = setter!(x, :recording_duration, F64Sec, val)

# TODO: should this return an enumerable
"""
    recording_type(x) -> String

Defines whether the recording is "continuous", "discontinuous" or "epoched";
this latter limited to time windows about events of interest (e.g., stimulus
presentations, subject responses etc.)

Corresponds to BIDS "RecordingType".
"""
recording_type(x) = getter(x, :recording_type, String, i -> "")
recording_type!(x, val) = setter!(x, :recording_type, String, val)

"""
    repetition_time(x) -> F64Sec

Returns the time in seconds between the beginning of an acquisition of one volume
and the beginning of acquisition of the volume following it (TR). Please note that
this definition includes time between scans (when no data has been acquired) in
case of sparse acquisition schemes. This value needs to be consistent with the
pixdim[4] field (after accounting for units stored in xyzt_units field) in the
NIfTI header. This field is mutually exclusive with volume_timing and is derived
from DICOM Tag 0018, 0080 and converted to seconds.

Corresponds to BIDS "RepetitionTime".
"""
repetition_time(x) = getter(x, :repetition_time, F64Sec, i -> 1.0u"s")
repetition_time!(x, val) = setter!(x, :repetition_time, F64Sec, val)

"""
    sampling_frequency -> F64Hz

Sampling frequency (in Hz) of all the data in the recording, regardless of their
type (e.g., 2400).

Corresponds to BIDS "SamplingFrequency".
"""
sampling_frequency(x) = getter(x, :sampling_frequency, F64Hz, i -> 1.0u"Hz")
sampling_frequency!(x, val) = setter!(x, :sampling_frequency, F64Hz, val)

"""
    scan_options(x)

Parameters of scanning_sequence. Corresponds to DICOM  `Scan Options`.

Corresponds to BIDS "ScanOptions".
"""
scan_options(x) = getter(x, :scan_options, Dict{Symbol,Any}, i -> Dict{Symbol,Any}())
scan_options!(x, val) = setter!(x, :scan_options, Dict{Symbol,Any}, val)

# TODO enumerable for bolus ingredient
"""
    scanning_sequence(x) -> String

Returns the description of the type of sequence data acquired.

Corresponds to BIDS "ScanningSequence".
"""
scanning_sequence(x) = getter(x, :scanning_sequence, String, i -> "")
scanning_sequence!(x, val) = setter!(x, :scanning_sequence, String, val)

"""
    sequence_name(x) -> String

Returns the manufacturer’s designation of the sequence name.

Corresponds to BIDS "SequenceName".
"""
sequence_name(x) = getter(x, :sequence_name, String, i -> "")
sequence_name!(x, val) = setter!(x, :sequence_name, String, val)

"""
    sequence_varient(x) -> String

Returns the variant of the `scanning_sequence` property.

Corresponds to BIDS "SequenceVarient".
"""
sequence_varient(x) = getter(x, :sequence_varient, String, i -> "")
sequence_varient!(x, val) = setter!(x, :sequence_varient, String, val)

"""
    slice_encoding_direction(x) -> EncodingDirection

Possible values: `i`, `j`, `k`, `ineg, `jneg`, `kneg` (the axis of the NIfTI data along which
slices were acquired, and the direction in which slice_timing is defined with
respect to). `i`, `j`, `k` identifiers correspond to the first, second and third axis
of the data in the NIfTI file. `*neg` indicates that the contents of
slice_timing are defined in reverse order - that is, the first entry corresponds
to the slice with the largest index, and the final entry corresponds to slice
index zero. When present, the axis defined by slice_encoding_direction needs to be
consistent with the ‘slicedim’ field in the NIfTI header. When absent, the
entries in slice_timing must be in the order of increasing slice index as defined
by the NIfTI header.

Corresponds to BIDS "SliceEncodingDirection".
"""
slice_encoding_direction(x) = EncodingDirection(slicedim(x))
slice_encoding_direction!(x, val) = slicedim!(x, val)

"""
    slice_timing(x) -> Vector{F64Sec}

The time at which each slice was acquired within each volume (frame) of the
acquisition. Slice timing is not slice order -- rather, it is a list of times
(in JSON format) containing the time (in seconds) of each slice acquisition in
relation to the beginning of volume acquisition. The list goes through the
slices along the slice axis in the slice encoding dimension (see below). Note
that to ensure the proper interpretation of the `slice_timing` field, it is
important to check if the OPTIONAL `slice_encoding_direction` exists. In
particular, if `slice_encoding_direction` is negative, the entries in
`slice_timing` are defined in reverse order with respect to the slice axis
(i.e., the final entry in the `slice_timing` list is the time of acquisition of
slice 0). This parameter is REQUIRED for sparse sequences that do not have the
`delay_time` field set. In addition without this parameter slice time correction
will not be possible.

Corresponds to BIDS "SliceTiming".
"""
slice_timing(x) = getter(x, :slice_timing, Vector{F64Sec}, i -> F64Sec[])
slice_timing!(x, val) = setter!(x, :slice_timing, Vector{F64Sec}, val)

"""
    software_versions(x) -> String

The manufacturer’s designation of software version of the equipment
that produced the composite instances.

Corresponds to BIDS "SoftwareVersions".
"""
software_versions(x) = getter(x, :software_versions, String, i -> "")
software_versions!(x, val) = setter!(x, :software_versions, String, val)

"""
    start_time(x) -> F64Sec

Start time in seconds in relation to the start of acquisition of the first data
sample in the corresponding neural dataset (negative values are allowed).

Corresponds to BIDS "StartTime".
"""
start_time(x) = first(timaxis(x))

"""
    station_name(x) -> String

Institution defined name of the machine that produced the composite instances.

Corresponds to BIDS "StationName".
"""
station_name(x) = getter(x, :station_name, String, i -> "")
station_name!(x, val) = setter!(x, :station_name, String, val)

"""
    subject_artefact_description(x) -> String

Freeform description of the observed subject artefact and its possible cause (e.g.,
"door open", "nurse walked into room at 2 min", "seizure at 10 min"). If this
field is left empty, it will be interpreted as absence of artifacts.

Corresponds to BIDS "SubjectArtefactDescription".
"""
subject_artefact_description(x) = getter(x, :subject_artefact_description, String, i -> "")
subject_artefact_description!(x, val) = setter!(x, :subject_artefact_description, String, val)

"""
    task_description(x) -> String

Longer description of the task.

Corresponds to BIDS "TaskDescription".
"""
task_description(x) = getter(x, :task_description, String, i -> "")
task_description!(x, val) = setter!(x, :task_description, String, val)

"""
    task_name(x) -> String

Name of the task. No two tasks should have the same name. Task label (task-) included
in the file name is derived from this field by removing all non alphanumeric
([a-zA-Z0-9]) characters. For example task name faces n-back will corresponds to
task label facesnback. A RECOMMENDED convention is to name resting state task using
labels beginning with rest.

Corresponds to BIDS "TaskName".
"""
task_name(x) = getter(x, :task_name, String, i -> "")
task_name!(x, val) = setter!(x, :task_name, String, val)

"""
    total_readout_time(x) -> F64Sec

The total readout time. This is actually the "effective" total readout time ,
defined as the readout duration, specified in seconds, that would have generated
data with the given level of distortion. It is NOT the actual, physical duration
of the readout train. If `effective_echo_spacing` has been properly computed, it
is just `effective_echo_spacing * (ReconMatrixPE - 1)`.3 .

* This parameter is
REQUIRED if corresponding "field/distortion" maps acquired with opposing phase
encoding directions are present.

3We use the "FSL definition", i.e, the time between the center of the
first "effective" echo and the center of the last "effective" echo.

Corresponds to BIDS "TotalReadoutTime".
"""
total_readout_time(x) = getter(x, :total_readout_time, F64Sec, i -> 1.0u"s")
total_readout_time!(x, val) = setter!(x, :total_readout_time, F64Sec, val)

"""
    volume_timing(x) -> Vector{F64Sec}

Returns the time at which each volume was acquired during the acquisition. It is
described using a list of times (in JSON format) referring to the onset of each
volume in the BOLD series. The list must have the same length as the BOLD
series, and the values must be non-negative and monotonically increasing. This
field is mutually exclusive with repetition_time and delay_time. If defined, this
requires acquisition time (TA) be defined via either slice_timing or
acquisition_duration be defined.

Corresponds to BIDS "VolumeTiming".
"""
volume_timing(x) = getter(x, :volume_timing, Vector{F64Sec}, i -> F64Sec[])
volume_timing!(x, val) = setter!(x, :volume_timing, Vector{F64Sec}, val)


text:
##################################
# Cross Plaform tests for spawn. #
##################################

# Assumes the following are available in the
#  - GNU coreutils (or equivalent)
#  - perl

#TODO:
# - Windows:
#   - Add a test whether coreutils are available and skip tests if not
yes = `perl -le 'while (1) {print STDOUT "y"}'`

#### Examples used in the manual ####

@test readall(`echo hello | sort`) == "hello | sort\n"
@test readall(`echo hello` |> `sort`) == "hello\n"
@test length(spawn(`echo hello` |> `sort`).processes) == 2

out = readall(`echo hello` & `echo world`)
@test search(out,"world") != (0,0)
@test search(out,"hello") != (0,0)
@test readall((`echo hello` & `echo world`) |> `sort`)=="hello\nworld\n"

@test (run(`printf "       \033[34m[stdio passthrough ok]\033[0m\n"`); true)

# Test for SIGPIPE being treated as normal termination (throws an error if broken)
@unix_only @test (run(yes|>`head`|>DevNull); true)

begin
    a = Base.Condition()
    @schedule begin
        p = spawn(yes|>DevNull)
        Base.notify(a,p)
        @test !success(p)
    end
    p = wait(a)
    kill(p)
end

@test_throws Base.UVError run(`foo_is_not_a_valid_command`)

if false
    prefixer(prefix, sleep) = `perl -nle '$|=1; print "'$prefix' ", $_; sleep '$sleep';'`
    @test success(`perl -le '$|=1; for(0..2){ print; sleep 1 }'` |>
                  prefixer("A",2) & prefixer("B",2))
    @test success(`perl -le '$|=1; for(0..2){ print; sleep 1 }'` |>
                  prefixer("X",3) & prefixer("Y",3) & prefixer("Z",3) |>
                  prefixer("A",2) & prefixer("B",2))
end

@test  success(`true`)
@test !success(`false`)
@test success(`true` |> `true`)
if false
    @test  success(ignorestatus(`false`))
    @test  success(ignorestatus(`false`) |> `true`)
    @test !success(ignorestatus(`false`) |> `false`)
    @test !success(ignorestatus(`false`) & `false`)
    @test  success(ignorestatus(`false` |> `false`))
    @test  success(ignorestatus(`false` & `false`))
end

# STDIN Redirection
file = tempname()
run(`echo hello world` |> file)
@test readall(file |> `cat`) == "hello world\n"
@test open(readall, file |> `cat`, "r") == "hello world\n"
rm(file)

# Stream Redirection
@unix_only begin
    @async begin
        server = listen(2326)
        client = accept(server)
        @test readall(client |> `cat`) == "hello world\n"
        close(server)
    end
    @async begin
        sock = connect(2326)
        run(`echo hello world` |> sock)
        close(sock)
    end
end

readall(setenv(`sh -c "echo \$TEST"`,["TEST=Hello World"])) == "Hello World\n"
readall(setenv(`sh -c "echo \$TEST"`,["TEST"=>"Hello World"])) == "Hello World\n"
readall(setenv(`sh -c "pwd"`;dir="/")) == readall(setenv(`sh -c "cd / && pwd"`))

# Here we test that if we close a stream with pending writes, we don't lose the writes.
str = ""
for i=1:1000
  str = "$str\n $(randstring(10))"
end
stdout, stdin, proc = readandwrite(`cat -`)
write(stdin, str)
close(stdin)
str2 = readall(stdout)
@test str2 == str

# This test hangs if the end of run walk across uv streams calls shutdown on a stream that is shutting down.
file = tempname()
open(`cat -` |> file, "w") do io
    write(io, str)
end
rm(file)

# issue #3373
# fixing up Conditions after interruptions
r = RemoteRef()
t = @async begin
    try
        wait(r)
    end
    p = spawn(`sleep 1`); wait(p)
    @test p.exitcode == 0
end
yield()
schedule(t, InterruptException(), error=true)
yield()
put!(r,11)
yield()

# Test marking of AsyncStream

@async begin
    server = listen(2327)
    client = accept(server)
    write(client, "Hello, world!\n")
    write(client, "Goodbye, world...\n")
    close(server)
end
sleep(0.1)
sock = connect(2327)
mark(sock)
@test ismarked(sock)
@test readline(sock) == "Hello, world!\n"
@test readline(sock) == "Goodbye, world...\n"
@test reset(sock) == 0
@test !ismarked(sock)
mark(sock)
@test ismarked(sock)
@test readline(sock) == "Hello, world!\n"
unmark(sock)
@test !ismarked(sock)
@test_throws ErrorException reset(sock)
@test !unmark(sock)
@test readline(sock) == "Goodbye, world...\n"
#@test eof(sock) ## doesn't work...
close(sock)

# issue #4535
exename=joinpath(JULIA_HOME,(ccall(:jl_is_debugbuild,Cint,())==0?"julia":"julia-debug"))
@test readall(`$exename -f -e 'println(STDERR,"Hello World")'` .> `cat`) == "Hello World\n"

# issue #6310
@test readall(`echo "2+2"` |> `$exename -f`) == "4\n"

# issue #5904
@test run(ignorestatus(`false`) |> `true`) === nothing

# issue #6010
# TODO: should create separate set of task tests
ducer = @async for i=1:100; produce(i); end
yield()
@test consume(ducer) == 1
@test consume(ducer) == 2

# redirect_*
OLD_STDOUT = STDOUT
fname = tempname()
f = open(fname,"w")
redirect_stdout(f)
println("Hello World")
redirect_stdout(OLD_STDOUT)
close(f)
@test "Hello World\n" == readall(fname)
@test is(OLD_STDOUT,STDOUT)

max_stars_repo_name: Aerlinger/julia

max_issues_repo_name: Aerlinger/julia

max_forks_repo_name: Aerlinger/julia


text:
function op_H( pw, Vpot, psi, ik )
  return op_K( pw, psi, ik ) + op_Vpot( pw, Vpot, psi, ik )
end


text:
using CUSPARSE
using CUDArt
using Base.Test

m = 25
n = 35
k = 10
blockdim = 5

@testset "sctr" begin
    @testset for elty in [Float32,Float64,Complex64,Complex128]
        x = sparsevec(rand(1:m,k), rand(elty,k), m)
        y = zeros(elty,m)
        @testset "sctr!" begin
            d_x = CudaSparseVector(x)
            d_y = CudaArray(y)
            d_y = CUSPARSE.sctr!(d_x,d_y,'O')
            h_y = to_host(d_y)
            y[x.nzind]  += x.nzval
            @test h_y ≈ y
        end
        y = zeros(elty,m)

        @testset "sctr" begin
            d_x = CudaSparseVector(x)
            d_y = CUSPARSE.sctr(d_x,'O')
            h_y = to_host(d_y)
            y = zeros(elty,m)
            y[x.nzind]  += x.nzval
            @test h_y ≈ y
        end
    end
end


text:
"""
Reserved Julia identifiers will be prepended with "_"
"""
const RESERVED_WORDS = ["begin", "while", "if", "for", "try", "return", "break", "continue",
                        "function", "macro", "quote", "let", "local", "global", "const", "do",
                        "struct", "module", "baremodule", "using", "import", "export", "end",
                        "else", "elseif", "catch", "finally", "true", "false"]

"""
    name_safe(name::AbstractString)
Return a valid Julia variable name, prefixed with "_" if the `name` is conflict with Julia's
reserved words.
"""
name_safe(name::AbstractString) = (name in RESERVED_WORDS) ? "_"*name : name

"""
    symbol_safe(name::AbstractString)
Same as [`name_safe`](@ref), but return a Symbol.
"""
symbol_safe(name::AbstractString) = Symbol(name_safe(name))

## dumpobj
printind(ind::Int, st...) = println(join([repeat(" ", 2*ind), st...]))

dumpobj(cursor::CLCursor) = dumpobj(0, cursor)
dumpobj(t::CLType) = join(typeof(t), " ", spelling(t))
dumpobj(t::CLInt) = t
dumpobj(t::CLPointer) = pointee_type(t)
dumpobj(ind::Int, t::CLType) = printind(ind, dumpobj(t))

function dumpobj(ind::Int, cursor::Union{CLFieldDecl, CLParmDecl})
    printind(ind+1, typeof(cursor), " ", dumpobj(type(cursor)), " ", name(cursor))
end

function dumpobj(ind::Int, node::Union{CLCursor,CLStructDecl,CLCompoundStmt,CLFunctionDecl,CLBinaryOperator})
    printind(ind, " ", typeof(node), " ", name(node))
    for c in children(node)
        dumpobj(ind + 1, c)
    end
end

"""
    copydeps(dst)
Copy dependencies to `dst`.
"""
function copydeps(dst)
    cp(joinpath(@__DIR__, "ctypes.jl"), joinpath(dst, "ctypes.jl"), force=true)
end

function print_template(path, libname="LibXXX")
    open(path, "w") do io
        print(io, """
        module $libname

        import Libdl

        # Load in `deps.jl`, complaining if it does not exist
        const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
        if !isfile(depsjl_path)
            error("$libname was not build properly. Please run Pkg.build(\\"$libname\\").")
        end
        include(depsjl_path)
        # Module initialization function
        function __init__()
            check_deps()
        end

        using CEnum

        include("ctypes.jl")
        export Ctm, Ctime_t, Cclock_t

        include(joinpath(@__DIR__, "..", "gen", "$(lowercasefirst(libname))_common.jl"))
        include(joinpath(@__DIR__, "..", "gen", "$(lowercasefirst(libname))_api.jl"))

        # export everything
        #foreach(names(@__MODULE__, all=true)) do s
        #    if startswith(string(s), "SOME_PREFIX")
        #        @eval export \$s
        #    end
        #end

        end # module
        """)
    end
end

function find_std_headers()
    headers = String[]
    @static if Sys.isapple()
        cmd_path = strip(read(`xcode-select --print-path`, String))
        sdk_path = joinpath(cmd_path, "SDKs", "MacOSX.sdk")
        xcode_path = cmd_path
        occursin("Xcode", xcode_path) && (xcode_path *= "/Toolchains/XcodeDefault.xctoolchain/")
        didfind = false
        lib = joinpath(xcode_path, "usr", "lib", "c++", "v1")
        inc = joinpath(xcode_path, "usr", "include", "c++", "v1")
        isdir(lib) && (push!(headers, lib); didfind = true;)
        isdir(inc) && (push!(headers, inc); didfind = true;)
        if isdir("/usr/include")
            push!(headers, "/usr/include")
        else isdir(joinpath(sdk_path, "usr", "include"))
            push!(headers, joinpath(sdk_path, "usr", "include"))
        end
        didfind || error("Could not find C++ standard library. Is XCode or CommandLineTools installed?")
    end
    return headers
end


text:
# Synchronization (B.6)

export sync_threads, sync_warp
export threadfence, threadfence_block, threadfence_system

"""
    sync_threads()

Waits until all threads in the thread block have reached this point and all global and
shared memory accesses made by these threads prior to `sync_threads()` are visible to all
threads in the block.
"""
@inline sync_threads() = ccall("llvm.nvvm.barrier0", llvmcall, Cvoid, ())

"""
    sync_warp(mask::Integer=0xffffffff)

Waits threads in the warp, selected by means of the bitmask `mask`, have reached this point
and all global and shared memory accesses made by these threads prior to `sync_warp()` are
visible to those threads in the warp. The default value for `mask` selects all threads in
the warp.

!!! note
   Requires CUDA >= 9.0 and sm_6.2
"""
sync_warp

if cuda_driver_version >= v"9.0" && v"6.0" in ptx_support
    @inline function sync_warp(mask::Integer=0xffffffff)
        @asmcall("bar.warp.sync \$0;", "r", true,
                 Cvoid, Tuple{UInt32}, convert(UInt32, mask))
    end
else
    @inline sync_warp(mask::Integer=0xffffffff) = nothing
end

"""
    threadfence_block()

A memory fence that ensures that:
- All writes to all memory made by the calling thread before the call to `threadfence_block()`
  are observed by all threads in the block of the calling thread as occurring before all writes
  to all memory made by the calling thread after the call to `threadfence_block()`
- All reads from all memory made by the calling thread before the call to `threadfence_block()`
  are ordered before all reads from all memory made by the calling thread after the call to `threadfence_block()`.
"""
@inline threadfence_block() = ccall("llvm.nvvm.membar.cta", llvmcall, Cvoid, ())

"""
    threadfence()

A memory fence that acts as [`threadfence_block`](@ref) for all threads in the block of the
calling thread and also ensures that no writes to all memory made by the calling thread after
the call to `threadfence()` are observed by any thread in the device as occurring before any
write to all memory made by the calling thread before the call to `threadfence()`.

Note that for this ordering guarantee to be true, the observing threads must truly observe the
memory and not cached versions of it; this is requires the use of volatile loads and stores,
which is not available from Julia right now.
"""
@inline threadfence() = ccall("llvm.nvvm.membar.gl", llvmcall, Cvoid, ())

"""
    threadfence_system()

A memory fence that acts as [`threadfence_block`](@ref) for all threads in the block of the
calling thread and also ensures that all writes to all memory made by the calling thread
before the call to `threadfence_system()` are observed by all threads in the device,
host threads, and all threads in peer devices as occurring before all writes to all
memory made by the calling thread after the call to `threadfence_system()`.
"""
@inline threadfence_system() = ccall("llvm.nvvm.membar.sys", llvmcall, Cvoid, ())


text:
import UUIDs: UUID, uuid1
import .ExpressionExplorer: SymbolsState, FunctionNameSignaturePair, FunctionName
import .Configuration
import .PkgCompat: PkgCompat, PkgContext
import Pkg
import TOML

const DEFAULT_NOTEBOOK_METADATA = Dict{String, Any}()

mutable struct BondValue
    value::Any
end
function Base.convert(::Type{BondValue}, dict::AbstractDict)
    BondValue(dict["value"])
end

const ProcessStatus = (
    ready="ready",
    starting="starting",
    no_process="no_process",
    waiting_to_restart="waiting_to_restart",
)

"Like a [`Diary`](@ref) but more serious. 📓"
Base.@kwdef mutable struct Notebook
    "Cells are ordered in a `Notebook`, and this order can be changed by the user. Cells will always have a constant UUID."
    cells_dict::Dict{UUID,Cell}
    cell_order::Vector{UUID}

    path::String
    notebook_id::UUID=uuid1()
    topology::NotebookTopology
    _cached_topological_order::Union{Nothing,TopologicalOrder}=nothing
    _cached_cell_dependencies_source::Union{Nothing,NotebookTopology}=nothing

    # buffer will contain all unfetched updates - must be big enough
    # We can keep 1024 updates pending. After this, any put! calls (i.e. calls that push an update to the notebook) will simply block, which is fine.
    # This does mean that the Notebook can't be used if nothing is clearing the update channel.
    pendingupdates::Channel=Channel(1024)

    executetoken::Token=Token()

    # per notebook compiler options
    # nothing means to use global session compiler options
    compiler_options::Union{Nothing,Configuration.CompilerOptions}=nothing
    # nbpkg_ctx::Union{Nothing,PkgContext}=nothing
    nbpkg_ctx::Union{Nothing,PkgContext}=PkgCompat.create_empty_ctx()
    nbpkg_ctx_instantiated::Bool=false
    nbpkg_restart_recommended_msg::Union{Nothing,String}=nothing
    nbpkg_restart_required_msg::Union{Nothing,String}=nothing
    nbpkg_terminal_outputs::Dict{String,String}=Dict{String,String}()
    nbpkg_busy_packages::Vector{String}=String[]
    nbpkg_installed_versions_cache::Dict{String,String}=Dict{String,String}()

    process_status::String=ProcessStatus.starting
    wants_to_interrupt::Bool=false
    last_save_time::typeof(time())=time()
    last_hot_reload_time::typeof(time())=zero(time())

    bonds::Dict{Symbol,BondValue}=Dict{Symbol,BondValue}()

    metadata::Dict{String, Any}=copy(DEFAULT_NOTEBOOK_METADATA)
end

_collect_cells(cells_dict::Dict{UUID,Cell}, cells_order::Vector{UUID}) =
    map(i -> cells_dict[i], cells_order)
_initial_topology(cells_dict::Dict{UUID,Cell}, cells_order::Vector{UUID}) =
    NotebookTopology(;
        cell_order=ImmutableVector(_collect_cells(cells_dict, cells_order)),
    )

function Notebook(cells::Vector{Cell}, @nospecialize(path::AbstractString), notebook_id::UUID)
    cells_dict=Dict(map(cells) do cell
        (cell.cell_id, cell)
    end)
    cell_order=map(x -> x.cell_id, cells)
    Notebook(;
        cells_dict,
        cell_order,
        topology=_initial_topology(cells_dict, cell_order),
        path,
        notebook_id
    )
end

Notebook(cells::Vector{Cell}, path::AbstractString=numbered_until_new(joinpath(new_notebooks_directory(), cutename()))) = Notebook(cells, path, uuid1())

function Base.getproperty(notebook::Notebook, property::Symbol)
    if property == :cells
        _collect_cells(notebook.cells_dict, notebook.cell_order)
    elseif property == :cell_inputs
        notebook.cells_dict
    else
        getfield(notebook, property)
    end
end

emptynotebook(args...) = Notebook([Cell()], args...)

function sample_notebook(name::String)
    file = project_relative_path("sample", name * ".jl")
    nb = load_notebook_nobackup(file)
    nb.path = tempname() * ".jl"
    nb
end

create_cell_metadata(metadata::Dict{String,<:Any}) = merge(DEFAULT_CELL_METADATA, metadata)
create_notebook_metadata(metadata::Dict{String,<:Any}) = merge(DEFAULT_NOTEBOOK_METADATA, metadata)
get_metadata(cell::Cell)::Dict{String,Any} = cell.metadata
get_metadata(notebook::Notebook)::Dict{String,Any} = notebook.metadata
get_metadata_no_default(cell::Cell)::Dict{String,Any} = Dict{String,Any}(setdiff(pairs(cell.metadata), pairs(DEFAULT_CELL_METADATA)))
get_metadata_no_default(notebook::Notebook)::Dict{String,Any} = Dict{String,Any}(setdiff(pairs(notebook.metadata), pairs(DEFAULT_NOTEBOOK_METADATA)))


text:
using DataDeps

register(DataDep(
  "letter-recognition",
  "http://archive.ics.uci.edu/ml/datasets/Letter+Recognition",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/letter-recognition/letter-recognition.data",
  "2b89f3602cf768d3c8355267d2f13f2417809e101fc2b5ceee10db19a60de6e2",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=1,
        feature_indices=2:17,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
# This file was generated, do not modify it. # hide
using PyPlot

figure(figsize=(8,6))
plot(X.Horsepower, y, ls="none", marker="o")

xlabel("Horsepower", fontsize=14)
xticks(50:50:250, fontsize=12)
yticks(10:10:50, fontsize=12)
ylabel("MPG", fontsize=14)

savefig(joinpath(@OUTPUT, "ISL-lab-5-g1.svg")) # hide

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
using Revise
using ADCME
using ADSeismic
using PyPlot
using DelimitedFiles

num_repeat = 3

open("elastic-time-$(has_gpu() ? "gpu" : "cpu").txt", "a") do io
    write(io, "===============================================\n")
end

for (k, scale) in enumerate(101:100:1001)

    param = ElasticPropagatorParams(NX=scale, NY=scale, NSTEP=100, DELTAT=1e-4,  DELTAX=1.0, DELTAY=1.0)

    source = Gauss(param, 1/pi, 100, 1e6) + 1
    srci = [div(param.NX,2)]
    srcj = [div(param.NY,2)]
    srctype = [0]
    srcv = reshape(source, :, 1)
    src = ElasticSource(srci, srcj, srctype, srcv)

    vp = 3300.
    vs = 3300. / 1.732
    rho = 2800.
    λ, ρ, μ = compute_lame_parameters(param.NX, param.NY, vp, vs, rho)
    model = ElasticPropagatorSolver(param, src, ρ, λ, μ)

    sess = Session(); init(sess)
    vs = run(sess, model.vx)

    time = 0
    for i = 1:num_repeat
      init(sess)
      time += (@timed run(sess, model.vx))[2]
    end
    time /= num_repeat

    @info time, scale
    open("elastic-time-$(has_gpu() ? "gpu" : "cpu").txt", "a") do io
        writedlm(io, [scale*scale time])
    end

end