text:
using Flux
using Flux: @functor

abstract type AbstractTransformer end

struct PwFFN
    din::Dense
    dout::Dense
end

@functor PwFFN

"just a wrapper for two dense layer."
PwFFN(size::Int, h::Int, act = relu) = PwFFN(
    Dense(size, h, act),
    Dense(h, size)
)

function (pw::PwFFN)(x)::AbstractMatrix
    # size(x) == (dims, seq_len)
    pw.dout(pw.din(x))
end

struct Transformer <: AbstractTransformer
    mh::MultiheadAttention
    mhn::LayerNorm
    pw::PwFFN
    pwn::LayerNorm
    drop::Dropout
end

@functor Transformer

"""
    Transformer(size::Int, head::Int, ps::Int;
                future::Bool = true, act = relu, pdrop = 0.1)
    Transformer(size::Int, head::Int, hs::Int, ps::Int;
                future::Bool = true, act = relu, pdrop = 0.1)

Transformer layer.

`size` is the input size. if `hs` is not specify, use `div(size, head)` as the hidden size of multi-head attention.
`ps` is the hidden size & `act` is the activation function of the positionwise feedforward layer.
When `future` is `false`, the k-th token can't see the j-th tokens where j > k. `pdrop` is the dropout rate.
"""
function Transformer(size::Int, head::Int, ps::Int; future::Bool = true, act = relu, pdrop = 0.1)
    rem(size, head) != 0 && error("size not divisible by head")
    Transformer(size, head, div(size, head), ps;future=future, act=act, pdrop=pdrop)
end

Transformer(size::Int, head::Int, hs::Int, ps::Int; future::Bool = true, act = relu, pdrop = 0.1) = Transformer(
    MultiheadAttention(head, size, hs, size; future=future, pdrop=pdrop),
    LayerNorm(size),
    PwFFN(size, ps, act),
    LayerNorm(size),
    Dropout(pdrop),
)

function (t::Transformer)(x::AbstractArray{T, N}, mask=nothing) where {T, N}
    a = t.mh(x, x, x; mask=mask)
    a = t.drop(a)
    res_a = x .+ a
    if N == 3
        insize = size(res_a)
        res_a = reshape(res_a, insize[1], :)
    end
    res_a = t.mhn(res_a)
    pwffn = t.pw(res_a)
    pwffn = t.drop(pwffn)
    res_pwffn = res_a .+ pwffn
    res_pwffn = t.pwn(res_pwffn)
    if N == 3
        res_pwffn = reshape(res_pwffn, :, Base.tail(insize)...)
    end
    res_pwffn
end

function Base.show(io::IO, t::Transformer)
    hs = div(size(t.mh.iqproj.W)[1], t.mh.head)
    h, ps = size(t.pw.dout.W)

    print(io, "Transformer(")
    print(io, "head=$(t.mh.head), ")
    print(io, "head_size=$(hs), ")
    print(io, "pwffn_size=$(ps), ")
    print(io, "size=$(h)")
    if Flux.istraining()
        print(io, ", dropout=$(t.drop.p))")
    else
        print(io, ")")
    end
end

struct TransformerDecoder <: AbstractTransformer
    mh::MultiheadAttention
    mhn::LayerNorm
    imh::MultiheadAttention
    imhn::LayerNorm
    pw::PwFFN
    pwn::LayerNorm
    drop::Dropout
end

@functor TransformerDecoder

"""
    TransformerDecoder(size::Int, head::Int, ps::Int; act = relu, pdrop = 0.1)
    TransformerDecoder(size::Int, head::Int, hs::Int, ps::Int; act = relu, pdrop = 0.1)

TransformerDecoder layer. Decode the value from a Encoder.

`size` is the input size. if `hs` is not specify, use `div(size, head)` as the hidden size of multi-head attention.
`ps` is the hidden size & `act` is the activation function of the positionwise feedforward layer.
`pdrop` is the dropout rate.
"""
function TransformerDecoder(size::Int, head::Int, ps::Int; act = relu, pdrop = 0.1)
    rem(size, head) != 0 && error("size not divisible by head")
    TransformerDecoder(size, head, div(size, head), ps; act=act, pdrop=pdrop)
end

TransformerDecoder(size::Int, head::Int, hs::Int, ps::Int; act = relu, pdrop = 0.1) = TransformerDecoder(
    MultiheadAttention(head, size, hs, size; future=false, pdrop=pdrop),
    LayerNorm(size),
    MultiheadAttention(head, size, hs, size; future=true, pdrop=pdrop),
    LayerNorm(size),
    PwFFN(size, ps, act),
    LayerNorm(size),
    Dropout(pdrop),
)

function (td::TransformerDecoder)(x::AbstractArray{T,N}, m, mask=nothing) where {T,N}
    a = td.mh(x,x,x)
    a = td.drop(a)
    res_a = x .+ a
    res_a = N == 3 ? @toNd(td.mhn(res_a)) : td.mhn(res_a)

    ia = td.imh(res_a, m, m, mask=mask)
    ia = td.drop(ia)
    res_ia = res_a .+ ia
    if N == 3
        insize = size(res_ia)
        res_ia = reshape(res_ia, insize[1], :)
    end
    res_ia = td.imhn(res_ia)
    pwffn = td.pw(res_ia)
    pwffn = td.drop(pwffn)
    res_pwffn = res_ia .+ pwffn
    res_pwffn = td.pwn(res_pwffn)
    if N == 3
        res_pwffn = reshape(res_pwffn, :, Base.tail(insize)...)
    end

    res_pwffn
end

function Base.show(io::IO, td::TransformerDecoder)
    hs = div(size(td.imh.iqproj.W)[1], td.imh.head)
    h, ps = size(td.pw.dout.W)

    print(io, "TransformerDecoder(")
    print(io, "head=$(td.mh.head), ")
    print(io, "head_size=$(hs), ")
    print(io, "pwffn_size=$(ps), ")
    print(io, "size=$(h)")
    if Flux.istraining()
        print(io, ", dropout=$(td.drop.p))")
    else
        print(io, ")")
    end
end


text:
include(joinpath("Duration", "predByMean.jl"))
include(joinpath("Duration", "linReg.jl"))
include(joinpath("Duration", "gammaReg.jl"))
include(joinpath("Duration", "decisionTree.jl"))
include(joinpath("Duration", "randomForest.jl"))

ext: jl

max_stars_repo_name: AlexLsn/CSO

max_issues_repo_name: AlexLsn/CSO

max_forks_repo_name: AlexLsn/CSO


text:
module ExampleServer
    using Immortal

    include("./Handlers.jl")
    import .Handlers

    function run()
        config_path::String = "/home/entropy/Проекты/ExampleServer/config.yml"

        connection, broker_config, channels_config = Immortal.init(config_path)
        channels_dicts = Immortal.declares(connection, broker_config, channels_config)
        Adapter = Immortal.get_adapter(broker_config)
    end
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Data

A domain together with a table of values for each element of the domain.
The i-th row of the table is a vector of features associated with the
i-th element of the domain. If the domain is a mesh, then the table
stores all the properties of the elements.
"""
abstract type Data end

"""
    domain(data)

Return underlying domain of the `data`.
"""
function domain end

"""
    values(data)

Return the values of `data` as a table.
"""
values(data::Data)

"""
    constructor(D)

Return the constructor of the data type `D` as a function.
The function takes a `domain` and `table` as inputs and
combines them into an instance of the data type.
"""
function constructor end

# ----------
# FALLBACKS
# ----------

function (D::Type{<:Data})(stable)
  # build domain from geometry column
  ctable = Tables.columns(stable)
  elms   = Tables.getcolumn(ctable, :geometry)
  domain = Collection(elms)

  # build table from remaining columns
  vars = setdiff(Tables.columnnames(ctable), (:geometry,))
  cols = map(vars) do var
    var => Tables.getcolumn(ctable, var)
  end
  table = (; cols...)

  # combine the two with constructor
  constructor(D)(domain, table)
end

==(data₁::Data, data₂::Data) =
  domain(data₁) == domain(data₂) && values(data₁) == values(data₂)

# implement Domain traits for convenience
embeddim(data::Data) = embeddim(domain(data))
coordtype(data::Data) = coordtype(domain(data))
nelements(data::Data) = nelements(domain(data))
centroid(data::Data, ind) = centroid(domain(data), ind)

# -----------------
# TABLES INTERFACE
# -----------------

Tables.istable(::Type{<:Data}) = true

Tables.rowaccess(::Type{<:Data}) = true

Tables.rows(data::Data) = DataRows(domain(data), Tables.rows(values(data)))

# wrapper type for rows of the data table
# so that we can easily inform the schema
struct DataRows{𝒟,𝒯}
  domain::𝒟
  rtable::𝒯
end

function Base.iterate(rows::DataRows, state=1)
  if state > nelements(rows.domain)
    nothing
  else
    row, _ = iterate(rows.rtable, state)
    elm, _ = iterate(rows.domain, state)
    (; NamedTuple(row)..., geometry=elm), state + 1
  end
end

Base.length(rows::DataRows) = nelements(rows.domain)

function Tables.schema(rows::DataRows)
  geomtype = eltype(rows.domain)
  schema = Tables.schema(rows.rtable)
  names, types = schema.names, schema.types
  Tables.Schema((names..., :geometry), (types..., geomtype))
end

# data table is compatible with the Queryverse
TableTraits.isiterabletable(data::Data) = true
IIE.getiterator(data::Data) = Tables.datavaluerows(Tables.rows(data))
IIE.isiterable(data::Data) = true

Tables.materializer(D::Type{<:Data}) = D

# -------------------
# VARIABLE INTERFACE
# -------------------

"""
    getindex(data, var)

Returns the data for the variable `var` in `data` as a column vector.
"""
Base.getindex(data::Data, var::Symbol) =
  Tables.getcolumn(values(data), var)

Base.getindex(data::Data, var::String) =
  getindex(data, Symbol(var))

function variables(data::Data)
  s = Tables.schema(values(data))
  @. Variable(s.names, nonmissingtype(s.types))
end

# ----------
# UTILITIES
# ----------

"""
    asarray(data, var)

Returns the data for the variable `var` in `data` as a Julia array
with size equal to the size of the underlying domain if the size is
defined, otherwise returns a vector.
"""
function asarray(data::Data, var::Symbol)
  D = domain(data)
  hassize = hasmethod(size, (typeof(D),))
  hassize ? reshape(data[var], size(D)) : data[var]
end

asarray(data::Data, var::String) =
  asarray(data, Symbol(var))

# -----------
# IO METHODS
# -----------

function Base.show(io::IO, data::Data)
  name = nameof(typeof(data))
  nelm = nelements(domain(data))
  Dim  = embeddim(domain(data))
  T    = coordtype(domain(data))
  print(io, "$nelm $name{$Dim,$T}")
end

function Base.show(io::IO, ::MIME"text/plain", data::Data)
  𝒟 = domain(data)
  𝒯 = values(data)
  s = Tables.schema(𝒯)
  vars = zip(s.names, s.types)
  println(io, data)
  println(io, "  variables")
  varlines = ["    └─$var ($V)" for (var,V) in vars]
  println(io, join(sort(varlines), "\n"))
  print(  io, "  domain: ", 𝒟)
end


text:
#!/usr/bin/env julia

using Eigenbroetler

function high_freq_sinusoid(x::Number, y::Number)
    # sin(100 * r)
    r = sqrt(x * x + y * y)
    return sin(100 * r)
end

eb = Eigenbrot(high_freq_sinusoid, 512, 512)
save("high_freq_sinusoid.png", eb, ImageSetting(RealPart, LinearScale))
save("high_freq_sinusoid_dft.png", fft(eb), ImageSetting(Magn, LogScale))


text:
using LinearAlgebra

using DelimitedFiles
fdata = "C:\\Users\\jcai.AUTH\\Downloads\\SCODE-master\\SCODE-master\\data2\\data.txt" # ARGS[1]
ftime = "C:\\Users\\jcai.AUTH\\Downloads\\SCODE-master\\SCODE-master\\data2\\time.txt" # ARGS[2]
X=readdlm(fdata,'\t')
pseudotime=readdlm(ftime,'\t')[:,2]
pseudotime = pseudotime/maximum(pseudotime)

cnum=length(pseudotime)
tfnum = 100   #parse(Int32,ARGS[4])
pnum = 4      #parse(Int32,ARGS[5])
maxite = 10

maxB = 2.0
minB = -10.0

W = zeros(tfnum, pnum)
Z = zeros(pnum, cnum)
WZ = zeros(tfnum, cnum)

new_B = zeros(pnum)
old_B = zeros(pnum)

#initialization
RSS = Inf
for i in 1:pnum
    new_B[i] = rand()*(maxB-minB)+minB
    old_B[i] = new_B[i]
end

function sample_Z()
    for i in 1:pnum
        for j in 1:cnum
            Z[i,j] = exp(new_B[i]*pseudotime[j]) + (rand()*0.002-0.001)
        end
    end
end

#optimization
for ite in 1:maxite
    #sampling B
    target = rand(1:pnum)
    global new_B[target] = rand()*(maxB-minB)+minB

    #for last calc
    if ite==maxite
        global new_B = copy(old_B)
    end

    #sample Z from new B
    sample_Z()

    #regression
    for i in 1:tfnum
        W[i,:] = (Z * Z')^-1 * (Z * X[i,:])
    end
    #prediction
    WZ = W * Z

    #RSS
    tmp_RSS = 0
    for i in 1:tfnum
        tmp_RSS += sum((X[i,:]-WZ[i,:]).*(X[i,:]-WZ[i,:]))
    end

    if tmp_RSS < RSS
        global RSS = tmp_RSS
        global old_B[target] = new_B[target]
    else
        global new_B[target] = old_B[target]
    end
end

#open("RSS.txt", "w") do f
#    println(f, RSS)
#end

#infer A
B = zeros(pnum, pnum)
for i in 1:pnum
    B[i,i] = new_B[i]
end
invW = pinv(W)
A = W * B * invW

#writedlm("A.txt",A)
#writedlm("B.txt",B)


text:
function get_eop_history(obj::CIRFProvider)
    return jcall(obj, "getEOPHistory", EOPHistory, ())
end

function get_non_interpolating_provider(obj::CIRFProvider)
    return jcall(obj, "getNonInterpolatingProvider", CIRFProvider, ())
end

function get_transform(obj::CIRFProvider, arg0::AbsoluteDate)
    return jcall(obj, "getTransform", Transform, (AbsoluteDate,), arg0)
end

function get_transform(obj::CIRFProvider, arg0::FieldAbsoluteDate)
    return jcall(obj, "getTransform", FieldTransform, (FieldAbsoluteDate,), arg0)
end


text:
using Genie.Router, Genie.Requests, Genie.Renderer.Json
using JSON3

route("/") do
    serve_static_file("index.html")
end

route("/graph", method = GET) do
    JSON3.write(NodeGraphEditor.get_nodegraph())
end

route("/addnode/:add_type::String", method = POST) do
    types = NodeGraphEditor.get_param_types()

    add_type = Symbol(payload(:add_type))

    if !(add_type in types)
        return Genie.Router.error(404, "Invalid node type $add_type", MIME"text/html")
    end

    data = jsonpayload()
    new_pos = (Float64(data["x"]), Float64(data["y"]))

    ng = NodeGraphEditor.get_nodegraph()
    nn_id = NodeGraphEditor.get_new_node_id()

    add_node!(ng, add_type, nn_id, new_pos)

    nn_id += 1

    NodeGraphEditor.set_nodegraph(ng)
    NodeGraphEditor.set_new_node_id(nn_id)

    return JSON3.write(ng)
end

route("/updatenode/:id::Int/:key::String", method = POST) do
    ng = NodeGraphEditor.get_nodegraph()
    id = payload(:id)
    index = id_to_index(id, ng)
    isnothing(index) && return Genie.Router.error(404, "Invalid node_id: $id", MIME"text/html")

    key = Symbol(payload(:key))
    value = jsonpayload()["value"]

    if key == :name
        update_node_name!(ng, index, value)
    elseif key == :pos
        update_node_position!(ng, index, Float64(value["x"]), Float64(value["y"]))
    else
        return Genie.Router.error(400, "Invalid key: $key", MIME"text/html")
    end

    NodeGraphEditor.set_nodegraph(ng)
    return JSON3.write(ng)
end

route("/updateparams/:id::Int/:key::String", method = POST) do
    ng = NodeGraphEditor.get_nodegraph()

    id = payload(:id)
    index = findnext(n -> n.id == id, ng.nodes, 1)
    isnothing(index) && return Genie.Router.error(404, "Invalid node_id: $id", MIME"text/html")

    key = Symbol(payload(:key))
    value = rawpayload()

    update_node_params!(ng, index, key, value)

    NodeGraphEditor.set_nodegraph(ng)
    return JSON3.write(ng.nodes[index].params)
end

route("/getparams/:id::Int", method = GET) do
    ng = NodeGraphEditor.get_nodegraph()

    id = payload(:id)
    index = findnext(n -> n.id == id, ng.nodes, 1)
    isnothing(index) && return Genie.Router.error(404, "Invalid node_id: $id", MIME"text/html")

    return JSON3.write(ng.nodes[index].params)
end

route("/addedge", method = POST) do
    ng = NodeGraphEditor.get_nodegraph()

    edge = JSON3.read(rawpayload(), Edge)

    isvalid = add_edge!(ng, edge)

    isvalid || return Genie.Router.error(400, "Invalid connection", MIME"text/html")

    NodeGraphEditor.set_nodegraph(ng)
    return JSON3.write(ng)
end

route("/deleteedge", method = POST) do
    ng = NodeGraphEditor.get_nodegraph()

    edge = JSON3.read(rawpayload(), Edge)

    delete_edge!(ng, edge)

    NodeGraphEditor.set_nodegraph(ng)

    return JSON3.write(ng)
end

route("/deletenode/:id::Int", method = POST) do
    ng = NodeGraphEditor.get_nodegraph()

    id = payload(:id)
    ind = id_to_index(id, ng)
    isnothing(ind) && return Genie.Router.error(404, "Invalid ID: $id", MIME"text/html")

    delete_node!(ng, ind)

    NodeGraphEditor.set_nodegraph(ng)

    return JSON3.write(ng)
end

route("/types", method = GET) do
    types = NodeGraphEditor.get_param_types()
    return JSON3.write(Dict(:types => types))
end


text:
abstract type FactorHistogramQuery <: Query end

abstract type FactorHistogramQueries <: Queries end

mutable struct Factor
    attributes::Array{Int, 1}
    histogram::Histogram
end

"""
    FactorHistogram

