text:
"""
parse_tuple(ps, ret)

`ret` is followed by a comma so tries to parse the rest of the
tuple.
"""
function parse_tuple end

iserrorwrapped(x, head) = x.head === head || (x.head === :errortoken && length(x.args) > 0 && iserrorwrapped(x.args[1], head))

function parse_tuple(ps::ParseState, ret::EXPR)
    op = EXPR(next(ps))
    if istuple(ret) && !(length(ret.trivia) > 1 && iserrorwrapped(last(ret.trivia), :RPAREN))
        if (isassignmentop(ps.nt) && kindof(ps.nt) != Tokens.APPROX)
            pushtotrivia!(ret, op)
        elseif closer(ps)
            @static if VERSION > v"1.1-"
                pushtotrivia!(ret, mErrorToken(ps, op, Unknown))
            else
                pushtotrivia!(ret, op)
            end
        else
            nextarg = @closer ps :tuple parse_expression(ps)
            if !(is_lparen(first(ret.trivia)))
                pushtotrivia!(ret, op)
                push!(ret, nextarg)
            else
                ret = EXPR(:tuple, EXPR[ret, nextarg], EXPR[op])
            end
        end
    else
        if (isassignmentop(ps.nt) && kindof(ps.nt) != Tokens.APPROX)
            ret = EXPR(:tuple, EXPR[ret], EXPR[op])
        elseif closer(ps)
            @static if VERSION > v"1.1-"
                ret = mErrorToken(ps, EXPR(:tuple, EXPR[ret], EXPR[op]), Unknown)
            else
                ret = EXPR(:tuple, EXPR[ret], EXPR[op])
            end
        else
            nextarg = @closer ps :tuple parse_expression(ps)
            ret = EXPR(:tuple, EXPR[ret, nextarg], EXPR[op])
        end
    end
    return ret
end

# XXX: Avert thine eyes.
function count_semicolons(ps, check_newline = true)
    dims = 0
    has_newline = false
    old_pos = position(ps.l.io)
    seek(ps.l.io, ps.ws.startbyte)
    while true
        c = readchar(ps.l.io)
        if c == ';'
            dims += 1
        elseif check_newline && c == '\n'
            # technically, only trailing newlines are allowed; we're a bit more lenient here
            has_newline = true
        else
            dims == 0 || break
        end
    end
    seek(ps.l.io, old_pos)
    return dims == 2 && has_newline ? 0 : dims
end

"""
    parse_array(ps)
Having hit '[' return either:
+ A vect
+ A vcat
+ A ncat
+ A comprehension
+ An array (vcat of hcats)
"""
function parse_array(ps::ParseState, isref = false)
    args = EXPR[]
    trivia = EXPR[EXPR(ps)]
    # [] or [;;;]
    if kindof(ps.nt) === Tokens.RSQUARE
        dims = 0
        if kindof(ps.ws) == SemiColonWS
            dims = count_semicolons(ps, false)
        end
        push!(trivia, accept_rsquare(ps))
        if dims > 0
            ret = EXPR(:ncat, args, trivia)
            pushfirst!(ret, EXPR(Symbol(dims), 0, 0, ""))
            return ret
        else
            return EXPR(:vect, args, trivia)
        end
    else
        ret = parse_array_outer(ps::ParseState, trivia, isref)
        pushtotrivia!(ret, accept_rsquare(ps))

        return ret
    end
end

binding_power(ps) =
    if kindof(ps.ws) == SemiColonWS
        -count_semicolons(ps)
    elseif kindof(ps.ws) == NewLineWS
        -1
    elseif kindof(ps.ws) == WS
        0
    else
        1
    end

function parse_array_outer(ps::ParseState, trivia, isref)
    args_list = EXPR[]
    trivia_bp = Int[]
    min_bp = 0
    is_start = true
    while kindof(ps.nt) !== Tokens.RSQUARE && kindof(ps.nt) !== Tokens.ENDMARKER
        a = @nocloser ps :newline @nocloser ps :newline @closesquare ps @closer ps :insquare @closer ps :ws @closer ps :wsop @closer ps :comma parse_expression(ps)

        if is_start
            args = EXPR[]
            if isref && _do_kw_convert(ps, a)
                a = _kw_convert(a)
            end
            if kindof(ps.nt) === Tokens.RSQUARE
                if headof(a) === :generator || headof(a) === :flatten
                    if isbinarycall(a.args[1]) && is_pairarrow(a.args[1].args[2])
                        return EXPR(:dict_comprehension, EXPR[a], trivia)
                    else
                        return EXPR(:comprehension, EXPR[a], trivia)
                    end
                elseif kindof(ps.ws) == SemiColonWS
                    dims = count_semicolons(ps)
                    push!(args, a)
                    if dims > 1
                        ret = EXPR(:ncat, args, trivia)
                        pushfirst!(ret, EXPR(Symbol(dims), 0, 0, ""))
                        return ret
                    else
                        return EXPR(:vcat, args, trivia)
                    end
                else
                    push!(args, a)
                    return EXPR(:vect, args, trivia)
                end
            elseif iscomma(ps.nt)
                args = EXPR[]
                push!(args, a)
                push!(trivia, accept_comma(ps))
                @closesquare ps parse_comma_sep(ps, args, trivia, isref, insert_params_at = 1)
                return EXPR(:vect, args, trivia)
            end
        end
        is_start = false

        push!(args_list, a)

        bp = binding_power(ps)
        bp < min_bp && (min_bp = bp)
        if bp <= 0
            push!(trivia_bp, bp)
        end
    end

    if length(trivia_bp) + 1 < length(args_list)
        return EXPR(:hcat, EXPR[EXPR(:errortoken, args_list, nothing)], trivia)
    end

    ret = _process_inner_array(args_list, trivia_bp)
    for t in trivia
        pushtotrivia!(ret, t)
    end

    if ret.head === :nrow
        if min_bp == 0
            popfirst!(ret.args)
            ret.head = :hcat
        elseif min_bp == -1
            popfirst!(ret.args)
            ret.head = :vcat
        else
            ret.head = :ncat
        end
    elseif min_bp == 0
        ret.head = :hcat
    elseif min_bp == -1
        ret.head = :vcat
    end

    return ret
end

function _process_inner_array(args_list, trivia_bp)
    if isempty(trivia_bp)
        if length(args_list) == 1
            arg = first(args_list)
            arg.args = something(arg.args, EXPR[])
            arg.trivia = something(arg.trivia, EXPR[])
            return arg
        else
            return EXPR(:errortoken, EXPR[], EXPR[])
        end
    end
    bp = minimum(trivia_bp)
    if all(==(bp), trivia_bp)
        if bp == 0
            ret = EXPR(:row, EXPR[], EXPR[])
        else
            ret = EXPR(:nrow, copy(args_list), EXPR[])
            pushfirst!(ret, EXPR(Symbol(-bp), 0, 0, ""))
        end
    end
    i = 1
    i0 = 1
    if bp == 0
        ret = EXPR(:row, EXPR[], EXPR[])
    else
        ret = EXPR(:nrow, EXPR[], EXPR[])
        push!(ret, EXPR(Symbol(-bp), 0, 0, ""))
    end
    while true
        i = findnext(==(bp), trivia_bp, i0)
        i === nothing && (i = length(args_list))
        if !checkbounds(Bool, args_list, i0:i) || !checkbounds(Bool, trivia_bp, i0:(i-1))
            return EXPR(:errortoken, args_list, EXPR[])
        end
        inner_args = args_list[i0:i]
        inner_trivia = trivia_bp[i0:(i-1)]
        i0 = i + 1

        push!(ret, _process_inner_array(inner_args, inner_trivia))
        i >= length(args_list) && break
    end
    return ret
end

"""
    parse_ref(ps, ret)

Handles cases where an expression - `ret` - is followed by
`[`. Parses the following bracketed expression and modifies it's
`.head` appropriately.
"""
function parse_ref(ps::ParseState, ret::EXPR)
    next(ps)
    ref = @closer ps :inref @nocloser ps :inwhere parse_array(ps, true)
    if headof(ref) === :vect
        args = EXPR[ret]
        for a in ref.args
            push!(args, a)
        end
        return EXPR(:ref, args, ref.trivia)
    elseif headof(ref) === :hcat
        args = EXPR[ret]
        for a in ref.args
            push!(args, a)
        end
        return EXPR(:typed_hcat, args, ref.trivia)
    elseif headof(ref) === :vcat
        args = EXPR[ret]
        for a in ref.args
            push!(args, a)
        end
        return EXPR(:typed_vcat, args, ref.trivia)
    elseif headof(ref) === :ncat
        args = EXPR[ret]
        for a in ref.args
            push!(args, a)
        end
        return EXPR(:typed_ncat, args, ref.trivia)
    else
        args = EXPR[ret]
        for a in ref.args
            push!(args, a)
        end
        return EXPR(:typed_comprehension, args, ref.trivia)
    end
end

"""
parse_curly(ps, ret)

Parses the juxtaposition of `ret` with an opening brace. Parses a comma
seperated list.
"""
function parse_curly(ps::ParseState, ret::EXPR)
    args = EXPR[ret]
    trivia = EXPR[EXPR(next(ps))]
    parse_comma_sep(ps, args, trivia, true)
    accept_rbrace(ps, trivia)
    return EXPR(:curly, args, trivia)
end

function parse_braces(ps::ParseState)
    return @default ps @nocloser ps :inwhere parse_barray(ps)
end

function parse_barray(ps::ParseState)
    args = EXPR[]
    trivia = EXPR[EXPR(ps)]

    if kindof(ps.nt) === Tokens.RBRACE
        accept_rbrace(ps, trivia)
        ret = EXPR(:braces, args, trivia)
    else
        first_arg = @nocloser ps :newline @closebrace ps  @closer ps :ws @closer ps :wsop @closer ps :comma parse_expression(ps)
        if kindof(ps.nt) === Tokens.RBRACE
            push!(args, first_arg)
            if kindof(ps.ws) == SemiColonWS
                pushfirst!(args, EXPR(:parameters, EXPR[], nothing))
            end
            accept_rbrace(ps, trivia)
            ret = EXPR(:braces, args, trivia)
        elseif iscomma(ps.nt)
            push!(args, first_arg)
            accept_comma(ps, trivia)
            @closebrace ps parse_comma_sep(ps, args, trivia, true, insert_params_at = 1)
            accept_rbrace(ps, trivia)
            return EXPR(:braces, args, trivia)
        elseif kindof(ps.ws) == NewLineWS
            ret = EXPR(:bracescat, args, trivia)
            push!(ret, first_arg)
            prevpos = position(ps)
            while kindof(ps.nt) != Tokens.RBRACE && kindof(ps.nt) !== Tokens.ENDMARKER
                a = @closebrace ps  parse_expression(ps)
                push!(ret, a)
                prevpos = loop_check(ps, prevpos)
            end
            pushtotrivia!(ret, accept_rbrace(ps))
            return ret
        elseif kindof(ps.ws) == WS || kindof(ps.ws) == SemiColonWS
            first_row = EXPR(:row, EXPR[first_arg])
            prevpos = position(ps)
            while kindof(ps.nt) !== Tokens.RBRACE && kindof(ps.ws) !== NewLineWS && kindof(ps.ws) !== SemiColonWS && kindof(ps.nt) !== Tokens.ENDMARKER
                a = @closebrace ps @closer ps :ws @closer ps :wsop parse_expression(ps)
                push!(first_row, a)
                prevpos = loop_check(ps, prevpos)
            end
            if kindof(ps.nt) === Tokens.RBRACE && kindof(ps.ws) != SemiColonWS
                if length(first_row.args) == 1
                    first_row = EXPR(:bracescat, first_row.args)
                end
                push!(args, first_row)
                push!(trivia, INSTANCE(next(ps)))
                return EXPR(:bracescat, args, trivia)
            else
                if length(first_row.args) == 1
                    first_row = first_row.args[1]
                else
                    first_row = EXPR(:row, first_row.args)
                end
                ret = EXPR(:bracescat, EXPR[first_row], trivia)
                prevpos = position(ps)
                while kindof(ps.nt) != Tokens.RBRACE
                    if kindof(ps.nt) === Tokens.ENDMARKER
                        break
                    end
                    first_arg = @closebrace ps @closer ps :ws @closer ps :wsop parse_expression(ps)
                    push!(ret, EXPR(:row, EXPR[first_arg]))
                    prevpos1 = position(ps)
                    while kindof(ps.nt) !== Tokens.RBRACE && kindof(ps.ws) !== NewLineWS && kindof(ps.ws) !== SemiColonWS && kindof(ps.nt) !== Tokens.ENDMARKER
                        a = @closebrace ps @closer ps :ws @closer ps :wsop parse_expression(ps)
                        push!(last(ret.args), a)
                        prevpos1 = loop_check(ps, prevpos1)
                    end
                    # if only one entry dont use :row
                    if length(last(ret.args).args) == 1
                        ret.args[end] = setparent!(ret.args[end].args[1], ret)
                    end
                    update_span!(ret)
                    prevpos = loop_check(ps, prevpos)
                end
                pushtotrivia!(ret, accept_rbrace(ps))
                update_span!(ret)
                return ret
            end
        else
            ret = EXPR(:braces, args, trivia)
            push!(ret, first_arg)
            pushtotrivia!(ret, accept_rbrace(ps))
        end
    end
    return ret
end


text:
export bissect_nwt
function bissect_nwt(h       :: AbstractNLPModel,
                          nlpstop :: AbstractStopping;
                          t1 = h.meta.x0, t2 = 100.0,
                          verbose = false)

  (tₐ, tᵦ) = trouve_intervalle(h, t1, 2.5)
  γ = 0.8; t = tᵦ; tₚ = tₐ;
  hₖ = [0.0]; hplus = [0.0]
  iter = 0

  gₖ = grad(h, t)

  OK = update_and_start!(nlpstop, x = t, gx = gₖ, g0 = copy(gₖ))

  verbose &&
    @printf(" iter        tₚ        t         dN         gₖ          ")
  verbose &&
    @printf("gplus        \n")
  verbose &&
    @printf(" %7.2e %7.2e %7.2e  %7.2e  %7.2e  %7.2e\n",
                iter, tₚ[1], t[1], 0.0, gₖ[1], 0.0)
  while !OK

      kₖ = hess(h, t)
      dN = vec(-gₖ ./ kₖ) #direction de Newton
      if (true in ((tₚ .- t) .* dN .> 0.0)) && (true in (dN ./ (tₚ - t) .< γ))
        tplus = t .+ dN
        gplus = grad(h, tplus)
      else
        tplus = (t .+ tₚ) ./ 2
        gplus = grad(h, tplus)
      end

      if true in (t .> tₚ)
        if true in (gplus .< 0.0)
          tₚ = t;  t = tplus
        else
           t = tplus
        end
      else
        if true in (gplus .> 0.0)
          tₚ = t;  t = tplus
        else
          t = tplus
        end
      end

      # updating values
      hₖ = hplus; gₖ = gplus
      OK = update_and_stop!(nlpstop, x = tₚ, gx = gₖ)
      iter = iter + 1
      verbose &&
        @printf(" %7.2e %7.2e %7.2e  %7.2e  %7.2e  %7.2e\n",
                iter, tₚ[1], t[1], dN[1], gₖ[1], gplus[1])
    end

  optimal = nlpstop.meta.optimal
  return optimal, nlpstop
