text:
using DiffEqSensitivity, OrdinaryDiffEq, ForwardDiff, Calculus
using Test
function fb(du,u,p,t)
  du[1] = dx = p[1]*u[1] - p[2]*u[1]*u[2]
  du[2] = dy = -t*p[3]*u[2] + t*u[1]*u[2]
end
function jac(J,u,p,t)
  (x, y, a, b, c) = (u[1], u[2], p[1], p[2], p[3])
  J[1,1] = a + y * b * -1
  J[2,1] = t * y
  J[1,2] = b * x * -1
  J[2,2] = t * c * -1 + t * x
end
function paramjac(pJ,u,p,t)
  (x, y, a, b, c) = (u[1], u[2], p[1], p[2], p[3])
  pJ[1,1] = x
  pJ[2,1] = 0.0
  pJ[1,2] = - x * y
  pJ[2,2] = 0.0
  pJ[1,3] = 0.0
  pJ[2,3] = - t * y
end

f = ODEFunction(fb,jac=jac,paramjac=paramjac)
p = [1.5,1.0,3.0]
prob = ODEForwardSensitivityProblem(f,[1.0;1.0],(0.0,10.0),p)
probInpl = ODEForwardSensitivityProblem(fb,[1.0;1.0],(0.0,10.0),p)
probnoad = ODEForwardSensitivityProblem(fb,[1.0;1.0],(0.0,10.0),p,
                                        ForwardSensitivity(autodiff=false))
probnoad2 = ODEForwardSensitivityProblem(f,[1.0;1.0],(0.0,10.0),p,
                                        ForwardSensitivity(autodiff=false))
probvecmat = ODEForwardSensitivityProblem(fb,[1.0;1.0],(0.0,10.0),p,
                                          ForwardSensitivity(autojacvec=false,autojacmat=true))
sol = solve(prob,Tsit5(),abstol=1e-14,reltol=1e-14)
@test_broken solve(probInpl,KenCarp4(),abstol=1e-14,reltol=1e-14).retcode == :Success
solInpl = solve(probInpl,KenCarp4(autodiff=false),abstol=1e-14,reltol=1e-14)
solInpl2 = solve(probInpl,Rodas4(autodiff=false),abstol=1e-10,reltol=1e-10)
solnoad = solve(probnoad,KenCarp4(autodiff=false),abstol=1e-14,reltol=1e-14)
solnoad2 = solve(probnoad,KenCarp4(autodiff=false),abstol=1e-14,reltol=1e-14)
solvecmat = solve(probvecmat,Tsit5(),abstol=1e-14,reltol=1e-14)

x = sol[1:sol.prob.f.numindvar,:]

@test sol(5.0) ≈ solnoad(5.0)
@test sol(5.0) ≈ solnoad2(5.0)
@test sol(5.0) ≈ solInpl(5.0)
@test isapprox(solInpl(5.0), solInpl2(5.0),rtol=1e-5)
@test sol(5.0) ≈ solvecmat(5.0)

# Get the sensitivities

da = sol[sol.prob.f.numindvar+1:sol.prob.f.numindvar*2,:]
db = sol[sol.prob.f.numindvar*2+1:sol.prob.f.numindvar*3,:]
dc = sol[sol.prob.f.numindvar*3+1:sol.prob.f.numindvar*4,:]

sense_res1 = [da[:,end] db[:,end] dc[:,end]]

prob = ODEForwardSensitivityProblem(f.f,[1.0;1.0],(0.0,10.0),p,
                                    ForwardSensitivity(autojacvec=true))
sol = solve(prob,Tsit5(),abstol=1e-14,reltol=1e-14,saveat=0.01)
x = sol[1:sol.prob.f.numindvar,:]

# Get the sensitivities

res = sol[1:sol.prob.f.numindvar,:]
da = sol[sol.prob.f.numindvar+1:sol.prob.f.numindvar*2,:]
db = sol[sol.prob.f.numindvar*2+1:sol.prob.f.numindvar*3,:]
dc = sol[sol.prob.f.numindvar*3+1:sol.prob.f.numindvar*4,:]

sense_res2 = [da[:,end] db[:,end] dc[:,end]]

function test_f(p)
  prob = ODEProblem(f,eltype(p).([1.0,1.0]),(0.0,10.0),p)
  solve(prob,Tsit5(),abstol=1e-14,reltol=1e-14,save_everystep=false)[end]
end

p = [1.5,1.0,3.0]
fd_res = ForwardDiff.jacobian(test_f,p)
calc_res = Calculus.finite_difference_jacobian(test_f,p)

@test sense_res1 ≈ sense_res2 ≈ fd_res
@test sense_res1 ≈ sense_res2 ≈ calc_res

## Check extraction

xall, dpall = extract_local_sensitivities(sol)
@test xall == res
@test dpall[1] == da

_,dpall_matrix = extract_local_sensitivities(sol,Val(true))
@test mapreduce(x->x[:, 2], hcat, dpall) == dpall_matrix[2]

x, dp = extract_local_sensitivities(sol,length(sol.t))
sense_res2 = reduce(hcat,dp)
@test sense_res1 == sense_res2

@test extract_local_sensitivities(sol,sol.t[3]) == extract_local_sensitivities(sol,3)

tmp = similar(sol[1])
@test extract_local_sensitivities(tmp,sol,sol.t[3]) == extract_local_sensitivities(sol,3)

# asmatrix=true
@test extract_local_sensitivities(sol, length(sol), true) == (x, sense_res2)
@test extract_local_sensitivities(sol, sol.t[end], true) == (x, sense_res2)
@test extract_local_sensitivities(tmp, sol, sol.t[end], true) == (x, sense_res2)

# Return type inferred
@inferred extract_local_sensitivities(sol, 1)
@inferred extract_local_sensitivities(sol, 1, Val(true))
@inferred extract_local_sensitivities(sol, sol.t[3])
@inferred extract_local_sensitivities(sol, sol.t[3], Val(true))
@inferred extract_local_sensitivities(tmp, sol, sol.t[3])
@inferred extract_local_sensitivities(tmp, sol, sol.t[3], Val(true))

### ForwardDiff version

prob = ODEForwardSensitivityProblem(f.f,[1.0;1.0],(0.0,10.0),p,
                                    ForwardDiffSensitivity())
sol = solve(prob,Tsit5(),abstol=1e-14,reltol=1e-14,saveat=0.01)

xall, dpall = extract_local_sensitivities(sol)
@test xall ≈ res
@test dpall[1] ≈ da atol=1e-9

_,dpall_matrix = extract_local_sensitivities(sol,Val(true))
@test mapreduce(x->x[:, 2], hcat, dpall) == dpall_matrix[2]

x, dp = extract_local_sensitivities(sol,length(sol.t))
sense_res2 = reduce(hcat,dp)
@test fd_res == sense_res2

@test extract_local_sensitivities(sol,sol.t[3]) == extract_local_sensitivities(sol,3)

tmp = similar(sol[1])
@test extract_local_sensitivities(tmp,sol,sol.t[3]) == extract_local_sensitivities(sol,3)

# asmatrix=true
@test extract_local_sensitivities(sol, length(sol), true) == (x, sense_res2)
@test extract_local_sensitivities(sol, sol.t[end], true) == (x, sense_res2)
@test extract_local_sensitivities(tmp, sol, sol.t[end], true) == (x, sense_res2)

# Return type inferred
@inferred extract_local_sensitivities(sol, 1)
@inferred extract_local_sensitivities(sol, 1, Val(true))
@inferred extract_local_sensitivities(sol, sol.t[3])
@inferred extract_local_sensitivities(sol, sol.t[3], Val(true))
@inferred extract_local_sensitivities(tmp, sol, sol.t[3])
@inferred extract_local_sensitivities(tmp, sol, sol.t[3], Val(true))

# Test Float32

function f32(du,u,p,t)
  du[1] = dx = p[1]*u[1] - p[2]*u[1]*u[2]
  du[2] = dy = -p[3]*u[2] + u[1]*u[2]
end
p = [1.5f0,1.0f0,3.0f0]
prob = ODEForwardSensitivityProblem(f32,[1.0f0;1.0f0],(0.0f0,10.0f0),p)
sol = solve(prob,Tsit5())


text:
"""
The string part of JSCode.
"""
struct JSString
    source::String
end

"""
Javascript code that supports interpolation of Julia Objects.
Construction of JSCode via string macro:
```julia
jsc = js"console.log(\$(some_julia_variable))"
```
This will decompose into:
```julia
jsc.source == [JSString("console.log("), some_julia_variable, JSString("\"")]
```
"""
struct JSCode
    source::Vector{Union{JSString, Any}}
end

"""
Represent an asset stored at an URL.
We try to always have online & local files for assets
If one gives an online resource, it will be downloaded, to host it locally.
"""
struct Asset
    media_type::Symbol
    # We try to always have online & local files for assets
    # If you only give an online resource, we will download it
    # to also be able to host it locally
    online_path::String
    local_path::String
    onload::Union{Nothing, JSCode}
end

"""
Encapsulates frontend dependencies. Can be used in the following way:

```Julia
const noUiSlider = Dependency(
    :noUiSlider,
    # js & css dependencies are supported
    [
        "https://cdn.jsdelivr.net/gh/leongersen/noUiSlider/distribute/nouislider.min.js",
        "https://cdn.jsdelivr.net/gh/leongersen/noUiSlider/distribute/nouislider.min.css"
    ]
)
# use the dependency on the frontend:
evaljs(session, js"\$(noUiSlider).some_function(...)")
```
jsrender will make sure that all dependencies get loaded.
"""
struct Dependency
    name::Symbol # The JS Module name that will get loaded
    assets::Vector{Asset}
    # The global -> Function name, JSCode -> the actual function code!
    functions::Dict{Symbol, JSCode}
end

"""
A web session with a user
"""
struct Session
    # indicates, whether session is in fuse mode, so it won't
    # send any messages, until fuse ends and send them all at once
    fusing::Base.RefValue{Bool}
    connections::Vector{WebSocket}
    observables::Dict{String, Tuple{Bool, Observable}} # Bool -> if already registered with Frontend
    # We need to keep track of any objects that is used on the js side, to not free it!
    objects::Vector{Any}
    message_queue::Vector{Dict{Symbol, Any}}
    dependencies::Set{Asset}
    on_document_load::Vector{JSCode}
    id::String
    js_fully_loaded::Channel{Bool}
    on_websocket_ready::Any
end

abstract type AbstractJSObject end

const GLOBAL_UUID = Base.Threads.Atomic{UInt}(0)

"""
References objects stored in Javascript.
Maps the following expressions to actual calls on the Javascript side:
```julia
jso = JSObject(name, scope, typ)
# getfield:
x = jso.property # returns a new JSObject
# setfield
jso.property = "newval" # Works with JSObjects or Julia objects as newval
# call:
result = jso.func(args...) # works with julia objects and other JSObjects as args

# constructors are wrapped this way:
scene = jso.new.Scene() # same as in JS: scene = new jso.Scene()
```
"""
mutable struct JSObject <: AbstractJSObject
    # fields are private and not accessible via jsobject.field
    name::Symbol
    session::Session
    typ::Symbol
    # transporting the UUID allows us to have a uuid different from the objectid
    # which will help to better capture === equivalence on the js side.
    uuid::UInt

    function JSObject(name::Symbol, scope::Session, typ::Symbol)
        obj = new(name, scope, typ)
        uuid = GLOBAL_UUID[] += 1
        setfield!(obj, :uuid, uuid)
        finalizer(remove_js_reference, obj)
        return obj
    end
end

struct Routes
    table::Vector{Pair{Any, Any}}
end

function Routes(pairs::Pair...)
    return Routes([pairs...])
end

pattern_priority(x::Pair) = pattern_priority(x[1])
pattern_priority(x::String) = 1
pattern_priority(x::Tuple) = 2
pattern_priority(x::Regex) = 3

function Base.setindex!(routes::Routes, f, pattern)
    idx = findfirst(pair-> pair[1] == pattern, routes.table)
    if idx !== nothing
        routes.table[idx] = pattern => f
    else
        push!(routes.table, pattern => f)
    end
    # Sort for priority so that exact string matches come first
    sort!(routes.table, by = pattern_priority)
    return
end

function apply_handler(f, args...)
    f(args...)
end

function apply_handler(chain::Tuple, context, args...)
    f = first(chain)
    result = f(args...)
    apply_handler(Base.tail(chain), context, result...)
end

function delegate(routes::Routes, application, request::Request, args...)
    for (pattern, f) in routes.table
        match = match_request(pattern, request)
        if match !== nothing
            context = (
                routes = routes,
                application = application,
                request = request,
                match = match
            )
            return apply_handler(f, context, args...)
        end
    end
    # If no route is found we have a classic case of 404!
    # What a classic this response!
    return response_404("Didn't find route for $(request.target)")
end

function match_request(pattern::String, request)
    return request.target == pattern ? pattern : nothing
end

function match_request(pattern::Regex, request)
    return match(pattern, request.target)
end

function match_request(pattern::Tuple, request)
    return Matcha.matchat(pattern, request.target)
end

"""
The application one serves
"""
struct Application
    url::String
    port::Int
    sessions::Dict{String, Dict{String, Session}}
    server_task::Ref{Task}
    server_connection::Ref{TCPServer}
    routes::Routes
    websocket_routes::Routes
end
local_url(application::Application, url) = string("http://", application.url, ":", application.port, url)

WebSockets.getrawstream(io::IO) = io