FactorHistogram represents the data as a product of histograms on disjoint
attributes. This corresponds to a distribution where variables in different
factors are independent.
"""
struct FactorHistogram <: Data
    factors::Array{Factor, 1}
    lookup::Dict{Int, Int}
end

function attributes(q::FactorHistogramQuery)
    throw("`attributes` not implemented for `$(typeof(q))`.")
end

function default_value(q::FactorHistogramQuery)
    throw("`default_value` not implemented for `$(typeof(q))`.")
end

function restrict(q::FactorHistogramQuery, attributes::Array{Int, 1})
    throw("`restrict` not implemented for `$(typeof(q))`.")
end

function get_parent(d::FactorHistogram, attributes::Array{Int, 1})
    factor_idxs = [d.lookup[i] for i in attributes]
    all(i -> i == factor_idxs[1], factor_idxs) ? factor_idxs[1] : 0
end

function merge_factors(f::Factor, g::Factor)
    Factor(append!(f.attributes, g.attributes),
           Histogram(kron(f.histogram.weights, g.histogram.weights)))
end

function merge_components!(q::FactorHistogramQuery, d::FactorHistogram)
    parent_index = d.lookup[attributes(q)[1]]
    for i in 2:length(attributes(q))
        current_index = d.lookup[attributes(q)[i]]
        if current_index != parent_index
            d.factors[parent_index] = merge_factors(d.factors[parent_index],
                                                    d.factors[current_index])
            for a in d.factors[current_index].attributes
                d.lookup[a] = parent_index
            end
            d.factors[current_index].attributes = Int[]
            d.factors[current_index].histogram = Histogram([1.0])
        end
    end
end

function evaluate(q::FactorHistogramQuery, d::FactorHistogram)
    parent_index = get_parent(d, attributes(q))
    if parent_index > 0
        c = d.factors[parent_index]
        qc = restrict(q, c.attributes)
        evaluate(qc, c.histogram)
    else
        default_value(q)
    end
end

function update!(q::FactorHistogramQuery, d::FactorHistogram, error::Float64)
    parent_index = get_parent(d, attributes(q))
    if parent_index > 0
        c = d.factors[parent_index]
        qc = restrict(q, c.attributes)
        update!(qc, c.histogram, error)
        normalize!(c.histogram)
    else
        merge_components!(q, d)
        update!(q, d, error)
    end
end

# normalize happens inside update!
normalize!(h::FactorHistogram) = h

function initialize(queries::FactorHistogramQueries, table::Tabular, ps::MWParameters)
    d, n = size(table.data)
    components = [Factor([i], Histogram([0.5, 0.5])) for i = 1:d]
    lookup = Dict{Int, Int}()
    for i = 1:d
        lookup[i] = i
    end
    synthetic = FactorHistogram(components, lookup)
    real_answers = evaluate(queries, table)
    scale = 2.0/(ps.epsilon*n)
    MWState(table, synthetic, queries, real_answers, Dict{Int, Float64}(), scale)
end


text:
import Base: match, replace

export match

const Σ = Symbol
>ₑ(a::Σ, b::Σ) = a > b

abstract type AbstractTheory end
abstract type AbstractTerm end
abstract type AbstractMatcher end
abstract type AbstractSubproblem end
abstract type AbstractRewriter end

include("variable.jl")
include("match.jl")
include("rewriter.jl")


text:
using DistQuads, Distributions
using Test
import LinearAlgebra: norm

@testset "analytic vs quadrature" begin
    @testset "Beta" begin
        for a = 0.1:0.1:3.0, b = 0.1:0.1:3.0
            bd = Beta(a, b)
            dq = DistQuad(bd)
            # Test that mean function didn't accidentally regress
            @test mean(dq) == E(identity, dq)
            # Test nodes and weights against Distributions
            @test norm(mean(dq)-mean(bd), Inf) < 1e-8
            @test norm(var(dq)-var(bd), Inf) < 1e-8
            # Test nodes and weights against know mean (Distributions could in principle regress)
            @test norm(mean(dq)-a/(a+b), Inf) < 1e-8
            @test norm(var(dq)-a*b/((a+b+1)*(a+b)^2), Inf) < 1e-8
        end
    end
    @testset "Gamma" begin
        for a = 0.1:0.1:3.0, b = 0.1:0.1:3.0
            bd = Gamma(a, b)
            dq = DistQuad(bd)
            # Test that mean function didn't accidentally regress
            @test mean(dq) == E(identity, dq)
            # Test nodes and weights against Distributions
            @test norm(mean(dq)-mean(bd), Inf) < 1e-8
            @test norm(var(dq)-var(bd), Inf) < 1e-8
            # Test nodes and weights against know mean (Distributions could in principle regress)
            @test norm(mean(dq)-a*b, Inf) < 1e-8
            @test norm(var(dq)-a*b^2, Inf) < 1e-8
        end
    end
    @testset "Normal" begin
        for μ = 0.1:0.1:3.0, σ = 0.1:0.1:3.0
            bd = Normal(μ, σ)
            dq = DistQuad(bd)
            # Test that mean function didn't accidentally regress
            @test mean(dq) == E(identity, dq)
            # Test nodes and weights against Distributions
            @test norm(mean(dq)-mean(bd), Inf) < 1e-8
            @test norm(var(dq)-var(bd), Inf) < 1e-8
            # Test nodes and weights against know mean (Distributions could in principle regress)
            @test norm(mean(dq)-μ, Inf) < 1e-8
            @test norm(var(dq)-σ^2, Inf) < 1e-8
        end
    end
    @testset "LogNormal" begin
        for μ = 0.1:0.1:3.0, σ = 0.1:0.1:2.0
            bd = LogNormal(μ, σ)
            dq = DistQuad(bd)
            # Test that mean function didn't accidentally regress
            @test mean(dq) == E(identity, dq)
            # Test nodes and weights against Distributions
            @test norm(mean(dq)-mean(bd), Inf) < 1e-8
            @test norm(var(dq)-var(bd), Inf) < 1e-8 # not too precise!
            # Test nodes and weights against know mean (Distributions could in principle regress)
            @test norm(mean(dq)-exp(μ+σ^2/2), Inf) < 1e-8
            @test norm(var(dq)-(exp(σ^2)-1)*exp(2μ+σ^2), Inf) < 1e-8 # not too precise!
        end
    end
end


text:
using MLBase, Random

@testset "cross validation" begin

datapath = joinpath(dirname(@__FILE__), "data")

data = readcsvmat(joinpath(datapath,"gamlr.gaussian.data.csv"))
y = data[:,1]
X = data[:,2:end]
offset = fill(0.001,size(y))

path = fit(LassoPath, X, y; offset=offset)
β = coef(path; select=:all)

coefsAICc = coef(path;select=:AICc)
segminAICc = minAICc(path)
@test segminAICc == 71
@test coefsAICc == β[:,segminAICc]

coefsBIC = coef(path;select=:BIC)
segminBIC = Lasso.minBIC(path)
@test segminBIC == 55
@test coefsBIC == β[:,segminBIC]

coefsAIC = coef(path;select=:AIC)
segminAIC = Lasso.minAIC(path)
@test segminAIC == 71
@test coefsAIC == β[:,segminAIC]

Random.seed!(13)
gen = Kfold(length(y),10)
segCVmin = cross_validate_path(path;gen=gen)
coefsCVmin = coef(path;select=:CVmin)
@test segCVmin == 71
@test coefsCVmin == β[:,segCVmin]

Random.seed!(13)
gen = Kfold(length(y),10)
segCVmin = cross_validate_path(path,X,y; gen=gen, offset=offset)
coefsCVmin = coef(path;select=:CVmin)
@test segCVmin == 71
@test coefsCVmin == β[:,segCVmin]

Random.seed!(13)
coefsCV1se = coef(path;select=:CV1se,nCVfolds=20)
segCV1se = cross_validate_path(path,X,y;select=:CV1se,gen=gen,offset=offset)
@test segCV1se == 42
@test coefsCV1se == β[:,segCV1se]

end


text:
# Julia wrapper for header: /usr/include/scip/pub_event.h
# Automatically generated using Clang.jl wrap_c

function SCIPeventhdlrGetName(eventhdlr)
    ccall((:SCIPeventhdlrGetName, libscip), Cstring, (Ptr{SCIP_EVENTHDLR},), eventhdlr)
end

function SCIPeventhdlrGetData(eventhdlr)
    ccall((:SCIPeventhdlrGetData, libscip), Ptr{SCIP_EVENTHDLRDATA}, (Ptr{SCIP_EVENTHDLR},), eventhdlr)
end

function SCIPeventhdlrSetData(eventhdlr, eventhdlrdata)
    ccall((:SCIPeventhdlrSetData, libscip), Cvoid, (Ptr{SCIP_EVENTHDLR}, Ptr{SCIP_EVENTHDLRDATA}), eventhdlr, eventhdlrdata)
end

function SCIPeventhdlrIsInitialized(eventhdlr)
    ccall((:SCIPeventhdlrIsInitialized, libscip), UInt32, (Ptr{SCIP_EVENTHDLR},), eventhdlr)
end

function SCIPeventhdlrGetSetupTime(eventhdlr)
    ccall((:SCIPeventhdlrGetSetupTime, libscip), Cdouble, (Ptr{SCIP_EVENTHDLR},), eventhdlr)
end

function SCIPeventhdlrGetTime(eventhdlr)
    ccall((:SCIPeventhdlrGetTime, libscip), Cdouble, (Ptr{SCIP_EVENTHDLR},), eventhdlr)
end

function SCIPeventGetType(event)
    ccall((:SCIPeventGetType, libscip), SCIP_EVENTTYPE, (Ptr{SCIP_EVENT},), event)
end

function SCIPeventGetVar(event)
    ccall((:SCIPeventGetVar, libscip), Ptr{SCIP_VAR}, (Ptr{SCIP_EVENT},), event)
end

function SCIPeventGetOldobj(event)
    ccall((:SCIPeventGetOldobj, libscip), Cdouble, (Ptr{SCIP_EVENT},), event)
end

function SCIPeventGetNewobj(event)
    ccall((:SCIPeventGetNewobj, libscip), Cdouble, (Ptr{SCIP_EVENT},), event)
end

function SCIPeventGetOldbound(event)
    ccall((:SCIPeventGetOldbound, libscip), Cdouble, (Ptr{SCIP_EVENT},), event)
end

function SCIPeventGetNewbound(event)
    ccall((:SCIPeventGetNewbound, libscip), Cdouble, (Ptr{SCIP_EVENT},), event)
end

function SCIPeventGetNode(event)
    ccall((:SCIPeventGetNode, libscip), Ptr{SCIP_NODE}, (Ptr{SCIP_EVENT},), event)
end

function SCIPeventGetSol(event)
    ccall((:SCIPeventGetSol, libscip), Ptr{SCIP_SOL}, (Ptr{SCIP_EVENT},), event)
end

function SCIPeventGetHoleLeft(event)
    ccall((:SCIPeventGetHoleLeft, libscip), Cdouble, (Ptr{SCIP_EVENT},), event)
end

function SCIPeventGetHoleRight(event)
    ccall((:SCIPeventGetHoleRight, libscip), Cdouble, (Ptr{SCIP_EVENT},), event)
end

function SCIPeventGetRow(event)
    ccall((:SCIPeventGetRow, libscip), Ptr{SCIP_ROW}, (Ptr{SCIP_EVENT},), event)
end

function SCIPeventGetRowCol(event)
    ccall((:SCIPeventGetRowCol, libscip), Ptr{SCIP_COL}, (Ptr{SCIP_EVENT},), event)
end

function SCIPeventGetRowOldCoefVal(event)
    ccall((:SCIPeventGetRowOldCoefVal, libscip), Cdouble, (Ptr{SCIP_EVENT},), event)
end

function SCIPeventGetRowNewCoefVal(event)
    ccall((:SCIPeventGetRowNewCoefVal, libscip), Cdouble, (Ptr{SCIP_EVENT},), event)
end

function SCIPeventGetRowOldConstVal(event)
    ccall((:SCIPeventGetRowOldConstVal, libscip), Cdouble, (Ptr{SCIP_EVENT},), event)
end

function SCIPeventGetRowNewConstVal(event)
    ccall((:SCIPeventGetRowNewConstVal, libscip), Cdouble, (Ptr{SCIP_EVENT},), event)
end

function SCIPeventGetRowSide(event)
    ccall((:SCIPeventGetRowSide, libscip), SCIP_SIDETYPE, (Ptr{SCIP_EVENT},), event)
end

function SCIPeventGetRowOldSideVal(event)
    ccall((:SCIPeventGetRowOldSideVal, libscip), Cdouble, (Ptr{SCIP_EVENT},), event)
end

function SCIPeventGetRowNewSideVal(event)
    ccall((:SCIPeventGetRowNewSideVal, libscip), Cdouble, (Ptr{SCIP_EVENT},), event)
end


text:
using WingBiology

using LinearAlgebra

using DelimitedFiles

include("airfoils.jl")
include("structural_properties.jl")

root_afl = afl350
root_afl_Arel = afl350_A

tip_afl = afl1
tip_afl_Arel = afl1_A

acft = Aircraft()

# geometry
b = 3.0
cr = 0.5
ct = 0.25
incidence = 8.0
twist = -5.0
sweep = 15.0

R_tube = 10e-3
t_tube = 2e-3
t_skin = 1e-3
e = -0.05

# root and tip structural properties
EAt, EIt, GJt, mt, Ixxt, mxt = structural_tube(
    R_tube, R_tube + t_tube, ct; e = e
) .+ structural_skin(
    tip_afl_Arel,
    t_skin,
    ct
)
EAr, EIr, GJr, mr, Ixxr, mxr = structural_tube(
    R_tube, R_tube + t_tube, cr; e = e
) .+ structural_skin(
    root_afl_Arel,
    t_skin,
    cr
)

s1 = WingSectInfo(
    [tand(sweep) * b / 2, - b / 2, 0.0],
    ct;
    EIy = EIt, # flexural rigidity. EIz can also be used for rigidity around the vertical axis
    GJ = GJt, # local torsional rigidity
    m = mt, # local mass per unit span
    xCG = mxt / mt, # chord perc. of CG from quarter chord
    Ixx = Ixxt, # inertia for rotation around the wing axis
    e = e, # distance from quarter chord to elastic axis. I didn't calculate this, but you get the idea
    incidence = incidence + twist
)
s2 = WingSectInfo(
    [0.0, 0.0, 0.0],
    cr;
    EIy = EIr,
    GJ = GJr,
    m = mr,
    xCG = mxr / mr,
    Ixx = Ixxr,
    e = e,
    incidence = incidence
)
s3 = WingSectInfo(
    [tand(sweep) * b / 2, b / 2, 0.0],
    ct;
    EIy = EIt,
    GJ = GJt,
    m = mt,
    xCG = mxt / mt,
    Ixx = Ixxt,
    e = e,
    incidence = incidence + twist
)

interp_points, interp_wing_strips = interpolate!(acft, 20, s1, s2) # discretize with 20 strips
interpolate!(acft, 20, s2, s3) # discretize with 20 strips

add_link!(acft, interp_points[end])

U∞ = 25.0

q = get_state(acft, U∞)
q0 = copy(q)

@info "Using Arnoldi method to solve for static aerostructural solution"

A, eig = Arnoldi_solve!(acft, q, U∞)

plot_aircraft(acft, q)
qd, fcon = state_space(acft, q, U∞)

n_eig = 6

#=
@info "Obtaining assumed modes from structural matrices"
ϕ, ωs = get_eigenmodes(acft, n_eig)

for (i, ω) in enumerate(ωs)
    @show i, ω

    plot_aircraft(acft, ϕ[:, i] .+ q0)
end
=#

# I haven't yet corrected the flutter prediction feature!

@info "Using assumed modes method to deduce assumed mode eigenvalues"

A, eig = assumed_modes_solve(
    acft, q, U∞;
    n_eig = n_eig, structural_damping = true, transient_aerodynamics = true
)
for (i, λ) in enumerate(eig.λ)
    @show i, λ
end

#=
@show get_̇w(A[:, 3])
@show get_̇w(eig.JA[:, 3])
@show get_̇θy(A[:, 3])
@show get_̇θy(eig.JA[:, 3])
@show get_̇u(A[:, 3]) ⋅ get_̇u(eig.JA[:, 3]), get_̇v(A[:, 3]) ⋅ get_̇v(eig.JA[:, 3]), get_̇w(A[:, 3]) ⋅ get_̇w(eig.JA[:, 3]), get_̇θx(A[:, 3]) ⋅ get_̇θx(eig.JA[:, 3]), get_̇θy(A[:, 3]) ⋅ get_̇θy(eig.JA[:, 3])
=#

writedlm("eig_basis.dat", [q0 real.(A) imag.(A)])

writedlm("assumed_eigs.dat", [real.(eig.ϕ) imag.(eig.ϕ) real.(eig.λ) imag.(eig.λ)])

@info "Performing dynamic simulation"

A, ts, qs = backward_Euler(acft, q0, U∞; transient_aerodynamics = true, dt = 1e-3, t = 0.2)

writedlm("Arnoldi_basis.dat", [q0 A])
writedlm("results.dat", [ts qs'])

#=
@info "Using full Jacobian to deduce eigenvalues"

eig = get_eigen(acft, q0, U∞) # ; fixed_points = fixed_points)

for (i, λ) in enumerate(eig.values)
    if abs(λ) < 40.0 * (2 * π)
        @show i, λ
    end
end
=#

@info "Calculating masses and inertias"

@show get_massic_properties(acft)

max_stars_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_issues_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_forks_repo_name: Equipe-AeroDesign-ITA/WingBiology


text:
using Test, LinearAlgebra
using WaveProp
using WaveProp.Nystrom
using WaveProp.Geometry
using WaveProp.Integration
using WaveProp.Mesh

@testset "Basic tests" begin
    Geometry.clear!()
    Ω,M   = WaveProp.IO.gmsh_sphere(dim=2,h=0.05)
    Γ     = boundary(Ω)
    mesh  = view(M,boundary(Ω))
    nmesh = NystromMesh(mesh,order=2)
    @test isapprox(sum(qweights(nmesh)),π,atol=0.1)
    Geometry.clear!()
    geo   = Circle(radius=0.5)
    Ω     = Domain(geo)
    M     = meshgen(Ω;h=(0.1))
    mesh  = view(M,boundary(Ω))
    nmesh = NystromMesh(mesh,order=2)
    @test isapprox(sum(qweights(nmesh)),π,atol=0.1)
end

@testset "Mesh quadrature" begin
    @testset "Cube" begin
        # generate a mesh
        Geometry.clear!()
        (lx,ly,lz) = widths = (1.,1.,2.)
        Ω, M  = WaveProp.IO.gmsh_box(;widths=widths)
        ∂Ω = boundary(Ω)
        mesh  = NystromMesh(view(M,∂Ω),order=1)
        A     = 2*(lx*ly + lx*lz + ly*lz)
        @test A ≈ sum(qweights(mesh))
        # generate a Nystrom mesh for volume. Note that since we are creating a
        # volume Nystrom mesh, we pass the kwarg compute_normal=false since
        # normals do not make sense in this case
        mesh  = NystromMesh(view(M,Ω),order=1,compute_normal=false)
        V     = prod(widths)
        # sum only weights corresponding to tetras
        @test V ≈ sum(qweights(mesh))
        # finally generate a Nystrom mesh for the surface AND volume
        Ω₁ = union(Ω,∂Ω)
        mesh  = NystromMesh(view(M,Ω₁),order=1,compute_normal=false)
        @test V+A ≈ sum(qweights(mesh))
    end
    @testset "Sphere" begin
        r = 0.5
        Geometry.clear!()
        Ω, M = WaveProp.IO.gmsh_sphere(;radius=r)
        Γ = boundary(Ω)
        A = 4π*r^2
        mesh = NystromMesh(view(M,Γ),order=2)
        # the coarse tolerance below is because we use flat elements to
        # approximate the surface area and volume of a sphere
        @test isapprox(A,sum(qweights(mesh));atol=0.1)
        # test volume
        mesh = NystromMesh(view(M,Ω),order=2;compute_normal=false)
        A = 4π*r^2
        V = (4/3)*π*r^3
        @test isapprox(V,sum(qweights(mesh));atol=0.1)
    end
    @testset "Disk" begin
        r = rx = ry = 0.5
        Geometry.clear!()
        Ω, M = WaveProp.IO.gmsh_disk(;rx,ry)
        Γ = boundary(Ω)
        mesh = NystromMesh(view(M,Ω),order=2,compute_normal=false)
        A = π*r^2
        # the coarse tolerance below is because we use flat elements to
        # approximate the surface area and volume of a sphere
        @test isapprox(A,sum(qweights(mesh));atol=0.1)
        # test perimeter
        mesh = NystromMesh(view(M,Γ),order=2)
        P = 2π*r
        @test isapprox(P,sum(qweights(mesh));atol=0.1)
    end
    @testset "Append" begin
        Geometry.clear!()
        Ω, M = WaveProp.IO.gmsh_disk()
        Γ = boundary(Ω)
        mesh = NystromMesh(view(M,Ω),order=2,compute_normal=false)
        Σ_mesh = NystromMesh{2,Float64}()
        append!(Σ_mesh,mesh)
        append!(Σ_mesh,mesh)
    end
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
struct At end
struct After end
struct Now end
struct Every end

const after = After()
const at = At()
const now = Now()
const every = Every()

abstract type AbstractEvent end

"""
    Event <: AbstractEvent