end

ext: jl


text:
"""
    methodlearner(method, data, backbone; <kwargs>...) -> Learner

Create a `Learner` to train a model for learning method `method` using
`data`. See also [`Learner`](#).
"""
function methodlearner(
        method::LearningMethod,
        data,
        backbone,
        callbacks...;
        isbackbone = true,
        pctgval = 0.2,
        batchsize = 16,
        validbsfactor = 2,
        optimizer = ADAM(),
        lossfn = methodlossfn(method),
        dlkwargs = (;),
    )
    model = isbackbone ? methodmodel(method, backbone) : backbone
    dls = methoddataloaders(data, method, batchsize;
        validbsfactor = validbsfactor, pctgval = pctgval)
    return Learner(model, dls, optimizer, lossfn, callbacks...)
end


text:
if pwd()[end - 2:end] != "src"
    cd("src")
end

module Arrhenius
    using YAML
    using NPZ
    using LinearAlgebra
    using SparseArrays

    include("Constants.jl")
    include("DataStructure.jl")
    include("Solution.jl")
    include("Magic.jl")
    include("Thermo.jl")
    include("Kinetics.jl")
    include("Transport.jl")
end

```
Following codes are used during development phase only.
```
using LinearAlgebra
using DifferentialEquations
using Sundials
using Profile
using ForwardDiff
using SparseArrays
using PyCall
using Test

R = Arrhenius.R
one_atm = Arrhenius.one_atm

mech = "../mechanism/h2o2.yaml"
gas = Arrhenius.CreateSolution(mech);
ns = gas.n_species

Y0 = zeros(ns);
Y0[Arrhenius.species_index(gas, "H2")] = 0.1;
Y0[Arrhenius.species_index(gas, "O2")] = 0.2;
Y0[Arrhenius.species_index(gas, "AR")] = 0.7;

P = Arrhenius.one_atm
T = 1200.0
Y = Y0;
mean_MW = 1. / dot(Y, 1 ./ gas.MW)
ρ_mass = P / R / T * mean_MW
X =  Arrhenius.Y2X(gas, Y, mean_MW);
C =  Arrhenius.Y2C(gas, Y, ρ_mass);

η_mix, λ_mix, Dkm = Arrhenius.mix_trans(gas, P, T, X, mean_MW)

## Call Cantera
ct = pyimport("cantera")
ct_gas = ct.Solution(mech)

ct_gas.TPY = T, P, "H2:0.1,O2:0.2,AR:0.7"

@test η_mix ≈ ct_gas.viscosity rtol = 1.e-4

@test λ_mix ≈ ct_gas.thermal_conductivity rtol = 1.e-4

@test Dkm ≈ ct_gas.mix_diff_coeffs rtol = 1.e-3

println("Dkm from Arrhenius.jl | Cantera")
hcat(Dkm, ct_gas.mix_diff_coeffs)


text:
# file for unit tests

using DrWatson
quickactivate(@__DIR__)

using PolarizationFramework
using Test

using DifferentialEquations

@testset "BinaryAttributes" begin
    g = 5
    b = BinaryAttributes(g)
    nodes = 3
    attr = get_attributes(b, nodes)

    @test size(attr) == (nodes, g) #test rozmiaru
    @test all(abs.(attr) .== 1) #test wartosci powinny byc +-1

    attr = ones(nodes, g)
    attr[1, :] .= -1
    attr[2, 1:2] .= -1
    #attr is: [-1, -1, -1, -1, -1;
    #          -1, -1, 1, 1, 1;
    #           1, 1, 1, 1, 1]
    #This should give following weights:
    # [0.0  -0.2  -1.0
    #  0.0   0.0   0.2
    #  0.0   0.0   0.0]
    exp_weights = zeros(3, 3)
    exp_weights[1, 2] = -0.2
    exp_weights[1, 3] = -1.0
    exp_weights[2, 3] = 0.2

    weights = get_attribute_layer_weights(b, attr)

    @test size(weights) == (nodes, nodes)
    @test weights == exp_weights
end

@testset "OrderedAttributes" begin
    g = 5
    threshold = 0.55
    v = 3
    b = OrderedAttributes(g, threshold, v)
    nodes = 3
    attr = get_attributes(b, nodes)

    @test min(attr...) == 1 #test rozmiaru
    @test max(attr...) == v #test wartosci powinny byc +-1

    attr = ones(Int, nodes, g)
    attr[1, :] .= 3
    attr[2, 1:2] .= 3
    attr[2, 3:5] .= 2
    #attr is:
    # [3  3  3  3  3
    #  3  3  2  2  2
    #  1  1  1  1  1]
    #This should give following weights:
    # [0.0  0.5  -0.9
    #  0.0  0.0  -0.3
    #  0.0  0.0   0.0]
    exp_weights = zeros(3, 3)
    exp_weights[1, 2] = 0.5
    exp_weights[1, 3] = -0.9
    exp_weights[2, 3] = -0.3

    weights = get_attribute_layer_weights(b, attr)

    @test size(weights) == (nodes, nodes)
    @test weights ≈ exp_weights
end

@testset "UnorderedAttributes" begin
    g = 5
    threshold = 0.55
    v = 3
    b = UnorderedAttributes(g, threshold, v)
    nodes = 3
    attr = get_attributes(b, nodes)

    @test min(attr...) == 1 #test rozmiaru
    @test max(attr...) == v #test wartosci powinny byc +-1

    attr = ones(Int, nodes, g)
    attr[1, 1:4] .= 3
    attr[2, 1:3] .= 3
    attr[2, 4:5] .= 2
    #attr is:
    # [3  3  3  3  1
    #  3  3  3  2  2
    #  1  1  1  1  1]
    #This should give following weights:
    # [0.0  0.1  -0.7
    #  0.0  0.0  -1.0
    #  0.0  0.0   0.0]
    exp_weights = zeros(3, 3)
    exp_weights[1, 2] = 0.1
    exp_weights[1, 3] = -0.7
    exp_weights[2, 3] = -1.1

    weights = get_attribute_layer_weights(b, attr)

    @test size(weights) == (nodes, nodes)
    @test weights ≈ exp_weights
end

@testset "Simple simulation run" begin
    g = 5
    attr = BinaryAttributes(g)
    (ishb_sim, t, u, u0, xy_attr, sol) = calc_heider_attr(
        45,
        attr,
        0.0,
        1000.0,
        Heider7!,
        AutoTsit5(Rodas5(autodiff = false)),
        true,
    )
end

@testset "Simulation with initial conditions" begin
    n = 15
    g = 5
    gamma = 2.7
    attr = BinaryAttributes(g)

    rl_bal = init_random_balanced_relations(n)

    val0_attr = get_attributes(attr, n)
    al_w = get_attribute_layer_weights(attr, val0_attr)

    (ishb_sim, t, u, u0, xy_attr, sol) = calc_heider_attr(
        n,
        attr,
        0.0,
        1000.0,
        Heider7!,
        AutoTsit5(Rodas5(autodiff = false)),
        false,
    )
    (ishb_sim, t, u, u0, xy_attr, sol) = calc_heider_attr(
        n,
        attr,
        gamma,
        1000.0,
        Heider7!,
        AutoTsit5(Rodas5(autodiff = false)),
        false,
    )
    (ishb_sim, t, u, u0, xy_attr, sol) = calc_heider_attr(
        n,
        attr,
        gamma,
        1000.0,
        Heider7!,
        AutoTsit5(Rodas5(autodiff = false)),
        false,
        rl_bal,
    )
    (ishb_sim, t, u, u0, xy_attr, sol) = calc_heider_attr(
        n,
        attr,
        gamma,
        1000.0,
        Heider7!,
        AutoTsit5(Rodas5(autodiff = false)),
        false,
        (),
        al_w,
    )
    (ishb_sim, t, u, u0, xy_attr, sol) = calc_heider_attr(
        n,
        attr,
        gamma,
        1000.0,
        Heider7!,
        AutoTsit5(Rodas5(autodiff = false)),
        false,
        rl_bal,
        al_w,
    )
end

@testset "Series of simulations" begin
    n = 5
    g = 5
    attr = BinaryAttributes(g)
    gammas = [0.5, 1.5, 3.5]
    zmax = 10
    maxtime = 1000.0
    ode_fun_name = "Heider7!"
    disp_each = 0.5
    disp_more_every = 60
    save_each = 60
    files_folder = ["test"]
    filename_prefix = "SoSNumerTest"

    using_heider_attr(
        n,
        attr,
        gammas,
        zmax,
        maxtime,
        ode_fun_name;
        disp_each = disp_each,
        disp_more_every = disp_more_every,
        save_each = save_each,
        files_folder = files_folder,
        filename_prefix = filename_prefix,
    )

    filename_prefix = "SoSDestabTest"
    larger_size = 3

    using_heider_attr_destab(
        n,
        attr,
        gammas,
        larger_size,
        zmax,
        maxtime,
        ode_fun_name;
        disp_each = disp_each,
        disp_more_every = disp_more_every,
        save_each = save_each,
        files_folder = files_folder,
        filename_prefix = filename_prefix,
    )
end

@testset "Longer series of simulations" begin
    n = 9
    g = 5
    attr = BinaryAttributes(g)
    gammas = [0.5, 1.5, 2.0, 3.5]
    zmax = 1000
    maxtime = 1000.0
    ode_fun_name = "Heider7!"
    disp_each = 0.5
    disp_more_every = 60
    save_each = 60
    files_folder = ["test"]
    filename_prefix = "LSoSNumerTest"

    using_heider_attr(
        n,
        attr,
        gammas,
        zmax,
        maxtime,
        ode_fun_name;
        disp_each = disp_each,
        disp_more_every = disp_more_every,
        save_each = save_each,
        files_folder = files_folder,
        filename_prefix = filename_prefix,
    )

    filename_prefix = "LSoSDestabTest"
    larger_size = Int((n + 1) / 2 + 2)
    zmax = 300

    using_heider_attr_destab(
        n,
        attr,
        gammas,
        larger_size,
        zmax,
        maxtime,
        ode_fun_name;
        disp_each = disp_each,
        disp_more_every = disp_more_every,
        save_each = save_each,
        files_folder = files_folder,
        filename_prefix = filename_prefix,
    )
end

max_stars_repo_name: pjgorski/PolarizationFramework

max_issues_repo_name: pjgorski/PolarizationFramework

max_forks_repo_name: pjgorski/PolarizationFramework


text:
## Preloads
module eqws

    using FFTW
    using LinearAlgebra
    using SparseArrays
    using Statistics
    using Test
    using DSP
    using Interpolations
    using StatsBase
    using NNlib
    using SparseArrays

    import CUDA
    using CUDA.CUSPARSE

    export
        fink_filter_bank,
        fink_filter_bank_p,
        fink_filter_list,
        fink_filter_hash,
        fink_filter_hash_p,
        fink_filter_bank_3dizer,
        wst_S1_deriv,
        wst_S20_deriv,
        wst_S20_deriv_sum,
        eqws_compute,
        eqws_compute_p,
        eqws_compute_3d,
        eqws_compute_RGB,
        eqws_compute_wrapper,
        eqws_compute_convmap,
        apodizer,
        wind_2d,
        transformMaker,
        S1_iso_matrix,
        S2_iso_matrix,
        S1_equiv_matrix,
        S2_equiv_matrix,
        S1_iso_matrix3d,
        S2_iso_matrix3d,
        pre_subdivide,
        post_subdivide,
        sub_strip,
        renorm

    include("compute.jl")
    include("filter_bank.jl")
    include("derivatives.jl")
    include("post_process.jl")
    include("pre_process.jl")

end # of module

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
using GeomDAE.Fields: Field, Q, P, V, F

"Parameters for right-hand side function of implicit partitioned Runge-Kutta methods."
type NonlinearFunctionParametersIPRK{T} <: NonlinearFunctionParameters{T}
    fields::Tuple{Vararg{Field}}
    functs::Tuple{Vararg{Field}}

    fq::Function
    fp::Function
    fv::Function
    ff::Function

    Δt::T
    t::T

    d::Int
    s::Int
    a_q::Matrix{T}
    a_p::Matrix{T}
    c_q::Vector{T}
    c_p::Vector{T}

    q::Vector{T}
    p::Vector{T}
    y::Vector{T}
    z::Vector{T}

    Q::Matrix{T}
    V::Matrix{T}
    P::Matrix{T}
    F::Matrix{T}

    Y::Array{T,3}

    tQ::Vector{T}
    tV::Vector{T}
    tP::Vector{T}
    tF::Vector{T}

    function NonlinearFunctionParametersIPRK(fields, functs, fq, fp, fv, ff, Δt, d, s, a_q, a_p, c_q, c_p)
        # create solution vectors
        q = zeros(T,d)
        p = zeros(T,d)
        y = zeros(T,d)
        z = zeros(T,d)

        # create internal stage vectors
        Q = zeros(T,d,s)
        V = zeros(T,d,s)
        P = zeros(T,d,s)
        F = zeros(T,d,s)

        # create temporary vectors
        Y  = zeros(T,2,d,s)

        tQ = zeros(T,d)
        tV = zeros(T,d)
        tP = zeros(T,d)
        tF = zeros(T,d)

        new(fields, functs, fq, fp, fv, ff, Δt, 0, d, s, a_q, a_p, c_q, c_p, q, p, y, z, Q, V, P, F, Y, tQ, tV, tP, tF)
    end
end