function websocket_request()
    headers = [
        "Host" => "127.0.0.1",
        "User-Agent" => "Mozilla/5.0 (Windows NT 10.0; Win64; x64; rv:68.0) Gecko/20100101 Firefox/68.0",
        "Accept" => "*/*",
        "Accept-Encoding" => "gzip, deflate, br",
        "Accept-Language" => "de,en-US;q=0.7,en;q=0.3",
        "Cache-Control" => "no-cache",
        "Connection" => "keep-alive, Upgrade",
        "Dnt" => "1",
        "Origin" => "https://localhost",
        "Pragma" => "no-cache",
        "Sec-Websocket-Extensions" => "permessage-deflate",
        "Sec-Websocket-Key" => "BL3d8I8KC5faPjubRM0riA==",
        "Sec-Websocket-Version" => "13",
        "Upgrade" => "websocket",
    ]
    msg = HTTP.Request(
        "GET",
        "/",
        headers,
        UInt8[],
        parent = nothing,
        version = v"1.1.0"
    )
    return Stream(msg, IOBuffer())
end

"""
warmup(application::Application)

Warms up the application, by sending a couple of request.
"""
function warmup(application::Application)
    yield() # yield to server task to give it a chance to get started
    task = application.server_task[]
    if Base.istaskdone(task)
        error("Webserver doesn't serve! Error: $(fetch(task))")
    end
    # Make a websocket request
    stream = websocket_request()
    try
        @async stream_handler(application, stream)
        write(stream, "blaaa")

    catch e
        # TODO make it not error so we can test this properly
        # This will error, since its not a propper websocket request
        @debug "Error in stream_handler" exception=e
    end
    target = AssetRegistry.register(JSCallLibLocal.local_path) # http target part
    asset_url = local_url(application, target)
    request = Request("GET", target)

    delegate(application.routes, application, request)

    if Base.istaskdone(task)
        error("Webserver doesn't serve! Error: $(fetch(task))")
    end
    resp = WebSockets.HTTP.get(asset_url, readtimeout=500, retries=1)
    if resp.status != 200
        error("Webserver didn't start succesfully")
    end
    return
end

function stream_handler(application::Application, stream::Stream)
    if WebSockets.is_upgrade(stream)
        WebSockets.upgrade(stream) do request, websocket
            delegate(
                application.websocket_routes, application, request, websocket
            )
        end
        return
    end
    f = HTTP.RequestHandlerFunction() do request
        delegate(
            application.routes, application, request,
        )
    end
    try
        HTTP.handle(f, stream)
    catch e
        # we expect the IOError to happen, if either the page gets closed
        # or we close the server!
        if !(e isa Base.IOError && e.msg == "stream is closed or unusable")
            rethrow(e)
        end
    end
end

const MATCH_HEX = r"[\da-f]"
const MATCH_UUID4 = MATCH_HEX^8 * r"-" * (MATCH_HEX^4 * r"-")^3 * MATCH_HEX^12
const MATCH_SESSION_ID = MATCH_UUID4 * r"/" * MATCH_HEX^4

function serve_dom(context, dom)
    application = context.application
    session_id = string(uuid4())
    session = Session()
    application.sessions[session_id] = Dict("base" => session)
    html_dom = Base.invokelatest(dom, session, context.request)
    return html(dom2html(session, session_id, html_dom))
end

"""
Application(
        dom, url::String, port::Int;
        verbose = false
    )

Creates an application that manages the global server state!
"""
function Application(
        dom, url::String, port::Int;
        verbose = false,
        routes = Routes(
            "/" => ctx-> serve_dom(ctx, dom),
            r"/assetserver/" * MATCH_HEX^40 * r"-.*" => file_server,
            r".*" => (context)-> response_404()
        ),
        websocket_routes = Routes(
            r"/" * MATCH_SESSION_ID => websocket_handler
        )
    )

    application = Application(
        url, port, Dict{String, Dict{String, Session}}(),
        Ref{Task}(), Ref{TCPServer}(),
        routes,
        websocket_routes
    )
    try
        start(application; verbose=verbose)
        # warmup server!
        warmup(application)
    catch e
        close(application)
        rethrow(e)
    end

    return application
end

function isrunning(application::Application)
    return (isassigned(application.server_task) &&
        isassigned(application.server_connection) &&
        !istaskdone(application.server_task[]) &&
        isopen(application.server_connection[]))
end

function Base.close(application::Application)
    # Closing the io connection should shut down the HTTP listen loop
    for (id, clients) in application.sessions
        for (id, session) in clients
            close(session)
        end
    end

    if isassigned(application.server_connection)
        close(application.server_connection[])
        @assert !isopen(application.server_connection[])
    end
    # For good measures, wait until the task finishes!
    if isassigned(application.server_task)
        try
            wait(application.server_task[])
            @assert !isdone(application.server_connection[])
        catch e
            @debug "Server task failed with an (expected) exception on close" exception=e
        end
    end
    # Sometimes, the first request after closing will still go through
    # see: https://github.com/JuliaWeb/HTTP.jl/pull/494
    # We need to make sure, that we are the ones making this request,
    # So that a newly opened connection won't get a faulty response from this server!
    try
        app_url = JSServe.local_url(application, "/")
        while true
            x = HTTP.get(app_url, readtimeout=1, retries=1)
            x.status != 200 && break # we got a bad request, maining server is closed!
        end
    catch e
        # This is expected to fail!
        @debug "Failed get request successfully after closing server!" exception=e
    end
    @assert !isrunning(application)
end

function start(application::Application; verbose=false)
    isrunning(application) && return
    address = Sockets.InetAddr(parse(Sockets.IPAddr, application.url), application.port)
    ioserver = Sockets.listen(address)
    application.server_connection[] = ioserver
    # pass tcp connection to listen, so that we can close the server
    application.server_task[] = @async HTTP.listen(
            application.url, application.port; server=ioserver, verbose=verbose
        ) do stream::Stream
        Base.invokelatest(stream_handler, application, stream)
    end
    return
end

function route!(application::Application, pattern_f::Pair)
    application.routes[pattern_f[1]] = pattern_f[2]
end

function route!(f, application::Application, pattern)
    route!(application, pattern => f)
end

function websocket_route!(application::Application, pattern_f::Pair)
    application.webscoket_routes[pattern_f[1]] = pattern_f[2]
end

ext: jl


text:
using JSON
using Markdown

Markdown.htmlesc(io::IO, str::AbstractString) = print(io, str)

file_config = open("config.json")
config = JSON.parse(file_config)
html_head = open("include/head.html")
head = read(html_head, String)
head = replace(head, r"\{\{[\w\-\.\/]+?\}\}" => s -> config[s])
html_header = open("include/header.html")
header = read(html_header, String)
header = replace(header, r"\{\{[\w\-\.\/]+?\}\}" => s -> config[s])
html_avatar = open("include/avatar.html")
avatar = read(html_avatar, String)
avatar = replace(avatar, r"\{\{[\w\-\.\/]+?\}\}" => s -> config[s])
html_footer = open("include/footer.html")
footer = read(html_footer, String)
footer = replace(footer, r"\{\{[\w\-\.\/]+?\}\}" => s -> config[s])

cd("_markdown/")
files = readdir()
for markdown_file in files
    file = open(markdown_file)
    filep = read(file, String)
    filep = replace(filep, r"\{\{[\w\-\.\/]+?\}\}" => s -> config[s])
    md = Markdown.parse(filep)
    md_html = Markdown.html(md)
    md_title = match(r"
(.+?)<\/h1>", md_html).captures[1]
    full_html = """<!DOCTYPE html>

<title>$md_title | NYTxx</title>
$head

    $header

        $md_html

        $avatar

    $footer

"""
    open("../blog/$(match(r"[\w\-\/]+", markdown_file).match)/index.html", "w") do target
        write(target, full_html)
    end
    println("file ($markdown_file | $md_title) loaded")
    close(file)
end

close(html_head)
close(file_config)


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

using Documenter
using Arrow

makedocs(;
    modules=[Arrow],
    repo="https://github.com/JuliaData/Arrow.jl/blob/{commit}{path}#L{line}",
    sitename="Arrow.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaData.github.io/Arrow.jl",
        assets=String[],
    ),
    pages = [
        "Home" => "index.md",
        "User Manual" => "manual.md",
        "API Reference" => "reference.md"
    ]
)

deploydocs(;
    repo="github.com/JuliaData/Arrow.jl",
    devbranch = "main"
)

max_stars_repo_name: jonalm/arrow-julia

max_issues_repo_name: jonalm/arrow-julia

max_forks_repo_name: jonalm/arrow-julia


text:
mutable struct PlateGenerateState{U,V}
    score::Float64
    weight::Float64
    is_empty::Bool
    args::U
    nodes::V
end

function generate_new_trace!(gen::Plate{T,U}, key::Int, state::PlateGenerateState) where {T,U}
    kernel_args = get_args_for_key(state.args, key)
    local subtrace::U
    if haskey(state.nodes, key)
        node = state.nodes[key]
        (subtrace, kernel_weight) = generate(gen.kernel, kernel_args, node)
        state.weight += kernel_weight
    else
        subtrace = simulate(gen.kernel, kernel_args)
    end
    state.is_empty = state.is_empty && !has_choices(subtrace)
    call::CallRecord = get_call_record(subtrace)
    state.score += call.score
    (subtrace, call.retval::T)
end

function generate_process_constraints!(nodes, key::Int, node, len)
    if key <= len
        nodes[key] = node
    else
        error("Constraints under key $key not part of trace")
    end
end

function generate_process_constraints!(nodes, key, node, len)
    error("Constraints under key $key not part of trace")
end

function generate(gen::Plate{T,U}, args, constraints) where {T,U}
    len = length(args[1])

    # collect constraints, indexed by key, and error if any are not in the trace
    nodes = Dict{Int, Any}()
    for (key, node) in get_internal_nodes(constraints)
        generate_process_constraints!(nodes, key, node, len)
    end

    # collect initial state
    state = PlateGenerateState(0., 0., true, args, nodes)
    subtraces = Vector{U}(undef, len)
    retvals = Vector{T}(undef, len)

    # visit each new application and run generate (or simulate) on it
    for key=1:len
        (subtraces[key], retvals[key]) = generate_new_trace!(gen, key, state)
    end

    call = CallRecord(state.score, PersistentVector{T}(retvals), args)
    trace = VectorTrace{T,U}(PersistentVector{U}(subtraces), call, state.is_empty)
    (trace, state.weight)
end

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
# glm.jl

typealias Hyper Dict{Symbol,Float64}

immutable Param_GLM
  β::Vector{Float64}
  θ::Hyper
end

function θ_vec(θ::Hyper, keys::Vector{Symbol})
  return [θ[k] for k in keys]
end

function toθ(v::Vector{Float64}, keys::Vector{Symbol})
  assert(length(v) == length(keys))
  d = Hyper()
  for i in 1:length(keys)
    d[keys[i]] = v[i]
  end
  return d
end

immutable GLM
  post_params::Vector{Param_GLM}
  y::Vector{Float64}
  X::Matrix{Float64}
end

"""
invlink(xi'β) should be the mean of distribution F
loglike(y::Vector{Float64}, Xb::Vector{Float64}, θ::Hyper) should return the loglikelihood
"""

lp_θ_default(θ::Hyper) = 0.0
β_logprior_default(β::Vector{Float64}) = 0.0
function glm(y::Vector{Float64}, X::Matrix{Float64}, # include your own intercept!
             Σ_β::Matrix{Float64}, loglike;
             β_logprior = β_logprior_default,
             Σ_θ::Matrix{Float64}=eye(1)*1.,
             θ_bounds::Matrix{Float64}=[0. Inf],
             θ_names::Vector{Symbol}=[:empty],
             θ_logprior=lp_θ_default,
             B::Int=10000, burn::Int=1000, printFreq::Int=0)

  assert(size(θ_bounds,2) == 2)

  const empty = θ_names[1] == :empty
  const θ_lower = θ_bounds[:,1]
  const θ_upper = θ_bounds[:,2]

  lp(β::Vector{Float64}, θ::Hyper) = β_logprior(β) + θ_logprior(θ)
  ll(β::Vector{Float64}, θ::Hyper) = loglike(y, X*β, θ)

  ll_plus_lp(β::Vector{Float64}, θ::Hyper) = ll(β, θ) + lp(β, θ)

  function update(p::Param_GLM)
    # Update Coefficients (β)
    β_cand = rand(MvNormal(p.β,Σ_β))
    if ll_plus_lp(β_cand, p.θ) - ll_plus_lp(p.β, p.θ) > log(rand())
      β_new = β_cand
    else
      β_new = copy(p.β)
    end

    # Update Hyper Params (θ)
    if !empty
      θ_curr = θ_vec(p.θ,θ_names)
      θ_cand = rand(MvNormal(θ_curr,Σ_θ))
      θ_cand_s = toθ(θ_cand,θ_names)

      if all(θ_lower .< θ_cand .< θ_upper) && ll_plus_lp(β_new, θ_cand_s) - ll_plus_lp(β_new, p.θ) > log(rand())
        θ_new = θ_cand_s
      else
        θ_new = copy(p.θ)
      end
    else
      θ_new = p.θ
    end

    return Param_GLM(β_new,θ_new)
  end #update

  const β₀ = zeros(Float64,size(Σ_β,1))
  const θ₀ = toθ(ones(Float64,size(Σ_θ,1)), θ_names)
  const init = Param_GLM(β₀,θ₀)

  const post = MCMC.gibbs(init, update, B, burn, printFreq=printFreq)
  return GLM(post,y,X)
end #glm

function stats(v::Vector{Float64}; α::Float64=.05)
  (mean(v), std(v), quantile(v,[α/2, 1-α/2]))
end