- `fun`: a function
- `time`: time of function call
- `id`: an event id
- `type`: type of event
- `description`: a description of the event
"""
mutable struct Event{F<:Function} <: AbstractEvent
    fun::F
    time::Float64
    id::String
    type::String
    description::String
end

abstract type AbstractScheduler end

"""
Scheduler <: AbstractScheduler
- `events`: a priority queue of events
- `time`: current time of the system
- `running`: simulation can run if true
- `trace`: will print out the events if true
- `store`: will store a vector of completed events if true
- `complete_events`: an optional vector of completed events
"""
mutable struct Scheduler{PQ<:PriorityQueue,E} <: AbstractScheduler
    events::PQ
    time::Float64
    running::Bool
    trace::Bool
    store::Bool
    complete_events::E
end

"""
Constructor for Scheduler with default keyword values:

```julia
Scheduler(;event=Event, time=0.0, running=true, trace=false, store=false)
```
"""
function Scheduler(;event=Event, time=0.0, running=true, trace=false, store=false)
    events = PriorityQueue{event,Float64}()
    return Scheduler(events, time, running, trace, store, Vector{event}())
end


text:
#
# Returns a single parameter and its weight - for simulation
#
function generate_parameters(
        priors::AbstractArray{D,1},
        ) where {
        D<:ContinuousUnivariateDistribution
        }
    n_dims = length(priors)
    parameters = zeros(n_dims)

    weight = 1.
    @inbounds for i in 1:n_dims
        parameters[i] = rand(priors[i])
        weight *= pdf(priors[i], parameters[i])
    end

    return reshape(parameters, (1, n_dims)), weight
end

#
# Returns a set parameter and its weight with size batch_size - for emulation
#
function generate_parameters(
        priors::AbstractArray{D,1},
        batch_size::Int,
        ) where {
        D<:ContinuousUnivariateDistribution
        }

    n_dims = length(priors)
    parameters = zeros(batch_size, n_dims)
    priors = reshape(priors, 1, n_dims)
    parameters .= rand.(priors)
    weights = prod(pdf.(priors, parameters), dims=2)
    weights = reshape(weights, batch_size) # So that this returns a 1D array, like the simulation version

    return parameters, weights
end

function ABCrejection(input::SimulatedABCRejectionInput;
    write_progress::Bool = true,
    progress_every::Int = 1000)

	checkABCInput(input)
    if write_progress
        @info "GpABC rejection simulation. ϵ = $(input.threshold)."
    end

	# initialise
    n_tries = 0
    n_accepted = 0
    accepted_parameters = zeros(input.n_particles, input.n_params)
    accepted_distances = zeros(input.n_particles)
    weight_values = zeros(input.n_particles)
    distance = 0.0

    # simulate
    while n_accepted < input.n_particles && n_tries < input.max_iter

        parameters, weight_value = generate_parameters(input.priors)
        try
            distance = simulate_distance(parameters, input.distance_simulation_input)
        catch e
            if isa(e, DimensionMismatch)
                # This prevents the whole code from failing if there is a problem
                # solving the differential equation(s). The exception is thrown by the
                # distance function
                @warn "The summarised simulated data does not have the same size as the summarised reference data. If this is not happening at every iteration it may be due to the behaviour of DifferentialEquations::solve - please check for related warnings. Continuing to the next iteration."
                n_tries += 1
                continue
            else
                throw(e)
            end
        end

        n_tries += 1

        if distance[1] <= input.threshold
            n_accepted += 1
            accepted_parameters[n_accepted,:] = parameters
            accepted_distances[n_accepted] = distance[1]
            weight_values[n_accepted] = weight_value
        end

        if write_progress && (n_tries % progress_every == 0)
            @info "GpABC rejection simulation. Accepted $(n_accepted)/$(n_tries) particles."
        end
    end

    if n_accepted < input.n_particles
        @warn "Simulation reached maximum iterations $(input.max_iter) before finding $(input.n_particles) particles - will return $n_accepted"
        accepted_parameters = accepted_parameters[1:n_accepted, :]
        accepted_distances = accepted_distances[1:n_accepted]
        weight_values = weight_values[1:n_accepted]
    end

    # output
    output = SimulatedABCRejectionOutput(input.n_params,
                                n_accepted,
                                n_tries,
                                input.threshold,
                                accepted_parameters,
                                accepted_distances,
                                StatsBase.Weights(weight_values ./ sum(weight_values)) # normalise weights
                                )

    return output

end

function ABCrejection(input::EmulatedABCRejectionInput;
    write_progress = true,
    progress_every = 1000)

    checkABCInput(input)

    if write_progress
        @info "GpABC rejection emulation. ϵ = $(input.threshold)."
    end
    # initialise
    n_accepted = 0
    n_tries = 0
    batch_no = 1
    accepted_parameters = zeros(input.n_particles, input.n_params)
    accepted_distances = zeros(input.n_particles)
    weights = ones(input.n_particles)

    # emulate
    emulator = abc_train_emulator(input.priors, input.emulator_training_input)
    while n_accepted < input.n_particles && batch_no <= input.max_iter
        parameter_batch, weight_batch = generate_parameters(input.priors, input.batch_size)

        # distances, vars = gp_regression(parameter_batch, emulator)
        n_tries += input.batch_size
        #
        # Check which parameter indices were accepted
        #
        distances, accepted_batch_idxs = abc_select_emulated_particles(emulator, parameter_batch, input.threshold, input.selection)

        n_accepted_batch = length(accepted_batch_idxs)

        #
        # If some parameters were accepted, store their values, (predicted) distances and weights
        #
        if n_accepted_batch > 0
            # Check that we won't accept too many parameters - throw away the extra parameters if so
            if n_accepted + n_accepted_batch > input.n_particles
                accepted_batch_idxs = accepted_batch_idxs[1:input.n_particles - n_accepted]
                n_accepted_batch = length(accepted_batch_idxs)
                distances = distances[1:n_accepted_batch]
            end

            accepted_parameters[n_accepted+1:n_accepted + n_accepted_batch,:] = parameter_batch[accepted_batch_idxs,:]
            accepted_distances[n_accepted+1:n_accepted + n_accepted_batch] = distances
            weights[n_accepted+1:n_accepted + n_accepted_batch] = weight_batch[accepted_batch_idxs]
            n_accepted += n_accepted_batch

        end

        if write_progress
            @info "GpABC rejection emulation. Accepted $(n_accepted)/$(n_tries) particles ($(batch_no) batches of size $(input.batch_size))."
        end

        batch_no += 1
    end

    if n_accepted < input.n_particles
        @warn "Emulation reached maximum $(input.max_iter) iterations before finding $(input.n_particles) particles - will return $n_accepted"
        accepted_parameters = accepted_parameters[1:n_accepted, :]
        accepted_distances = accepted_distances[1:n_accepted]
        weights = weights[1:n_accepted]
    end

    # output
    output = EmulatedABCRejectionOutput(input.n_params,
                                n_accepted,
                                n_tries,
                                input.threshold,
                                accepted_parameters,
                                accepted_distances,
                                StatsBase.Weights(weights ./ sum(weights)),
                                emulator
                                )

    return output
end


text:
module clock

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
@testset "Coupling constraints" begin
    cp = convert(CoreModelCoupled, test_LP())
    @test size(cp.lm.S) == (4, 3)
    @test size(stoichiometry(convert(CoreModel, cp))) == (4, 3)
    new_cp = add_coupling_constraints(cp, stoichiometry(cp)[end, :], -1.0, 1.0)
    @test n_coupling_constraints(cp) + 1 == n_coupling_constraints(new_cp)

    new_cp =
        add_coupling_constraints(cp, stoichiometry(cp)[1:2, :], [-1.0; -1.0], [1.0; 1.0])
    @test n_coupling_constraints(cp) + 2 == n_coupling_constraints(new_cp)

    n_c = n_coupling_constraints(cp)
    add_coupling_constraints!(cp, stoichiometry(cp)[end, :], -1.0, 1.0)
    @test n_c + 1 == n_coupling_constraints(cp)
    add_coupling_constraints!(cp, stoichiometry(cp)[1:2, :], [-1.0; -1.0], [1.0; 1.0])
    @test n_c + 3 == n_coupling_constraints(cp)

    n_c = n_coupling_constraints(cp)
    remove_coupling_constraints!(cp, 1)
    @test n_c - 1 == n_coupling_constraints(cp)
    remove_coupling_constraints!(cp, [1, 2])
    @test n_c - 3 == n_coupling_constraints(cp)
    @test n_coupling_constraints(cp) == 0

    cp = test_coupledLP()
    n_c = n_coupling_constraints(cp)
    new_cp = remove_coupling_constraints(cp, 1)
    @test size(coupling(cp)) == (n_c, n_reactions(cp))
    @test n_c - 1 == n_coupling_constraints(new_cp)
    @test n_coupling_constraints(cp) == n_c
    new_cp = remove_coupling_constraints(cp, [1, 2])
    @test n_c - 2 == n_coupling_constraints(new_cp)
    new_cp = remove_coupling_constraints(cp, Vector(1:n_coupling_constraints(cp)))
    @test n_coupling_constraints(new_cp) == 0
    @test n_coupling_constraints(cp) == n_c

    cp = test_coupledLP()
    change_coupling_bounds!(cp, [3, 1], cl = [-10.0, -20], cu = [10.0, 20])
    cl, cu = coupling_bounds(cp)
    @test cl[[1, 3]] == [-20, -10]
    @test cu[[1, 3]] == [20, 10]
    change_coupling_bounds!(cp, [1000, 1001], cl = [-50.0, -60.0])
    cl, cu = coupling_bounds(cp)
    @test cl[[1000, 1001]] == [-50.0, -60.0]
end


text:
@testset "979.distribute-coins-in-binary-tree.jl" begin
    @test distribute_coins(TreeNode(3, TreeNode(0), TreeNode(0))) == 2
    @test distribute_coins(TreeNode(0, TreeNode(3), TreeNode(0))) == 3
    @test distribute_coins(TreeNode(1, TreeNode(2), TreeNode(0))) == 2
end


text:
function check_hypers_size(gpem::GPModel,
        hp_lower::AbstractArray{Float64, 1}, hp_upper::AbstractArray{Float64, 1})
    expected_hypers_size = get_hyperparameters_size(gpem)
    if length(hp_lower) > 0 && length(hp_lower) != expected_hypers_size
        error("Incorrect size of hyperparameters vector lower bound for ",
            "$(typeof(gpem.kernel)): $(length(hp_lower)). ",
            "Expected $(expected_hypers_size).")
    end
    if length(hp_upper) > 0 && length(hp_upper) != expected_hypers_size
        error("Incorrect size of hyperparameters vector upper bound for ",
            "$(typeof(gpem.kernel)): $(length(hp_upper)). ",
            "Expected $(expected_hypers_size).")
    end
end

"""
    gp_train(gpm::GPModel; <optional keyword arguments>)