"Compute stages of implicit partitioned Runge-Kutta methods."
@generated function function_stages!{T}(y::Vector{T}, b::Vector{T},
                             params::NonlinearFunctionParametersIPRK{T},
                             fields::Field...)

    local fargs = []

     if Q ∈ fields
         push!(fargs, :(params.tQ))
     end

     if V ∈ fields
         push!(fargs, :(params.tV))
     end

     if P ∈ fields
         push!(fargs, :(params.tP))
     end

     if F ∈ fields
         push!(fargs, :(params.tF))
     end

     # TODO Below, only the case (fp, ff) is implemented. All other combinations
     #      need to be implemented as well.

     code = quote
         local tqᵢ::T
         local tpᵢ::T

         # loop through stages
         for i in 1:params.s
             for k in 1:params.d
                 # copy y to Y
                 params.Y[1,k,i] = y[2*(params.d*(i-1)+k-1)+1]
                 params.Y[2,k,i] = y[2*(params.d*(i-1)+k-1)+2]

                 # compute Q and V
                 params.Q[k,i] = params.q[k] + params.Δt * params.Y[1,k,i]
                 params.V[k,i] = params.Y[2,k,i]
             end

             # compute f(X)
             tqᵢ = params.t + params.Δt * params.c_q[i]
             tpᵢ = params.t + params.Δt * params.c_p[i]
             simd_copy_xy_first!(params.tQ, params.Q, i)
             simd_copy_xy_first!(params.tV, params.V, i)
             params.fp(tpᵢ, params.tQ, params.tV, params.tP)
             params.ff(tpᵢ, params.tQ, params.tV, params.tF)
             simd_copy_yx_first!(params.tP, params.P, i)
             simd_copy_yx_first!(params.tF, params.F, i)
         end

         # compute b = - [(Y-AV), (P-AF)]
         for i in 1:params.s
             for k in 1:params.d
                 b[2*(params.d*(i-1)+k-1)+1] = - params.Y[1,k,i]
                 b[2*(params.d*(i-1)+k-1)+2] = - params.P[k,i] + params.p[k]
                 for j in 1:params.s
                     b[2*(params.d*(i-1)+k-1)+1] += params.a_q[i,j] * params.V[k,j]
                     b[2*(params.d*(i-1)+k-1)+2] += params.a_p[i,j] * params.F[k,j] * params.Δt
                 end
             end
         end
    end

    return code
end

"Compute stages of implicit partitioned Runge-Kutta methods."
function function_stages!{T}(y::Vector{T}, b::Vector{T},
                             params::NonlinearFunctionParametersIPRK{T})
    function_stages!(y, b, params, params.fields...)
end

"Implicit partitioned Runge-Kutta integrator."
immutable IntegratorIPRK{T, ST} <: Integrator{T}
    equation::SODE{T}
    tableau::TableauIPRK{T}
    Δt::T

    solver::ST

    q::Array{T,1}
    p::Array{T,1}
    y::Array{T,1}
    z::Array{T,1}
    Q::Array{T,2}
    V::Array{T,2}
    P::Array{T,2}
    F::Array{T,2}
end

function IntegratorIPRK{T}(equation::SODE{T}, tableau::TableauIPRK{T}, Δt::T; nonlinear_solver=QuasiNewtonSolver)
    D = equation.d
    S = tableau.s

    # create solution vector for internal stages / nonlinear solver
    z = zeros(T, 2*D*S)

    # create params
    params = NonlinearFunctionParametersIPRK{T}(equation.fields, equation.functs,
                                                equation.q, equation.p, equation.v, equation.f,
                                                Δt, D, S,
                                                tableau.a_q, tableau.a_p, tableau.c_q, tableau.c_p)

    # create solver
    solver = nonlinear_solver(z, params)

    # create integrator
    IntegratorIPRK{T, typeof(solver)}(equation, tableau, Δt, solver,
                                      params.q, params.p, params.y, params.z,
                                      params.Q, params.V, params.P, params.F)
end

"Integrate ODE with special partitioned Runge-Kutta integrator."
function integrate!(int::IntegratorIPRK, sol::SolutionPODE)
    local nt::Int

    # copy initial conditions from solution
    simd_copy_xy_first!(int.q, sol.q, 1)
    simd_copy_xy_first!(int.p, sol.p, 1)

    for n in 1:sol.ntime
        # compute initial guess
        # TODO
        for i in 1:int.tableau.s
            for k in 1:int.equation.d
                int.solver.x[2*(sol.nd*(i-1)+k-1)+1] = int.q[k]
                # TODO initial guess for velocity
                # int.solver.x[2*(params.d*(i-1)+k-1)+1] = 0.
            end
        end

        # call nonlinear solver
        solve!(int.solver)

        if !solverStatusOK(int.solver.status, int.solver.params)
            println(int.solver.status)
        end

        # compute final update
        simd_mult!(int.y, int.V, int.tableau.b_q)
        simd_mult!(int.z, int.F, int.tableau.b_p)
        simd_axpy!(int.Δt, int.y, int.q)
        simd_axpy!(int.Δt, int.z, int.p)

        # copy to solution
        if mod(n, sol.nsave) == 0
            nt = div(n, sol.nsave)+1
            simd_copy_yx_first!(int.q, sol.q, nt)
            simd_copy_yx_first!(int.p, sol.p, nt)
        end
    end
    nothing
end


text:
abstract type AbstractWatson <: AbstractFloat end

struct DrWatson64{f} <: AbstractWatson
    val::Float64
end

# conversions back from Sherlog64
Base.Float64(x::DrWatson64) = x.val
Base.Float32(x::DrWatson64) = Float32(x.val)
Base.Float16(x::DrWatson64) = Float16(x.val)
Base.Int64(x::DrWatson64) = Int64(x.val)
Base.Int32(x::DrWatson64) = Int32(x.val)
Base.Int16(x::DrWatson64) = Int16(x.val)

# conversion among DrWatsons
DrWatson64{f}(x::DrWatson64) where f = DrWatson64{f}(x.val)
DrWatson64(x::DrWatson64{f}) where {f} = DrWatson64{f}(x.val)

Base.promote_rule(::Type{Int64},::Type{DrWatson64{f}}) where f = DrWatson64{f}
Base.promote_rule(::Type{Int32},::Type{DrWatson64{f}}) where f = DrWatson64{f}
Base.promote_rule(::Type{Int16},::Type{DrWatson64{f}}) where f = DrWatson64{f}

Base.promote_rule(::Type{Float64},::Type{DrWatson64{f}}) where f = DrWatson64{f}
Base.promote_rule(::Type{Float32},::Type{DrWatson64{f}}) where f = DrWatson64{f}
Base.promote_rule(::Type{Float16},::Type{DrWatson64{f}}) where f = DrWatson64{f}

Base.bitstring(x::DrWatson64) = bitstring(x.val)

Base.eps(::Type{DrWatson64{f}}) where f = eps(Float64)
Base.eps(x::DrWatson64) = eps(x.val)

Base.typemin(::DrWatson64{f}) where f = DrWatson64{f}(typemin(Float64))
Base.typemax(::DrWatson64{f}) where f = DrWatson64{f}(typemax(Float64))
Base.floatmin(::DrWatson64{f}) where f = DrWatson64{f}(floatmin(Float64))
Base.floatmax(::DrWatson64{f}) where f = DrWatson64{f}(floatmax(Float64))
Base.precision(::DrWatson64{f}) where f = DrWatson64{f}(precision(Float64))

function +(x::DrWatson64{f},y::DrWatson64{f}) where f
    r = x.val + y.val
    scent(f,r)
    return DrWatson64{f}(r)
end

function -(x::DrWatson64{f},y::DrWatson64{f}) where f
    r = x.val - y.val
    scent(f,r)
    return DrWatson64{f}(r)
end

function *(x::DrWatson64{f},y::DrWatson64{f}) where f
    r = x.val * y.val
    scent(f,r)
    return DrWatson64{f}(r)
end

function /(x::DrWatson64{f},y::DrWatson64{f}) where f
    r = x.val / y.val
    scent(f,r)
    return DrWatson64{f}(r)
end

function ^(x::DrWatson64{f},y::DrWatson64{f}) where f
    r = x.val ^ y.val
    scent(f,r)
    return DrWatson64{f}(r)
end

for O in ( :(-), :(+),
           :sign,
           :prevfloat, :nextfloat,
           :round, :trunc, :ceil, :floor,
           :inv, :abs, :sqrt, :cbrt,
           :exp, :expm1, :exp2, :exp10,
           :log, :log1p, :log2, :log10,
           :rad2deg, :deg2rad, :mod2pi, :rem2pi,
           :sin, :cos, :tan, :csc, :sec, :cot,
           :asin, :acos, :atan, :acsc, :asec, :acot,
           :sinh, :cosh, :tanh, :csch, :sech, :coth,
           :asinh, :acosh, :atanh, :acsch, :asech, :acoth,
           :sinc, :sinpi, :cospi,
           :sind, :cosd, :tand, :cscd, :secd, :cotd,
           :asind, :acosd, :atand, :acscd, :asecd, :acotd
          )
    @eval begin
        function Base.$O(x::DrWatson64{f}) where f
            r = $O(x.val)
            scent(f,r)
            return DrWatson64{f}(r)
        end
    end
end

for O in ( :(<), :(<=))     # do not trigger scent for those
    @eval begin
        Base.$O(x::DrWatson64, y::DrWatson64) = $O(x.val, y.val)
    end
end


text:
using MLKernels, MLLabelUtils

@testset "oc_classifier" begin
    dummy_data = [1.0 2.0 3.0; 4.0 5.0 6.0]
    pools = fill(:U, size(dummy_data, 2))

    init_strategy = SVDD.FixedParameterInitialization(MLKernels.GaussianKernel(0.5), 0.5)

    models = []
    svdd_neg = SVDD.SVDDneg(dummy_data, pools)
    SVDD.initialize!(svdd_neg, init_strategy)
    SVDD.fit!(svdd_neg, TEST_SOLVER)
    push!(models, svdd_neg)

    vanilla_svdd = SVDD.VanillaSVDD(dummy_data)
    SVDD.initialize!(vanilla_svdd, init_strategy)
    SVDD.fit!(vanilla_svdd, TEST_SOLVER)
    push!(models, vanilla_svdd)

    ssad = SVDD.SSAD(dummy_data, pools)
    SVDD.initialize!(ssad, init_strategy)
    SVDD.fit!(ssad, TEST_SOLVER)
    push!(models, ssad)

    update_data = [3.0 2.0 3.0; 4.0 5.0 6.0]
    pools_updated = copy(pools)
    pools_updated[1] = :Lout

    @testset "instantiate" begin
        actual = SVDD.instantiate(SVDD.VanillaSVDD, dummy_data, pools, Dict(:C => 0.5))
        @test actual.C == 0.5
        actual = SVDD.instantiate(SVDD.SVDDneg, dummy_data, pools, Dict(:C1 => 0.5))
        @test actual.C1 == 0.5
        actual = SVDD.instantiate(SVDD.SSAD, dummy_data, pools, Dict(:C1 => 0.5))
        @test actual.C1 == 0.5
    end

    @testset "set_param!" begin
        actual = deepcopy(vanilla_svdd)
        @test_throws ArgumentError SVDD.set_param!(actual, Dict(:C => "0.4", :Ψ => "ipd"))
        @test_throws ArgumentError SVDD.set_param!(actual, Dict(:C => "0.4"))
        SVDD.set_param!(actual, Dict(:C => 0.5))
        @test actual.state == SVDD.model_fitted
        SVDD.set_param!(actual, Dict(:C => 0.4))
        @test actual.state == SVDD.model_initialized
        @test actual.C == 0.4
        @test_throws ArgumentError SVDD.set_param!(actual, Dict(:C => 1.1))
    end

    for m in models
        @testset "get_kernel $(typeof(m))" begin
            actual = SVDD.get_kernel(m)
            if isdefined(m, :kernel_fct)
                @test m.kernel_fct == actual
            end
        end
    end

    for m in models
        @testset "set_adjust_kernel $(typeof(m))" begin
            actual = deepcopy(m)
            @assert actual.state == SVDD.model_fitted
            @assert !SVDD.is_K_adjusted(actual)
            SVDD.set_adjust_K!(actual, false)
            @test actual.state == SVDD.model_fitted
            @test !isdefined(actual, :K_adjusted)
            SVDD.set_adjust_K!(actual, true)
            @test isdefined(actual, :K_adjusted)
            @test size(actual.K) == size(actual.K_adjusted)
            @test actual.state == SVDD.model_initialized
        end
    end

    for m in models
        @testset "model updates $(typeof(m))" begin
            @testset "no change" begin
                actual = deepcopy(m)
                @assert actual.state == SVDD.model_fitted
                SVDD.set_pools!(actual, pools)
                SVDD.set_pools!(actual, labelmap(pools))
                @test actual.state == SVDD.model_fitted
                SVDD.update_K!(actual)
                @test actual.state == SVDD.model_fitted
                SVDD.set_data!(actual, actual.data)
                @test actual.state == SVDD.model_fitted
                SVDD.set_kernel!(actual, actual.kernel_fct)
                @test actual.state == SVDD.model_fitted
            end

            @testset "invalidate model after change" begin
                actual = deepcopy(m)
                SVDD.set_data!(actual, update_data)
                @test actual.state == SVDD.model_initialized
                @test actual.data == update_data

                if (typeof(actual) <: SVDD.SSAD || typeof(actual) <: SVDD.SVDDneg)
                    actual = deepcopy(m)
                    SVDD.set_pools!(actual, pools_updated)
                    @test actual.state == SVDD.model_initialized

                    actual = deepcopy(m)
                    SVDD.set_pools!(actual, labelmap(pools_updated))
                    @test actual.state == SVDD.model_initialized
                end

                actual = deepcopy(m)
                SVDD.set_kernel!(actual, LinearKernel())
                @test actual.state == SVDD.model_initialized
                @test actual.K != m.K
                @test actual.kernel_fct == LinearKernel()
            end
        end
    end
end


text:
module ITSE

export itse

"""
    ```julia
    itse(error, time)
    ```
    Integral of Time-Weighted Squared Error (ITSE) fromm error vector.
    Inputs:
    error -> error vector.
    time  -> time vector.
"""
function itse(error::Vector{T}, time::Vector{T}) where T
    sum(time .* (error.^2))
end

"""
    ```julia
    itse(error, time)
    ```
    Integral of Time-Weighted Squared Error (ITSE) from vector of error vector.
    Inputs:
    error -> vector of error vector.
    time  -> time vector.
"""
function itse(error::Vector{Vector{T}}, time::Vector{T}) where T
    map(x->sum(time .* (x.^2)), error)
end

"""
    ```julia
    itse(x, t, desired)
    ```
    Integral of Time-Weighted Squared Error (ITSE) from error betwwen desired value and controler output.
    Inputs:
    x -> controler output vector;
    t -> time vector;
    desired -> desired value function.
"""
function itse(x::Vector{T}, t::Vector{T}, desired::Function) where T
    x = x .- map(desired, t)
    sum( t .* (x.^2))
end

end