function stats(m::Matrix{Float64}; α::Float64=.05)
  const P = size(m,2)
  return (vec(mean(m,1)), vec(std(m,1)),
          hcat([quantile(m[:,p], [α/2, 1-α/2]) for p in 1:P]...)')
end

immutable Coef
  mean::Vector{Float64}
  std::Vector{Float64}
  q::Matrix{Float64}
  ne0::Vector{Bool}
  acc::Float64
end

function summary(model::GLM; α::Float64=.05)
  const β = hcat(map(o -> o.β, model.post_params)...)'
  const β_stats = stats(β)
  const ne0 = !(β_stats[3][:,1] .<= 0 .<= β_stats[3][:,2])
  const β_acc = size(unique(β,1),1) / size(β,1)

  const coef = Coef(β_stats[1], β_stats[2], β_stats[3], ne0, β_acc)
  const hyper = map(m -> m.θ, model.post_params)
  return  Summary_glm(coef,hyper)
end

function show(io::IO, coef::Coef)
  const P = length(coef.mean)
  @printf "%5s%10s%10s%10s%10s\n" "" "mean" "std" "CI_lower" "CI_upper"
  for k in 1:P
    @printf "%5s%10.4f%10.4f%10.4f%10.4f%3s\n" string("β",k-1) coef.mean[k] coef.std[k] coef.q[k,1] coef.q[k,2] (coef.ne0[k] ? "*" : "")
  end
end

immutable Summary_glm
  coef::Coef
  hyper::Vector{Hyper}
end

function show(io::IO, hyper::Vector{Hyper})
  for k in keys(hyper[1])
    s = stats(map(h -> h[k], hyper))
    @printf "%5s%10.4f%10.4f%10.4f%10.4f\n" string(k) s[1] s[2] s[3][1] s[3][2]
  end
end

function show(io::IO, sglm::Summary_glm)
  const empty = collect(keys(sglm.hyper[1]))[1] == :empty

  show(io, sglm.coef)
  if !empty
    show(io, sglm.hyper)
  end

  k = collect(keys(sglm.hyper[1]))[1]
  some_hyper = map(h -> h[k], sglm.hyper)
  acc_hyper = length(unique(some_hyper)) / length(some_hyper)

  println()

  @printf "%5s%10.4f\n" "accβ" sglm.coef.acc
  if !empty
    @printf "%5s%10.4f\n" "accθ" acc_hyper
  end
end

function latex(sglm::Summary_glm; digits::Int=4)
  const empty = collect(keys(sglm.hyper[1]))[1] == :empty
  const coef = sglm.coef
  const hyper = sglm.hyper

  const P = length(coef.mean)
  println("\\begin{center}\\begin{tabular}{crrrrc}")
  println("& mean & sd & CI 2.5\\% & CI 97.5\\% \\\\")
  for k in 1:P
    @printf "%s&%10.4f&%10.4f&%10.4f&%10.4f&%3s\\\\\n" string("\$\\beta_{",k-1,"}\$") coef.mean[k] coef.std[k] coef.q[k,1] coef.q[k,2] (coef.ne0[k] ? "*" : "")
  end

  if !empty
    for k in keys(hyper[1])
      s = stats(map(h -> h[k], hyper))
      @printf "%s&%10.4f&%10.4f&%10.4f&%10.4f&\\\\\n" string(k) s[1] s[2] s[3][1] s[3][2]
    end
  end

  k = collect(keys(sglm.hyper[1]))[1]
  some_hyper = map(h -> h[k], sglm.hyper)
  acc_hyper = length(unique(some_hyper)) / length(some_hyper)

  println("\\hline\\\\")
  @printf "%s&%10.4f\\\\\n" "Acceptance \$\\beta\$" sglm.coef.acc
  if !empty
    @printf "%s&%10.4f\\\\\n" "Acceptance Hyper" acc_hyper
  end
  println("\\end{tabular}\\end{center}")
end

#function dic(g::GLM)
#end


text:
"""
    (newids, removedids) = split_trio_from_ped("pedigree", "id")

Take input id and remove that id and any parents of that id
from a pedigree file.  Write a cleaned pedigree file and a
pedigree file for the removed ids to disk.  Return strings
containing the ids for the new files.
"""
function split_trio_from_ped(pedigree, id)

    newids = String[]
    dtrioids = String[]
    ped = String[]
    badids = Dict{String,String}()

    open(pedigree) do file

        replace_parent = String[]

        for i in eachline(file)
            l = split(i, '\t')
            if l[2] == id
                badids[l[2]] = "bad"
                badids[l[3]] = "bad"
                badids[l[4]] = "bad"
                replace_parent = [l[1], l[2], "0", "0", l[5], l[6]]
            end
            push!(ped, i)
        end

        OUT = open("$pedigree.fam", "w")
        OUT2 = open("$pedigree.trio", "w")

        for i in ped
            l = split(i, '\t')
            x = l[2]
            if haskey(badids, x)
                push!(dtrioids, x)
                println(OUT2, i)
            else
                push!(newids, x)
                println(OUT, i)
            end
        end

        println(OUT, join(replace_parent, '\t'))
        push!(newids, id)

        close(OUT)
        close(OUT2)

    end

    newids = join(newids, ",")
    dtrioids = join(dtrioids, ",")
    return (newids, dtrioids)

end

max_stars_repo_name: ScottWatkins/STRdiff

max_issues_repo_name: ScottWatkins/STRdiff

max_forks_repo_name: ScottWatkins/STRdiff


text:
@defstruct BinaryAccuracyLayer Layer (
  name :: AbstractString = "binary-accuracy",
  report_error :: Bool = false,
  (threshold :: Float64 = 0.0, isfinite(threshold)),
  (bottoms :: Vector{Symbol} = Symbol[], length(bottoms) == 2),
)
@characterize_layer(BinaryAccuracyLayer,
  is_sink    => true,
  has_stats  => true,
)

type BinaryAccuracyLayerState <: LayerState
  layer :: BinaryAccuracyLayer

  accuracy :: Float64
  n_wrong :: Int
  n_accum  :: Int

  etc :: Any
end

function setup_etc(backend::CPUBackend, layer::BinaryAccuracyLayer, inputs)
  nothing
end

function setup(backend::Backend, layer::BinaryAccuracyLayer, inputs::Vector{Blob}, diffs::Vector{Blob})
  etc = setup_etc(backend, layer, inputs)

  return BinaryAccuracyLayerState(layer, 0.0, 0, 0, etc)
end
function shutdown(backend::CPUBackend, state::BinaryAccuracyLayerState)
end

function reset_statistics(state::BinaryAccuracyLayerState)
  state.n_accum = 0
  state.n_wrong = 0
  state.accuracy = 0.0
end

function dump_statistics(storage, state::BinaryAccuracyLayerState, show::Bool)
  update_statistics(storage, "$(state.layer.name)-accuracy", state.accuracy)
  if state.layer.report_error
    update_statistics(storage, "$(state.layer.name)-error", 1-state.accuracy)
  end

  if show
    accuracy = @sprintf("%.4f%%", state.accuracy*100)
    @info("  Accuracy (avg over $(state.n_accum)) = $accuracy")
  end
end

function forward(backend::CPUBackend, state::BinaryAccuracyLayerState, inputs::Vector{Blob})
  pred = inputs[1].data
  label = inputs[2].data

  T = eltype(pred)
  threshold = convert(T, state.layer.threshold)

  for i=1:length(pred)
    if (pred[i]>threshold) $ (label[i]>threshold)
      state.n_wrong += 1
    end
  end
  state.n_accum += length(pred)
  state.accuracy = (state.n_accum-state.n_wrong) / state.n_accum
end

function backward(backend::Backend, state::BinaryAccuracyLayerState, inputs::Vector{Blob}, diffs::Vector{Blob})
end


text:
#=
    MadNLP utils
=#

@kernel function _copy_diag!(dest, src)
    i = @index(Global)
    dest[i] = src[i, i]
end

function MadNLP.diag!(dest::CuVector{T}, src::CuMatrix{T}) where T
    @assert length(dest) == size(src, 1)
    ev = _copy_diag!(CUDADevice())(dest, src, ndrange=length(dest))
    wait(ev)
end

@kernel function _add_diagonal!(dest, src1, src2)
    i = @index(Global)
    dest[i, i] = src1[i] + src2[i]
end

function MadNLP.diag_add!(dest::CuMatrix, src1::CuVector, src2::CuVector)
    ev = _add_diagonal!(CUDADevice())(dest, src1, src2, ndrange=size(dest, 1))
    wait(ev)
end

#=
   Contiguous views do not dispatch to the correct copyto! kernel
   in CUDA.jl. To avoid fallback to the (slow) implementation in Julia Base,
   we overload the copyto! operator locally
=#
_copyto!(dest::CuArray, src::Array) = copyto!(dest, src)
function _copyto!(dest::CuArray, src::SubArray)
    @assert src.stride1 == 1 # src array should be one-strided
    n = length(dest)
    offset = src.offset1
    p_src = parent(src)
    copyto!(dest, 1, p_src, offset+1, n)
end

#=
    MadNLP kernels
=#

# Overload MadNLP.is_valid to avoid fallback to default is_valid, slow on GPU
MadNLP.is_valid(src::CuArray) = true

# Constraint scaling
function MadNLP.scale_constraints!(
    nlp::AbstractNLPModel,
    con_scale::AbstractVector,
    jac::CuMatrix;
    max_gradient=1e-8,
)
    # Compute reduction on the GPU with built-in CUDA.jl function
    d_con_scale = maximum(abs, jac, dims=2)
    copyto!(con_scale, d_con_scale)
    con_scale .= min.(1.0, max_gradient ./ con_scale)
end

@kernel function _treat_fixed_variable_kernell!(dest, ind_fixed)
    k, j = @index(Global, NTuple)
    i = ind_fixed[k]

    if i == j
        dest[i, i] = 1.0
    else
        dest[i, j] = 0.0
        dest[j, i] = 0.0
    end
end

function MadNLP.treat_fixed_variable!(kkt::MadNLP.AbstractKKTSystem{T, MT}) where {T, MT<:CuMatrix{T}}
    length(kkt.ind_fixed) == 0 && return
    aug = kkt.aug_com
    d_ind_fixed = kkt.ind_fixed |> CuVector # TODO: allocate ind_fixed directly on the GPU
    ndrange = (length(d_ind_fixed), size(aug, 1))
    ev = _treat_fixed_variable_kernell!(CUDADevice())(aug, d_ind_fixed, ndrange=ndrange)
    wait(ev)
end

#=
    AbstractDenseKKTSystem
=#

function MadNLP.jtprod!(y::AbstractVector, kkt::MadNLP.AbstractDenseKKTSystem{T, VT, MT}, x::AbstractVector) where {T, VT<:CuVector{T}, MT<:CuMatrix{T}}
    # Load buffers
    m = size(kkt.jac, 1)
    nx = size(kkt.jac, 2)
    ns = length(kkt.ind_ineq)
    haskey(kkt.etc, :jac_w1) || (kkt.etc[:jac_w1] = CuVector{T}(undef, m))
    haskey(kkt.etc, :jac_w2) || (kkt.etc[:jac_w2] = CuVector{T}(undef, nx))
    haskey(kkt.etc, :jac_w3) || (kkt.etc[:jac_w3] = CuVector{T}(undef, ns))

    d_x = kkt.etc[:jac_w1]::VT
    d_yx = kkt.etc[:jac_w2]::VT
    d_ys = kkt.etc[:jac_w3]::VT

    # x and y can be host arrays. Copy them on the device to avoid side effect.
    _copyto!(d_x, x)

    # / x
    LinearAlgebra.mul!(d_yx, kkt.jac', d_x)
    copyto!(parent(y), 1, d_yx, 1, nx)

    # / s
    d_ys .= -d_x[kkt.ind_ineq] .* kkt.constraint_scaling[kkt.ind_ineq]
    copyto!(parent(y), nx+1, d_ys, 1, ns)
    return
end

function MadNLP.set_aug_diagonal!(kkt::MadNLP.AbstractDenseKKTSystem{T, VT, MT}, ips::MadNLP.InteriorPointSolver) where {T, VT<:CuVector{T}, MT<:CuMatrix{T}}
    haskey(kkt.etc, :pr_diag_host) || (kkt.etc[:pr_diag_host] = Vector{T}(undef, length(kkt.pr_diag)))
    pr_diag_h = kkt.etc[:pr_diag_host]::Vector{T}
    # Broadcast is not working as MadNLP array are allocated on the CPU,
    # whereas pr_diag is allocated on the GPU
    pr_diag_h .= ips.zl./(ips.x.-ips.xl) .+ ips.zu./(ips.xu.-ips.x)
    copyto!(kkt.pr_diag, pr_diag_h)
    fill!(kkt.du_diag, 0.0)
end

#=
    DenseKKTSystem kernels
=#

function LinearAlgebra.mul!(y::AbstractVector, kkt::MadNLP.DenseKKTSystem{T, VT, MT}, x::AbstractVector) where {T, VT<:CuVector{T}, MT<:CuMatrix{T}}
    # Load buffers
    haskey(kkt.etc, :hess_w1) || (kkt.etc[:hess_w1] = CuVector{T}(undef, size(kkt.aug_com, 1)))
    haskey(kkt.etc, :hess_w2) || (kkt.etc[:hess_w2] = CuVector{T}(undef, size(kkt.aug_com, 1)))

    d_x = kkt.etc[:hess_w1]::VT
    d_y = kkt.etc[:hess_w2]::VT

    # x and y can be host arrays. Copy them on the device to avoid side effect.
    copyto!(d_x, x)
    LinearAlgebra.mul!(d_y, kkt.aug_com, d_x)
    copyto!(y, d_y)
end
function LinearAlgebra.mul!(y::MadNLP.ReducedKKTVector, kkt::MadNLP.DenseKKTSystem{T, VT, MT}, x::MadNLP.ReducedKKTVector) where {T, VT<:CuVector{T}, MT<:CuMatrix{T}}
    LinearAlgebra.mul!(MadNLP.full(y), kkt, MadNLP.full(x))
end

@kernel function _build_dense_kkt_system_kernel!(
    dest, hess, jac, pr_diag, du_diag, diag_hess, ind_ineq, con_scale, n, m, ns
)
    i, j = @index(Global, NTuple)
    if (i <= n)
        # Transfer Hessian
        if (i == j)
            dest[i, i] = pr_diag[i] + diag_hess[i]
        else
            dest[i, j] = hess[i, j]
        end
    elseif i <= n + ns
        # Transfer slack diagonal
        dest[i, i] = pr_diag[i]
        # Transfer Jacobian wrt slack
        js = i - n
        is = ind_ineq[js]
        dest[is + n + ns, is + n] = - con_scale[is]
        dest[is + n, is + n + ns] = - con_scale[is]
    elseif i <= n + ns + m
        # Transfer Jacobian wrt variable x
        i_ = i - n - ns
        dest[i, j] = jac[i_, j]
        dest[j, i] = jac[i_, j]
        # Transfer dual regularization
        dest[i, i] = du_diag[i_]
    end
end

function MadNLP._build_dense_kkt_system!(
    dest::CuMatrix, hess::CuMatrix, jac::CuMatrix,
    pr_diag::CuVector, du_diag::CuVector, diag_hess::CuVector, ind_ineq, con_scale, n, m, ns
)
    ind_ineq_gpu = ind_ineq |> CuArray
    ndrange = (n+m+ns, n)
    ev = _build_dense_kkt_system_kernel!(CUDADevice())(
        dest, hess, jac, pr_diag, du_diag, diag_hess, ind_ineq_gpu, con_scale, n, m, ns,
        ndrange=ndrange
    )
    wait(ev)
end

#=
    DenseCondensedKKTSystem
=#
function MadNLP.get_slack_regularization(kkt::MadNLP.DenseCondensedKKTSystem{T, VT, MT}) where {T, VT<:CuVector{T}, MT<:CuMatrix{T}}
    n, ns = MadNLP.num_variables(kkt), kkt.n_ineq
    return view(kkt.pr_diag, n+1:n+ns) |> Array
end
function MadNLP.get_scaling_inequalities(kkt::MadNLP.DenseCondensedKKTSystem{T, VT, MT}) where {T, VT<:CuVector{T}, MT<:CuMatrix{T}}
    return kkt.constraint_scaling[kkt.ind_ineq] |> Array
end

@kernel function _build_jacobian_condensed_kernel!(
    dest, jac, pr_diag, ind_ineq, con_scale, n, m_ineq,
)
    i, j = @index(Global, NTuple)
    is = ind_ineq[i]
    @inbounds dest[i, j] = jac[is, j] * sqrt(pr_diag[n+i]) / con_scale[is]
end

function MadNLP._build_ineq_jac!(
    dest::CuMatrix, jac::CuMatrix, pr_diag::CuVector,
    ind_ineq::AbstractVector, ind_fixed::AbstractVector, con_scale::CuVector, n, m_ineq,
)
    (m_ineq == 0) && return # nothing to do if no ineq. constraints
    ind_ineq_gpu = ind_ineq |> CuArray
    ndrange = (m_ineq, n)
    ev = _build_jacobian_condensed_kernel!(CUDADevice())(
        dest, jac, pr_diag, ind_ineq_gpu, con_scale, n, m_ineq,
        ndrange=ndrange, dependencies=Event(CUDADevice()),
    )
    wait(ev)
    # need to zero the fixed components
    dest[:, ind_fixed] .= 0.0
    return
end

@kernel function _build_condensed_kkt_system_kernel!(
    dest, hess, jac, pr_diag, du_diag, ind_eq, n, m_eq,
)
    i, j = @index(Global, NTuple)

    # Transfer Hessian
    if i <= n
        if i == j
            @inbounds dest[i, i] += pr_diag[i] + hess[i, i]
        else
            @inbounds dest[i, j] += hess[i, j]
        end
    elseif i <= n + m_eq
        i_ = i - n
        @inbounds is = ind_eq[i_]
        # Jacobian / equality
        @inbounds dest[i_ + n, j] = jac[is, j]
        @inbounds dest[j, i_ + n] = jac[is, j]
        # Transfer dual regularization
        @inbounds dest[i_ + n, i_ + n] = du_diag[is]
    end
end

function MadNLP._build_condensed_kkt_system!(
    dest::CuMatrix, hess::CuMatrix, jac::CuMatrix,
    pr_diag::CuVector, du_diag::CuVector, ind_eq::AbstractVector, n, m_eq,
)
    ind_eq_gpu = ind_eq |> CuArray
    ndrange = (n + m_eq, n)
    ev = _build_condensed_kkt_system_kernel!(CUDADevice())(
        dest, hess, jac, pr_diag, du_diag, ind_eq_gpu, n, m_eq,
        ndrange=ndrange, dependencies=Event(CUDADevice()),
    )
    wait(ev)
end

function LinearAlgebra.mul!(y::AbstractVector, kkt::MadNLP.DenseCondensedKKTSystem{T, VT, MT}, x::AbstractVector) where {T, VT<:CuVector{T}, MT<:CuMatrix{T}}
    if length(y) == length(x) == size(kkt.aug_com, 1)
        # Load buffers
        haskey(kkt.etc, :hess_w1) || (kkt.etc[:hess_w1] = CuVector{T}(undef, size(kkt.aug_com, 1)))
        haskey(kkt.etc, :hess_w2) || (kkt.etc[:hess_w2] = CuVector{T}(undef, size(kkt.aug_com, 1)))

        d_x = kkt.etc[:hess_w1]::VT
        d_y = kkt.etc[:hess_w2]::VT

        # Call parent() as CUDA does not dispatch on proper copyto! when passed a view
        copyto!(d_x, 1, parent(x), 1, length(x))
        LinearAlgebra.mul!(d_y, kkt.aug_com, d_x)
        copyto!(y, d_y)
    else
        # Load buffers
        haskey(kkt.etc, :hess_w3) || (kkt.etc[:hess_w3] = CuVector{T}(undef, length(x)))
        haskey(kkt.etc, :hess_w4) || (kkt.etc[:hess_w4] = CuVector{T}(undef, length(y)))

        d_x = kkt.etc[:hess_w3]::VT
        d_y = kkt.etc[:hess_w4]::VT

        # Call parent() as CUDA does not dispatch on proper copyto! when passed a view
        copyto!(d_x, 1, parent(x), 1, length(x))
        MadNLP._mul_expanded!(d_y, kkt, d_x)
        copyto!(y, d_y)
    end
end
function LinearAlgebra.mul!(y::MadNLP.ReducedKKTVector, kkt::MadNLP.DenseCondensedKKTSystem{T, VT, MT}, x::MadNLP.ReducedKKTVector) where {T, VT<:CuVector{T}, MT<:CuMatrix{T}}
    LinearAlgebra.mul!(MadNLP.full(y), kkt, MadNLP.full(x))
end

function MadNLP.jprod_ineq!(y::AbstractVector, kkt::MadNLP.DenseCondensedKKTSystem{T, VT, MT}, x::AbstractVector) where {T, VT<:CuVector{T}, MT<:CuMatrix{T}}
    # Create buffers
    haskey(kkt.etc, :jac_ineq_w1) || (kkt.etc[:jac_ineq_w1] = CuVector{T}(undef, kkt.n_ineq))
    haskey(kkt.etc, :jac_ineq_w2) || (kkt.etc[:jac_ineq_w2] = CuVector{T}(undef, size(kkt.jac_ineq, 2)))

    y_d = kkt.etc[:jac_ineq_w1]::VT
    x_d = kkt.etc[:jac_ineq_w2]::VT

    # Call parent() as CUDA does not dispatch on proper copyto! when passed a view
    copyto!(x_d, 1, parent(x), 1, length(x))
    LinearAlgebra.mul!(y_d, kkt.jac_ineq, x_d)
    copyto!(parent(y), 1, y_d, 1, length(y))
end


text:
###############################################################################
## LASSO Problem
###############################################################################
function get_randomlasso(n, m, sparsity; reg=regularizer_l1, seed=1234)
    # A is drawn from the standard normal distribution
    # b = Ax0 + e where e is taken from the normal distribution with standard deviation 0.001.
    # We set λ1 so that the original sparsity is ultimately recovered.

    Random.seed!(seed)
    A = rand(Normal(), m, n)

    Random.seed!(seed)
    x0 = rand(Normal(), n)

    if reg == regularizer_l1
        Random.seed!(seed)
        inds_nz = Set(randperm(n)[1:sparsity])

        for i in 1:n
            !(i in inds_nz) && (x0[i] = 0)
        end

    elseif reg <: regularizer_l12
        T = reg.parameters[1]
        ngroups = Int(ceil(n/T))

        Random.seed!(seed)
        inds_nz = Set(randperm(ngroups)[1:sparsity])

        for i in 1:ngroups-1
            !(i in inds_nz) && (x0[T*i-3:T*i] .= 0)
        end
        !(ngroups in inds_nz) && (x0[T*ngroups-3:end] .= 0)

    elseif reg == regularizer_lnuclear
        n_mat = Int(sqrt(n))

        Random.seed!(seed)
        basisvecs = rand(Normal(), n_mat, sparsity)

        x0 = zeros(n)
        for sp in 1:sparsity
            x0 += vec(basisvecs[:, sp] * transpose(basisvecs[:, sp]))
        end

    elseif reg <: regularizer_distball
        p = reg.parameters[1]

        Random.seed!(seed)
        x0 = rand(Normal(), n); x0 /= norm(x0, p)

    else
        @error "hunhandled regularizer $reg"
    end

    delta = 0.01
    Random.seed!(seed)
    e = rand(Normal(0, delta^2), m)

    y = A*x0+e
    pb = LassoPb{reg}(A, y, delta, n, x0)

    return pb
end

function get_IM_testcase(; n=80, m=130, sparsity=0.1)
    # A is drawn from the standard normal distribution
    # b = Ax0 + e where x0 is a 10% sparse vector taken from the normal distribution,
    # e is taken from the normal distribution with standard deviation 0.001.
    # We set λ1 so that the original sparsity is ultimately recovered.

    Random.seed!(23)
    A = rand(Normal(), m, n)

    Random.seed!(23)
    x0 = rand(Normal(), n)
    for i in 1:n
        rand() < sparsity && (x0[i]=0)
    end

    Random.seed!(23)
    e = rand(Normal(0, 0.001), m)

    y = A*x0+e

    return LassoPb{regularizer_l1}(A, y, 1., n, x0)
end

function get_randomlassopb_TV(; n, p, λ=1, TVsparsity=0.05)
    Random.seed!(23)
    A = rand(Float64, p, n)

	y_orig = zeros(n)
	y = zeros(n)
	y_orig[1] = 10

	## Ground truth vector
	Random.seed!(23); sparsity = floor.(rand(n-1) .+ TVsparsity)
	Random.seed!(23);
	aux = rand(Normal(0, 4.0), n-1) .* sparsity
	for i in 2:n
		y_orig[i] = y_orig[i-1] + aux[i-1]
	end

	## Noisy observation
	Random.seed!(23); pert = rand(Normal(0, 1.0), n)
    y = y_orig .+ pert

    y_obs = A * y

    return LassoPb{regularizer_TV}(A, y_obs, λ, n, y_orig)
end

function build_2dlasso_normdsol(;p=2, θ = π*8/(9*2))
    xcenter = [cos(θ), sin(θ)]
    xcenter /= norm(xcenter, p)

    diagscaling = [1, 10]
    A = Diagonal(diagscaling)
    y = diagscaling .* xcenter
    λ = 1

    return LassoPb{regularizer_distball{p}}(A, y, λ, 2, xcenter)
end

###############################################################################
## Logistic Problem
###############################################################################

## ionosphere problem
function get_ionopb()
    @assert isfile("data/ionosphere.data")

    rawdata = readdlm("data/ionosphere.data", ',')
    A = Matrix{Float64}(rawdata[:, 1:end-1])
    n = size(A, 2)

    y = map(x-> x=="g" ? 1 : -1, rawdata[:, end])

    return LogisticPb{regularizer_l1}(A, y, 0.1, n, zeros(1))
end

## Random sparse logit
function get_logit(;n=80, m=85, sparsity=0.5)

    A = rand(m, n)*10
    y = Vector([rand()>sparsity for i in 1:m])

    return LogisticPb{regularizer_l1}(A, y, 0.1, n, zeros(1))
end


text:
using GeometryTypes, ColorTypes
using GeometryTypes: contains
using Test
using Test: @inferred

@test !isnan(Point3f0(1))
@test isnan(Point3f0(NaN))
@test unit(Point3f0, 1) == Point3f0(1, 0, 0)
@test unit(Point3f0, 2) == Point3f0(0, 1, 0)
@test unit(Point3f0, 3) == Point3f0(0, 0, 1)

@testset "GeometryTypes" begin
    include("baseutils.jl")
    include("convexhulls.jl")
    include("polygons.jl")
    include("hyperrectangles.jl")
    include("faces.jl")
    include("meshes.jl")
    include("distancefields.jl")
    include("primitives.jl")
    include("decompose.jl")
    include("simplerectangle.jl")
    include("hypersphere.jl")
    include("typeutils.jl")
    include("simplices.jl")
    include("gjk.jl")
    include("lines.jl")
    include("polygons.jl")
    include("cylinder.jl")
    include("capsule.jl")
    include("algorithms.jl")
end


text:
"""
    DB.jl

This is a Julia port of a MATLAB implementation of batch and incremental
Davies-Bouldin (DB) Cluster Validity Index.

Authors:
MATLAB implementation: Leonardo Enzo Brito da Silva
Julia port: Sasha Petrenko <sap625@mst.edu>

REFERENCES
[1] D. L. Davies and D. W. Bouldin, "A cluster separation measure,"
IEEE Transaction on Pattern Analysis and Machine Intelligence, vol. 1,
no. 2, pp. 224–227, Feb. 1979.
[2] M. Moshtaghi, J. C. Bezdek, S. M. Erfani, C. Leckie, and J. Bailey,
"Online Cluster Validity Indices for Streaming Data," ArXiv e-prints, 2018,
arXiv:1801.02937v1 [stat.ML]. [Online].
[3] M. Moshtaghi, J. C. Bezdek, S. M. Erfani, C. Leckie, J. Bailey, "Online
cluster validity indices for performance monitoring of streaming data clustering,"
Int. J. Intell. Syst., pp. 1–23, 2018.
"""

"""
    DB

The stateful information of the Davies-Bouldin (DB) Cluster Validity Index.

# References
1. D. L. Davies and D. W. Bouldin, "A cluster separation measure," IEEE Transaction on Pattern Analysis and Machine Intelligence, vol. 1, no. 2, pp. 224–227, Feb. 1979.
2. M. Moshtaghi, J. C. Bezdek, S. M. Erfani, C. Leckie, and J. Bailey, "Online Cluster Validity Indices for Streaming Data," ArXiv e-prints, 2018, arXiv:1801.02937v1 [stat.ML]. [Online].
3. M. Moshtaghi, J. C. Bezdek, S. M. Erfani, C. Leckie, J. Bailey, "Online cluster validity indices for performance monitoring of streaming data clustering," Int. J. Intell. Syst., pp. 1–23, 2018.
"""
mutable struct DB <: CVI
    label_map::LabelMap
    dim::Integer
    n_samples::Integer
    mu_data::RealVector         # dim
    n::IntegerVector            # dim
    v::RealMatrix               # dim x n_clusters
    CP::RealVector              # dim
    S::RealVector               # dim
    R::RealMatrix               # dim x n_clusters
    G::RealMatrix               # dim x n_clusters
    D::RealMatrix               # n_clusters x n_clusters
    n_clusters::Integer
    criterion_value::Float
end # DB <: CVI

"""
    DB()

Default constructor for the Davies-Bouldin (DB) Cluster Validity Index.
"""
function DB()
    DB(
        LabelMap(),                     # label_map
        0,                              # dim
        0,                              # n_samples
        Array{Float, 1}(undef, 0),      # mu_data
        Array{Integer, 1}(undef, 0),    # n
        Array{Float, 2}(undef, 0, 0),   # v
        Array{Float, 1}(undef, 0),      # CP
        Array{Float, 1}(undef, 0),      # S
        Array{Float, 2}(undef, 0, 0),   # R
        Array{Float, 2}(undef, 0, 0),   # G
        Array{Float, 2}(undef, 0, 0),   # D
        0,                              # n_clusters
        0.0                             # criterion_value
    )
end # DB()

"""
    setup!(cvi::DB, sample::Vector{T}) where {T<:RealFP}
"""
function setup!(cvi::DB, sample::Vector{T}) where {T<:RealFP}
    # Get the feature dimension
    cvi.dim = length(sample)
    # Initialize the augmenting 2-D arrays with the correct feature dimension
    # NOTE: R is emptied and calculated in evaluate!, so it is not defined here
    cvi.v = Array{T, 2}(undef, cvi.dim, 0)
    cvi.G = Array{T, 2}(undef, cvi.dim, 0)
end # setup!(cvi::DB, sample::Vector{T}) where {T<:RealFP}

function param_inc!(cvi::DB, sample::RealVector, label::Integer)
    # Get the internal label
    i_label = get_internal_label!(cvi.label_map, label)

    n_samples_new = cvi.n_samples + 1
    if isempty(cvi.mu_data)
        mu_data_new = sample
        setup!(cvi, sample)
    else
        mu_data_new = (1 - 1/n_samples_new) .* cvi.mu_data + (1/n_samples_new) .* sample
    end

    if i_label > cvi.n_clusters
        n_new = 1
        v_new = sample
        CP_new = 0.0
        G_new = zeros(cvi.dim)
        S_new = 0.0
        if cvi.n_clusters == 0
            D_new = zeros(1,1)
        else
            D_new = zeros(cvi.n_clusters + 1, cvi.n_clusters + 1)
            D_new[1:cvi.n_clusters, 1:cvi.n_clusters] = cvi.D
            d_column_new = zeros(cvi.n_clusters + 1)
            for jx = 1:cvi.n_clusters
                d_column_new[jx] = sum((v_new - cvi.v[:, jx]).^2)
            end
            D_new[:, i_label] = d_column_new
            D_new[i_label, :] = transpose(d_column_new)
        end
        # Update 1-D parameters with a push
        cvi.n_clusters += 1
        push!(cvi.CP, CP_new)
        push!(cvi.n, n_new)
        push!(cvi.S, S_new)
        # Update 2-D parameters with appending and reassignment
        cvi.v = [cvi.v v_new]
        cvi.G = [cvi.G G_new]
        cvi.D = D_new
    else
        n_new = cvi.n[i_label] + 1
        v_new = (1 - 1/n_new) .* cvi.v[:, i_label] + (1/n_new) .* sample
        delta_v = cvi.v[:, i_label] - v_new
        diff_x_v = sample .- v_new
        CP_new = cvi.CP[i_label] + transpose(diff_x_v)*diff_x_v + cvi.n[i_label]*transpose(delta_v)*delta_v + 2*transpose(delta_v)*cvi.G[:, i_label]
        G_new = cvi.G[:, i_label] + diff_x_v + cvi.n[i_label].*delta_v
        S_new = CP_new / n_new
        d_column_new = zeros(cvi.n_clusters)
        for jx = 1:cvi.n_clusters
            # Skip the current i_label index
            if jx == i_label
                continue
            end
            d_column_new[jx] = sum((v_new - cvi.v[:, jx]).^2)
        end
        # Update parameters
        cvi.n[i_label] = n_new
        cvi.v[:, i_label] = v_new
        cvi.CP[i_label] = CP_new
        cvi.G[:, i_label] = G_new
        cvi.S[i_label] = S_new
        cvi.D[:, i_label] = d_column_new
        cvi.D[i_label, :] = transpose(d_column_new)
    end
    cvi.n_samples = n_samples_new
    cvi.mu_data = mu_data_new
end # param_inc!(cvi::DB, sample::RealVector, label::Integer)

function param_batch!(cvi::DB, data::RealMatrix, labels::IntegerVector)
    cvi.dim, cvi.n_samples = size(data)
    # Take the average across all samples, but cast to 1-D vector
    cvi.mu_data = mean(data, dims=2)[:]
    # u = findfirst.(isequal.(unique(labels)), [labels])
    u = unique(labels)
    cvi.n_clusters = length(u)
    cvi.n = zeros(Integer, cvi.n_clusters)
    cvi.v = zeros(cvi.dim, cvi.n_clusters)
    cvi.CP = zeros(cvi.n_clusters)
    cvi.D = zeros(cvi.n_clusters, cvi.n_clusters)
    cvi.S = zeros(cvi.n_clusters)
    for ix = 1:cvi.n_clusters
        subset = data[:, findall(x->x==u[ix], labels)]
        cvi.n[ix] = size(subset, 2)
        cvi.v[1:cvi.dim, ix] = mean(subset, dims=2)
        diff_x_v = subset - cvi.v[:, ix] * ones(1, cvi.n[ix])
        cvi.CP[ix] = sum(diff_x_v.^2)
        cvi.S[ix] = cvi.CP[ix] / cvi.n[ix]
    end
    for ix = 1 : (cvi.n_clusters - 1)
        for jx = ix + 1 : cvi.n_clusters
            cvi.D[jx, ix] = sum((cvi.v[:, ix] - cvi.v[:, jx]).^2)
        end
    end
    cvi.D = cvi.D + transpose(cvi.D)
end # function param_batch!(cvi::DB, data::RealMatrix, labels::IntegerVector)

function evaluate!(cvi::DB)
    cvi.R = zeros(cvi.n_clusters, cvi.n_clusters)
    for ix = 1:(cvi.n_clusters - 1)
        for jx = ix + 1 : cvi.n_clusters
            cvi.R[ix, jx] = (cvi.S[ix] + cvi.S[jx]) / cvi.D[ix, jx]
        end
    end
    cvi.R = cvi.R + transpose(cvi.R)
    cvi.criterion_value = sum(maximum(cvi.R, dims=2)) / cvi.n_clusters
end # evaluate(cvi::DB)


text:
@testset "MSA Annotations" begin

    pfam = read(joinpath(pwd(), "data", "PF09645_full.stockholm"), Stockholm,
        generatemapping=true, useidcoordinates=true)

    F112_SSV1 = collect(string(".....QTLNSYKMAEIMYKILEKKGELTLEDILAQFEISVPSAYNIQRALKAIC",
        "ERHPDECEVQYKNRKTTFKWIKQEQKEEQKQEQTQDNIAKIFDAQPANFEQTDQGFIKAKQ....."))

    fasta = read(joinpath(pwd(), "data", "Gaoetal2011.fasta"), FASTA, generatemapping=true)

    @testset "Mapping" begin

        @test minimum(getsequencemapping(pfam, 4)) == 3
        @test maximum(getsequencemapping(pfam, "F112_SSV1/3-112")) == 112
        @test getcolumnmapping(pfam) == filter(i->F112_SSV1[i] !== '.', 1:length(F112_SSV1))
        @test length(pfam.annotations.sequences) == 9

        @test getsequencemapping(fasta, 1) == [1, 2, 3, 4, 5, 6]
        @test getsequencemapping(fasta, "SEQ2") == [1, 2, 3, 4, 5, 6]
        @test getcolumnmapping(fasta) == [1, 2, 3, 4, 5, 6]
    end

    @testset "Modifications" begin

        buffer = IOBuffer()

        print(buffer, annotations(pfam))
        printed = split(String(take!(buffer)), '\n')
        @test length(printed) == 17
        @test printed[1] == string("#=GF NCol\t120")
        @test printed[2] ==
            string("#=GF ColMap\t6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22",
                  ",23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,",
                  "43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,",
                  "63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,",
                  "84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,",
                  "103,104,105,106,107,108,109,110,111,112,113,114,115")
        @test occursin(r"MIToS_", printed[3])
        @test occursin(r"MIToS_", printed[4])
    end

    @testset "Delete modification annotations" begin

        buffer = IOBuffer()

        delete_annotated_modifications!(pfam)
        print(buffer, annotations(pfam))
        printed = split(String(take!(buffer)), '\n')
        @test length(printed) == 17 - 2 # 2 MIToS annotations
        @test !occursin(r"MIToS_", printed[3])
        @test !occursin(r"MIToS_", printed[4])
    end
end


text:
"""
Elliptic Curve
==============

References
----------

- https://en.wikipedia.org/wiki/Elliptic_curve
- "An Introduction to the Theory of Elliptic Curves", Joseph H. Silverman
- https://andrea.corbellini.name/2015/05/17/elliptic-curve-cryptography-a-gentle-introduction/
"""
module EllipticCurveCrypto

include("field.jl")

export
    Point, EllipticCurve, is_on_curve,
    ec_zero, ec_add, ec_neg, ec_scalar_mul, ec_scalar_mul_ml

struct Point{T}
    x::T
    y::T
end

# `T` might be mutable (e.g. `BigInt`) and compare using object identity by default.
function Base.:(==)(a::Point{T}, b::Point{T}) where T
    a.x == b.x && a.y == b.y
end

function Base.:(≈)(a::Point{T}, b::Point{T}) where T
    a.x ≈ b.x && a.y ≈ b.y
end

"""
Elliptic curve (with characteristic >= 5) `y ^ 2 = x ^ 3 + a * x + b`.
"""
struct EllipticCurve{T}
    field::Field{T}

    a::T
    b::T

    function EllipticCurve(f::F, a::T, b::T) where {T, F<:Field{T}}
        a = f_from(f, a)
        b = f_from(f, b)
        if f_add(f,
            f_mul(f, f_pow(f, a, 3), f_from(f, T(4))),
            f_mul(f, f_pow(f, b, 2), f_from(f, T(27)))) == f_zero(f)
            error("The curve is singular")
        end
        new{T}(f, a, b)
    end
end

function is_on_curve(ec::EllipticCurve{T}, p::Point{T})::Bool where T
    if p == ec_zero(ec)
        return true
    end

    f = ec.field
    a = ec.a
    b = ec.b
    x = p.x
    y = p.y
    lhs = f_pow(f, y, 2)
    rhs = f_add(f,
        f_add(f, f_pow(f, x, 3), f_mul(f, x, a)),
        b)
    lhs == rhs
end

# Zero, also the infinity.
function ec_zero(ec::EllipticCurve{T})::Point{T} where T
    z = f_zero(ec.field)
    Point(z, z)
end

function ec_add(ec::EllipticCurve{T}, p1::Point{T}, p2::Point{T})::Point{T} where T
    f = ec.field
    z = ec_zero(ec)

    s = if p1 == z
        return p2
    elseif p2 == z
        return p1
    elseif p1.x == p2.x
        if f_add(f, p1.y, p2.y) == f_zero(f)
            return z
        else
            # Doubling
            (x, y) = (p1.x, p1.y)
            numer = f_add(f, f_mul(f, f_from(f, T(3)), f_pow(f, x, 2)), ec.a)
            denom = f_mul(f, f_from(f, T(2)), y)
            f_div(f, numer, denom)
        end
    else
        f_div(f, f_sub(f, p2.y, p1.y), f_sub(f, p2.x, p1.x))
    end

    x = f_sub(f, f_mul(f, s, s), f_add(f, p1.x, p2.x))
    y = f_sub(f, f_mul(f, s, f_sub(f, p1.x, x)), p1.y)
    Point(x, y)
end

function ec_neg(ec::EllipticCurve{T}, p::Point{T})::Point{T} where T
    Point(p.x, f_neg(ec.field, p.y))
end

# Double-and-add algorithm
function ec_scalar_mul(ec::EllipticCurve{T}, p::Point{T}, s::Int)::Point{T} where T
    acc = ec_zero(ec)
    while s != 0
        if s & 1 != 0
            acc = ec_add(ec, acc, p)
        end
        p = ec_add(ec, p, p)
        s >>= 1
    end
    acc
end

# Side-channel attack resistent multiplication.
"Point multiplication using Montgomery ladder."
function ec_scalar_mul_ml(ec::EllipticCurve{T}, p::Point{T}, s::Int)::Point{T} where T
    # p * d = \sum_{i = 0..m} p * (d_i 2^i)          ...0-based
    #       = \sum_{i = 1..m} p * (d_i 2^(i - 1))    ...1-based
    d = digits(s, base=2)
    m = length(d)

    # Invariants:
    #
    #     r0 = p * t0 = p * (s >> (i - 1))
    #     r1 = p * t1 = p * (t0 + 1)
    r0 = ec_zero(ec)
    r1 = p
    for i in m:-1:1
        if d[i] == 0
            r1 = ec_add(ec, r0, r1)  # t0 + (t0 + 1) => (2 * t0) + 1
            r0 = ec_add(ec, r0, r0)  # 2 * t0
        else
            r0 = ec_add(ec, r0, r1)  # t0 + (t0 + 1) => 2 * t0 + 1
            r1 = ec_add(ec, r1, r1)  # 2 * (t0 + 1) => (2 * t0 + 1) + 1
        end
    end
    r0
end

end # module

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
type HullNode
    m::Float64
    b::Float64
    left::Float64
    right::Float64
    pr::Float64

    #empty constructor
    function HullNode()
        return new()
    end
end


text:
# DioOption represents a set of options that is used in DIO,
# a fundamental frequency analysis.
immutable DioOption
    f0floor::Cdouble
    f0ceil::Cdouble
    channels_in_octave::Cdouble
    period::Cdouble # ms
    speed::Cint
    allowed_range::Cdouble # added in v0.2.1-2 (WORLD 0.2.0_2)

    function DioOption(f0floor, f0ceil, channels_in_octave, period, speed,
                       allowed_range = 0.1)
        f0floor > f0ceil && throw(ArgumentError("F0floor must be larger than F0ceil"))
        f0floor < 0 && throw(ArgumentError("f0floor must be positive"))
        if channels_in_octave < 0
            throw(ArgumentError("channels_in_octave must be positive"))
        end
        period <= 0 && throw(ArgumentError("period must be positive"))
        (1 <= speed <= 12) || throw(ArgumentError("1 ≤ speed ≤ 12 is supprted"))
        allowed_range >= 0 || throw(ArgumentError("allowed_range >= 0 is supported"))
        new(f0floor, f0ceil, channels_in_octave, period, speed, allowed_range)
    end
end

immutable CheapTrickOption
    q1::Cdouble
    CheapTrickOption(q1=-0.09) = new(q1)
end

immutable D4COption
    dummy::Cdouble
    D4COption(dummy=0.0) = new(dummy)
end

# Note that the default options assume that the sampling frequency of a input
# speech signal is 44.1 kHz.
function DioOption(;
                   f0floor::Float64=71.0,
                   f0ceil::Float64=800.0,
                   channels_in_octave::Float64=2.0,
                   period::Float64=5.0,
                   speed::Integer=1,
                   allowed_range::Float64=0.1
                   )
    DioOption(f0floor, f0ceil, channels_in_octave, period, speed, allowed_range)
end

function get_samples_for_dio(fs::Real, len::Integer, period::Real)
    ccall((:GetSamplesForDIO, libworld), Cint,
          (Cint, Cint, Cdouble), fs, len, period)
end

function dio(x::StridedVector{Cdouble}, fs::Real, opt::DioOption=DioOption())
    expectedlen = get_samples_for_dio(fs, length(x), opt.period)
    f0 = Array(Cdouble, expectedlen)
    timeaxis = Array(Cdouble, expectedlen)
    # Note that value passinig of julia-type to C-struct doesn't work.
    ccall((:DioByOptPtr, libworld),  Void,
          (Ptr{Cdouble}, Cint, Cint, Ptr{DioOption}, Ptr{Cdouble}, Ptr{Cdouble}),
          x, length(x), fs, &opt, timeaxis, f0)
    f0, timeaxis
end

function dio(x::StridedVector{Cdouble}, fs::Real; opt::DioOption=DioOption())
    dio(x, fs, opt)
end

function stonemask(x::StridedVector{Cdouble}, fs::Integer,
                   timeaxis::StridedVector{Cdouble},
                   f0::StridedVector{Cdouble})
    refinedF0 = Array(Cdouble, length(f0))
    ccall((:StoneMask, libworld),  Void,
          (Ptr{Cdouble}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cint,
           Ptr{Cdouble}),
          x, length(x), fs, timeaxis, f0, length(f0), refinedF0)
    refinedF0
end

# Array{T,2} -> Array{Ptr{T}}
function ptrarray2d!{T<:Real}(dst::Array{Ptr{T},1}, src::Array{T,2})
    for i=1:size(src, 2)
        @inbounds dst[i] = pointer(sub(src, 1:size(src, 1), i), 1)
    end
end

function get_fftsize_for_cheaptrick(fs::Integer)
    fftsize = ccall((:GetFFTSizeForCheapTrick, libworld), Cint, (Cint,), fs)
    convert(Int, fftsize)
end

function cheaptrick(x::StridedVector{Cdouble}, fs::Integer,
                    timeaxis::StridedVector{Cdouble},
                    f0::StridedVector{Cdouble};
                    opt::CheapTrickOption=CheapTrickOption()
    )
    freqbins = get_fftsize_for_cheaptrick(fs)>>1 + 1
    spectrogram = Array(Cdouble, freqbins, length(f0))

    # Array{Cdouble,2} -> Array{Ptr{Cdouble}}
    cspectrogram = Array(Ptr{Cdouble}, size(spectrogram, 2))
    ptrarray2d!(cspectrogram, spectrogram)

    ccall((:CheapTrick, libworld), Void,
          (Ptr{Cdouble}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cint,
           Ptr{CheapTrickOption}, Ptr{Ptr{Cdouble}}),
          x, length(x), fs, timeaxis, f0, length(f0), &opt, cspectrogram)

    # Array{Cdouble,2} <- Array{Ptr{Cdouble}}
    for i=1:length(f0), j=1:freqbins
        @inbounds spectrogram[j,i] = unsafe_load(cspectrogram[i], j)
    end

    spectrogram
end

function d4c(x::StridedVector{Cdouble}, fs::Integer,
             timeaxis::StridedVector{Cdouble},
             f0::StridedVector{Cdouble};
             opt::D4COption=D4COption())
    fftsize = get_fftsize_for_cheaptrick(fs)
    freqbins = fftsize>>1 + 1
    aperiodicity = zeros(Cdouble, freqbins, length(f0))

    # Array{Cdouble,2} -> Array{Ptr{Cdouble}}
    caperiodicity = Array(Ptr{Cdouble}, size(aperiodicity, 2))
    ptrarray2d!(caperiodicity, aperiodicity)

    ccall((:D4C, libworld), Void,
          (Ptr{Cdouble}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Cint,
           Ptr{D4COption}, Ptr{Ptr{Cdouble}}),
          x, length(x), fs, timeaxis, f0, length(f0), fftsize, &opt,
          caperiodicity)

    # Array{Cdouble,2} <- Array{Ptr{Cdouble}}
    for i=1:length(f0), j=1:freqbins
        aperiodicity[j,i] = unsafe_load(caperiodicity[i], j)
    end

    aperiodicity
end

function synthesis(f0::StridedVector{Cdouble},
                   spectrogram::StridedMatrix{Cdouble},
                   aperiodicity::StridedMatrix{Cdouble},
                   period::Real, fs::Integer, len::Integer)
    fftsize = get_fftsize_for_cheaptrick(fs)

    # Array{Cdouble,2} -> Array{Ptr{Cdouble}}
    cspectrogram = Array(Ptr{Cdouble}, size(spectrogram, 2))
    ptrarray2d!(cspectrogram, spectrogram)

    caperiodicity = Array(Ptr{Cdouble}, size(aperiodicity, 2))
    ptrarray2d!(caperiodicity, aperiodicity)

    synthesized = Array(Cdouble, len)
    ccall((:Synthesis, libworld), Void,
          (Ptr{Cdouble}, Cint, Ptr{Ptr{Cdouble}}, Ptr{Ptr{Cdouble}},
           Cint, Cdouble, Cint, Cint, Ptr{Cdouble}),
          f0, length(f0), cspectrogram, caperiodicity, fftsize, period, fs, len,
          synthesized)

    synthesized
end

# matlabfunctions

function interp1!(x::StridedVector{Cdouble},
                  y::StridedVector{Cdouble},
                  xi::StridedVector{Cdouble},
                  yi::StridedVector{Cdouble})
    @assert length(x) == length(y)
    @assert length(xi) == length(yi)
    ccall((:interp1, libworld), Void,
          (Ptr{Cdouble}, Ptr{Cdouble}, Cint, Ptr{Cdouble}, Cint, Ptr{Cdouble}),
          x, y, length(x), xi, length(xi), yi)
    yi
end

function interp1(x::StridedVector{Cdouble},
                 y::StridedVector{Cdouble},
                 xi::StridedVector{Cdouble})
    yi = similar(xi)
    interp1!(x, y, xi, yi)
end


text:
# Program to Swap two numbers without temp

print("Enter first value a : ")
a = parse(Int64,readline())
print("Enter second value b: ")
b = parse(Int64,readline())

print("\n")
print("Before swapping : a = $a and b = $b")

a = a-b
b = a+b
a = b-a

print("\n")
print("After swapping : a = $a and b = $b")

max_stars_repo_name: cashewcodes/Basic-Julia-Projects

max_issues_repo_name: cashewcodes/Basic-Julia-Projects

max_forks_repo_name: cashewcodes/Basic-Julia-Projects


text:
import NCDatasets
using Test

ncfile1 = tempname()
ncfile2 = tempname()
jlfile = tempname()
#jlfile = "/tmp/out.jl"

ds = NCDatasets.NCDataset(ncfile1,"c")
ds.dim["lon"] = 3;
ds.dim["unlimited"] = Inf;
nclon = NCDatasets.defVar(ds,"variable with space", Float32, ("lon",))
nclon.attrib["string"] = "degrees_east";
nclon.attrib["float32"] = Float32(1.)
nclon.attrib["float64"] = 1.
nclon.attrib["float32_vector"] = Float32[1.,2.,3.]
nclon.attrib["float64_vector"] = [1.,2.,3.]
nclon.attrib["int32_vector"] = Int32[1,2,3]

ds.attrib["dollar"] = "a dollar \$ stop";
ds.attrib["backslash"] = "a backslash \\ stop";
ds.attrib["doublequote"] = "a doublequote \" stop";

close(ds)

#ncfile1 = "/tmp/sresa1b_ncar_ccsm3-example.nc"

NCDatasets.ncgen(ncfile1,jlfile; newfname = ncfile2)
include(jlfile)

buf1 = IOBuffer()
buf2 = IOBuffer()

NCDatasets.ncgen(buf1,ncfile1)
NCDatasets.ncgen(buf2,ncfile2)

@test String(take!(buf1)) == String(take!(buf2))

rm(ncfile1)
rm(ncfile2)
rm(jlfile)


text:
module GEMPIC

import Base.Threads: @sync, @spawn, nthreads, threadid

# Utilities
include("mesh.jl")
include("low_level_bsplines.jl")
include("splinepp.jl")
include("distributions.jl")
include("distributions_spin.jl")

# Maxwell solvers
include("maxwell_1d_fem.jl")
include("linear_solver_spline_mass_2d.jl")
include("poisson_2d_fem.jl")
include("maxwell_2d_fem.jl")

# Particle Groups
include("particle_group.jl")

# Particle-Mesh coupling
abstract type AbstractParticleMeshCoupling end
include("particle_mesh_coupling_1d.jl")
include("particle_mesh_coupling_2d.jl")
include("particle_mesh_coupling_spin.jl")

# Particle sampling
include("particle_sampling.jl")
include("particle_sampling_spin.jl")
include("landau_damping.jl")

# Splittings
include("hamiltonian_splitting.jl")
include("hamiltonian_splitting_boris.jl")
include("hamiltonian_splitting_spin.jl")

# Diagnostics
include("diagnostics.jl")
include("diagnostics_spin.jl")

end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct FixedDCSource <: DCSource
        voltage::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int64
        internal::InfrastructureSystemsInternal
    end

Parameters of a Fixed DC Source that returns a fixed DC voltage

# Arguments
- `voltage::Float64`: rated VA, validation range: (0, nothing)
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`
- `n_states::Int64`
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct FixedDCSource <: DCSource
    "rated VA"
    voltage::Float64
    ext::Dict{String, Any}
    states::Vector{Symbol}
    n_states::Int64
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function FixedDCSource(voltage, ext=Dict{String, Any}(), )
    FixedDCSource(voltage, ext, Vector{Symbol}(), 0, InfrastructureSystemsInternal(), )
end

function FixedDCSource(; voltage, ext=Dict{String, Any}(), )
    FixedDCSource(voltage, ext, )
end

# Constructor for demo purposes; non-functional.
function FixedDCSource(::Nothing)
    FixedDCSource(;
        voltage=0,
        ext=Dict{String, Any}(),
    )
end

"""Get FixedDCSource voltage."""
get_voltage(value::FixedDCSource) = value.voltage
"""Get FixedDCSource ext."""
get_ext(value::FixedDCSource) = value.ext
"""Get FixedDCSource states."""
get_states(value::FixedDCSource) = value.states
"""Get FixedDCSource n_states."""
get_n_states(value::FixedDCSource) = value.n_states
"""Get FixedDCSource internal."""
get_internal(value::FixedDCSource) = value.internal

"""Set FixedDCSource voltage."""
set_voltage!(value::FixedDCSource, val::Float64) = value.voltage = val
"""Set FixedDCSource ext."""
set_ext!(value::FixedDCSource, val::Dict{String, Any}) = value.ext = val
"""Set FixedDCSource states."""
set_states!(value::FixedDCSource, val::Vector{Symbol}) = value.states = val
"""Set FixedDCSource n_states."""
set_n_states!(value::FixedDCSource, val::Int64) = value.n_states = val
"""Set FixedDCSource internal."""
set_internal!(value::FixedDCSource, val::InfrastructureSystemsInternal) = value.internal = val


text:
const defaultparser = TensorParser()

"""
    @notensor(block)

Marks a block which should be ignored within an `@tensor` environment. Has no effect outside of `@tensor`.
"""
macro notensor(ex::Expr)
    return esc(ex)
end

"""
    @tensor(block)

Specify one or more tensor operations using Einstein's index notation. Indices can
be chosen to be arbitrary Julia variable names, or integers. When contracting several
tensors together, this will be evaluated as pairwise contractions in left to right
order, unless the so-called NCON style is used (positive integers for contracted
indices and negative indices for open indices).

A second argument to the `@tensor` macro can be provided of the form `order=(...)`, where
the list specifies the contraction indices in the order in which they will be contracted.
"""
macro tensor(ex::Expr)
    return esc(defaultparser(ex))
end

macro tensor(ex::Expr, orderex::Expr)
    parser = TensorParser()
    if !(orderex.head == :(=) && orderex.args[1] == :order &&
            orderex.args[2] isa Expr && orderex.args[2].head == :tuple)
        throw(ArgumentError("unkown first argument in @tensor, should be `order = (...,)`"))
    end
    indexorder = map(normalizeindex, orderex.args[2].args)
    parser = TensorParser()
    parser.contractiontreebuilder = network->indexordertree(network, indexorder)
    return esc(parser(ex))
end

"""
    @tensoropt(optex, block)
    @tensoropt(block)

Specify one or more tensor operations using Einstein's index notation. Indices can
be chosen to be arbitrary Julia variable names, or integers. When contracting several
tensors together, the macro will determine (at compile time) the optimal contraction
order depending on the cost associated to the individual indices. If no `optex` is
provided, all indices are assumed to have an abstract scaling `χ` which is optimized
in the asympotic limit of large `χ`.

The cost can be specified in the following ways:

```julia
@tensoropt (a=>χ,b=>χ^2,c=>2*χ,e=>5) C[a,b,c,d] := A[a,e,c,f,h]*B[f,g,e,b]*C[g,d,h]
# asymptotic cost as specified for listed indices, unlisted indices have cost 1 (any symbol for χ can be used)
@tensoropt (a,b,c,e) C[a,b,c,d] := A[a,e,c,f,h]*B[f,g,e,b]*C[g,d,h]
# asymptotic cost χ for indices a,b,c,e, other indices (d,f) have cost 1
@tensoropt !(a,b,c,e) C[a,b,c,d] := A[a,e,c,f,h]*B[f,g,e,b]*C[g,d,h]
# cost 1 for indices a,b,c,e; other indices (d,f) have asymptotic cost χ
@tensoropt C[a,b,c,d] := A[a,e,c,f,h]*B[f,g,e,b]*C[g,d,h]
# asymptotic cost χ for all indices (a,b,c,d,e,f)
```

Note that `@tensoropt` will optimize any tensor contraction sequence it encounters
in the (block of) expressions. It will however not break apart expressions that have
been explicitly grouped with parenthesis, i.e. in
```julia
@tensoroptC[a,b,c,d] := A[a,e,c,f,h]*(B[f,g,e,b]*C[g,d,h])
```
it will always contract `B` and `C` first. For a single tensor contraction sequence,
the optimal contraction order and associated (asymptotic) cost can be obtained using
`@optimalcontractiontree`.
"""
macro tensoropt(expressions...)
    if length(expressions) == 1
        ex = expressions[1]
        optdict = optdata(ex)
    elseif length(expressions) == 2
        optex = expressions[1]
        ex = expressions[2]
        optdict = optdata(optex, ex)
    end

    parser = TensorParser()
    parser.contractiontreebuilder = network->optimaltree(network, optdict)[1]
    return esc(parser(ex))
end

"""
    @tensoropt_verbose(optex, block)
    @tensoropt_verbose(block)

Similar to `@tensoropt`, but prints information details regarding the optimization process to the standard output.
"""
macro tensoropt_verbose(expressions...)
    if length(expressions) == 1
        ex = expressions[1]
        optdict = optdata(ex)
    elseif length(expressions) == 2
        optex = expressions[1]
        ex = expressions[2]
        optdict = optdata(optex, ex)
    end

    parser = TensorParser()
    parser.contractiontreebuilder = network->optimaltree(network, optdict; verbose = true)[1]
    return esc(parser(ex))
end

macro optimalcontractiontree(expressions...)
    if length(expressions) == 1
        ex = expressions[1]
        optdict = optdata(ex)
    elseif length(expressions) == 2
        optex = expressions[1]
        ex = expressions[2]
        optdict = optdata(optex, ex)
    end

    if isassignment(ex) || isdefinition(ex)
        ex = getrhs(ex)
    end
    if !(ex.head == :call && ex.args[1] == :*)
        error("cannot compute optimal contraction tree for this expression")
    end
    network = [getindices(ex.args[k]) for k = 2:length(ex.args)]
    tree, cost = optimaltree(network, optdict)
    return tree, cost
end

"""
    @ncon(tensorlist, indexlist; order = ..., output = ...)

Contract the tensors in `tensorlist` (of type `Vector` or `Tuple`) according to the network
as specified by `indexlist`. Here, `indexlist` is a list (i.e. a `Vector` or `Tuple`) with
the same length as `tensorlist` whose entries are themselves lists (preferably
`Vector{Int}`) where every integer entry provides a label for corresponding index/dimension
of the corresponding tensor in `tensorlist`. Positive integers are used to label indices
that need to be contracted, and such thus appear in two different entries within
`indexlist`, whereas negative integers are used to label indices of the output tensor, and
should appear only once.

By default, contractions are performed in the order such that the indices being contracted
over are labelled by increasing integers, i.e. first the contraction corresponding to label
`1` is performed. The output tensor had an index order corresponding to decreasing
(negative, so increasing in absolute value) index labels. The keyword arguments `order` and
`output` allow to change these defaults.

The advantage of the macro `@ncon` over the function call `ncon` is that the former automatically generates a unique symbol that hooks into the cache. Furthermore, if `tensorlist` is not just some variable but an actual list (as a tuple with parentheses or a vector with square brackets) at the call site, the `@ncon` macro will scan for conjugation calls, e.g. `conj(A)`, and replace this with just `A` but build a matching list of conjugation flags to be specified to `ncon`. This makes it more convenient to specify tensor conjugation, without paying the cost of actively performing the conjugation beforehand.

See also the function [`ncon`](@ref).
"""
macro ncon(args...)
    if length(args) == 2
        return _nconmacro(args[1], args[2])
    else
        return _nconmacro(args[2], args[3], args[1])
    end
end
function _nconmacro(tensors, indices, kwargs = nothing)
    if !(tensors isa Expr) # there is not much that we can do
        if kwargs === nothing
            ex = Expr(:call, :ncon, tensors, indices,
                        Expr(:call, :fill, false, Expr(:call, :length, tensors)),
                        QuoteNode(gensym()))
        else
            ex = Expr(:call, :ncon, kwargs, tensors, indices,
                        Expr(:call, :fill, false, Expr(:call, :length, tensors)),
                        QuoteNode(gensym()))
        end
        return esc(ex)
    end
    if tensors.head == :vect || tensors.head == :tuple
        tensorargs = tensors.args
    elseif tensors.head == :ref
        tensorargs = tensors.args[2:end]
    else
        throw(ArgumentError("invalid @ncon syntax"))
    end
    conjlist = fill(false, length(tensorargs))
    for i = 1:length(tensorargs)
        if tensorargs[i] isa Expr
            if tensorargs[i].head == :call && tensorargs[i].args[1] == :conj
                tensorargs[i] = tensorargs[i].args[2]
                conjlist[i] = true
            end
        end
    end
    if tensors.head == :ref
        tensorex = Expr(:ref, tensors.args[1], tensorargs...)
    else
        tensorex = Expr(:ref, :Any, tensorargs...)
    end
    if kwargs === nothing
        ex = Expr(:call, :ncon, tensorex, indices, conjlist, QuoteNode(gensym()))
    else
        ex = Expr(:call, :ncon, kwargs, tensorex, indices, conjlist, QuoteNode(gensym()))
    end
    return esc(ex)
end


text:
function makejkl(n)

m = (n-2)*(n-1)*n/6;
indj = zeros(Int64,m);
indk = zeros(Int64,m);
indl = zeros(Int64,m);
count = 0;
for j = 1:(n-2)
    for k = (j+1):(n-1)
        for l = (k+1):n
            count = count+1;
            indj[count] = j;
            indk[count] = k;
            indl[count] = l;
        end
    end
end

return indj, indk, indl

end

function neighquadletlist(pairlist, pairnumsum, ilist)

pairnum = pairnumsum[2:end]-pairnumsum[1:end-1];
quadletnum  = (pairnum-2).*(pairnum-1).*pairnum/6;

inum = length(ilist);
quadletlistj = [];
quadletlistk = [];
quadletlistl = [];

for ii = 1:inum
    n1 = pairnumsum[ii];
    n2 = pairnumsum[ii+1];
    g = pairlist[(n1+1):n2]; # a list of neighbors around atom i

    indj, indk, indl = makejkl(n);
    quadletlistj = [quadletlistj; g[indj]];
    quadletlistk = [quadletlistk; g[indk]];
    quadletlistl = [quadletlistl; g[indl]];
end

quadletlist = [quadletlistj[:] quadletlistk[:] quadletlistl[:]];
quadletnum = [0; cumsum(quadletnum)];

quadletlist = Int64.(quadletlist)
quadletnum  = Int64.(quadletnum)

return quadletlist, quadletnum

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
using Semirings
using Test

include("implementationtest.jl")
include("interfacetest.jl")
include("coretest.jl")


text:
# # Deep Kernel Learning with Flux

## Background

# This example trains a GP whose inputs are passed through a neural network.
# This kind of model has been considered previously [^Calandra] [^Wilson], although it has been shown that some care is needed to avoid substantial overfitting [^Ober].
# In this example we make use of the `FunctionTransform` from [KernelFunctions.jl](github.com/JuliaGaussianProcesses/KernelFunctions.jl/) to put a simple Multi-Layer Perceptron built using Flux.jl inside a standard kernel.

# [^Calandra]: Calandra, R., Peters, J., Rasmussen, C. E., & Deisenroth, M. P. (2016, July). [Manifold Gaussian processes for regression.](https://ieeexplore.ieee.org/abstract/document/7727626) In 2016 International Joint Conference on Neural Networks (IJCNN) (pp. 3338-3345). IEEE.

# [^Wilson]: Wilson, A. G., Hu, Z., Salakhutdinov, R. R., & Xing, E. P. (2016). [Stochastic variational deep kernel learning.](https://proceedings.neurips.cc/paper/2016/hash/bcc0d400288793e8bdcd7c19a8ac0c2b-Abstract.html) Advances in Neural Information Processing Systems, 29.

# [^Ober]: Ober, S. W., Rasmussen, C. E., & van der Wilk, M. (2021, December). [The promises and pitfalls of deep kernel learning.](https://proceedings.mlr.press/v161/ober21a.html) In Uncertainty in Artificial Intelligence (pp. 1206-1216). PMLR.

# ### Package loading
# We use a couple of useful packages to plot and optimize
# the different hyper-parameters
using AbstractGPs
using Distributions
using Flux
using KernelFunctions
using LinearAlgebra
using Plots
default(; legendfontsize=15.0, linewidth=3.0);

# ## Data creation
# We create a simple 1D Problem with very different variations

xmin, xmax = (-3, 3)  # Limits
N = 150
noise_std = 0.01
x_train_vec = rand(Uniform(xmin, xmax), N) # Training dataset
x_train = collect(eachrow(x_train_vec)) # vector-of-vectors for Flux compatibility
target_f(x) = sinc(abs(x)^abs(x)) # We use sinc with a highly varying value
y_train = target_f.(x_train_vec) + randn(N) * noise_std
x_test_vec = range(xmin, xmax; length=200) # Testing dataset
x_test = collect(eachrow(x_test_vec)) # vector-of-vectors for Flux compatibility

plot(xmin:0.01:xmax, target_f; label="ground truth")
scatter!(x_train_vec, y_train; label="training data")

# ## Model definition
# We create a neural net with 2 layers and 10 units each.
# The data is passed through the NN before being used in the kernel.
neuralnet = Chain(Dense(1, 20), Dense(20, 30), Dense(30, 5))

# We use the Squared Exponential Kernel:
k = SqExponentialKernel() ∘ FunctionTransform(neuralnet)

# We now define our model:
gpprior = GP(k)  # GP Prior
fx = AbstractGPs.FiniteGP(gpprior, x_train, noise_std^2)  # Prior at the observations
fp = posterior(fx, y_train)  # Posterior of f given the observations

# This computes the negative log evidence of `y` (the negative log marginal likelihood of
# the neural network parameters), which is going to be used as the objective:
loss(y) = -logpdf(fx, y)

@info "Initial loss = $(loss(y_train))"

# Flux will automatically extract all the parameters of the kernel
ps = Flux.params(k)

# We show the initial prediction with the untrained model
p_init = plot(; title="Loss = $(round(loss(y_train); sigdigits=6))")
plot!(vcat(x_test...), target_f; label="true f")
scatter!(vcat(x_train...), y_train; label="data")
pred_init = marginals(fp(x_test))
plot!(vcat(x_test...), mean.(pred_init); ribbon=std.(pred_init), label="Prediction")

# ## Training
nmax = 200
opt = Flux.ADAM(0.1)

anim = Animation()
for i in 1:nmax
    grads = gradient(ps) do
        loss(y_train)
    end
    Flux.Optimise.update!(opt, ps, grads)

    if i % 10 == 0
        L = loss(y_train)
        @info "iteration $i/$nmax: loss = $L"

        p = plot(; title="Loss[$i/$nmax] = $(round(L; sigdigits=6))")
        plot!(vcat(x_test...), target_f; label="true f")
        scatter!(vcat(x_train...), y_train; label="data")
        pred = marginals(posterior(fx, y_train)(x_test))
        plot!(vcat(x_test...), mean.(pred); ribbon=std.(pred), label="Prediction")
        frame(anim)
        display(p)
    end
end
gif(anim, "train-dkl.gif"; fps=3)
nothing #hide

# ![](train-dkl.gif)

max_stars_repo_name: JuliaGaussianProcesses/AbstractGP

max_issues_repo_name: JuliaGaussianProcesses/AbstractGP

max_forks_repo_name: JuliaGaussianProcesses/AbstractGP


text:
import Base: -, +, *, /, \, ∘, zero, conj, real, imag, abs, abs2

# Identity operator

struct Identity end

(*)(::Identity,s::GridData) = s

### On scalar grid data ####

# Set it to negative of itself
function (-)(p_in::GridData)
  p = deepcopy(p_in)
  p.data .= -p.data
  return p
end

function (-)(p1::T,p2::T) where {T <: GridData}
   return T(p1.data .- p2.data)
 end

function (+)(p1::T,p2::T) where {T <: GridData}
  return T(p1.data .+ p2.data)
end

# Multiply and divide by a constant
function (*)(p::T,c::Number) where {T<:GridData}
  return T(c*p.data)
end

(*)(c::Number,p::T) where {T<:GridData} = *(p,c)

function (/)(p::T,c::Number) where {T<:GridData}
  return T(p.data ./ c)
end

"""
    product!(out::GridData,p::GridData,q::GridData)

Compute the Hadamard (i.e. element by element) product of grid data
`p` and `q` (of the same type) and return the result in `out`.

# Example

```jldoctest
julia> q = Edges(Dual,(8,6));

julia> out = p = deepcopy(q);

julia> q.u[3,2] = 0.3;

julia> p.u[3,2] = 0.2;

julia> product!(out,p,q)
Edges{Dual,8,6,Float64} data
u (in grid orientation)
6×7 Array{Float64,2}:
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
 0.0  0.0  0.06  0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
v (in grid orientation)
5×8 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
```
"""
@inline product!(out::T, p::T, q::T) where {T <: GridData} = (out .= p.*q)

"""
    product(p::Edges/Nodes,q::Edges/Nodes) --> Edges/Nodes

Compute the Hadamard product of edge or nodal (primal or dual) data `p` and `q` and return
the result. This operation can also be carried out with the `∘` operator:

# Example

```jldoctest
julia> q = Edges(Dual,(8,6));

julia> p = deepcopy(q);

julia> q.u[3,2] = 0.3;

julia> p.u[3,2] = 0.2;

julia> p∘q
Edges{Dual,8,6,Float64} data
u (in grid orientation)
6×7 Array{Float64,2}:
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
 0.0  0.0  0.06  0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
v (in grid orientation)
5×8 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
```
"""
@inline product(p::T, q::T) where {T <: GridData} = product!(T(), p, q)

@inline function (∘)(p::GridData{T,NX,NY}, q::GridData) where {T,NX,NY}
  product!(typeof(p)(), p, q)
end

zero(::Type{T}) where {T <: GridData} = T()

#### ON COMPLEX GRID DATA

for f in (:conj,)
    @eval function $f(A::GridData{NX,NY,T}) where {NX,NY,T <: ComplexF64}
        Acopy = deepcopy(A)
        Acopy .= broadcast($f,Acopy)
        return Acopy
    end
end

for f in (:real, :imag, :abs, :abs2)
  @eval function $f(A::GridData{NX,NY,T}) where {NX,NY,T <: ComplexF64}
      Acopy = similar(A,element_type=Float64)
      Acopy .= broadcast($f,A)
      return Acopy
  end
end


text:
# Planet, Reluplex and Reluval

using NeuralVerification
using Test

macro no_error(ex)
    quote
        try $(esc(ex))
            true
        catch
            false
        end
    end
end

at = @__DIR__

# Problem type - input:Hyperrectangle, output:Hyperrectangle
print("###### Problem type - input:Hyperrectangle, output:Hyperrectangle ######\n")
# Small MNIST Network

mnist_small = read_nnet("$at/../examples/networks/mnist_small.nnet")

# entry 23 in MNIST datset
input_center = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,121,254,136,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,13,230,253,248,99,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,118,253,253,225,42,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,61,253,253,253,74,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,32,206,253,253,186,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,211,253,253,239,69,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,254,253,253,133,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,142,255,253,186,8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,149,229,254,207,21,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,54,229,253,254,105,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,152,254,254,213,26,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,112,251,253,253,26,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,29,212,253,250,149,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,36,214,253,253,137,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,75,253,253,253,59,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,93,253,253,189,17,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,224,253,253,84,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,43,235,253,126,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,99,248,253,119,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,225,235,49,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
output_center = [-1311.1257826380004,4633.767704436501,-654.0718535670002,-1325.349417307,1175.2361184373997,-1897.8607293569007,-470.3405972940001,830.8337987382,-377.7467076115001,572.3674015264198]

in_epsilon = 1 #0-255
out_epsilon = 10 #logit domain

input_low = input_center .- in_epsilon
input_high = input_center .+ in_epsilon

output_low = output_center .- out_epsilon
output_high = output_center .+ out_epsilon

inputSet = Hyperrectangle(low=input_low, high=input_high)
outputSet = Hyperrectangle(low=output_low, high=output_high)

print("\n\n################## Small ##################\n")

problem_hyperrect_small = Problem(mnist_small, inputSet, outputSet)
solver_reluVal = ReluVal(max_iter = 2)
print("\nReluval - Small")
timed_result = @timed solve(solver_reluVal, problem_hyperrect_small)
t = timed_result[2]
out = timed_result[1]
print(" - Time: $t s - Output: $out")

print("\nPlanet - Small")
optimizer = GLPKSolverMIP()
solver = Planet(optimizer)
timed_result = @timed solve(solver, problem_hyperrect_small)
t = timed_result[2]
out = timed_result[1]
print(" - Time: $t s - Output: $out")

print("\nReluplex - Small")
solver=Reluplex()
timed_result = @timed solve(solver, problem_hyperrect_small)
t = timed_result[2]
out = timed_result[1]
print(" - Time: $t s - Output: $out")

# Deep MNIST network
mnist_large = read_nnet("$at/../examples/networks/mnist_large.nnet")

# entry 23 in MNIST datset
input_center = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,121,254,136,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,13,230,253,248,99,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,118,253,253,225,42,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,61,253,253,253,74,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,32,206,253,253,186,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,211,253,253,239,69,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,254,253,253,133,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,142,255,253,186,8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,149,229,254,207,21,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,54,229,253,254,105,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,152,254,254,213,26,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,112,251,253,253,26,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,29,212,253,250,149,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,36,214,253,253,137,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,75,253,253,253,59,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,93,253,253,189,17,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,224,253,253,84,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,43,235,253,126,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,99,248,253,119,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,225,235,49,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
output_center = [131.8781755,134.8987015,141.6166255,158.34307,129.8803525,104.8286425,98.64196,133.6358395,131.1716215,105.10621]

in_epsilon = 1 #0-255
out_epsilon = 10 #logit domain

input_low = input_center .- in_epsilon
input_high = input_center .+ in_epsilon

output_low = output_center .- out_epsilon
output_high = output_center .+ out_epsilon

inputSet = Hyperrectangle(low=input_low, high=input_high)
outputSet = Hyperrectangle(low=output_low, high=output_high)

print("\n\n################## Deep ##################\n")
problem_hyperrect_deep = Problem(mnist_large, inputSet, outputSet)
solver_reluVal = ReluVal(max_iter = 2)
print("\nReluval - Deep")
timed_result = @timed solve(solver_reluVal, problem_hyperrect_deep)
t = timed_result[2]
out = timed_result[1]
print(" - Time: $t s - Output: $out")

print("\nPlanet - Deep")
optimizer = GLPKSolverMIP()
solver = Planet(optimizer)
timed_result = @timed solve(solver, problem_hyperrect_deep)
t = timed_result[2]
out = timed_result[1]
print(" - Time: $t s - Output: $out")

print("\nReluplex - Deep")
solver=Reluplex()
#timed_result = @timed solve(solver, problem_hyperrect_deep)
#t = timed_result[2]
#out = timed_result[1]
#print(" - Time: $t s - Output: $out")

# Wide MNIST network

mnist_wide = read_nnet("$at/../examples/networks/mnist-1-100.nnet")
input_center = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,121,254,136,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,13,230,253,248,99,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,118,253,253,225,42,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,61,253,253,253,74,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,32,206,253,253,186,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,211,253,253,239,69,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,254,253,253,133,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,142,255,253,186,8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,149,229,254,207,21,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,54,229,253,254,105,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,152,254,254,213,26,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,112,251,253,253,26,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,29,212,253,250,149,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,36,214,253,253,137,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,75,253,253,253,59,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,93,253,253,189,17,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,224,253,253,84,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,43,235,253,126,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,99,248,253,119,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,225,235,49,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
output_center = [-1219.6781446362043,2713.019467032775,-3826.9579099907123,-3306.084501236369,2641.9250628243003,-1829.0475478735848,-2477.0765438475523,84.03007306734946,-2072.196780418711,145.0876643197917]

in_epsilon = 1 #0-255
out_epsilon = 10 #logit domain

input_low = input_center .- in_epsilon
input_high = input_center .+ in_epsilon

output_low = output_center .- out_epsilon
output_high = output_center .+ out_epsilon

inputSet = Hyperrectangle(low=input_low, high=input_high)
outputSet = Hyperrectangle(low=output_low, high=output_high)

print("\n\n################## Wide ##################\n")

problem_hyperrect_wide = Problem(mnist_wide, inputSet, outputSet)
solver_reluVal = ReluVal(max_iter = 2)
print("\nReluval - Wide")
timed_result = @timed solve(solver_reluVal, problem_hyperrect_wide)
t = timed_result[2]
out = timed_result[1]
print(" - Time: $t s - Output: $out")

print("\nPlanet - Wide")
optimizer = GLPKSolverMIP()
solver = Planet(optimizer)
timed_result = @timed solve(solver, problem_hyperrect_wide)
t = timed_result[2]
out = timed_result[1]
print(" - Time: $t s - Output: $out")

print("\nReluplex - Wide")
solver=Reluplex()
#timed_result = @timed solve(solver, problem_hyperrect_wide)
#t = timed_result[2]
#out = timed_result[1]
#print(" - Time: $t s - Output: $out")

print("\n\n################## Acas ##################\n")

acas_nnet = read_nnet("$at/../examples/networks/ACASXU_run2a_1_1_tiny_4.nnet")

b_upper = [0.58819589, 0.4999999 , -0.4999999, 0.52838384, 0.4]
b_lower = [0.21466922, 0.11140846, -0.4999999, 0.52838384, 0.4]

inputSet = Hyperrectangle(low=b_lower, high=b_upper)

out_lower = [-0.02088955, -0.01863226,  0.01521364, -0.01918539, 0.01563704]
out_upper = [-0.02077696, -0.01852366,  0.01523927, -0.01912529, 0.015789]

outputSet = Hyperrectangle(low=out_lower, high=out_upper)

# TODO: compute output region - need to do some forward passes in the network
# to get a sense of the output ranges.

problem_hyperrectangle_hyperrectangle_acas = Problem(acas_nnet, inputSet, outputSet)
solver_reluVal = ReluVal(max_iter = 2)
print("\nReluval - Acas")
timed_result = @timed solve(solver_reluVal, problem_hyperrectangle_hyperrectangle_acas)
t = timed_result[2]
out = timed_result[1]
print(" - Time: $t s - Output: $out")

print("\nPlanet - Acas")
optimizer = GLPKSolverMIP()
solver = Planet(optimizer)
timed_result = @timed solve(solver, problem_hyperrectangle_hyperrectangle_acas)
t = timed_result[2]
out = timed_result[1]
print(" - Time: $t s - Output: $out")

print("\nReluplex - Acas")
solver=Reluplex()
timed_result = @timed solve(solver, problem_hyperrectangle_hyperrectangle_acas)
t = timed_result[2]
out = timed_result[1]
print(" - Time: $t s - Output: $out")