Find Maximum Likelihood Estimate of Gaussian Process hyperparameters by maximising
[`gp_loglikelihood`](@ref), using [`Optim`](http://julianlsolvers.github.io/Optim.jl/stable/) package. The optimisation target is [`gp_loglikelihood_log`](@ref), with gradient computed by [`gp_loglikelihood_grad`](@ref). Internally, this function optimises the MLE with respect to logarithms of hyperparameters. This is done for numerical stability. Logarithmisation and exponentiation is performed by this funtion, i.e. real hyperparameters, not logarithms, are taken in and returned back.

By default,
[Conjugate Gradient](http://julianlsolvers.github.io/Optim.jl/stable/algo/cg/) bounded box optimisation is used, as long as the gradient
with respect to hyperparameters ([`covariance_grad`](@ref)) is implemented for the kernel function. If the gradient
implementation is not provided, [Nelder Mead](http://julianlsolvers.github.io/Optim.jl/stable/algo/nelder_mead/) optimiser is used by default.

# Mandatory argument
- `gpm`: the [`GPModel`](@ref), that contains the training data (x and y),
  the kernel and the starting hyperparameters that will be used for optimisation.

# Optional keyword arguments
- `optimiser::Type{<:Optim.AbstractOptimizer}`: the solver to use.
  If not given, then `ConjugateGradient` will be used for kernels that have gradient
  implementation, and `NelderMead` will be used for those that don't.
- `hp_lower::AbstractArray{Float64, 1}`: the lower boundary for box optimisation.
  Defaults to ``e^{-10}`` for all hyperparameters.
- `hp_upper::AbstractArray{Float64, 1}`: the upper boundary for box optimisation.
  Defaults to ``e^{10}`` for all hyperparameters.
- `log_level::Int`: log level. Default is `0`, which is no logging at all. `1`
  makes `gp_train` print basic information to standard output. `2` switches `Optim`
  logging on, in addition to `1`.

# Return
The list of all hyperparameters, including the standard deviation of the measurement
noise ``\\sigma_n``. Note that after this function returns, the hyperparameters of `gpm`
will be set to the optimised value, and there is no need to call [`set_hyperparameters`](@ref)
once again.
"""
function gp_train(gpem::GPModel;
        optimiser::Union{Nothing, TOpt}=nothing,
        hp_lower::AbstractArray{Float64, 1}=zeros(0),
        hp_upper::AbstractArray{Float64, 1}=zeros(0),
        log_level::Int=0) where {TOpt<:Optim.AbstractOptimizer}
    check_hypers_size(gpem, hp_lower, hp_upper)
    hypers = log.(gpem.gp_hyperparameters)
    kernel_hypers = hypers[1:end-1]
    test_grad_result = covariance_grad(gpem.kernel,
        kernel_hypers, gpem.gp_training_x,
        zeros(size(gpem.gp_training_x, 1), size(gpem.gp_training_x, 1)))
    show_trace = log_level >= 2
    opt_res = nothing
    if test_grad_result === :covariance_gradient_not_implemented
        if optimiser === nothing
            optimiser = NelderMead()
        end
        if log_level > 0
            println("Unbound optimisation using $(typeof(optimiser)). ",
            "No gradient provided. Start point: $(hypers).")
        end
        try
            opt_res = optimize(theta->-gp_loglikelihood_log(theta, gpem),
                hypers,
                optimiser)
        catch ex
            last_hypers = gpem.cache.theta
            gpem.gp_hyperparameters = exp.(last_hypers)
            println("Exception ", ex, " for hyperparameters $(gpem.gp_hyperparameters)")
            rethrow(ex)
        end
    else
        if optimiser === nothing
            optimiser = ConjugateGradient()
        end
        expected_hypers_size = get_hyperparameters_size(gpem)
        if length(hp_lower) == 0
            hp_lower = fill(-10.0, expected_hypers_size)
        else
            hp_lower = log.(hp_lower)
        end
        if length(hp_upper) == 0
            hp_upper = fill(10.0, expected_hypers_size)
        else
            hp_upper = log.(hp_upper)
        end

        if log_level > 0
            println("Bound optimisation using $(typeof(optimiser)). ",
            "Gradient provided. Start point: $(gpem.gp_hyperparameters); ",
            "lower bound: $(hp_lower); upper bound: $(hp_upper)")
        end
        try
            opt_res = optimize(theta->-gp_loglikelihood_log(theta, gpem),
                (storage, theta) -> storage[:] = -gp_loglikelihood_grad(theta, gpem),
                hp_lower, hp_upper,
                hypers,
                Fminbox(optimiser),
                Optim.Options(show_trace = show_trace))
        catch ex
            last_hypers = gpem.cache.theta
            gpem.gp_hyperparameters = exp.(last_hypers)
            println("Exception ", ex, " for hyperparameters $(last_hypers)")
            rethrow(ex)
        end
    end
    if log_level > 0
         println(opt_res)
         println("Optimized log hyperparameters: ", opt_res.minimizer)
         println("Optimized hyperparameters: ", exp.(opt_res.minimizer))
    end
    gpem.gp_hyperparameters = exp.(opt_res.minimizer)
    return exp.(opt_res.minimizer)
end;


text:
"""
    icrs2g(dr, dd)

Transformation from ICRS to Galactic Coordinates.

### Given ###

- `dr`: ICRS right ascension (radians)
- `dd`: ICRS declination (radians)

### Returned ###

- `dl`: Galactic longitude (radians)
- `db`: Galactic latitude (radians)

### Notes ###

1. The IAU 1958 system of Galactic coordinates was defined with
   respect to the now obsolete reference system FK4 B1950.0.  When
   interpreting the system in a modern context, several factors have
   to be taken into account:

   . The inclusion in FK4 positions of the E-terms of aberration.

   . The distortion of the FK4 proper motion system by differential
     Galactic rotation.

   . The use of the B1950.0 equinox rather than the now-standard
     J2000.0.

   . The frame bias between ICRS and the J2000.0 mean place system.

   The Hipparcos Catalogue (Perryman & ESA 1997) provides a rotation
   matrix that transforms directly between ICRS and Galactic
   coordinates with the above factors taken into account.  The
   matrix is derived from three angles, namely the ICRS coordinates
   of the Galactic pole and the longitude of the ascending node of
   the galactic equator on the ICRS equator.  They are given in
   degrees to five decimal places and for canonical purposes are
   regarded as exact.  In the Hipparcos Catalogue the matrix
   elements are given to 10 decimal places (about 20 microarcsec).
   In the present ERFA function the matrix elements have been
   recomputed from the canonical three angles and are given to 30
   decimal places.

2. The inverse transformation is performed by the function eraG2icrs.

### Called ###

- `eraAnp`: normalize angle into range 0 to 2pi
- `eraAnpm`: normalize angle into range +/- pi
- `eraS2c`: spherical coordinates to unit vector
- `eraRxp`: product of r-matrix and p-vector
- `eraC2s`: p-vector to spherical

### Reference ###

- Perryman M.A.C. & ESA, 1997, ESA SP-1200, The Hipparcos and Tycho
    catalogues.  Astrometric and photometric star catalogues
    derived from the ESA Hipparcos Space Astrometry Mission.  ESA
    Publications Division, Noordwijk, Netherlands.

"""
function icrs2g(a, b)
    r1 = Ref(0.0)
    r2 = Ref(0.0)
    ccall((:eraIcrs2g, liberfa), Cvoid,
            (Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
            a, b, r1, r2)
    r1[], r2[]
end


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

module AbsTest
using ParallelAccelerator

#ParallelAccelerator.DomainIR.set_debug_level(4)
#ParallelAccelerator.ParallelIR.set_debug_level(4)
#ParallelAccelerator.CGen.set_debug_level(4)
#ParallelAccelerator.set_debug_level(4)

@acc function example(x)
    abs(x)
end

function test1()
    return example(-3)
end

function test2()
    return example(-3.0)
end

function test3()
    return example(3.0)
end

function test4()
    A = zeros(10, 10).-1
    return example(A)
end

function test5()
    A = zeros(10, 10).-1.1
    return example(A)
end

function test6()
    A = zeros(10, 10).+1.0
    return example(A)
end

end

using Base.Test

println("Testing abs()...")
@test AbsTest.test1() == 3
@test AbsTest.test2() == 3.0
@test AbsTest.test3() == 3.0
@test AbsTest.test4() == ones(10, 10)
@test AbsTest.test5() == ones(10, 10).+0.1
@test AbsTest.test6() == ones(10, 10)
println("Done testing abs().")


text:
module VectorizationBase

using LinearAlgebra, Libdl
const LLVM_SHOULD_WORK = isone(length(filter(lib->occursin(r"LLVM\b", basename(lib)), Libdl.dllist())))

# isfile(joinpath(@__DIR__, "cpu_info.jl")) || throw("File $(joinpath(@__DIR__, "cpu_info.jl")) does not exist. Please run `using Pkg; Pkg.build()`.")

# using Base: llvmcall
# using Base: llvmcall
@inline llvmcall(s::String, args...) = Base.llvmcall(s, args...)
@inline llvmcall(s::Tuple{String,String}, args...) = Base.llvmcall(s, args...)

export Vec, VE, SVec, Mask, _MM,
    firstval, gep, gesp,
    extract_data,
    pick_vector_width,
    pick_vector_width_shift,
    stridedpointer,
    PackedStridedPointer, RowMajorStridedPointer,
    StaticStridedPointer, StaticStridedStruct,
    vload, vstore!, vbroadcast, Static, mask, masktable

# @static if VERSION < v"1.4"
#     # I think this is worth using, and simple enough that I may as well.
#     # I'll uncomment when I find a place to use it.
#     function only(x)
#         @boundscheck length(x) == 0 && throw(ArgumentError("Collection is empty, must contain exactly 1 element"))
#         @boundscheck length(x) > 1 && throw(ArgumentError("Collection has multiple elements, must contain exactly 1 element"))
#         @inbounds x[1]
#     end
#     export only
# end
const IntTypes = Union{Int8, Int16, Int32, Int64, Int128}
const UIntTypes = Union{UInt8, UInt16, UInt32, UInt64, UInt128}
const IntegerTypes = Union{IntTypes, UIntTypes, Ptr, Bool}
const FloatingTypes = Union{Float16, Float32, Float64}
const ScalarTypes = Union{IntegerTypes, FloatingTypes}

const VE{T} = Core.VecElement{T}
const Vec{W,T<:Number} = NTuple{W,VE{T}}
const _Vec{W,T<:Number} = Tuple{VE{T},Vararg{VE{T},W}}

abstract type AbstractStructVec{W,T<:Number} end
struct SVec{W,T} <: AbstractStructVec{W,T}
    data::Vec{W,T}
    # SVec{N,T}(v) where {N,T} = new(v)
end
struct Mask{W,U<:Unsigned} <: AbstractStructVec{W,Bool}
    u::U
end
const AbstractMask{W} = Union{Mask{W}, SVec{W,Bool}}
@inline Mask{W}(u::U) where {W,U<:Unsigned} = Mask{W,U}(u)
# Const prop is good enough; added an @inferred test to make sure.
@inline Mask(u::U) where {U<:Unsigned} = Mask{sizeof(u)<<3,U}(u)

@inline Base.broadcastable(v::AbstractStructVec) = Ref(v)

# SVec{N,T}(x) where {N,T} = SVec(ntuple(i -> VE(T(x)), Val(N)))
# @inline function SVec{N,T}(x::Number) where {N,T}
    # SVec(ntuple(i -> VE(T(x)), Val(N)))
# end
# @inline function SVec{N,T}(x::Vararg{<:Number,N}) where {N,T}
    # SVec(ntuple(i -> VE(T(x[i])), Val(N)))
# end
# @inline function SVec(v::Vec{N,T}) where {N,T}
    # SVec{N,T}(v)
# end
@inline SVec(u::Unsigned) = u # Unsigned integers are treated as vectors of bools
@inline SVec{W}(u::U) where {W,U<:Unsigned} = Mask{W,U}(u) # Unsigned integers are treated as vectors of bools
@inline SVec(v::SVec{W,T}) where {W,T} = v
@inline SVec{W}(v::SVec{W,T}) where {W,T} = v
@inline SVec{W,T}(v::SVec{W,T}) where {W,T} = v
@inline SVec{W}(v::Vec{W,T}) where {W,T} = SVec{W,T}(v)
@inline vbroadcast(::Val, b::Bool) = b
@generated function vbroadcast(::Type{_Vec{_W,Ptr{T}}}, s::Ptr{T}) where {_W, T}
    W = _W + 1
    typ = "i$(8sizeof(Int))"
    vtyp = "<$W x $typ>"
    instrs = String[]
    push!(instrs, "%ie = insertelement $vtyp undef, $typ %0, i32 0")
    push!(instrs, "%v = shufflevector $vtyp %ie, $vtyp undef, <$W x i32> zeroinitializer")
    push!(instrs, "ret $vtyp %v")
    quote
        $(Expr(:meta,:inline))
        llvmcall( $(join(instrs,"\n")), Vec{$W,Ptr{$T}}, Tuple{Ptr{$T}}, s )
    end
end
@generated function vbroadcast(::Type{_Vec{_W,T}}, s::T) where {_W, T <: Integer}
    W = _W + 1
    typ = "i$(8sizeof(T))"
    vtyp = "<$W x $typ>"
    instrs = String[]
    push!(instrs, "%ie = insertelement $vtyp undef, $typ %0, i32 0")
    push!(instrs, "%v = shufflevector $vtyp %ie, $vtyp undef, <$W x i32> zeroinitializer")
    push!(instrs, "ret $vtyp %v")
    quote
        $(Expr(:meta,:inline))
        llvmcall( $(join(instrs,"\n")), Vec{$W,$T}, Tuple{$T}, s )
    end
end
@generated function vbroadcast(::Type{_Vec{_W,T}}, s::T) where {_W, T <: Union{Float16,Float32,Float64}}
    W = _W + 1
    typ = llvmtype(T)
    vtyp = "<$W x $typ>"
    instrs = String[]
    push!(instrs, "%ie = insertelement $vtyp undef, $typ %0, i32 0")
    push!(instrs, "%v = shufflevector $vtyp %ie, $vtyp undef, <$W x i32> zeroinitializer")
    push!(instrs, "ret $vtyp %v")
    quote
        $(Expr(:meta,:inline))
        llvmcall( $(join(instrs,"\n")), Vec{$W,$T}, Tuple{$T}, s )
    end
end
@generated function vbroadcast(::Type{_Vec{_W,T}}, ptr::Ptr{T}) where {_W, T}
    W = _W + 1
    typ = llvmtype(T)
    ptyp = JuliaPointerType
    vtyp = "<$W x $typ>"
    instrs = String[]
    alignment = Base.datatype_alignment(T)
    push!(instrs, "%ptr = inttoptr $ptyp %0 to $typ*")
    push!(instrs, "%res = load $typ, $typ* %ptr, align $alignment")
    push!(instrs, "%ie = insertelement $vtyp undef, $typ %res, i32 0")
    push!(instrs, "%v = shufflevector $vtyp %ie, $vtyp undef, <$W x i32> zeroinitializer")
    push!(instrs, "ret $vtyp %v")
    quote
        $(Expr(:meta,:inline))
        llvmcall( $(join(instrs,"\n")), Vec{$W,$T}, Tuple{Ptr{$T}}, ptr )
    end
end
@generated function vzero(::Type{_Vec{_W,T}}) where {_W,T}
    W = _W + 1
    typ = llvmtype(T)
    instrs = "ret <$W x $typ> zeroinitializer"
    quote
        $(Expr(:meta,:inline))
        llvmcall($instrs, Vec{$W,$T}, Tuple{}, )
    end
end
@generated function vzero(::Val{W}, ::Type{T}) where {W,T}
    typ = llvmtype(T)
    instrs = "ret <$W x $typ> zeroinitializer"
    quote
        $(Expr(:meta,:inline))
        SVec(llvmcall($instrs, Vec{$W,$T}, Tuple{}, ))
    end
end

# @inline vzero(::Type{Vec{W,T}}) where {W,T} = vzero(Val{W}(), T)
@inline vbroadcast(::Val{W}, s::T) where {W,T} = SVec(vbroadcast(Vec{W,T}, s))
@inline vbroadcast(::Val{W}, ptr::Ptr{T}) where {W,T} = SVec(vbroadcast(Vec{W,T}, ptr))
@inline vbroadcast(::Type{_Vec{_W,T1}}, s) where {_W,T1} = vbroadcast(_Vec{_W,T1}, convert(T1,s))
@inline vbroadcast(::Type{_Vec{_W,T1}}, s::T2) where {_W,T1<:Integer,T2<:Integer} = vbroadcast(_Vec{_W,T1}, s % T1)
@inline vbroadcast(::Type{_Vec{_W,T}}, ptr::Ptr) where {_W,T} = vbroadcast(_Vec{_W,T}, Base.unsafe_convert(Ptr{T},ptr))
@inline vbroadcast(::Type{SVec{W,T}}, s) where {W,T} = SVec(vbroadcast(Vec{W,T}, s))
@inline vbroadcast(::Type{_Vec{_W,T}}, v::_Vec{_W,T}) where {_W,T} = v
@inline vbroadcast(::Type{SVec{W,T}}, v::SVec{W,T}) where {W,T} = v
@inline vbroadcast(::Type{SVec{W,T}}, v::Vec{W,T}) where {W,T} = SVec(v)

@inline vone(::Type{_Vec{_W,T}}) where {_W,T} = vbroadcast(_Vec{_W,T}, one(T))
# @inline vzero(::Type{Vec{W,T}}) where {W,T} = vbroadcast(Vec{W,T}, zero(T))
@inline vone(::Type{SVec{W,T}}) where {W,T} = SVec(vbroadcast(Vec{W,T}, one(T)))
@inline vzero(::Type{SVec{W,T}}) where {W,T} = SVec(vzero(Vec{W,T}))
@inline vone(::Type{T}) where {T} = one(T)
@inline vzero(::Type{T}) where {T<:Number} = zero(T)
@inline vzero() = SVec(vzero(pick_vector_width_val(Float64)))
@inline sveczero(::Type{T}) where {T} = Svec(vzero(pick_vector_width_val(T)))
@inline sveczero() = Svec(vzero(pick_vector_width_val(Float64)))
@inline VectorizationBase.SVec{W}(s::T) where {W,T<:Number} = SVec(vbroadcast(Vec{W,T}, s))
@inline VectorizationBase.SVec{W,T}(s::T) where {W,T<:Number} = SVec(vbroadcast(Vec{W,T}, s))
@inline VectorizationBase.SVec{W,T}(s::Number) where {W,T} = SVec(vbroadcast(Vec{W,T}, convert(T, s)))

@inline VectorizationBase.SVec{W}(v::Vec{W,T}) where {W,T<:Number} = SVec(v)
@inline VectorizationBase.SVec{W}(v::SVec{W,T}) where {W,T<:Number} = v

@inline SVec{W,T}(v::Vararg{T,W}) where {W,T} = @inbounds SVec{W,T}(ntuple(Val(W)) do w Core.VecElement(v[w]) end)
@inline SVec{W,T1}(v::Vararg{T2,W}) where {W,T1<:Number,T2<:Number} = @inbounds SVec{W,T1}(ntuple(Val(W)) do w Core.VecElement(convert(T1,v[w])) end)
@inline SVec{1,T1}(v::Vararg{T2,1}) where {T1<:Number,T2<:Number} = @inbounds SVec{1,T1}((Core.VecElement(convert(T1,v[1])),))

@inline Base.length(::AbstractStructVec{N}) where N = N
@inline Base.size(::AbstractStructVec{N}) where N = (N,)
@inline Base.eltype(::AbstractStructVec{N,T}) where {N,T} = T
@inline Base.conj(v::AbstractStructVec) = v # so that things like dot products work.
@inline Base.adjoint(v::AbstractStructVec) = v # so that things like dot products work.
@inline Base.transpose(v::AbstractStructVec) = v # so that things like dot products work.
@inline Base.getindex(v::SVec, i::Integer) = v.data[i].value

# @inline function SVec{N,T}(v::SVec{N,T2}) where {N,T,T2}
    # @inbounds SVec(ntuple(n -> Core.VecElement{T}(T(v[n])), Val(N)))
# end

@inline Base.one(::Type{<:AbstractStructVec{W,T}}) where {W,T} = SVec(vbroadcast(Vec{W,T}, one(T)))
@inline Base.one(::AbstractStructVec{W,T}) where {W,T} = SVec(vbroadcast(Vec{W,T}, one(T)))
@inline Base.zero(::Type{<:AbstractStructVec{W,T}}) where {W,T} = SVec(vbroadcast(Vec{W,T}, zero(T)))
@inline Base.zero(::AbstractStructVec{W,T}) where {W,T} = SVec(vbroadcast(Vec{W,T}, zero(T)))

# Use with care in function signatures; try to avoid the `T` to stay clean on Test.detect_unbound_args
const AbstractSIMDVector{W,T} = Union{Vec{W,T},AbstractStructVec{W,T}}

@inline extract_data(v) = v
@inline extract_data(v::SVec) = v.data
@inline extract_data(v::AbstractStructVec) = v.data
@inline extract_value(v::Vec, i) = v[i].value
@inline extract_value(v::SVec, i) = v.data[i].value

@inline firstval(x::Vec) = first(x).value
@inline firstval(x::SVec) = first(extract_data(x)).value
@inline firstval(x) = first(x)

function Base.show(io::IO, v::SVec{W,T}) where {W,T}
    print(io, "SVec{$W,$T}<")
    for w ∈ 1:W
        print(io, repr(v[w]))
        w < W && print(io, ", ")
    end
    print(io, ">")
end
Base.bitstring(m::Mask{W}) where {W} = bitstring(extract_data(m))[end-W+1:end]
function Base.show(io::IO, m::Mask{W}) where {W}
    bits = bitstring(m)
    bitv = split(bits, "")
    print(io, "Mask{$W,Bool}<")
    for w ∈ 0:W-1
        print(io, bitv[W-w])
        w < W-1 && print(io, ", ")
    end
    print(io, ">")
end

struct _MM{W,I<:Number}
    i::I
    @inline _MM{W}(i::T) where {W,T} = new{W,T}(i)
end

include("cartesianvindex.jl")
include("static.jl")
include("vectorizable.jl")
include("strideprodcsestridedpointers.jl")
@static if Sys.ARCH === :x86_64 || Sys.ARCH === :i686
    @static if Base.libllvm_version >= v"8" && LLVM_SHOULD_WORK
        include("cpu_info_x86_llvm.jl")
    else
        include("cpu_info_x86_cpuid.jl")
    end
else
    include("cpu_info_generic.jl")
end

include("vector_width.jl")
include("number_vectors.jl")
include("masks.jl")
include("alignment.jl")
include("precompile.jl")
_precompile_()

end # module


text:
function jetrgb(x::Float64)
  y = 4x
  r = min(y - 1.5, -y + 4.5)
  g = min(y - 0.5, -y + 3.5)
  b = min(y + 0.5, -y + 2.5)
  [clamp(r,0.,1.), clamp(g,0.,1.), clamp(b,0.,1.)]
end

function grayrgb(x::Float64)
  y = clamp(x,0.,1.)
  [y,y,y]
end

jetgray(x::Float64) = x <= 0.01 ? grayrgb(100x) : jetrgb(x)

function redblue(x::Float64)
  cmap = colormap("RdBu")
  idx = round(Int,99x + 1)
  cmap[idx]
end

function cubehelix(x::Float64; start=0.3, rot=-0.5, gamma=1.0,
                   minSat=1.2, maxSat=1.2)
  angle = 2pi*(start/3.0 + rot*x + 1.)
  x = x.^gamma
  satar = x.*(maxSat - minSat) + minSat
  amp = satar.*x.*(1. - x) / 2.0
  red = x + amp.*(-0.14861*cos(angle) + 1.78277*sin(angle))
  grn = x + amp.*(-0.29227*cos(angle) - 0.90649*sin(angle))
  blu = x + amp.*(1.97294*cos(angle))
  red = clamp(red, 0., 1.)
  grn = clamp(grn, 0., 1.)
  blu = clamp(blu, 0., 1.)
  [red, grn, blu]
end

function overlay(front::Array{Float64,2}, back::Array{Float64,2}, mask::Array{Bool,2})
  @assert(length(front) == length(back))
  z = similar(front)
  scale= maximum(front)
  back[:] = scale*(back[:] - minimum(back)) / (maximum(back) - minimum(back))
  for k in 1:length(mask)
    z[k] = mask[k] ? max(front[k], nextfloat(nextfloat(0.01scale))) : 0.01back[k]
  end
  z
end

function oplot2(front::Array{Float64,2}, back::Array{Float64,2}, mask::Array{Bool,2})
  @assert length(front) == length(back)
  # normalize the images
  r = maximum(back) - minimum(back)
  if r == 0.0
    back[:] = 0.0
  else
    @. back = (back - minimum(back)) / r
  end
  n,m = size(front)
  img = zeros(n,m,3)
  for k in 1:m, j in 1:n
    if mask[j,k]
      s = front[j,k]
      s = clamp(s, 0.0, 1.0)
      cmidx = 100 - round(Int,99.0*s)
      img[j,k,:] = jetrgb(s)
    else
      s = back[j,k]
      if !isnan(s)
        cmidx = 100 - round(Int,99.0*s)
        img[j,k,1] = s
        img[j,k,2] = s
        img[j,k,3] = s
      end
    end
  end
  img
end

function makeplots(mat::Dict; outdir::AbstractString="results")
  isdir(outdir) || mkdir(outdir)
  R1map = mat["R10"]
  S0map = mat["S0"]
  modelmap = mat["modelmap"]
  Ct = mat["Ct"]
  Kt = mat["Kt"]
  ve = mat["ve"]
  vp = mat["vp"]
  resid = mat["resid"]
  q = quantile(S0map[:], 0.99)
  S0map[S0map .> q] .= q
  back = @. (S0map - minimum(S0map)) / (maximum(S0map) - minimum(S0map))
  mask = convert(Array{Bool,2}, mat["mask"])

  figure(figsize=(4.5,4.5))
  clf()
  plot(mat["t"], mat["Cp"], "ko-")
  xlabel("time (min)")
  yticks(collect(0:5))
  ylabel("[Gd-DTPA] (mM)")
  title("arterial input function, \$C_p\$")
  savefig("$outdir/aif.pdf")

  figure()
  clf()
  imshow(mat["SER"], interpolation="nearest", cmap="cubehelix", vmin=0, vmax=10)
  colorbar()
  title("signal enhancement ratio")
  savefig("$outdir/ser.pdf")

  figure()
  clf()
  imshow(mask, interpolation="nearest", cmap="gray")
  title("mask")
  savefig("$outdir/mask.pdf")

  figure()
  clf()
  imshow(R1map, interpolation="nearest", cmap="cubehelix", vmin=0, vmax=5)
  title("\$R_1\$ relaxation rate (s\$^{-1}\$)")
  colorbar(ticks=collect(0:5))
  savefig("$outdir/R1.pdf")

  figure()
  clf()
  imshow(S0map, interpolation="nearest", cmap="gray")
  title("\$S_0\$")
  colorbar()
  savefig("$outdir/S0.pdf")

  figure()
  clf()
  Ct = dropdims(maximum(Ct; dims=1); dims=1)
  x = oplot2(clamp.(Ct, 0.0, 5.0), back, mask)
  imshow(x, interpolation="nearest", cmap="jet", vmin=0, vmax=5)
  title("max tissue conc., \$C_t\$ (mmol)")
  colorbar()
  savefig("$outdir/Ct.pdf")

  figure()
  clf()
  x = oplot2(clamp.(Kt, 0.0, 1.0), back, mask)
  imshow(x, interpolation="nearest", cmap="jet", vmin=0, vmax=1)
  title("\$K^\\mathrm{trans}\$ (min\$^{-1}\$)")
  colorbar(ticks=collect(0:2:10)/10.0)
  savefig("$outdir/Kt.pdf")

  figure()
  clf()
  x = oplot2(clamp.(ve, 0.0, 1.0), back, mask)
  imshow(x, interpolation="nearest", cmap="jet", vmin=0, vmax=1)
  title("\$v_e\$")
  colorbar(ticks=[0,0.2,0.4,0.6,0.8,1])
  savefig("$outdir/ve.pdf")

  figure()
  clf()
  x = oplot2(clamp.(vp, 0.0, 1.0), back, mask)
  imshow(x, interpolation="nearest", cmap="jet", vmin=0, vmax=1)
  title("\$v_p\$")
  colorbar(ticks=[0,0.2,0.4,0.6,0.8,1])
  savefig("$outdir/vp.pdf")

  figure()
  clf()
  kep = clamp.(Kt./ve, 0.0, 10.0)
  x = oplot2(kep, back, mask)
  imshow(x, interpolation="nearest", cmap="jet", vmin=0, vmax=10)
  title("\$k_{ep}\$")
  colorbar(ticks=[0,2,4,6,8,10])
  savefig("$outdir/kep.pdf")

  figure()
  clf()
  imshow(modelmap, interpolation="nearest", cmap="cubehelix")
  title("model used")
  colorbar(ticks=[0,1,2,3])
  savefig("$outdir/modelmap.pdf")

  figure()
  clf()
  x = oplot2(clamp.(100*resid, 0.0, 1.0), back, mask)
  imshow(x, interpolation="nearest", cmap="jet", vmin=0, vmax=1.0)
  title("residual \$\\times\$ 100")
  #colorbar(ticks=[0:2:10]/10000.0)
  colorbar()
  savefig("$outdir/resid.pdf")
end

function makeplots(matfile::AbstractString)
  println("creating plots from the file ", matfile)
  mat = read(matopen(matfile))
  makeplots(mat)
end


text:
# Equation 33 of https://www.mdpi.com/2227-7390/6/2/16/htm

# Loaded Packages
using FdeSolver
using Plots
using SpecialFunctions

# Parameters
tSpan = [0, 5]         # Time Span
y0 = [1, 0.5, 0.3]     # Initial values
β = [0.5, 0.2, 0.6]    # Order of derivation

# Definition of the System
function F(t, y)

    F1 = 1 / sqrt(pi) * (((y[2] - 0.5) * (y[3] - 0.3))^(1 / 6) + t^(1 / 2))
    F2 = gamma(2.2) * (y[1] - 1)
    F3 = gamma(2.8) / gamma(2.2) * (y[2] - 0.5)

    return [F1, F2, F3]

end

function JF(t,y)
    # System equation
    J11 = 0
    J12 = (y[2]-0.5).^(-5/6).*(y[3]-0.3).^(1/6)/6/sqrt(pi)
    J13 = (y[2]-0.5).^(1/6).*(y[3]-0.3).^(-5/6)/6/sqrt(pi)
    J21 =  gamma(2.2)
    J22 =  0
    J23 =  0
    J31 =  0
    J32 =  gamma(2.8)/gamma(2.2)
    J33 =  0

    J = [J11 J12 J13
         J21 J22 J23
         J31 J32 J33]
    return J
end

# Numerical Solution
t, Yapp = FDEsolver(F, tSpan, y0, β, nc = 5)
t_J, Yapp_J = FDEsolver(F, tSpan, y0, β, JF = JF)

# Plot
plot(t, Yapp, linewidth = 5, title = "Solution of system 33",
     xaxis = "Time (t)", yaxis = "y(t)", label = "Approximation")
# plot!(t_J, Yapp_J, linewidth = 5, ls= :dot, label = "Approximation with Jacob")
plot!(t, t -> (t .+ 1), lw = 3, ls = :dash, color= "red", label = "Exact solution")
plot!(t, t -> (t.^1.2 .+ 0.5), lw = 3, ls = :dash, color= "cyan", label = "Exact solution")
plot!(t, t -> (t.^1.8 .+ 0.3), lw = 3, ls = :dash, color= "black" ,label = "Exact solution")


text:
function check_for_focus(cur_item)

  has_focus = check_for_test_param(cur_item, "focus")

  has_focus

end


text:
# Submodule Julia HTTP Server
# HTTP and WebSockets need to be loaded in the calling context.
# LOAD_PATH must include logutils
# Intended for accepting echoing clients, such as ws_jce.jl, and
# running echo tests with that client.
# The server stays open until close_hts or the websocket is closed.
module ws_hts
using ..HTTP
import HTTP.Header
using ..WebSockets
# We want to log to a separate file, so
# we use our own instance of logutils_ws here.
import logutils_ws: logto, clog, zlog, zflush, clog_notime
const SRCPATH = Base.source_dir() == nothing ? Pkg.dir("WebSockets", "benchmark") :Base.source_dir()
const SERVEFILE = "bce.html"
const PORT = 8000
const SERVER = "127.0.0.1"
const WSMAXTIME = Base.Dates.Second(600)
const WEBSOCKET = Vector{WebSockets.WebSocket}()
const TCPREF = Ref{HTTP.Sockets.TCPServer}()
"Run asyncronously or in separate process"
function listen_hts()
    id = "listen_hts"
    try
        clog(id,"listen_hts starts on ", SERVER, ":", PORT)
        zflush()
        HTTP.listen(SERVER, UInt16(PORT), tcpref = TCPREF) do http
            if WebSockets.is_upgrade(http.message)
                acceptholdws(http)
                clog(id, "Websocket closed, server stays open until ws_hts.close_hts()")
            else
                HTTP.Servers.handle_request(handlerequest, http)
            end
        end
    catch err
        clog(id, :red, err)
        clog_notime.(catch_stacktrace()[1:4])
        zflush()
    end
end

"Accepts an incoming connection, upgrades to websocket,
and waits for timeout or a closed socket.
Also stops the server from accepting more connections on exit."
function acceptholdws(http)
    id = "ws_hts.acceptholdws"
    zlog(id);zflush()
    # If the ugrade is successful, just hold the reference and thread
    # of execution. Other tasks may do useful things with it.
    WebSockets.upgrade(http) do ws
        if length(WEBSOCKET) > 0
            # unexpected behaviour.
            if !isopen(WEBSOCKET[1])
                pop!(WEBSOCKET)
            else
                msg = " A websocket is already open. Not accepting the attempt at opening more."
                clog(id, :red, msg);zflush()
                return
            end
        end
        push!(WEBSOCKET, ws)
        zlog(id, ws);zflush()
        t1 = now() + WSMAXTIME
        while isopen(ws) && now() < t1
            yield()
        end
        length(WEBSOCKET) > 0 && pop!(WEBSOCKET)
        zlog(id, " exiting");zflush()
    end
end
"Returns a websocket or a string"
function getws_hts()
    id = "getws_hts"
    if length(WEBSOCKET) > 0
        if isopen(WEBSOCKET[1])
            zlog(id, " return reference")
            return WEBSOCKET[1]
        else
            msg = " Websocket is referred but not open. Acceptholdws might not have been scheduled after is was closed."
            clog(id, msg)
            zflush()
            return msg
        end
    else
        if !isdefined(TCPREF, :x)
            msg = " No server running yet, run ws_hts.listen_hts() or wait"
        else
            msg = " No websocket has connected yet, run ws_hce.echowithdelay_jce() or wait"
        end
        zlog(id, msg)
        return id * msg
    end
end

"HTTP request -> HTTP response."
function handlerequest(request::HTTP.Request)
    id = "handlerequest"
    zlog(id, request)
    response = responseskeleton(request)
    try
        if request.method == "GET"
            response = resp_HTTP(request.target, response)
        elseif request.method == "HEAD"
            response = resp_HTTP(request.target, response)
            response.body = Array{UInt8,1}()
        else
            response = HTTP.Response(501, "Not implemented method: $(request.method), fix $id")
        end
    end
    zlog(id, response)
    response
end

"""
Tell browser about the methods this server supports.
"""
function responseskeleton(request::HTTP.Request)
    r = HTTP.Response()
    HTTP.Messages.appendheader(r, Header("Allow" => "GET,HEAD"))
    HTTP.Messages.appendheader(r, Header("Connection" => "close"))
    r #
end

"request.target -> HTTP.Response , building on a skeleton response"
function resp_HTTP(resource::String, resp::HTTP.Response)
    id = "resp_HTTP"
    if resource == "/favicon.ico"
        s = read(joinpath(SRCPATH, "favicon.ico"))
        push!(resp.headers, "Content-Type" => "image/x-icon")
    else
        s = read(joinpath(SRCPATH, SERVEFILE))
        push!(resp.headers, "Content-Type" => "text/html")
    end
    resp.body = s
    push!(resp.headers, "Content-Length" => string(length(s)))
    resp.status = 200
    resp
end

"Close the websocket, stop the server. TODO improve"
function close_hts()
    clog("ws_hts.close_hts")
    length(WEBSOCKET) >0 && isopen(WEBSOCKET[1]) && close(WEBSOCKET[1]) && sleep(0.5)
    isassigned(TCPREF) && close(TCPREF.x)
end

end # module
"""
For debugging:

import HTTP
using WebSockets
const SRCPATH = Base.source_dir() == nothing ? Pkg.dir("WebSockets", "benchmark") :Base.source_dir()
const LOGGINGPATH = realpath(joinpath(SRCPATH, "../logutils/"))
# for finding local modules
SRCPATH ∉ LOAD_PATH && push!(LOAD_PATH, SRCPATH)
LOGGINGPATH ∉ LOAD_PATH && push!(LOAD_PATH, LOGGINGPATH)
import ws_hts.listen_hts
tas = @schedule ws_hts.listen_hts()
sleep(7)
hts = ws_hts.getws_hts()
"""


text:
function _initialize_preamble()
    pre = config(_DEFAULT_PREAMBLE_KEY; info=false)
    if pre === nothing
        _DEFAULT_PREAMBLE[] = [
            "\\usepackage{amsmath,amsthm,amssymb}",
            "\\usepackage{color}"
        ]
    else
        _DEFAULT_PREAMBLE[] = pre
    end
    return nothing
end

function _initialize_latexengine()
    eng = config(_DEFAULT_ENGINE_KEY; info=false)
    if eng !== nothing
        if eng == "PDFLaTeX"
            texengine!(PDFLaTeX())
        elseif eng == "XeLaTeX"
            texengine!(XeLaTeX())
        elseif eng == "LuaLaTeX"
            texengine!(LuaLaTeX())
        else
            @warn "Unknown LaTeX engine $eng; please set a valid default LaTeX engine with `config!`."
            return nothing
        end
    else
        for (engname, ENG) in zip(("xelatex", "pdflatex", "lualatex"), (XeLaTeX(), PDFLaTeX(), LuaLaTeX()))
            @debug "Looking for $engname"
            if Sys.which(engname) !== nothing
                texengine!(ENG)
                @info "Using $engname. Use `texengine!` to change the LaTeX engine for this session, or `config!` to set the default permanently (which will also suppress this message.)"
                return nothing
            end
        end
        @warn "No LaTeX engine is found. Make sure that LaTeX is properly installed on your system, and `xelatex`, `pdflatex`, or `lualatex` are available in your PATH."
        return nothing
    end
end

function _initialize_latex()
    _initialize_latexengine()
    if Sys.which("dvisvgm") === nothing
        @warn "dvisvgm not found. Your LaTeX installation should come with dvisvgm; make sure it is in your PATH, since this package won't function without it."
    end
    _initialize_preamble()
end


text:
using ModelingToolkit, NonlinearSolve

@parameters t
@variables u1(t) u2(t) u3(t) u4(t) u5(t)

eqs = [
        0 ~ u1 - sin(u5),
        0 ~ u2 - cos(u1),
        0 ~ u3 - hypot(u1, u2),
        0 ~ u4 - hypot(u2, u3),
        0 ~ u5 - hypot(u4, u1)
]

sys = NonlinearSystem(eqs, [u1, u2, u3, u4, u5], [])

simple_sys = structural_simplify(sys)

equations(simple_sys)

prob = NonlinearProblem(simple_sys, [u5=>0.0])
sol = solve(prob, NewtonRaphson())

sol[u5]
sol[u1]

# stochastic differential equations
using ModelingToolkit, DifferentialEquations

@parameters t σ ρ β
@variables x(t) y(t) z(t)
D = Differential(t)

eqs = [
        D(x) ~ σ*(y-x),
        D(y) ~ x*(ρ-z)-y,
        D(z) ~ x*y - β*z
       ]

noise_eqs = [0.1*x, 0.1*y, 0.1*z]

de = SDESystem(eqs, noise_eqs, t, [x, y, z], [σ, ρ, β])

u0map = [x=> 1.0, y=> 0.0, z=>0.0]
parammap = [σ => 10.0, ρ => 28.0, β => 8/3]

prob = SDEProblem(de, u0map, (0.0, 100.0), parammap)

sol = solve(prob)
plot(sol, vars=(x,y,z))

# ------------ optimization problem ----------------------
using ModelingToolkit, GalacticOptim, Optim

@variables x y
@parameters a b
loss = (a-x)^2 + b * (y -x^2)^2

sys = OptimizationSystem(loss, [x,y], [a,b])

# initial guess
u0 = [x=>1.0, y=>2.0]
p = [a => 6.0, b=> 7.0]

prob = OptimizationProblem(sys, u0, grad=true, hess=true)
sol = solve(prob) # doesn't work yet

# ------------ dsl demo via Catalyst ---------------------
using DifferentialEquations, Catalyst, Latexify, ModelingToolkit, Plots

repressilator = @reaction_network begin
    hillr(P₃,α,K,n), ∅ --> m₁
    hillr(P₁,α,K,n), ∅ --> m₂
    hillr(P₂,α,K,n), ∅ --> m₃
    (δ,γ), m₁ ↔ ∅
    (δ,γ), m₂ ↔ ∅
    (δ,γ), m₃ ↔ ∅
    β, m₁ --> m₁ + P₁
    β, m₂ --> m₂ + P₂
    β, m₃ --> m₃ + P₃
    μ, P₁ --> ∅
    μ, P₂ --> ∅
    μ, P₃ --> ∅
end α K n δ γ β μ;


text:
using InventoryManagement
using Test

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


text:
module CharDisplay

using BaremetalPi
#import BaremetalPi.gpio_set, BaremetalPi.gpio_clear

include("./display.jl")
include("./gpio.jl")

end # module


text:
using PyPlot
using Roots

ntime = 100
time = collect(linspace(0,10,ntime))
gaussian=exp(-time.^2/2.)

function psd_func(x,omega)
psd = 0.
for k=1:nexp
  a1 = x[1+(k-1)*4]
  a2 = x[2+(k-1)*4]
  b  = x[3+(k-1)*4]
  om = x[4+(k-1)*4]
  psdt = (a2*om*(b^2 + om^2 - omega^2) + a1*b*(b^2 + om^2 + omega^2))/
         (b^4 + (om^2 - omega^2)^2 + 2*b^2*(om^2 + omega^2))
  psd += psdt
end
return psd
end

nexp = 4
x = [-0.5247762751816804,0.2828915570634826,1.3837656415003445,1.8323467099030555, 63.120917132667905,0.0,1.478130410600365,0.0, -32.17183668366479,0.0,1.5929285243364661,0.0, -29.423970266940803,0.0,1.38693874263686,0.0]

nfreq = 1000
omega = logspace(-1,1,nfreq)
psd = zeros(nfreq)
for i=1:nfreq
  psd[i]=psd_func(x,omega[i])
end

# model = gauss_model(atrial)
function dpsd_domega(omega)
# Let the kernel components be:
# a1*Exp[-b*t]*Cos[omega_0*t] + a2*Exp[-b*t]*Sin[omega_0*t]
# x = [{a1,a2,b,omega_0}_1; {a1,a2,b,omega_0}_2; {a1,a2,b,omega_0}_3; ... ]
#
dpsd = zero(Real)
for i=1:nexp
  a1 = x[(i-1)*4+1]
  a2 = x[(i-1)*4+2]
  b  = x[(i-1)*4+3]
  om0 = x[(i-1)*4+4]
  dpsd += (-2*sqrt(2/pi)*omega*(a1*b*(b^4 - 3*om0^4 + 2*om0^2*omega^2 + omega^4 - 2*b^2*(om0^2 - omega^2)) -
     a2*om0*(-3*b^4 + (om0^2 - omega^2)^2 - 2*b^2*(om0^2 + omega^2))))/
     ((b^2 + (om0 - omega)^2)^2*(b^2 + (om0 + omega)^2)^2)
end
return dpsd
end

# Now, solve for roots of PSD:

#omega_min = newton(dpsd_domega,5.0)
omega_min = fzero(dpsd_domega,[1e-4,1e4])
fmin = psd_func(x,omega_min)

plot(omega,psd)
scatter(omega_min,fmin)

max_stars_repo_name: dfm/ess

max_issues_repo_name: dfm/ess

max_forks_repo_name: dfm/ess


text:
# Test library utility functions
function t_libutils()

	# Test confusion matrix

	references = ["a","a","c","c","b","b"];
	predictions = ["a","b","c","c","a","a"];

	C = j4pr.confusionmatrix(predictions,references)
	Cref = [1.0  2.0  0.0; 1.0  0.0  0.0; 0.0  0.0  2.0]
	Test.@test C == Cref

	Cn = j4pr.confusionmatrix(predictions,references;normalize=true)
	Cref = [0.5  1.0  0.0; 0.5  0.0  0.0; 0.0  0.0  1.0]
	Test.@test Cn == Cref

	Cp = j4pr.confusionmatrix(predictions,references;positive="a")
	Cref = [1.0 2.0; 1.0  2.0]
	Test.@test Cp == Cref

	Cpn = j4pr.confusionmatrix(predictions,references;normalize=true,positive="a")
	Cref = [0.5  0.5; 0.5  0.5]
	Test.@test Cpn == Cref

	# the 'showmatrix' option is not tested.
	#confusionmatrix(predictions,references;normalize=true,positive="a",showmatrix=true);

	Test.@test sum(j4pr.pintgen(23,[0.3,0.7])) == 23
end


text:
# This script validate the linearized tofts models

using MAT
using DCEMRI

test_folder = @__DIR__
phantom_ids = (4, 6) # QIBA Phantom #'s
conditions = (:noiseless, :noisy)

# Create noisy data
DCEMRI.makeQibaNoisy(4)
DCEMRI.makeQibaNoisy(6)

for phantom_id in phantom_ids, condition in conditions
  if condition == :noiseless
    mat_filename = "qiba$(phantom_id).mat"
    results_foldername = "results_linear"
  elseif condition == :noisy
    mat_filename = "qiba$(phantom_id)noisy.mat"
    results_foldername = "results_linear_noisy"
  else
    error("Unknown condition: $(condition)")
  end
  qiba_file = joinpath(test_folder, "q$(phantom_id)", mat_filename)
  results_folder = joinpath(test_folder, "q$(phantom_id)", results_foldername)
  results_file = joinpath(results_folder, "results.mat")
  if !isdir(results_folder)
    mkdir(results_folder)
  end
  # The phantom contains 10x10 blocks. Noiseless phantom gobbles them into 1 block.
  voxel_reduction_factor = (condition == :noiseless) ? 10 : 1
  # [5] -> Linearized extended tofts model, [4] -> Linearized tofts model
  models = (phantom_id == 4) ? [5] : [4]

  results = fitdata(datafile=qiba_file, outfile=results_file, models=models)
  DCEMRI.analyze(phantom_id, results, results_folder; dx=voxel_reduction_factor)
end


text:
function run(s)
    lines = split(s, "\n")
    timestamp = parse(Int, lines[1])
    minutes = [parse(Int, el) for el in split(lines[2], ",") if el != "x"]
    chosen::Int = -1
    waiting_time::Int = -1
    for m in minutes
        wt = m - timestamp % m
        if (wt < waiting_time) || (waiting_time < 0)
            waiting_time = wt
            chosen = m
        end
    end
    return waiting_time * chosen
end

#########################################

function main()
    res, time, memory = @timed run(ARGS[1])
    println("_duration:$(time * 1000)")
    println(res)
end

main()


text:
mutable struct EvoDuplex <: AbstractDuplex
   duplex::RNADuplex
   alignment::Array{DNA,2}
   first::Int
   bracket::Vector{Char}
   score::Float64

   function EvoDuplex( rdup::RNADuplex, tree::PhyloTree, dsa::RNADuplexArray, fwd_index, rev_index )
      alignment = Array{DNA,2}(length(tree.order), sum(lengths(rdup.path)))
      first   = 1
      last    = size(alignment, 2)
      bracket = Vector{Char}(size(alignment, 2))
      i,j     = 0,0
      for k in 1:length(rdup.path)
         if isbulge(rdup[k]) && isfiveprime(rdup[k])
            i += 1
            col = DNA[dsa.fwd.species[s][i][fwd_index] for s in 1:length(tree.order)-1]
            unshift!(col, dsa.fwd.depth[i][fwd_index])
            alignment[:,first] = col
            bracket[first] = '.'
            first += 1
         elseif isbulge(rdup[k]) && !isfiveprime(rdup[k])
            j += 1
            col = DNA[dsa.rev.species[s][j][rev_index] for s in 1:length(tree.order)-1]
            unshift!(col, dsa.rev.depth[j][rev_index])
            alignment[:,last] = col
            bracket[last] = '.'
            last -= 1
          else
            i += 1
            j += 1
            const lcol = DNA[dsa.fwd.species[s][i][fwd_index] for s in 1:length(tree.order)-1]
            const rcol = DNA[dsa.rev.species[s][j][rev_index] for s in 1:length(tree.order)-1]
            unshift!(lcol, dsa.fwd.depth[i][fwd_index])
            unshift!(rcol, dsa.rev.depth[j][rev_index])
            #println("dsa.rev.depth $j $rev_index")
            alignment[:,first] = lcol
            alignment[:,last]  = rcol
            bracket[first] = '('
            bracket[last]  = ')'
            first += 1
            last  -= 1
         end
      end
      return new( deepcopy(rdup), alignment, first-1, bracket, 0.0 )
   end

   function EvoDuplex( rdup::RNADuplex, single::PhyloTree, paired::PhyloTree, dsa::RNADuplexArray, fwd_index, rev_index )
      evo = EvoDuplex( rdup, single, dsa, fwd_index, rev_index )
      score!( evo, single, paired )
      return evo
   end
end

path( evo::EvoDuplex )        = evo.duplex.path

energy( evo::EvoDuplex )      = energy( evo.duplex )
energies( evo::EvoDuplex )    = map( energy, [evo.duplex; evo.species] )

npairs( evo::EvoDuplex )      = npairs( evo.duplex )
nmismatches( evo::EvoDuplex ) = nmismatches( evo.duplex )
nbulges( evo::EvoDuplex )     = nbulges( evo.duplex )

strings( evo::EvoDuplex )     = strings( evo.duplex )
brackets( evo::EvoDuplex )    = brackets( evo.duplex )

Base.push!{NP <: NucleotidePair}( evo::EvoDuplex, pair::NP )         = push!( evo.duplex, pair )
Base.push!{NP <: NucleotidePair}( evo::EvoDuplex, path::Vector{NP} ) = push!( evo.duplex, path )

function join_duplex!( left::EvoDuplex, right::EvoDuplex, npairs, npairs_first, npairs_last )
   push!( left.duplex, path(right.duplex)[npairs+1:end] )
   last  = (size(left.alignment, 2) - left.first) - npairs_last - 1
   first = left.first - npairs_first
#   try
   left.alignment = hcat(left.alignment[:,1:first], right.alignment, left.alignment[:,(end-last:end)])
   left.bracket = vcat(left.bracket[1:first], right.bracket, left.bracket[end-last:end])
   left.first = right.first + first
#   catch
#      println(STDERR, "$left\n$right\n$npairs\n$npairs_first\n$npairs_last")
#   end
   left
end

mutable struct EvoScore
   unstr_cons::Float64
   struc_neut::Float64
   struc_cons::Float64
   mean_covar::Float64
   mean_cons::Float64
end

score( evo::EvoDuplex ) = evo.score

function score!( evo::EvoDuplex, single::PhyloTree, paired::PhyloTree )
   uns_like = 0.0
   str_like = 0.0
   first   = 1
   last    = size(evo.alignment, 2)
   for k in 1:length(evo.duplex.path)
      if isbulge(evo.duplex[k]) && isfiveprime(evo.duplex[k])
         like = likelihood( single, evo.alignment[:,first] )
         str_like += log2(like)
         uns_like += log2(like)
         first += 1
      elseif isbulge(evo.duplex[k]) && !isfiveprime(evo.duplex[k])
         like = likelihood( single, evo.alignment[:,last] )
         str_like += log2(like)
         uns_like += log2(like)
         last -= 1
      else
         str_like += log2(likelihood(paired, evo.alignment[:,first], evo.alignment[:,last]))
         uns_like += log2(likelihood(single, evo.alignment[:,first])) + log2(likelihood(single, evo.alignment[:,last]))
         first += 1
         last  -= 1
      end
   end
   evo.score = str_like - uns_like
   str_like
end

function score( evo::EvoDuplex, single::PhyloTree; gapdenom=1.0 )
   single_like = 0.0
   for k in 1:size(evo.alignment,2)
      single_like += log2(likelihood( single, evo.alignment[:,k], gapdenom=gapdenom ))
   end
   single_like
end

function pscore!( evo::EvoDuplex, single::PhyloTree, paired::PhyloTree )
   uns_like = 0.0
   str_like = 0.0
   first   = 1
   last    = size(evo.alignment, 2)
   for k in 1:length(evo.duplex.path)
      if isbulge(evo.duplex[k]) && isfiveprime(evo.duplex[k])
         like = likelihood( single, evo.alignment[:,first] )
         str_like += log2(like)
         uns_like += log2(like)
         first += 1
      elseif isbulge(evo.duplex[k]) && !isfiveprime(evo.duplex[k])
         like = likelihood( single, evo.alignment[:,last] )
         str_like += log2(like)
         uns_like += log2(like)
         last -= 1
      else
         str_like += log2(likelihood(paired, evo.alignment[:,first], evo.alignment[:,last]))
         uns_like += log2(likelihood(single, evo.alignment[:,first])) + log2(likelihood(single, evo.alignment[:,last]))
         first += 1
         last  -= 1
      end
   end
   evo.score = str_like - uns_like
   evo.score
end

function covariance( evo::EvoDuplex )
   first   = 1
   last    = size(evo.alignment, 2)
   vals    = 0.0
   npairs  = 0
   for k in 1:length(evo.duplex.path)
      if isbulge(evo.duplex[k]) && isfiveprime(evo.duplex[k])
         first += 1
      elseif isbulge(evo.duplex[k]) && !isfiveprime(evo.duplex[k])
         last -= 1
      else
         cov = covariance_score(evo.alignment[:,first], evo.alignment[:,last])
         #print(" $cov ")
         vals   += cov
         first  += 1
         last   -= 1
         npairs += 1
      end
   end
   #println( " $maxval \n" )
   vals / npairs
end

const PI_MATRIX = [0.0 0.0 0.0 1.0
                   0.0 0.0 1.0 0.0
                   0.0 1.0 0.0 1.0
                   1.0 0.0 1.0 0.0]

delta( a::DNA, b::DNA ) = a == b ? 1 : 0

function pi_matrix( a::DNA, b::DNA )
   if isgap( a ) || isgap( b )
      return 0.0
   else
      const ind_a = trailing_zeros(reinterpret(UInt8, a)) + 1
      const ind_b = trailing_zeros(reinterpret(UInt8, b)) + 1
      return PI_MATRIX[ ind_a, ind_b ]
   end
end

# Calculated as described by Hofacker et al. Journal of Molecular Biology 2002.
# https://www.ncbi.nlm.nih.gov/pubmed/12079347
function covariance_score( a::Vector{DNA}, b::Vector{DNA} )
   denom = binomial(length(a), 2)
   numer = 0.0
   pairs = 0
   for i in 1:length(a)-1
      pairs += pi_matrix( a[i], b[i] )
      for j in i+1:length(a)
         numer += (2 - delta(a[i], a[j]) - delta(b[i], b[j])) * pi_matrix( a[i], b[i] ) * pi_matrix( a[j], b[j] )
      end
   end
   pairs += pi_matrix( a[end], b[end] )
   return (numer / denom) * (pairs / length(a))
end


text:
using MATLAB

check(a::Basic{ClosedNest}) = checkbasic(a.path, ["nest", "feeder"])
check(a::Basic{TransferNest}) = checkbasic(a.path, ["south", "north", "feeder", "originalnest"])
check(a::Basic{Transfer}) = checkbasic(a.path, ["south", "north", "feeder"])
check(a::Basic{TransferNestBelen}) = checkbasic(a.path, ["southbefore", "northbefore", "feederbefore", "nestbefore", "north", "south"])
check(a::Basic{DawaySandpaper}) = checkbasic(a.path, ["nest", "feeder", "rightdowninitial", "leftdowninitial", "rightupinitial", "leftupinitial", "rightdownfinal", "leftdownfinal", "rightupfinal", "leftupfinal"])
check(a::Basic{Daway}) = checkbasic(a.path, ["nest", "feeder", "initialfeeder"])

check!(x, a::MATcalib) = _check!(x, a.file, ["cameraParams"])
check!(x, a::Calib) = _check!(x, a.file, ["fc", "cc", "alpha_c", "kc", "Tc_ext", "Rc_ext"])
check!(x, a::Calibration) = nothing

function _check!(x, filename::String, fields::Vector{String})
    file = MatFile(filename, "r")
    # println(filename)
    try
        variables = variable_names(file)
        for field in fields
            x["no $field in the calibration file"] = field ∉ variables
        end
        close(file)
    catch
        x["unknown error in the calibration file"] = true
    end
end

function checkbasic(path::String, fields::Vector{String})
    push!(fields, "track")
    @assert isdir(path) "Path does not exists"
    x = Dict{String, Bool}()
    x["no calibrations file"] = !isfile(joinpath(path, "calibrations.csv"))
    x["no factors file"] = !isfile(joinpath(path, "factors.csv"))
    x["no columns file"] = !isfile(joinpath(path, "columns.csv"))
    if !x["no calibrations file"]
        res2calib = CSV.read(joinpath(path, "calibrations.csv"), header=["resfile", "calibfile"], types=[String, String], transforms = Dict(1 => strip, 2 => strip))
        res2poi = CSV.read(joinpath(path, "columns.csv"), header=["resfile", "columnnumber", "poi"], types=[String, Int, String], transforms = Dict(3 => strip, 1 => strip))
        todo = join(res2calib, res2poi, on = :resfile)
        for f in fields
            x["no $f in columns file"] = f ∉ todo[:poi]
            x["no column for $f in columns file"] = ismissing(todo[todo[:poi] .== f, :columnnumber])
            x["no resfile for $f in columns file"] = ismissing(todo[todo[:poi] .== f, :resfile])
            x["no calibfile for $f in columns file"] = ismissing(todo[todo[:poi] .== f, :calibfile])
        end
        resfiles = unique(todo[:resfile])
        for file in resfiles
            x["no $file"] = !isfile(joinpath(path, file))
        end
        calibfiles = unique(todo[:calibfile])
        for file in calibfiles
            x["no $file"] = !isfile(joinpath(path, file))
            calib = get_calibration(path, file)
            check!(x, calib)
        end
        for r in eachrow(todo)
            if r[:poi] ∈ [fields; "pellet"]
                file = MatFile(joinpath(path, r[:resfile]), "r")
                X = get_variable(file, "xdata")
                XI = nonzeros(X[:, r[:columnnumber]])
                x["empty column $(r[:columnnumber]) for $(r[:poi]) in $(r[:resfile])"] = isempty(XI)
                close(file)
            end
        end
    end
    if !x["no factors file"]
        factors = CSV.read(joinpath(path, factor_filename), header=["Factor", "Level"], types=[String, Union{String, Missing}], transforms = Dict(1 => strip))
        i = findfirst(factors[:Factor], "check")
        if i ≠ 0
            x["check corrupt in factors file"] = !r"^(?:\d+(\.\d+)?)\s(?:mm|cm|m)$"(strip(factors[i, :Level]))
        end
        i = findfirst(factors[:Factor], "azimuth")
        if i ≠ 0
            x["azimuth corrupt in factors file"] = !r"^(?:\d+(\.\d+)?)\s(?:rad|°)$"(strip(factors[i, :Level]))
        end
        i = findfirst(factors[:Factor], "nest2feeder")
        if i ≠ 0
            x["nest2feeder corrupt in factors file"] = !r"^(?:\d+(\.\d+)?)\s(?:mm|cm|m)$"(strip(factors[i, :Level]))
        end
    end
    [k for (k,v) in x if v]
end


text:
using AmplNLWriter, Compat, FactCheck, JuMP
using Base.Test

include("nl_convert.jl")
include("nl_linearity.jl")
include("nl_write.jl")

solver = JuMP.UnsetSolver()
solvers = Any[]
push!(solvers, IpoptNLSolver(["print_level=0"]))
push!(solvers, BonminNLSolver(["bonmin.nlp_log_level=0";
                               "bonmin.bb_log_level=0"]))
push!(solvers, CouenneNLSolver(["bonmin.nlp_log_level=0";
                                "bonmin.bb_log_level=0"]))
# Add tolerance option for couenne
open("couenne.opt", "w") do f
  write(f, "allowable_fraction_gap 0.0001")
end

examples_path = joinpath(dirname(dirname(@__FILE__)), "examples")
for solver in solvers
    solvername = getsolvername(solver)
    facts("[examples] test solver $solvername") do
        for example in ["jump_nltrig.jl", "jump_nlexpr.jl"]
            include(joinpath(examples_path, example))
        end
        if solvername != "ipopt"
            for example in ["jump_pruning.jl", "jump_minlp.jl",
                            "jump_nonlinearbinary.jl"]
                include(joinpath(examples_path, example))
            end
        end
        if solvername == "ipopt"
            for example in ["jump_maxmin.jl"]
                include(joinpath(examples_path, example))
            end
        end
    end
end

include(Pkg.dir("JuMP","test","solvers.jl"))
include(Pkg.dir("JuMP","test","nonlinear.jl"))

rm("couenne.opt")

FactCheck.exitstatus()


text:
include("core.jl")
include("conversion.jl")
include("colordiff.jl")
include("din99.jl")
include("layout.jl")


text:
import Base.getindex, Base.setindex!

immutable ETContainer{T,R} <: VectorizedExpression{T,R}
    data::Array{T, R}
end

function ETContainer{T,R}(A::VectorizedExpression{T,R})
    dat = Array(T, size(A)...)
    L = length(dat)
    block_lim = 8*Int64(floor(L/8))
    for i = 1:8:block_lim
        @inbounds dat[i] = A[i]
        @inbounds dat[i+1] = A[i+1]
        @inbounds dat[i+2] = A[i+2]
        @inbounds dat[i+3] = A[i+3]
        @inbounds dat[i+4] = A[i+4]
        @inbounds dat[i+5] = A[i+5]
        @inbounds dat[i+6] = A[i+6]
        @inbounds dat[i+7] = A[i+7]
    end
    for i = (block_lim+1):L
        @inbounds dat[i] = A[i]
    end
    return ETContainer(dat)
end

@inline function getindex{T,R}(A::ETContainer{T,R}, i...)
    @inbounds return A.data[i...]
end

@inline function setindex!{T,R}(A::ETContainer{T,R}, x::T, i...)
    A.data[i...] = x
end

@inline function length{T,R}(A::ETContainer{T,R})
    return length(A.data)
end

@inline function size{T,R}(A::ETContainer{T,R})
    return size(A.data)
end

export ETContainer


text:
import SpecialFunctions
import SpecialFunctions: besselj

@testset "Bessel Roots" begin
    @test_throws DomainError besselroots(0.0, -1)

    tol = 1e-11

    # Check if besselj(nu, besselroots(nu, n) ) is small:

    for nu = 0.:0.1:5.
        n = 10
        @test norm( besselj.(nu, besselroots(nu, n) ), Inf ) < tol
    end
end


text:
module TrafficAssignment

using LightGraphs, Optim, BinDeps, DataFrames, OrderedCollections
using Distributed, Printf, LinearAlgebra, SparseArrays

TNTP_SHA = "f730be5e3366e910bb7e9ada4665d32e9cbc219b"

include("load_network.jl")
include("frank_wolfe.jl")

export
        load_ta_network, download_tntp, read_ta_network,
        summarize_ta_data, read_ta_summary, net_dataframe,
        ta_frank_wolfe,
        TA_Data

end # module