text:
"""
    lorenz_key(x_init, y_init, z_init, num_keys; α=10.0, ρ=28.0, β=2.667, dt=0.01, scaling_factor=10.0^16, upper_bound=256.0)

Generates 3 vectors of pseudo-random numbers using Lorenz system of differential
equations.

The equations -

\$\\frac{dx}{dt} = α * (y - x)\$
\$\\frac{dy}{dt} = x * (ρ - z) - y\$
\$\\frac{dz}{dt} = x * y - β * z\$

# Arguments
- `x_init::Float64`: Initial value of x.
- `y_init::Float64`: Initial value of y.
- `z_init::Float64:` Initial value of z.
- `num_keys::Int64`: Number of keys (in a single list) to be generated.
- `α::Float64`: Constant associated with Lorenz system of differential equations.
- `ρ::Float64`: Constant associated with Lorenz system of differential equations.
- `β::Float64`: Constant associated with Lorenz system of differential equations.
- `scaling_factor::Float64=10.0^16`: Factor to be multiplied to the generated value of pseudo-random
    number. Ideally, the factor should be > upper_bound.
- `upper_bound::Float64=256.0`: Upper bound of keys (not included). Use 256 for encrypting images
    as the RGB values of a pixel varies from 0 to 255.

# Returns
- `x::Vector{Int64}`: Generated pseudo-random keys corresponding to x values.
- `y::Vector{Int64}`: Generated pseudo-random keys corresponding to y values.
- `z::Vector{Int64}`: Generated pseudo-random keys corresponding to z values.

# Example
```jldoctest
julia> lorenz_key(0.01, 0.02, 0.03, 20)
([0, 0, 256, 24, 129, 42, 54, 134, 43, 179, 85, 19, 24, 44, 71, 210, 238, 152, 22, 27], [0, 0, 240, 55, 25, 163, 89, 243, 123, 5, 197, 64, 227, 54, 188, 226, 154, 134, 64, 69], [0, 0, 80, 227, 178, 204, 89, 33, 144, 139, 105, 208, 108, 155, 61, 254, 57, 102, 149, 47])
```
"""
function lorenz_key(
    x_init::Float64,
    y_init::Float64,
    z_init::Float64,
    num_keys::Int64;
    α::Float64=10.0,
    ρ::Float64=28.0,
    β::Float64=2.667,
    dt::Float64=0.01,
    scaling_factor::Float64=10.0^16,
    upper_bound::Float64=256.0
)
    # Initializing 3 empty lists
    x = zeros(Float64, num_keys)
    y = zeros(Float64, num_keys)
    z = zeros(Float64, num_keys)

    # Initializing initial values
    x[1], y[1], z[1] = x_init, y_init, z_init

    # System of equations
    for i = 1:num_keys - 1
        x[i + 1] = x[i] + (α * (y[i] - x[i]) * dt)
        y[i + 1] = y[i] + ((x[i] * (ρ - z[i]) - y[i]) * dt)
        z[i + 1] = z[i] + ((x[i] * y[i] - β * z[i]) * dt)
    end

    for i = 1:length(x)
        x[i] = x[i] * scaling_factor % upper_bound
        y[i] = y[i] * scaling_factor % upper_bound
        z[i] = z[i] * scaling_factor % upper_bound
    end

    x = round.(Int, x)
    y = round.(Int, y)
    z = round.(Int, z)

    return x, y, z
end


text:
module Utilities

using Logging
using MLBase
using LinearAlgebra
using PyCall; pygui(:qt); using PyPlot

using ..AReg

logit(x) = log(x/(1-x))
export logit

inv_logit(x) = 1/(1+exp(-x))
export inv_logit

link_R_in_0_1(x) = inv_logit(x)
export link_R_in_0_1

inv_link_R_in_0_1(x) = logit(x)
export inv_link_R_in_0_1

link_R_in_R_pos(x) = exp(x)
export link_R_in_R_pos

inv_link_R_in_R_pos(x) = log(x)
export inv_link_R_in_R_pos

function make_pos_def( M; negOrSmallTh = 100*eps(), smallVal = 100*eps(), warnTh = 0.01)
    eig = eigen(M)
    minEigenVal = minimum(eig.values)

    if minEigenVal < negOrSmallTh
        Logging.@info("negative eigenvalue $minEigenVal")
        absMinEigVal = abs(minEigenVal)

        deltaEigenVal = absMinEigVal + smallVal

        absMinEigVal > warnTh ? Logging.@warn("large negative eigenvalue $minEigenVal") : ()

        eig.values[eig.values.< smallVal] .+= deltaEigenVal

        M = Matrix(eig)
    end

    return M, minEigenVal
end
export make_pos_def

drop_nan_col(x::Matrix) = x[:, .!dropdims(any(isnan.(x), dims=1), dims=1) ]
export drop_nan_col

drop_nan_row(x::Matrix) = x[.!dropdims(any(isnan.(x), dims=2), dims=2), : ]
export drop_nan_row

#Funzioni che uso spesso
sumSq(array :: AbstractArray,Dim::Int) = dropdims(sum(array,dims = Dim),dims =Dim)
sumSq(array :: AbstractArray) = sumSq(array,1)
export sumSq

meanSq(array::AbstractArray,Dim::Int) = dropdims(mean(array,dims = Dim),dims =Dim)
export meanSq

splitVec(Vec::AbstractArray{<:Any,1}) = (Nhalf = round(Int,length(Vec)/2); (Nhalf,Vec[1:Nhalf],Vec[Nhalf+1:2Nhalf]) )
export splitVec

splitMat(Mat::AbstractArray{<:Any,2}) = (Nhalf = round(Int,length(Mat[1,:])/2); (Nhalf,Mat[:,1:Nhalf],Mat[:,Nhalf+1:2Nhalf]) )
export splitMat

"""put to zero diagonal of matrices"""
putZeroDiag!(Matr::Array{T,2} where T<: Real) =   for i=1:length(Matr[:,1]) Matr[i,i] = 0 end
export putZeroDiag!

#putZeroDiag!(Matr::Symmetric{T,Array{T,2}} where T<: Real) =   for i=1:length(Matr[:,1]) Matr[i,i] = 0 end
putZeroDiag(Matr::Array{T,2} where T<: Real) = (tmp = Matr; putZeroDiag!(tmp); tmp)
#putZeroDiag(Matr::Symmetric{T,Array{T,2}} where T<: Real) = (tmp = Matr; putZeroDiag!(tmp); tmp)
putZeroDiag!(Matr::BitArray{2}) =   for i=1:length(Matr[:,1]) Matr[i,i] = false end
putZeroDiag(Matr::BitArray{2}) = (tmp = Matr; putZeroDiag!(tmp); tmp)
export putZeroDiag

putZeroDiag_no_mut(Matr::Array{T,2} where T<: Real) = (Matr - Diagonal(Diagonal(Matr)))
export putZeroDiag_no_mut

function sortrowsUtilities(A::Matrix{<:Number},rowInd::Int;rev = true)
    tmp  = sortperm(A[:,rowInd],rev=rev)
    return A[tmp,:]
end
export sortrowsUtilities

#define an adjacency matrix that has falses on diagonals: central upper and lower
tmpAjacMatTridiagFalses(N::Int) = (tmp = trues(N,N);for i=1:N tmp[i,i] = false; i>1 ? tmp[i,i-1]= 0 : (); i<N ? tmp[i,i+1]=0 : ()  end;tmp)
export tmpAjacMatTridiagFalses

function tmpArrayAdjmatrEqual(N::Int,T::Int;matGenFun::Function = tmpAjacMatTridiagFalses)
    out = trues(N,N,T)
    mat = matGenFun(N)
    for t=1:T
        out[:,:,t] = mat
    end
    return out
end
export tmpArrayAdjmatrEqual

function distributeAinB!(A::Array{<:Real,1},B::Array{<:Real,1})
    NA = length(A)
    NB = length(B)
    GroupSize  = floor(NB/NA)
    for i =0:NA-1 B[1 + Int64( GroupSize *i):Int64(GroupSize*(i+1))] = A[i+1] end;
    B[Int64(GroupSize*(NA-1))+1:end] = A[NA];
end
export distributeAinB!

function distributeAinVecN(A::Array{<:Real,1},N::I where I<: Int64)
    if length(A) == 0
        return zeros(Int64,N);
    else
        vecN = zeros(typeof(A[1]),N);
        distributeAinB!(A,vecN);
        return vecN
    end
end
distributeAinVecN(A::UnitRange,N::I where I<: Int64) = distributeAinVecN(Array{Int,1}(A),N)
export distributeAinVecN

function rocCurve(realVals::BitArray{1},foreVals::Vector{Float64}; Th::Vector{Float64}=Vector(0:0.00001:1),plotFlag = true)
    #Given binary observations and estimated probabilityes plot the roc and return true positives and false positives
    rocOut = roc(realVals,foreVals,Th)
    Nth = length(Th)
    tpr = [rocOut[i].tp/rocOut[i].p for i=1:Nth]
    fpr = [rocOut[i].fp/rocOut[i].n for i=1:Nth]
    auc = sum(-Base.diff(fpr).*tpr[1:end-1])
    if plotFlag
        plot(fpr,tpr)
        xlim([0,1])
        ylim([0,1])
        title("AUC =  $(round(auc,3))")
        grid()
    end
    return tpr, fpr,auc
end
export rocCurve

function collapseArr3(arr)
     N1,N2,N3 = length(arr),length(arr[1]),length(arr[1][1])
     tmp = zeros(N1,N2,N3)
     for i=1:N1,j=1:N2,k=1:N3
         tmp[i,j,k] = arr[i][j][k]
     end
     tmp = permutedims(tmp,[2,1,3])
     return tmp
end
export collapseArr3

function randSteps(startVal,endVal,Nsteps::Int,T::Int;rand=true)
     out = zeros(T)
     if Nsteps>1
         heights = linspace(startVal,endVal,Nsteps)
         if rand
             heights = shuffle(heights)
         end
         Tstep = round(Int,T/Nsteps)
         last = 0
         for i=1:Nsteps-1
             out[1+last:last+Tstep] = heights[i]
             last +=Tstep
         end
         out[last:end] = heights[Nsteps]
     elseif Nsteps == 1
         out =startVal
     end
     return out
end
export randSteps

function dgpSin(minVal::Ty  where Ty <:Real  ,maxVal::Ty  where Ty <:Real ,Ncycles,T::Int; phase =0)
     out = zeros(T)
     maxVal<=minVal ? error() : ()
     medVal = (maxVal + minVal)./2
     amplitude = abs.(maxVal-minVal)
     out = (medVal .+ (amplitude/2) .* sin.(2π*Ncycles/T .*  Vector(1:T) .+ phase))
     return out
end
export dgpSin

function dgpAR(mu,B,sigma,T::Int;minMax=zeros(2),scaling = "uniform")
  ar1 = AReg.ARp([mu*(1-B) , B],sigma,[mu])
  path = AReg.simulate(ar1,T)

  if minMax !=zeros(2)
      min = minimum(minMax)
      max = maximum(minMax)
      if scaling == "uniform"
          minPath = minimum(path)
          maxPath = maximum(path)
          Δ = (max-min)/( maxPath - minPath)
          rescPath = min .+ (path .- minPath).*Δ
      elseif scaling == "nonlinear"
          rescPath = min + (max - min)* 1 ./(1 + exp.(path))
      end

  else
      rescPath = path
  end
  return rescPath
end
export dgpAR

function drop_bad_un_estimates(inputEstMat)

    nanInds = dropdims(any(isnan.(inputEstMat), dims=1), dims=1)
    infInds = dropdims(any(.!isfinite.(inputEstMat), dims=1), dims=1)
    integrInds = dropdims(any(inputEstMat[2:3:6,:] .> 20, dims=1), dims=1)
    constInds = dropdims(any(inputEstMat[3:3:6,:] .< log(0.0005), dims=1), dims=1)

    dropInds = nanInds .| infInds  .| integrInds .| constInds
    @show (sum(nanInds), sum(infInds), sum(integrInds), sum(constInds))

    inputEstMat[:, .!dropInds]
end
export drop_bad_un_estimates

""" return an edge list with nodes indexed starting from 1 to the number of nodes present in the list """

function edge_list_pres(list)
    pres = sort(unique(reshape(list, :)))
    list_pres =  [findfirst(isequal(n), pres) for n in list]
    sorted_list_pres = sortslices(list_pres, dims=1)
    return sorted_list_pres
end
export edge_list_pres

unconditional_mean_auto_regr(w, B) = w./(1 .- B)
export unconditional_mean_auto_regr

"""
Barrier function to handle single N and sequence of Ns
"""
get_N_t(N_seq::AbstractArray, t::Int) = N_seq[t]
get_N_t(N::Int, t) = N
export get_N_t

end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "AzStorage"
version = v"0.4.0"

# Collection of sources required to build AzStorage
sources = [
    GitSource(
        "https://github.com/ChevronETC/AzStorage.jl.git",
        "2d45d02ac9a7b36a1e35e3e46dd54c73895a2c74"
    )
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/AzStorage.jl/src

# We need to tell the makefile where to find libssh2 on windows
if [[ ${target} == *mingw* ]]; then
    export LDFLAGS="-L${libdir}"
fi

make

cp libAzStorage.so ${libdir}/libAzStorage.${dlext}
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)

# The products that we will ensure are always built
# TODO - add libgomp dependency
products = [
    LibraryProduct("libAzStorage", :libAzStorage)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("CompilerSupportLibraries_jll"),
    Dependency("LibCURL_jll", v"7.73.0"),
    # MbedTLS is only an indirect dependency (through LibCURL), but we want to
    # be sure to have the right version of MbedTLS for the corresponding version
    # of Julia.
    BuildDependency(PackageSpec(; name="MbedTLS_jll", version="2.24.0")),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using Measurements, SpecialFunctions, QuadGK, Calculus
using Test, LinearAlgebra, Statistics

import Base: isapprox
import Measurements: value, uncertainty

# A few tests before overloading `isapprox`.
@testset "isapprox" begin
    @test (5 ± 1) ≈ (5 ± 2)
    @test (nextfloat(5.0) ± 1) ≈ 5 ± 1
    @test 5.1 ± 1 ≉ 5 ± 1
end

isapprox(x::Measurement, y::Measurement; rest...) =
    isapprox(x.val, y.val; nans = true, rest...) &&
    isapprox(x.err, y.err; nans = true, rest...)
isapprox(x::Complex{Measurement{<:AbstractFloat}},
         y::Complex{Measurement{<:AbstractFloat}}; rest...) =
             isapprox(real(x), real(y); nans = true, rest...) &&
             isapprox(imag(x), imag(y); nans = true, rest...)

w = -0.5 ± 0.03
x = 3 ± 0.1
y = 4 ± 0.2
z = complex(x)

@testset "meaningful" begin
    @test meaningful(w) == -0.5 ± 0.03
    @test meaningful(x) == 3 ± 0.1
    @test meaningful(y) == 4 ± 0.2
    @test meaningful(z) == meaningful(x)

    # test rounding up of uncertainty
    @test meaningful(0.5 ± 0.00100001) == 0.5 ± 0.0011
    @test meaningful(0.5 ± 0.00109999) == 0.5 ± 0.0011
    @test meaningful(0.5 ± 0.009900001) == 0.5 ± 0.01
    @test meaningful(0.5 ± 10000.00100001) == 0.0 ± 11000
    @test meaningful(0.5 ± 10999) == 0.0 ± 11000
    @test meaningful(0.5 ± 99001) == 0.0 ± 100000

    # test rounding of value
    @test meaningful(0.5 ± 1.1) == 0.5 ± 1.1
    @test meaningful(0.5 ± 0.1) == 0.5 ± 0.1
    @test meaningful(0.5 ± 0.01) == 0.5 ± 0.01
    @test meaningful(0.5001 ± 11) == 1.0 ± 11
    # we don't make any assumption about behavious of 0.5 rounding
    @test meaningful(0.4999 ± 11) == 0.0 ± 11
    @test meaningful(0.5 ± 110) == 0.0 ± 110
    @test meaningful(0.9999 ± 99.1) == 0.0 ± 100

    # test complex rounding
    @test meaningful(0.5 ± 1.1 + im*(0.5 ± 0.1)) == 0.5 ± 1.1 + im*(0.5 ± 0.1)
end

@testset "Standard Score" begin
    @test stdscore(w, x.val) ≈ -350/3
    @test stdscore(x, y) ≈ -4.472135954999579
    @test stdscore(w, y) ≈ stdscore(w - y, 0)
    @test stdscore(y, 4.1) ≈ stdscore(y, 4.1 ± 0)
end

@testset "measurement" begin
    @test measurement(x) === x
    @test measurement.(1:2:5, 2:2:6) == [1±2, 3±4, 5±6]
    @test iszero(@inferred(measurement(1)).err)
    @test @inferred(measurement(pi)).tag === UInt64(0)
    @test length(@inferred(measurement(4//5)).der) == 0
    @test @inferred(measurement(ℯ, big(0))).tag === UInt64(0)
    @test length(@inferred(measurement(7//3, 0)).der) == 0
    @test length(@inferred(measurement(4.7f0, 0.3)).der) == 1
    @test typeof(@inferred(measurement(1, big(0)))) == Measurement{BigFloat}
    @test typeof(@inferred(Measurement{Float64}(big"3.14"))) == Measurement{Float64}
    @test typeof(@inferred(Measurement{BigFloat}(-5.43f2))) == Measurement{BigFloat}
end

@testset "Weighted Average" begin
    @test weightedmean((w, x, y)) ≈ measurement(-0.12584269662921355, 0.028442727788398632)
end

@testset "Derivative" begin
    c = 1 ± 0
    @test Measurements.derivative(3*x^2, (x.val, x.err, x.tag)) ≈ 18
    @test Measurements.derivative(3*x^2, x) ≈ 18
    @test Measurements.derivative.(2x + y - w, [x, y, w]) ≈ [2, 1, -1]
    @test Measurements.derivative.(x - x + y - y + w - w + c, [x, y, w, c]) ≈ [0, 0, 0, 0]
    @test length((x - x + y - y + w - w + c).der) == 0
end

@testset "Contributions to uncertainty" begin
    @test sort(collect(values(Measurements.uncertainty_components(w * x * y)))) ≈
        [0.2, 0.3, 0.36]
end

@testset "Conversion and Promotion" begin
    @test convert(Measurement{Float64}, pi) == pi ± 0
    @test convert(Measurement{Float64}, 1//2) == 0.5 ± 0
    @test convert(Measurement{Float64}, 3) == 3.0 ± 0.0
    @test convert(Measurement{Float64}, 3 ± 1) == 3.0 ± 1.0
    for a in (w, x, y)
        @test convert(typeof(a), a) === a
        @test convert(Measurement, a) === a
        @test float(a) === a
    end
    @test convert(Measurement, pi) == pi ± 0
    @test convert(Measurement, 1//2) == 0.5 ± 0
    @test convert(Measurement, 3) == 3 ± 0
    @test convert(Int, measurement(3, 0)) === 3
    @test_throws AssertionError convert(Int, measurement(3, 1))
    @test float(3 ± 1) == 3.0 ± 1.0
    @test promote(Measurement{Float32}(3 ± 0.5), Measurement{Float64}(4 ± 0)) ==
        (Measurement{Float64}(3 ± 0.5), Measurement{Float64}(4 ± 0))
    @test promote(x, complex(7)) == (complex(3.0 ± 0.1),
                                     complex(measurement(7.0)))
    @test promote(complex(0, 1 ± 0), 2.1 ± 0.2) ==
        (complex(0, 1.0 ± 0), complex(2.1 ± 0.2))
    @test promote_type(Measurement{Float32}, Float64) ==
        Measurement{Float64}
    @test promote_type(Measurement{BigFloat}, Measurement{Float64}) ==
        Measurement{BigFloat}
end

@testset "Comparisons and Tests" begin
    # isapprox is overloaded only for tests.  Make sure it works as expected.
    @test 5 ± 1 ≉ 5 ± 2
    @test (nextfloat(5.0) ± nextfloat(1.0)) ≈ 5 ± 1
    @test 5.1 ± 1 ≉ 5 ± 1
    @test z == x ≠ y
    @test x == x
    @test -2 < w <= x < y < 5
    @test big(1) < big(x) < big(4) <= big(y) <= big(5)
    @test -5//3 < w <= -1//2 <= x < 11//3 < y <= 24//6
    @test 3 == 3 ± 0 ≠ x == 3 ± 0.1
    @test 4 ± 0.2 == y ≠ 4 == 4 ± 0
    @test measurement(big"0.75", 0.01) == 3//4 ± 1//100
    @test measurement(big(π)) ≠ π
    @test ℯ ≠ measurement(Float32(ℯ))
    @test 3//4 == measurement(Float32(0.75), Float32(0)) ≠ 4//3
    @test isnan(x) == false
    @test isfinite(y) == true && isfinite(measurement(Inf)) == false
    @test isinf(measurement(Inf)) == true && isinf(x) == false
    @test isinteger(-3 ± 0)
    @test !isinteger(3 ± 1)
    @test !isinteger(pi ± 0)
    @test iszero(0 ± 0)
    @test iszero(zero(Measurement))
    @test !iszero(0 ± 1)
    @test !iszero(1 ± 0)
    @test isone(1 ± 0)
    @test isone(one(Measurement))
    @test !isone(1 ± 1)
    @test !isone(0 ± 0)
end

##### Mathematical Operations
@testset "Addition" begin
    @test @inferred(x + y) ≈ measurement(7, 0.22360679774997896)
    @test @inferred(2 + x) ≈ measurement(5, 0.1)
    @test @inferred(x + 5//2) ≈ measurement(5.5, 0.1)
    @test @inferred(x + true) ≈ measurement(4, 0.1)
    @test @inferred(x + 2) ≈ measurement(5, 0.1)
    @test @inferred(x + big(1)) == (x.val + big(1)) ± big(x.err)
    for a in (w, x, y)
        @test @inferred(+a) === a
        @test @inferred(a + zero(a)) ≈ a ≈ @inferred(measurement(0) + a) # Neutral element
        @test @inferred(a + a + a) ≈ 3a # Correlation
        for b in (w, x, y)
            @test @inferred(a + b) ≈ @inferred(b + a) # Commutativity
        end
    end
    # Test a weird operation just to be sure it doesn't throw an error.
    @test @inferred((x ± y) + (w ± x)) ≈ @inferred((x + w) ± hypot(y, x))
end

@testset "Subtraction" begin
    @test @inferred(-x) ≈ measurement(-3, 0.1)
    @test @inferred(x - y) ≈ measurement(-1, 0.22360679774997896)
    @test @inferred(2 - x) ≈ measurement(-1, 0.1)
    @test @inferred(x - 2) ≈ measurement(1, 0.1)
    # Test correlation
    for a in (w, x, y); @test @inferred(a - a) ≈ zero(a) ≈ measurement(0); end
end

@testset "Multiplication" begin
    @test @inferred(x * y) ≈ measurement(12, 0.7211102550927979)
    @test @inferred(2x) ≈ measurement(6, 0.2)
    @test @inferred(x * 3) ≈ measurement(9, 0.3)
    for a in (w, x, y)
        @test @inferred(a*(0 ± 1)) ≈ measurement(0, abs(a.val))
        @test @inferred(a * one(a)) ≈ a ≈ @inferred(measurement(1) * a) # Neutral element
        @test @inferred(a * a * a) ≈ a^3 # Correlation
        for b in (w, x, y)
            @test @inferred(a * b) ≈ @inferred(b * a) # Commutativity
        end
        @test @inferred(true * a) ≈ a
        @test @inferred(a * false) ≈ measurement(0)
        @test @inferred(a * 0) ≈ measurement(0)
    end
    for f in (muladd, fma)
        @test @inferred(f(w, x, y)) ≈ w*x + y
        @test @inferred(f(w, x, 1)) ≈ w*x + one(y)
        @test @inferred(f(w, 2, y)) ≈ w*2 + y
        @test @inferred(f(w, 2, 1)) ≈ w*2 + one(x)
        @test @inferred(f(3, x, y)) ≈ 3*x + y
        @test @inferred(f(3, x, 1)) ≈ 3*x + one(w)
        @test @inferred(f(3, 2, y)) ≈ 3*2 + y
    end
end

@testset "Division" begin
    @test @inferred(x / y) ≈ measurement(0.75, 0.04506939094329987)
    @test @inferred(x / 10) ≈ measurement(0.3, 0.01)
    @test @inferred(1 / y) ≈ measurement(0.25, 0.0125)
    @test @inferred(div(1.2*x, w)) ≈ measurement(-7)
    @test @inferred(div(x, 1.2)) ≈ measurement(2)
    @test @inferred(div(9.4, y)) ≈ measurement(2)
    @test @inferred(fld(1.2*x, w)) ≈ measurement(-8)
    @test @inferred(fld(x, 1.2)) ≈ measurement(2)
    @test @inferred(fld(9.4, y)) ≈ measurement(2)
    @test @inferred(cld(1.2*x, w)) ≈ measurement(-7)
    @test @inferred(cld(x, 1.2)) ≈ measurement(3)
    @test @inferred(cld(9.4, y)) ≈ measurement(3)
    for a in (w, x, y)
        @test @inferred(0 / a) ≈ measurement(0)
        @test @inferred((0 ± 1) / a) ≈ measurement(0, 1/abs(a.val))
        # Test correlation
        @test @inferred(a / a) ≈ one(a) ≈ measurement(1)
        # Test derivatives of "div", "fld", and "cld".  They're defined to be exactly 0.
        # Should you discover this is not correct, update the test accordingly.
        for f in (div, fld, cld), b in (w, x, y)
            @test Measurements.derivative.(f(a, b), [a, b]) == [0.0, 0.0]
        end
    end
end

@testset "Inverse" begin
    for a in (w, x, y); @test @inferred(inv(a)) ≈ 1/a; end
end

@testset "signbit" begin
    @test signbit(x) == false
    @test signbit(w) == true
end

@testset "Power" begin
    # Test numerical values
    @test @inferred(x^y) ≈ measurement(81, 20.818061515800505)
    @test @inferred(x^2) ≈ measurement(9, 0.6)
    @test @inferred(y^(1//2)) ≈ measurement(2, 0.05)
    @test @inferred(x^(pi)) ≈ measurement(31.54428070019754, 3.3033093503504967)
    @test @inferred(2^x) ≈ measurement(8, 0.5545177444479562)
    @test @inferred(pi^x) ≈ measurement(31.006276680299816, 3.5493811564854525)
    @test @inferred(z ^ 2.5) ≈ @inferred(x ^ 2.5)
    @test @inferred(z ^ 3) ≈ @inferred(x ^ 3)
    for a in (w, x, y)
        @test @inferred(x ^ a) ≈ @inferred(exp(a * log(x)))
        @test @inferred(abs(a) ^ w) ≈ @inferred(exp(w * log(abs(a))))
        @test @inferred(abs(a) ^ a) ≉ @inferred(exp(a * log(abs(a.val) ± a.err))) # correlation
        @test @inferred(a ^ (-1)) ≈ @inferred(inv(a))
        @test @inferred(a ^ 2) ≈ @inferred(a ^ 2.0)
        @test @inferred(a ^ (4//2)) ≈ @inferred(a * a) # correlation
        @test @inferred(2 ^ a) ≈ @inferred(2.0 ^ a)
        @test @inferred(ℯ ^ a) ≈ @inferred(exp(a))
        @test @inferred(exp2(a)) ≈ @inferred(2 ^ a)
    end
    # Make sure "p ± 0" behaves like "p", in particular with regard to the
    # uncertainty.
    for p in (-3, 0, 3); @test @inferred((0 ± 0.1)^(p ± 0)) ≈ (0 ± 0.1)^p; end
end

@testset "Deg and rad" begin
    @test @inferred(deg2rad(y)) ≈ 0.06981317007977318 ± 0.003490658503988659
    @test @inferred(rad2deg(x)) ≈ measurement(171.88733853924697, 5.729577951308232)
    for a in (w, x, y); @inferred(rad2deg(deg2rad(a))) ≈ a; end
end

@testset "Trig functions" begin
    @test @inferred(cos(x)) ≈ measurement(-0.9899924966004454, 0.014112000805986721)
    @test @inferred(cosd(x)) ≈ measurement(0.9986295347545738, 9.134347536190512e-5)
    @test @inferred(cosh(x)) ≈ measurement(10.067661995777765, 1.0017874927409902)
    @test @inferred(sin(y)) ≈ measurement(-0.7568024953079282, 0.13072872417272238)
    @test @inferred(sind(y)) ≈ measurement(0.0697564737441253, 0.0034821554353128255)
    @test @inferred(sinh(y)) ≈ measurement(27.28991719712775, 5.461646567203298)
    for a in (w, x, y)
        @test @inferred(cos(a) ^ 2 + sin(a) ^ 2) ≈ one(a)
        @test @inferred(tan(a))  ≈ sin(a)  / cos(a)
        @test @inferred(tand(a)) ≈ sind(a) / cosd(a)
        @test @inferred(tanh(a)) ≈ sinh(a) / cosh(a)
        @test @inferred(sinpi(a))≈ sin(pi * a) atol = 1e-15
        @test @inferred(cospi(a))≈ cos(pi * a) atol = 1e-15
        @test @inferred(sinc(a)) ≈ (sin(pi * a) / (pi * a)) atol = 1e-16
        @test @inferred(cosc(a)) ≈ ((cos(pi * a) - sin(pi * a) / (pi * a)) / a) atol = 1e-15
        if isdefined(Base, :sincos)
            # Check we got the sign of derivatives in `sincos` right.
            s, c = @inferred(sincos(a))
            @test s + c ≈ sin(a) + cos(a)
            @test s ^ 2 + c ^ 2 ≈ one(a)
        end
    end
    for c in (-1 ± 1, 0 ± 1, 1 ± 1)
        @test sinc(c) ≈ @uncertain(sinc(c)) rtol = 1e-7
        @test cosc(c) ≈ @uncertain(cosc(c)) rtol = 1e-7
        @test Measurements.value(sinc(c) ^ 2 + cosc(c) ^ 2) == 1
    end
end

@testset "Inverse trig functions" begin
    @test @inferred(acos(w)) ≈ measurement(2.0943951023931957, 0.034641016151377546)
    @test @inferred(acosd(w)) ≈ measurement(120, 1.9847840235184515)
    @test @inferred(acosh(x)) ≈ measurement(1.7627471740390859, 0.035355339059327376)
    @test @inferred(cos(acos(w))) ≈ w
    @test @inferred(cosd(acosd(w))) ≈ w
    @test @inferred(cosh(acosh(x))) ≈ x
    @test @inferred(asin(w)) ≈ measurement(-0.5235987755982989, 0.034641016151377546)
    @test @inferred(asind(w)) ≈ measurement(-30, 1.9847840235184515)
    @test @inferred(asinh(x)) ≈ measurement(1.8184464592320668, 0.0316227766016838)
    @test @inferred(sin(asin(w))) ≈ w
    @test @inferred(sind(asind(w))) ≈ w
    @test @inferred(sinh(asinh(x))) ≈ x
    @test @inferred(atan(w)) ≈ measurement(-0.4636476090008061, 0.024)
    @test @inferred(atand(w)) ≈ measurement(-26.56505117707799, 1.3750987083139758)
    @test @inferred(atanh(w)) ≈ measurement(-0.5493061443340548, 0.04)
    for a in (w, x, y); @test @inferred(tan(atan(a))) ≈ a; end
    for a in (w, x, y); @test @inferred(tand(atand(a))) ≈ a; end
    for a in (x, y)
        @test @inferred(sec(asec(a))) ≈ a
        @test @inferred(csc(acsc(a))) ≈ a
        @test @inferred(cot(acot(a))) ≈ a
        @test @inferred(coth(acoth(a))) ≈ a
    end
    @test @inferred(tanh(atanh(w))) ≈ w
    @test @inferred(sech(asech(-w))) ≈ -w
    @test @inferred(csch(acsch(w))) ≈ w
    @test @inferred(atan(x, y)) ≈ measurement(0.6435011087932844, 0.028844410203711916)
    @test @inferred(atan(x, 5)) ≈ measurement(0.5404195002705842, 0.014705882352941178)
    @test @inferred(atan(-3, y)) ≈ measurement(-0.6435011087932844, 0.024)
end

@testset "Reciprocal trig functions" begin
    @test @inferred(csc(y)) ≈ measurement(-1.3213487088109024, 0.228247438348944)
    @test @inferred(cscd(y)) ≈ measurement(14.335587026203676, 0.7156144742178738)
    @test @inferred(csch(y)) ≈ measurement(0.03664357032586561, 0.007333632734561779)
    @test @inferred(sec(w)) ≈ measurement(1.139493927324549, 0.018675251089778414)
    @test @inferred(secd(w)) ≈ measurement(1.000038078385737, 4.5695512845605615e-6)
    @test @inferred(sech(w)) ≈ measurement(0.886818883970074, 0.012294426649942352)
    @test @inferred(cot(x)) ≈ measurement(-7.015252551434534, 5.021376836040872)
    @test @inferred(cotd(x)) ≈ measurement(19.08113668772821, 0.6372018679183983)
    @test @inferred(coth(x)) ≈ measurement(1.0049698233136892, 0.000996434577114765)
    for a in (w, x, y)
        @test @inferred(csc(a)) ≈ 1/sin(a)
        @test @inferred(cscd(a)) ≈ 1/sin(deg2rad(a))
        @test @inferred(csch(a)) ≈ 1/sinh(a)
        @test @inferred(sec(a)) ≈  1/cos(a)
        @test @inferred(secd(a)) ≈ 1/cos(deg2rad(a))
        @test @inferred(sech(a)) ≈ 1/cosh(a)
        @test @inferred(cot(a)) ≈ 1/tan(a)
        @test @inferred(cotd(a)) ≈ 1/tan(deg2rad(a))
        @test @inferred(coth(a)) ≈ 1/tanh(a)
    end
end

@testset "Exponentials" begin
    @test @inferred(exp(x)) ≈ measurement(20.085536923187668, 2.008553692318767)
    for a in (w, 3//5*w, x/10, x, y/50, y)
        @test @inferred(expm1(a)) ≈ exp(a) - one(a)
        @test @inferred(exp10(a)) ≈ 10^a
        @test @inferred(ldexp(frexp(a)...)) ≈ a
        @test @inferred(log2(exp2(a))) ≈ a
        @test @inferred(log10(exp10(a))) ≈ a
        @test @inferred(log1p(expm1(a))) ≈ a
    end
end

@testset "Logarithm" begin
    @test @inferred(log(x, y)) ≈ measurement(1.261859507142915, 0.059474298734200806)
    @test @inferred(log(y)) ≈ measurement(1.3862943611198906, 0.05)
    @test @inferred(log(y, 4)) ≈ measurement(1, 0.03606737602222409)
    @test @inferred(log(pi, x)) ≈ measurement(0.9597131185693899, 0.029118950894341064)
    @test @inferred(log(big(3), x)) == big(1) ± (x.err / (log(big(3)) * x.val))
    for a in (abs(w), x, y)
        @test @inferred(log(ℯ, a)) ≈ log(a)
        @test @inferred(log(2, a)) ≈ log2(a)
        @test @inferred(log(10, a)) ≈ log10(a)
        @test @inferred(log1p(a)) ≈ log(1 + a)
        @test @inferred(log(exp(-a))) ≈ -a
        @test @inferred(exp(log(a))) ≈ a
    end
end

@testset "Hypotenuse" begin
    for a in (w, x, y)
        @test @inferred(hypot(a, 3)) ≈ @inferred(sqrt(abs2(a) + 9))
        @test @inferred(hypot(4, a)) ≈ @inferred(sqrt(16 + a * a))
        for b in (w, x, y)
            @test @inferred(hypot(a, b)) ≈ @inferred(sqrt(abs2(a) + b * b)) ≈
                @inferred(hypot(b, a)) # Commutativity
        end
    end
end

@testset "Square root" begin
    @test @inferred(sqrt(y)) ≈ measurement(2, 0.05)
    for a in (abs(w), x, y)
        @test @inferred(sqrt(a)) ≈ a ^ 0.5
        @test @inferred(sqrt(a) * sqrt(a)) ≈ a
        # Derivative of sqrt diverges in 0, but if the measurement is exact (like "a-a" is)
        # also the resulting quantity must have 0 uncertainty.
        @test @inferred(sqrt(a - a)) ≈ zero(a)
    end
end

@testset "Cube root" begin
    @test @inferred(cbrt(x)) ≈ measurement(1.4422495703074083, 0.01602499522563787)
    for a in (x, y)
        @test @inferred(cbrt(a)) ≈ a^(1/3)
        @test @inferred(cbrt(a)*cbrt(a)^2) ≈ a
    end
end

@testset "Absolute value" begin
    @test @inferred(abs(-x)) ≈ x
    @test @inferred(abs(measurement(0, 2))) ≈ measurement(0, 2)
    for a in (w, x, y); @test @inferred(abs2(a)) ≈ a*a; end
end

@testset "Sign" begin
    @test @inferred(sign(x)) ≈ measurement(1)
    @test @inferred(sign(-y)) ≈ measurement(-1)
    @test @inferred(sign(measurement(0, 5))) ≈ measurement(0)
    @test @inferred(copysign(x, -5)) == -x
    @test @inferred(copysign(x, w)) == -x
    @test @inferred(copysign(5, w)) == -5
    @test @inferred(copysign(1//2, w)) == -1//2
    @test @inferred(copysign(w, w)) == w
    @test @inferred(copysign(Float32(5), w)) === -Float32(5)
    @test @inferred(copysign(Float64(5), w)) === -Float64(5)
    @test copysign(pi, w) == -pi # This isn't inferrable in Julia
    @test @inferred(flipsign(x, -5)) == -x
    @test @inferred(flipsign(w, w)) == -w
    @test @inferred(flipsign(-5, w)) == 5
    @test @inferred(flipsign(Float32(-5), w)) === Float32(5)
    @test @inferred(flipsign(Float64(-5), w)) === Float64(5)
    @test flipsign(pi, w) == -pi # This isn't inferrable in Julia
end

@testset "One and zero" begin
    @test @inferred(one(y)) ≈ measurement(1)
    @test @inferred(zero(x)) ≈ measurement(0)
end

@testset "Error function" begin
    @test @inferred(erf(x)) ≈ measurement(0.9999779095030014, 1.3925305194674787e-5)
    @test @inferred(erfinv(erf(w))) ≈ w
    @test @inferred(erfc(w)) ≈ 1 - erf(w)
    @test @inferred(erfi(x)) ≈ 1629.9946226015657 ± 914.3351093102547
    for a in (w, x, y)
        @test @inferred(erfcinv(erfc(a))) ≈ a
        @test @inferred(erfcx(a)) ≈ exp(a^2)*erfc(a)
        @test @inferred(dawson(a)) ≈ 0.5*exp(-a^2)*erfi(a)*sqrt(pi)
    end
end

@testset "Cis" begin
    for a in (w, x, y); @test @inferred(cis(a) - exp(im * a)) ≈ 0; end
end

@testset "Factorial and gamma" begin
    @test @inferred(SpecialFunctions.factorial(x)) ≈ measurement(6, 0.7536706010590813)
    @test @inferred(digamma(y)) ≈ 1.256117668431802 ± 0.056764591147422994
    @test @inferred(polygamma(3, w)) ≈ 193.40909103400242 ± 0.10422749480000776
    for a in (w, x, y)
        @test @inferred(gamma(a)) ≈ SpecialFunctions.factorial(a - one(a))
        @test @inferred(gamma(a + one(a))) ≈ SpecialFunctions.factorial(a)
        @test @inferred(lgamma(abs(a))) ≈ log(gamma(abs(a)))
        @test @inferred(digamma(a)) ≈ polygamma(0, a)
        @test @inferred(digamma(invdigamma(a))) ≈ a + zero(a)
        @test @inferred(trigamma(a)) ≈ polygamma(1, a)
    end
end

@testset "Beta" begin
    for a in (w, x, y)
        @test beta(a, x) ≈ gamma(a)*gamma(x)/gamma(a + x)
        @test beta(a, pi) ≈ gamma(a)*gamma(pi)/gamma(a + pi)
        @test beta(ℯ, a) ≈ gamma(ℯ)*gamma(a)/gamma(ℯ + a)
        @test lbeta(abs(a), x) ≈ log(beta(abs(a), x))
        @test lbeta(abs(a), pi) ≈ log(beta(abs(a), pi))
        @test lbeta(ℯ, abs(a)) ≈ log(beta(ℯ, abs(a)))
    end
end

@testset "Airy" begin
    @test airyai(x) ≈ 0.006591139357460721 ± 0.0011912976705951322
    @test airyaiprime(x) ≈ -0.01191297670595132 ± 0.0019773418072382165
    @test airybi(x) ≈ 14.037328963730136 ± 2.2922214966382017
    @test airybiprime(x) ≈ 22.922214966382015 ± 4.211198689119041
end

@testset "Bessel" begin
    for a in (x, y); @test 2*besselj1(a)/a ≈ besselj0(a) + besselj(2, a); end
    for a in (x, y); @test besselj(-1/2, a) ≈ sqrt(2/(pi*a))*cos(a); end
    for a in (x, y); @test besselj(1/2, a) ≈ sqrt(2/(pi*a))*sin(a); end
    @test bessely0(x) ≈ 0.3768500100127904 ± 0.03246744247917999
    for a in (x, y); @test bessely1(a) ≈ -bessely(-1, a); end
    for a in (x, y); @test bessely(5/2, a) ≈
        (besselj(5/2, a)*cos(2.5pi) - besselj(-5/2, a))/sin(2.5pi); end
    for a in (x, y), k in (1, 2), nu in -1:1
        sgn = k == 1 ? +1 : -1
        @test besselh(nu, k, a) ≈
            besselj(nu, a) + sgn*im*bessely(nu, a)
    end
    @test besseli(5//2, y) ≈ 4.757626874823528 ± 1.0398232869843944
    for a in (x, y); @test besselix(3, a) ≈ besseli(3, a)*exp(-abs(a)); end
    @test besselk(7//3, x) ≈ 0.07521953258226349 ± 0.010340691203742959
    for a in (x, y); @test besselkx(3//4, a) ≈ besselk(3//4, a)*exp(a); end
end

@testset "Modulo" begin
    frac, int = @inferred(modf(x + w))
    @test frac ≈ x + w - int
    @test int ≈ @inferred(floor((x + w).val))
    for a in (w, x/2, y)
        @test @inferred(mod(a, -2.1)) ≈ @inferred(a + fld(a, -2.1) * 2.1)
        @test @inferred(mod(-5.8, a)) ≈ @inferred(-5.8 - fld(-5.8, a) * a)
        @test @inferred(rem(a, -3.7)) ≈ @inferred(a + div(a, -3.7) * 3.7)
        @test @inferred(rem(-4.9, a)) ≈ @inferred(-4.9 - div(-4.9, a) * a)
        for b in (10w, x, y)
            @test @inferred(mod(a, b)) ≈ a - fld(a, b)*b
            @test @inferred(rem(a, b)) ≈ a - div(a, b)*b
        end
    end
    @test @inferred(mod2pi(pi*x)) ≈ measurement(pi, 0.1*pi)
end

@testset "Machine precisionx" begin
    @test eps(Measurement{Float64}) ≈ eps(Float64)
    @test eps(x) ≈ eps(x.val)
    @test nextfloat(x) ≈ nextfloat(x.val)
    @test nextfloat(x, 3) ≈ nextfloat(x.val, 3)
    @test maxintfloat(Measurement{Float64}) ≈ maxintfloat(Float64)
end

@testset "Rounding" begin
    @test @inferred(round(w)) ≈ measurement(round(w.val), round(w.err))
    @test @inferred(round(w, digits=1)) ≈ measurement(round(w.val, digits=1),
                                                      round(w.err, digits=1))
    @test @inferred(round(Int, w)) ≈ round(Int, w.val)
    @test @inferred(floor(w)) ≈ measurement(floor(w.val))
    @test @inferred(floor(Int, w)) ≈ floor(Int, w.val)
    @test @inferred(ceil(w)) ≈ measurement(ceil(w.val))
    @test @inferred(ceil(Int, w)) ≈ ceil(Int, w.val)
    @test @inferred(trunc(w)) ≈ measurement(trunc(w.val))
    @test @inferred(trunc(Int, w)) ≈ trunc(Int, w.val)
end

@testset "widening" begin
    @test @inferred(widen(Measurement{Float32})) == Measurement{Float64}
    @test @inferred(widemul(x, y)) ≈ big(x)*big(y)
end

@testset "Type representation" begin
    @test repr("text/plain", [w 10x; 100y 1000w]) ==
        "2×2 Array{Measurement{Float64},2}:\n  -0.5±0.03    30.0±1.0 \n 400.0±20.0  -500.0±30.0"
    @test repr("text/plain", complex(x, w)) == "(3.0 ± 0.1) - (0.5 ± 0.03)im"
    @test repr("text/plain", complex(w, y)) == "(-0.5 ± 0.03) + (4.0 ± 0.2)im"
    @test repr("text/latex", x) == "\$3.0 \\pm 0.1\$"
    @test repr("text/x-tex", y) == repr("text/x-latex", y) == "4.0 \\pm 0.2"
    @test Base.alignment(devnull, x) == (5,4)
    # Make sure the printed representation of a Measurement object is correctly parsed as
    # the same number (note that the tag will be different, but that's not important here).
    for a in (w, x, y); @test eval(Meta.parse(repr(a))) == a; end
end

@testset "sum" begin
    @test @inferred(sum((x, w, x, y, -w))) ≈ 2x + y
    @test @inferred(sum([x, w, x, y, -w])) ≈ 2x + y
end

@testset "prod" begin
    @test @inferred(prod((w, x, y, -w))) ≈ -w^2 * x * y
    @test @inferred(prod([w, x, y, -w])) ≈ -w^2 * x * y
    @test @inferred(prod((w, x, y, -w, 0))) ≈ 0 ± 0
    @test @inferred(prod([w, x, y, -w, 0])) ≈ 0 ± 0
end

@testset "Miscellanea" begin
    @test mean((x, w, x, y, -w)) ≈ (2x + y)/5
    @test mean([x, w, x, y, -w]) ≈ (2x + y)/5
    @test min(w, x, y) === w
    @test max(w, x, y) === y
    @test extrema([w, x, y]) == (w, y)
    @test sort([y, w, x]) == [w, x, y]
end

@testset "Linear algebra" begin
    A = [(14 ± 0.1) (23 ± 0.2); (-12 ± 0.3) (24 ± 0.4)]
    c = [(7 ± 0.5), (-13 ± 0.6)]
    b = @inferred(A \ c)
    @test @inferred(A * b) ≈ c
    @test @inferred(b ⋅ c) ≈ 7.423202614379084 ± 0.5981875954418516
    @test @inferred(det(A)) ≈ 612 ± 9.51262319236918
    @test @inferred(A * inv(A)) ≈ I
end

@testset "NaNs" begin
    # NaN as nominal value
    @test @inferred(2*(NaN ± 3)) ≈ NaN ± 6
    # NaN as uncertainty
    @test @inferred(2*(3 ± NaN)) ≈ 6 ± NaN
    # Both
    @test @inferred(2*(NaN ± NaN)) ≈ NaN ± NaN
    # Make sure isapprox is working well here
    @test 2*(3 ± NaN) ≉ 6 ± 0
    @test 2*(NaN ± NaN) ≉ NaN ± 0
end

@testset "@uncertain" begin
    @test @uncertain(tan(x)) ≈ tan(x)
    @testset "non-Measurement arguments" begin
        @test @uncertain(sin(-7.8)) == measurement(sin(-7.8))
        @test @uncertain(log(5.4, 3.6)) == measurement(log(5.4, 3.6))
        @test @uncertain(atan(x, 1.2)) ≈ atan(x, 1.2)
        @test @uncertain(hypot(-2.9, y)) ≈ hypot(-2.9, y)
    end
    @testset "correlation" begin
        @test @uncertain((a -> a + a + a)(x)) ≈ 3x
        @test @uncertain(+(x, -x, y, log(y), -w, w^2)) ≈ + y + log(y) - w * (1 - w)
    end
    @testset "zeta function" begin
        @test @uncertain(zeta(x)) ≈
            measurement(1.2020569031595951, 0.019812624290876782)
        # Test a couple of identities of the zeta function:
        # http://functions.wolfram.com/ZetaFunctionsandPolylogarithms/Zeta/17/01/01/
        @test @uncertain(zeta(w)) ≈
            gamma(1 - w) * 2 ^ w * pi ^ (w - 1) * sinpi(w / 2) * @uncertain(zeta(1 - w))
        @test @uncertain(zeta(w)) ≈
            pi ^ (w - 1/2) * gamma((1 - w) / 2) / gamma(w / 2) * @uncertain(zeta(1 - w))
        # Test a few derivatives of the zeta functions which are known:
        # http://mathworld.wolfram.com/RiemannZetaFunction.html
        @test @uncertain(zeta(0 ± 1)).der.value ≈ -log(2pi) / 2
        @test @uncertain(zeta(0.5 ± 1)).der.value ≈ -3.92264613
        @test @uncertain(zeta(2 ± 1)).der.value ≈ -0.93754825431
    end
    @testset "2-arg functions: $f" for f in (log, hypot, atan)
        @test @uncertain(f(x, y)) ≈ f(x, y)
    end
    ## XXX: re-enable the following test when this bug is fixed:
    ## https://github.com/JuliaLang/julia/issues/23878
    # @testset "ccall" begin
    #     f(x) = x*x
    #     ptr = @cfunction(f, Cdouble, (Cdouble,))
    #     g(x) = ccall(ptr, Cdouble, (Cdouble,), x)*x
    #     @test @uncertain(g(x)) ≈ x^3
    # end
end

@testset "value and uncertainty" begin
    @test Measurements.value.([w, x, y]) == [-0.5, 3.0, 4.0]
    @test Measurements.value.([complex(w, x)]) == [complex(-0.5, 3.0)]
    @test Measurements.uncertainty.([w, x, y]) == [0.03, 0.1, 0.2]
    @test Measurements.uncertainty.([complex(w, x)]) == [complex(0.03, 0.1)]
    @test Measurements.uncertainty(5) === 0
    @test Measurements.uncertainty(7.4) === 0.0
    @test Measurements.uncertainty(-3.6 + 7.4im) === complex(0.0)
    @test Measurements.value(5) === 5
    @test Measurements.value(7.4) === 7.4
    @test Measurements.value(-3.6 + 7.4im) === complex(-3.6 + 7.4im)
end

@testset "Derivatives type" begin
    @test_throws KeyError getindex(x.der, (0, 0, -1))
    @test length((w + w + 2x + y).der) == 3
end

@testset "BigFloat's" begin
    @test big(Measurement) == Measurement{BigFloat}
    @test typeof(big(x)) == big(typeof(x))
    for a in (w, x, y); @test big(x) ^ 2 - x * x ≈ zero(x); end
    @test Measurements.derivative(big(x) / 3 + x, x) ≈ inv(big(3)) + 1
    @test Measurements.derivative(big(x) + y, y) ≈ 1
    @test typeof(big(z)) == big(typeof(z))
    a = big"3.00000001" ± big"1e-17"
    b = big"4.00000001" ± big"1e-17"
    c = big"5.00000001" ± big"1e-15"
    d = big"4.0000000100000001" ± big"1e-17"
    for X in (a, b, c, d)
        @test @inferred(sin(X*X + X*X)/cos(X*X + X*X)) ≈ @inferred(tan(2X^2))
    end
    @test @inferred(a * b) ≈
        big"12.0000000700000001" ±
        big"5.000000014000000000399999998880000003119999991353600023834879934652928178154734e-17"
    @test @inferred(a / x + w / b) ≈
        big"8.750000036458333325520833352864583284505208455403645528157552846272784550984708e-01" ±
        big"3.416666677095189699499391052146002403356891669706809079678073000173723744722607e-02"
    # The following test is present also in the examples in the manual.
    @test stdscore(a * d, big"12.00000007") > 7.9
end

@testset "QuadGK" begin
    @test QuadGK.quadgk(cos, x, y)[1] ≈ sin(y) - sin(x)
    @test QuadGK.quadgk(sin, -y, y)[1] ≈ cos(-y) - cos(y) atol = eps(Float64)
    @test QuadGK.quadgk(exp, 0.4, x)[1] ≈ exp(x) - exp(0.4)
    @test QuadGK.quadgk(sin, w, 2.7)[1] ≈ cos(w) - cos(2.7)
    @test QuadGK.quadgk(t -> cos(x - t), 0, 2pi)[1] ≈ measurement(0) atol = 7e-16
    @test QuadGK.quadgk(t -> exp(t / w), 0, 1)[1] ≈ w * (exp(1 / w) - one(x))
    for a in (w, x, y)
        @test QuadGK.quadgk(t -> 1 / abs2(t / a), 1, Inf)[1] ≈ a ^ 2
        @test QuadGK.quadgk(t -> exp(-abs2(t * a)), -Inf, Inf)[1] ≈ sqrt(pi) / abs(a)
    end
    @test QuadGK.quadgk(t -> cos(t - w), -w, 0)[1] + QuadGK.quadgk(t -> cos(t - w), 0, w)[1] ≈ sin(2w)
    @test QuadGK.quadgk(t -> exp(t / x), w, y)[1] ≈ x * (exp(y / x) - exp(w / x))
    @test QuadGK.quadgk(t -> sin(t - w), 0, y)[1] ≈ cos(w) - cos(y - w)
    @test QuadGK.quadgk(t -> log(y - t), w, -pi)[1] ≈
        (y - w)*log(y - w) - (y + pi)*log(y + pi) + w + pi
    # Compare some of the above integrals with results with "@uncertain" macro.
    @test QuadGK.quadgk(cos, x, y)[1] ≈
        @uncertain(((x,y) -> QuadGK.quadgk(cos, x, y)[1])(x, y))
    @test QuadGK.quadgk(sin, -y, y)[1] ≈
        @uncertain((x -> QuadGK.quadgk(sin, -x, x)[1])(y)) atol = 1e-10
    @test QuadGK.quadgk(exp, 0.4, x)[1] ≈
        @uncertain((x -> QuadGK.quadgk(exp, 0.4, x)[1])(x))
end

@testset "Calculus" begin
    for a in (w, x, y)
        @test Calculus.derivative(exp, a + w) ≈ exp(a + w)
        @test Calculus.derivative(t -> sin(t * x), a) ≈  x * cos(a * x)
        @test Calculus.derivative(t -> y * cos(t), a) ≈ -y * sin(a)
    end
end

@testset "String parsing" begin
    @test @inferred(measurement("  -123.4(56)  ")) ≈         -123.4 ± 5.6
    @test @inferred(measurement("  +1234(56)e-1  ")) ≈        123.4 ± 5.6
    @test @inferred(measurement(" ( +1.234 ± 0.056 )e2 ")) ≈  123.4 ± 5.6
    @test @inferred(measurement(" -12.34e1  +-  0.56e1 ")) ≈ -123.4 ± 5.6
    @test @inferred(measurement(" 1234e-1  +/-  5.6e0 ")) ≈   123.4 ± 5.6
    @test @inferred(measurement("  -1234e-1  ")) ≈ measurement(-1234e-1)
    @test @inferred(parse(Measurement{BigFloat}, "5.1 ± 3.3")) == big"5.1" ± big"3.3"
    @test @inferred(parse(Measurement{Float32}, "-7.6 ± 0.4")) == Float32(-7.6) ± Float32(0.4)
    @test @inferred(parse(Measurement{Float32}, "(3±0.1)e2")) == (Float32(3) ± Float32(0.1)) * Float32(100.0)
    @test @inferred(parse(Measurement{BigFloat}, "(3±0.1)e2")) == (big"3" ± big"0.1") * big"100.0"
    @test @inferred(parse(Measurement{Float16}, "+9.5 ± 2.8")) == Float16(+9.5) ± Float16(2.8)
    for a in (w, x, y); @test @inferred(parse(Measurement{Float64}, repr(a))) == a; end
    @test_throws ArgumentError measurement("abc")
    @test_throws ArgumentError measurement("(2±1")
    @test_throws ArgumentError measurement("2±1)e2")
    @test_throws ArgumentError measurement("(2)")
    @test_throws ArgumentError measurement("(2)e-2")
end

@testset "Fixed bugs" begin
    @test 1 * (big(3) ± 0) ≈ 3 ± 0
    @test Measurement{Float64}(2//1) === Measurement{Float64}(2.0)
end

@testset "Complex measurements" begin
    include("complex.jl")
end

@testset "Plot recipes" begin
    include("plots.jl")
end


text:
# ----------------------- CONSTANTS --------------------------------- #
const EARTHRADIUS = 6.37122e6  # meters
const WATERDENSITY = 1000.

# area
const METERSPERKM = 1000.
const METERSPERMILE = 1609.34
const METERS2PERACRE = 4046.856

# time
const HOURSPERDAY = 24.
const SECSPERHOUR = 3600.
const MINSPERHOUR = 60.
const MINSPERDAY = HOURSPERDAY * MINSPERHOUR
const SECSPERDAY = HOURSPERDAY * SECSPERHOUR
const MONTHSPERYEAR = 12

# length
const MMPERMETER = 1000.
const CMPERMETER = 100.


text:
import Dates
struct Clock
    hours::Int8
    minutes::Int8
    Clock(h,m) = new(rem(h + div(m,60,RoundDown),24,RoundDown),rem(m,60,RoundDown))
end

Base.:+(c::Clock,m::Dates.Minute) = Clock(c.hours, c.minutes + m.value)
Base.:-(c::Clock,m::Dates.Minute) = Clock(c.hours, c.minutes - m.value)
Base.:+(c::Clock,h::Dates.Hour) = Clock(c.hours + h.value, c.minutes)
Base.:-(c::Clock,h::Dates.Hour) = Clock(c.hours - h.value, c.minutes)

Base.show(io::IO, c::Clock) = print(io, '"', lpad(c.hours,2,'0'), ":", lpad(c.minutes,2,'0'), '"')


text:
export AbstractPreconditioner,
    ColumnwiseLUPreconditioner,
    NoPreconditioner,
    preconditioner_update!,
    preconditioner_solve!,
    preconditioner_counter_update!

"""
    Abstract base type for all preconditioners.
"""
abstract type AbstractPreconditioner end

"""
mutable struct NoPreconditioner
end

Do nothing
"""
mutable struct NoPreconditioner <: AbstractPreconditioner end

"""
Do nothing, when there is no preconditioner, preconditioner = Nothing
"""
function preconditioner_update!(
    op,
    dg,
    preconditioner::NoPreconditioner,
    args...,
) end

"""
Do nothing, when there is no preconditioner, preconditioner = Nothing
"""
function preconditioner_solve!(preconditioner::NoPreconditioner, Q) end

"""
Do nothing, when there is no preconditioner, preconditioner = Nothing
"""
function preconditioner_counter_update!(preconditioner::NoPreconditioner) end

"""
mutable struct ColumnwiseLUPreconditioner{AT}
    A::DGColumnBandedMatrix
    Q::AT
    PQ::AT
    counter::Int
    update_freq::Int
end

...
# Arguments
- `A`: the lu factor of the precondition (approximated Jacobian), in the DGColumnBandedMatrix format
- `Q`: MPIArray container, used to update A
- `PQ`: MPIArray container, used to update A
- `counter`: count the number of Newton, when counter > update_freq or counter < 0, update precondition
- `update_freq`: preconditioner update frequency
...
"""
mutable struct ColumnwiseLUPreconditioner{AT} <: AbstractPreconditioner
    A::DGColumnBandedMatrix
    Q::AT
    PQ::AT
    counter::Int
    update_freq::Int
end

"""
ColumnwiseLUPreconditioner constructor
build an empty ColumnwiseLUPreconditioner

...
# Arguments
- `dg`: DG model, use only the grid information
- `Q0`: MPIArray, use only its size
- `counter`: = -1, which indicates the preconditioner is empty
- `update_freq`: preconditioner update frequency
...
"""
function ColumnwiseLUPreconditioner(dg, Q0, update_freq = 100)
    single_column = false
    Q = similar(Q0)
    PQ = similar(Q0)

    A = empty_banded_matrix(dg, Q; single_column = single_column)

    # counter = -1, which indicates the preconditioner is empty
    ColumnwiseLUPreconditioner(A, Q, PQ, -1, update_freq)
end

"""
update the DGColumnBandedMatrix by the finite difference approximation
...
# Arguments
- `op`: operator used to compute the finte difference information
- `dg`: the DG model, use only the grid information
...
"""
function preconditioner_update!(
    op,
    dg,
    preconditioner::ColumnwiseLUPreconditioner,
    args...,
)

    # preconditioner.counter < 0, means newly constructed empty preconditioner
    if preconditioner.counter >= 0 &&
       (preconditioner.counter < preconditioner.update_freq)
        return
    end

    A = preconditioner.A
    Q = preconditioner.Q
    PQ = preconditioner.PQ

    update_banded_matrix!(A, op, dg, Q, PQ, args...)
    band_lu!(A)

    preconditioner.counter = 0
end

"""
Inplace applying the preconditioner

Q = P⁻¹ * Q
"""
function preconditioner_solve!(preconditioner::ColumnwiseLUPreconditioner, Q)
    A = preconditioner.A
    band_forward!(Q, A)
    band_back!(Q, A)

end

"""
Update the preconditioner counter, after each Newton iteration
"""
function preconditioner_counter_update!(
    preconditioner::ColumnwiseLUPreconditioner,
)
    preconditioner.counter += 1
end


text:
## A SortedDict is a wrapper around balancedTree with
## methods similiar to those of Julia container Dict.

type SortedDict{K, D, Ord <: Ordering} <: Associative{K,D}
    bt::BalancedTree23{K,D,Ord}

    ## Base constructors

    SortedDict(o::Ord) = new(BalancedTree23{K,D,Ord}(o))

    function SortedDict(o::Ord, kv)
        s = new(BalancedTree23{K,D,Ord}(o))

        if eltype(kv) <: Pair
            # It's (possibly?) more efficient to access the first and second
            # elements of Pairs directly, rather than destructure
            for p in kv
                s[p.first] = p.second
            end
        else
            for (k,v) in kv
                s[k] = v
            end
        end
        return s
    end

end

# Any-Any constructors
SortedDict() = SortedDict{Any,Any,ForwardOrdering}(Forward)
SortedDict{Ord <: Ordering}(o::Ord) = SortedDict{Any,Any,Ord}(o)

function not_iterator_of_pairs(kv)
    return any(x->isempty(methodswith(typeof(kv), x, true)),
               [start, next, done]) ||
           any(x->!isa(x, Union{Tuple,Pair}), kv)
end

# Construction from Pairs
# TODO: fix SortedDict(1=>1, 2=>2.0)
SortedDict(ps::Pair...) = SortedDict(Forward, ps)
SortedDict(o::Ordering, ps::Pair...) = SortedDict(o, ps)
@compat (::Type{SortedDict{K,D}}){K,D}(ps::Pair...) = SortedDict{K,D,ForwardOrdering}(Forward, ps)
@compat (::Type{SortedDict{K,D}}){K,D,Ord<:Ordering}(o::Ord, ps::Pair...) = SortedDict{K,D,Ord}(o, ps)

# Construction from Associatives
SortedDict{K,D,Ord<:Ordering}(o::Ord, d::Associative{K,D}) = SortedDict{K,D,Ord}(o, d)

## Construction from iteratables of Pairs/Tuples

# Construction specifying Key/Value types
# e.g., SortedDict{Int,Float64}([1=>1, 2=>2.0])
@compat (::Type{SortedDict{K,D}}){K,D}(kv) = SortedDict{K,D}(Forward, kv)
@compat function (::Type{SortedDict{K,D}}){K,D,Ord<:Ordering}(o::Ord, kv)
    try
        SortedDict{K,D,Ord}(o, kv)
    catch e
        if not_iterator_of_pairs(kv)
            throw(ArgumentError("SortedDict(kv): kv needs to be an iterator of tuples or pairs"))
        else
            rethrow(e)
        end
    end
end

# Construction inferring Key/Value types from input
# e.g. SortedDict{}
SortedDict(kv, o::Ordering=Forward) = SortedDict(o, kv)
function SortedDict(o::Ordering, kv)
    try
        _sorted_dict_with_eltype(o, kv, eltype(kv))
    catch e
        if not_iterator_of_pairs(kv)
            throw(ArgumentError("SortedDict(kv): kv needs to be an iterator of tuples or pairs"))
        else
            rethrow(e)
        end
    end
end

_sorted_dict_with_eltype{K,D,Ord}(o::Ord, ps, ::Type{Pair{K,D}}) = SortedDict{  K,  D,Ord}(o, ps)
_sorted_dict_with_eltype{K,D,Ord}(o::Ord, kv, ::Type{Tuple{K,D}}) = SortedDict{  K,  D,Ord}(o, kv)
_sorted_dict_with_eltype{K,  Ord}(o::Ord, ps, ::Type{Pair{K}}  ) = SortedDict{  K,Any,Ord}(o, ps)
_sorted_dict_with_eltype{    Ord}(o::Ord, kv, ::Type            ) = SortedDict{Any,Any,Ord}(o, kv)

## TODO: It seems impossible (or at least very challenging) to create the eltype below.
##       If deemed possible, please create a test and uncomment this definition.
# if VERSION < v"0.6.0-dev.2123"
#     _sorted_dict_with_eltype{  D,Ord}(o::Ord, ps, ::Type{Pair{TypeVar(:K),D}}) = SortedDict{Any,  D,Ord}(o, ps)
# else
#     eval(parse("_sorted_dict_with_eltype{  D,Ord}(o::Ord, ps, ::Type{Pair{K,D} where K}) = SortedDict{Any,  D,Ord}(o, ps)"))
# end

typealias SDSemiToken IntSemiToken

typealias SDToken Tuple{SortedDict,IntSemiToken}

## This function implements m[k]; it returns the
## data item associated with key k.

@inline function getindex(m::SortedDict, k_)
    i, exactfound = findkey(m.bt, convert(keytype(m),k_))
    !exactfound && throw(KeyError(k_))
    return m.bt.data[i].d
end

## This function implements m[k]=d; it sets the
## data item associated with key k equal to d.

@inline function setindex!{K, D, Ord <: Ordering}(m::SortedDict{K,D,Ord}, d_, k_)
    insert!(m.bt, convert(K,k_), convert(D,d_), false)
    m
end

## push! is an alternative to insert!; it returns the container.

@inline function push!{K,D}(m::SortedDict{K,D}, pr::Pair)
    insert!(m.bt, convert(K, pr[1]), convert(D, pr[2]), false)
    m
end

## This function looks up a key in the tree;
## if not found, then it returns a marker for the
## end of the tree.

@inline function find(m::SortedDict, k_)
    ll, exactfound = findkey(m.bt, convert(keytype(m),k_))
    IntSemiToken(exactfound? ll : 2)
end

## This function inserts an item into the tree.
## Unlike m[k]=d, it also returns a bool and a token.
## The bool is true if the inserted item is new.
## It is false if there was already an item
## with that key.
## The token points to the newly inserted item.

@inline function insert!{K,D, Ord <: Ordering}(m::SortedDict{K,D,Ord}, k_, d_)
    b, i = insert!(m.bt, convert(K,k_), convert(D,d_), false)
    b, IntSemiToken(i)
end

@inline eltype{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord}) =  Pair{K,D}
@inline eltype{K,D,Ord <: Ordering}(::Type{SortedDict{K,D,Ord}}) =  Pair{K,D}
@inline function in{K,D,Ord <: Ordering}(pr::Pair, m::SortedDict{K,D,Ord})
    i, exactfound = findkey(m.bt,convert(K,pr[1]))
    return exactfound && isequal(m.bt.data[i].d,convert(D,pr[2]))
end

@inline in(::Tuple{Any,Any}, ::SortedDict) =
    throw(ArgumentError("'(k,v) in sorteddict' not supported in Julia 0.4 or 0.5.  See documentation"))

@inline keytype{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord}) = K
@inline keytype{K,D,Ord <: Ordering}(::Type{SortedDict{K,D,Ord}}) = K
@inline valtype{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord}) = D
@inline valtype{K,D,Ord <: Ordering}(::Type{SortedDict{K,D,Ord}}) = D
@inline ordtype{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord}) = Ord
@inline ordtype{K,D,Ord <: Ordering}(::Type{SortedDict{K,D,Ord}}) = Ord

## First and last return the first and last (key,data) pairs
## in the SortedDict.  It is an error to invoke them on an
## empty SortedDict.

@inline function first(m::SortedDict)
    i = beginloc(m.bt)
    i == 2 && throw(BoundsError())
    return Pair(m.bt.data[i].k, m.bt.data[i].d)
end

@inline function last(m::SortedDict)
    i = endloc(m.bt)
    i == 1 && throw(BoundsError())
    return Pair(m.bt.data[i].k, m.bt.data[i].d)
end

@inline orderobject(m::SortedDict) = m.bt.ord

@inline function haskey(m::SortedDict, k_)
    i, exactfound = findkey(m.bt,convert(keytype(m),k_))
    exactfound
end

@inline function get{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord}, k_, default_)
    i, exactfound = findkey(m.bt, convert(K,k_))
   return exactfound? m.bt.data[i].d : convert(D,default_)
end

function get!{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord}, k_, default_)
    k = convert(K,k_)
    i, exactfound = findkey(m.bt, k)
    if exactfound
        return m.bt.data[i].d
    else
        default = convert(D,default_)
        insert!(m.bt,k, default, false)
        return default
    end
end

function getkey{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord}, k_, default_)
    i, exactfound = findkey(m.bt, convert(K,k_))
    exactfound? m.bt.data[i].k : convert(K, default_)
end

## Function delete! deletes an item at a given
## key

@inline function delete!(m::SortedDict, k_)
    i, exactfound = findkey(m.bt,convert(keytype(m),k_))
    !exactfound && throw(KeyError(k_))
    delete!(m.bt, i)
    m
end

@inline function pop!(m::SortedDict, k_)
    i, exactfound = findkey(m.bt,convert(keytype(m),k_))
    !exactfound && throw(KeyError(k_))
    d = m.bt.data[i].d
    delete!(m.bt, i)
    d
end

## Check if two SortedDicts are equal in the sense of containing
## the same (K,D) pairs.  This sense of equality does not mean
## that semitokens valid for one are also valid for the other.

function isequal(m1::SortedDict, m2::SortedDict)
    ord = orderobject(m1)
    if !isequal(ord, orderobject(m2)) || !isequal(eltype(m1), eltype(m2))
        throw(ArgumentError("Cannot use isequal for two SortedDicts unless their element types and ordering objects are equal"))
    end
    p1 = startof(m1)
    p2 = startof(m2)
    while true
        if p1 == pastendsemitoken(m1)
            return p2 == pastendsemitoken(m2)
        end
        if p2 == pastendsemitoken(m2)
            return false
        end
        k1,d1 = deref((m1,p1))
        k2,d2 = deref((m2,p2))
        if !eq(ord,k1,k2) || !isequal(d1,d2)
            return false
        end
        p1 = advance((m1,p1))
        p2 = advance((m2,p2))
    end
end

function mergetwo!{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord},
                                        m2::Associative{K,D})
    for (k,v) in m2
        m[convert(K,k)] = convert(D,v)
    end
end

function packcopy{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord})
    w = SortedDict(Dict{K,D}(),orderobject(m))
    mergetwo!(w,m)
    w
end

function packdeepcopy{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord})
    w = SortedDict(Dict{K,D}(),orderobject(m))
    for (k,v) in m
        newk = deepcopy(k)
        newv = deepcopy(v)
        w[newk] = newv
    end
    w
end

function merge!{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord},
                                     others::Associative{K,D}...)
    for o in others
        mergetwo!(m,o)
    end
end

function merge{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord},
                                    others::Associative{K,D}...)
    result = packcopy(m)
    merge!(result, others...)
    result
end

similar{K,D,Ord<:Ordering}(m::SortedDict{K,D,Ord}) =
    SortedDict{K,D,Ord}(orderobject(m))

isordered{T<:SortedDict}(::Type{T}) = true