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.

import Base: +

"""
    +(args...)
    .+(args...)

Elementwise summation of `SymbolicNode`.
"""
function +(x::SymbolicNode, ys::SymbolicNodeOrReal...)
  ret = x
  for y ∈ ys
    if y isa SymbolicNode
      ret = _plus(ret, y)
    else
      ret = _plus_scalar(ret, scalar=MX_float(y))
    end
  end
  ret
end

+(s::Real, x::SymbolicNode, ys::SymbolicNodeOrReal...) = +(x + s, ys...)

broadcasted(::typeof(+), x::SymbolicNode, ys::SymbolicNodeOrReal...) = +(x, ys...)
broadcasted(::typeof(+), s::Real, x::SymbolicNode, ys::SymbolicNodeOrReal...) = +(x + s, ys...)

import Base: -

"""
    -(x, y)
    .-(x, y)

Elementwise substraction of `SymbolicNode`.
Operating with `Real` is available.
"""
x::SymbolicNode - y::SymbolicNode = _minus(x, y)
x::SymbolicNode - s::Real         = _minus_scalar(x,  scalar=MX_float(s))
s::Real         - x::SymbolicNode = _rminus_scalar(x, scalar=MX_float(s))

-(x::SymbolicNode) = 0 - x

broadcasted(::typeof(-), x::SymbolicNode, y::SymbolicNodeOrReal) = x - y
broadcasted(::typeof(-), s::Real, x::SymbolicNode) = s - x

import Base: *

"""
    .*(x, y)

Elementwise multiplication of `SymbolicNode`.
"""
x::SymbolicNode * s::Real = _mul_scalar(x, scalar=MX_float(s))
s::Real * x::SymbolicNode = _mul_scalar(x, scalar=MX_float(s))

function broadcasted(::typeof(*), x::SymbolicNode, ys::SymbolicNodeOrReal...)
  ret = x
  for y in ys
    if y isa SymbolicNode
      ret = _mul(ret, y)
    else
      ret = _mul_scalar(ret, scalar=MX_float(y))
    end
  end
  ret
end

broadcasted(::typeof(*), s::Real, x::SymbolicNode, ys::SymbolicNodeOrReal...) =
  broadcasted(*, x * s, ys...)

import Base: /

"""
    ./(x, y)

* Elementwise dividing a `SymbolicNode` by a scalar or another `SymbolicNode`
of the same shape.

* Elementwise divide a scalar by an `SymbolicNode`.

* Matrix division (solving linear systems) is not implemented yet.
"""
x::SymbolicNode / s::Real = _DivScalar(x, scalar=MX_float(s))

broadcasted(::typeof(/), x::SymbolicNode, y::SymbolicNode) = _div(x, y)
broadcasted(::typeof(/), x::SymbolicNode, s::Real) = _div_scalar(x,  scalar=MX_float(s))
broadcasted(::typeof(/), s::Real, x::SymbolicNode) = _rdiv_scalar(x, scalar=MX_float(s))

import Base: ^

"""
    .^(x, y)

Elementwise power of `SymbolicNode` and `NDArray`.
Operating with `Real` is available.
"""
^

broadcasted(::typeof(^), x::SymbolicNode, y::SymbolicNode) = _power(x, y)
broadcasted(::typeof(^), x::SymbolicNode, s::Real) = _power_scalar(x,  scalar = s)
broadcasted(::typeof(^), s::Real, x::SymbolicNode) = _rpower_scalar(x, scalar = s)
broadcasted(::typeof(Base.literal_pow), ::typeof(^), x::SymbolicNode, ::Val{s}) where {s} =
  _power_scalar(x, scalar = s)

broadcasted(::typeof(^), ::Irrational{:ℯ}, x::SymbolicNode) = exp(x)
broadcasted(::typeof(^), x::SymbolicNode, s::Irrational) =
  _power_scalar(x, scalar=MX_float(s))
broadcasted(::typeof(^), s::Irrational, x::SymbolicNode) =
  _rpower_scalar(x, scalar=MX_float(s))


text:
__precompile__(true)

module GtkReactive

using Gtk, Colors, FixedPointNumbers, Reexport
@reexport using Reactive
using Graphics
using Graphics: set_coordinates, BoundingBox
using IntervalSets, RoundingIntegers
# There's a conflict for width, so we have to scope those calls
import Cairo

using Gtk: GtkWidget
# Constants for event analysis
using Gtk.GConstants.GdkModifierType: SHIFT, CONTROL, MOD1
using Gtk.GConstants.GdkScrollDirection: UP, DOWN, LEFT, RIGHT
using Gtk.GdkEventType: BUTTON_PRESS, DOUBLE_BUTTON_PRESS, BUTTON_RELEASE

# Re-exports
export set_coordinates, BoundingBox, SHIFT, CONTROL, MOD1, UP, DOWN, LEFT, RIGHT,
       BUTTON_PRESS, DOUBLE_BUTTON_PRESS, destroy

## Exports
export slider, button, checkbox, togglebutton, dropdown, textbox, textarea, spinbutton, cyclicspinbutton, progressbar
export label
export canvas, DeviceUnit, UserUnit, XY, MouseButton, MouseScroll, MouseHandler
export player, timewidget, datetimewidget
export signal, widget, frame
# Zoom/pan
export ZoomRegion, zoom, pan_x, pan_y, init_zoom_rubberband, init_zoom_scroll,
       init_pan_scroll, init_pan_drag

# The generic Widget interface
abstract type Widget end

# A widget that gives out a signal of type T
abstract type InputWidget{T}  <: Widget end

"""
    signal(w) -> s

Return the Reactive.jl Signal `s` associated with widget `w`.
"""
signal(w::Widget) = w.signal
signal(x::Signal) = x

"""
    widget(w) -> gtkw::GtkWidget

Return the GtkWidget `gtkw` associated with widget `w`.
"""
widget(w::Widget) = w.widget

Base.push!(w::Widget, val) = push!(signal(w), val)

Base.show(io::IO, w::Widget) = print(io, typeof(widget(w)), " with ", signal(w))
Gtk.destroy(w::Widget) = destroy(widget(w))
Reactive.value(w::Widget) = value(signal(w))
Base.map(f, w::Union{Widget,Signal}, ws::Union{Widget,Signal}...; kwargs...) = map(f, signal(w), map(signal, ws)...; kwargs...)
Base.foreach(f, w::Union{Widget,Signal}, ws::Union{Widget,Signal}...; kwargs...) = foreach(f, signal(w), map(signal, ws)...; kwargs...)

# Define specific widgets
include("widgets.jl")
include("extrawidgets.jl")
include("graphics_interaction.jl")
include("rubberband.jl")

## More convenience functions
# Containers
Gtk.GtkWindow(w::Union{Widget,Canvas}) = GtkWindow(widget(w))
Gtk.GtkFrame(w::Union{Widget,Canvas}) = GtkFrame(widget(w))
Gtk.GtkAspectFrame(w::Union{Widget,Canvas}, args...) =
    GtkAspectFrame(widget(w), args...)

Base.push!(container::Union{Gtk.GtkBin,GtkBox}, child::Widget) =
    push!(container, widget(child))
Base.push!(container::Union{Gtk.GtkBin,GtkBox}, child::Canvas) =
    push!(container, widget(child))

Base.:|>(parent::Gtk.GtkContainer, child::Union{Widget,Canvas}) = push!(parent, child)

widget(c::Canvas) = c.widget

Gtk.set_gtk_property!(w::Union{Widget,Canvas}, key, val) = set_gtk_property!(widget(w), key, val)
Gtk.get_gtk_property(w::Union{Widget,Canvas}, key) = get_gtk_property(widget(w), key)
Gtk.get_gtk_property(w::Union{Widget,Canvas}, key, ::Type{T}) where {T} = get_gtk_property(widget(w), key, T)

Base.unsafe_convert(::Type{Ptr{Gtk.GLib.GObject}}, w::Union{Widget,Canvas}) =
    Base.unsafe_convert(Ptr{Gtk.GLib.GObject}, widget(w))

Graphics.getgc(c::Canvas) = getgc(c.widget)
Graphics.width(c::Canvas) = Graphics.width(c.widget)
Graphics.height(c::Canvas) = height(c.widget)

Graphics.set_coordinates(c::Union{GtkCanvas,Canvas}, device::BoundingBox, user::BoundingBox) =
    set_coordinates(getgc(c), device, user)
Graphics.set_coordinates(c::Union{GtkCanvas,Canvas}, user::BoundingBox) =
    set_coordinates(c, BoundingBox(0, Graphics.width(c), 0, Graphics.height(c)), user)
function Graphics.set_coordinates(c::Union{GraphicsContext,Canvas,GtkCanvas}, zr::ZoomRegion)
    xy = zr.currentview
    bb = BoundingBox(xy)
    set_coordinates(c, bb)
end
function Graphics.set_coordinates(c::Union{Canvas,GtkCanvas}, inds::Tuple{AbstractUnitRange,AbstractUnitRange})
    y, x = inds
    bb = BoundingBox(first(x), last(x), first(y), last(y))
    set_coordinates(c, bb)
end
function Graphics.BoundingBox(xy::XY)
    BoundingBox(minimum(xy.x), maximum(xy.x), minimum(xy.y), maximum(xy.y))
end

function Base.push!(zr::Signal{ZoomRegion{T}}, cv::XY{ClosedInterval{S}}) where {T,S}
    fv = value(zr).fullview
    push!(zr, ZoomRegion{T}(fv, cv))
end

function Base.push!(zr::Signal{ZoomRegion{T}}, inds::Tuple{ClosedInterval,ClosedInterval}) where T
    push!(zr, XY{ClosedInterval{T}}(inds[2], inds[1]))
end

function Base.push!(zr::Signal{ZoomRegion{T}}, inds::Tuple{AbstractUnitRange,AbstractUnitRange}) where T
    push!(zr, ClosedInterval{T}.(inds))
end

Gtk.reveal(c::Canvas, args...) = reveal(c.widget, args...)

const _ref_dict = IdDict{Any, Any}()

"""
    gc_preserve(widget::GtkWidget, obj)

Preserve `obj` until `widget` has been [`destroy`](@ref)ed.
"""
function gc_preserve(widget::Union{GtkWidget,GtkCanvas}, obj)
    _ref_dict[obj] = true
    signal_connect(widget, :destroy) do w
        delete!(_ref_dict, obj)
    end
end

include("deprecations.jl")

end # module

ext: jl


text:
# Scatter plots

import .LibCImPlot: PlotScatter

function PlotScatter(x::AbstractArray{T}, y::AbstractArray{T};
                  count::Integer = min(length(x), length(y)),
                  offset::Integer = 0, stride::Integer = 1,
                  label_id::String = "") where {T<:ImPlotData}

    LibCImPlot.PlotScatter(label_id, x, y, count, offset, stride * sizeof(T))
end

function PlotScatter(x::AbstractArray{T1}, y::AbstractArray{T2}; kwargs...) where {T1<:Real,T2<:Real}
    PlotScatter(promote(x,y)...; kwargs...)
end

function PlotScatter(y::AbstractArray{T}; label_id::String="", count::Integer=length(y),
                  xscale::Real = 1.0, x0::Real = 0,
                  offset::Integer=0, stride::Integer=1) where {T<:ImPlotData}
    LibCImPlot.PlotScatter(label_id, y, count, xscale, x0, offset, stride * sizeof(T))
end

function PlotScatter(x::UnitRange{<:Integer}, y::AbstractArray{T};
                  xscale::Real = 1.0, x0::Real = 0,
                  label_id::String="") where {T<:ImPlotData}

    count::Cint = length(x)
    offset::Cint = x.start >= 1 ? x.start - 1 : throw("Range out of bounds")
    stride = Cint(sizeof(T))
    LibCImPlot.PlotScatter(label_id, y, count, xscale, x0, offset, stride)
end

function PlotScatter(x::StepRange, y::AbstractArray{T};
                     xscale::Real = 1.0, x0::Real = 0,
                     label_id::String="") where {T<:ImPlotData}

    x.stop < 1 && throw("Range out of bounds")
    count::Cint = length(x)
    offset::Cint = x.start >= 1 ? x.start - 1 : throw("Range out of bounds")
    stride = Cint(x.step * sizeof(T))
    LibCImPlot.PlotScatter(label_id, y, count, xscale, x0, offset, stride)
end


text:
# Handers for execute_request and related messages, which are
# the core of the Jupyter protocol: execution of Julia code and
# returning results.

import Base.Libc: flush_cstdio

using Compat
import Compat.String

const text_plain = MIME("text/plain")
const image_svg = MIME("image/svg+xml")
const image_png = MIME("image/png")
const image_jpeg = MIME("image/jpeg")
const text_markdown = MIME("text/markdown")
const text_html = MIME("text/html")
const text_latex = MIME("text/latex") # Jupyter expects this
const text_latex2 = MIME("application/x-latex") # but this is more standard?

include("magics.jl")

# return a String=>Any dictionary to attach as metadata
# in Jupyter display_data and pyout messages
metadata(x) = Dict()

# return a String=>String dictionary of mimetype=>data
# for passing to Jupyter display_data and execute_result messages.
function display_dict(x)
    data = Dict{String,String}("text/plain" => stringmime(text_plain, x))
    if mimewritable(image_svg, x)
        data[string(image_svg)] = stringmime(image_svg, x)
    end
    if mimewritable(image_png, x)
        data[string(image_png)] = stringmime(image_png, x)
    elseif mimewritable(image_jpeg, x) # don't send jpeg if we have png
        data[string(image_jpeg)] = stringmime(image_jpeg, x)
    end
    if mimewritable(text_markdown, x)
        data[string(text_markdown)] = stringmime(text_markdown, x)
    elseif mimewritable(text_html, x)
        data[string(text_html)] = stringmime(text_html, x)
    elseif mimewritable(text_latex, x)
        data[string(text_latex)] = stringmime(text_latex, x)
    elseif mimewritable(text_latex2, x)
        data[string(text_latex)] = stringmime(text_latex2, x)
    end
    return data
end

# queue of objects to display at end of cell execution
const displayqueue = Any[]

# remove x from the display queue
function undisplay(x)
    i = findfirst(displayqueue, x)
    if i > 0
        splice!(displayqueue, i)
    end
    return x
end

function show_bt(io::IO, top_func::Symbol, t, set)
    if VERSION >= v"0.5.0-pre+5636" # julia PR #17570
        # follow PR #17570 code in removing top_func from backtrace
        eval_ind = findlast(addr->Base.REPL.ip_matches_func(addr, top_func), t)
        eval_ind != 0 && (t = t[1:eval_ind-1])
        Base.show_backtrace(io, t)
    else
        Base.show_backtrace(io, top_func, t, set)
    end
end

# return the content of a pyerr message for exception e
function error_content(e, bt=catch_backtrace(); backtrace_top::Symbol=:include_string, msg::AbstractString="")
    tb = map(x->convert(Compat.UTF8String, x), split(sprint(show_bt,
                                        backtrace_top,
                                        bt, 1:typemax(Int)),
                                 "\n", keep=true))
    ename = string(typeof(e))
    evalue = try
        sprint((io, e, bt) -> showerror(io, e, bt, backtrace=false), e, bt)
    catch
        "SYSTEM: show(lasterr) caused an error"
    end
    unshift!(tb, evalue) # fperez says this needs to be in traceback too
    if !isempty(msg)
        unshift!(tb, msg)
    end
    Dict("ename" => ename, "evalue" => evalue,
                 "traceback" => tb)
end

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

# global variable so that display can be done in the correct Msg context
execute_msg = Msg(["julia"], Dict("username"=>"julia", "session"=>"????"), Dict())

function helpcode(code::AbstractString)
    code_ = strip(code)
    # as in base/REPL.jl, special-case keywords so that they parse
    if !haskey(Docs.keywords, Symbol(code_))
        return "Base.Docs.@repl $code_"
    elseif VERSION < v"0.5.0-dev+3831"
        return "eval(:(Base.Docs.@repl \$(symbol(\"$code_\"))))"
    else
        return "eval(:(Base.Docs.@repl \$(Symbol(\"$code_\"))))"
    end
end

# use a global array to accumulate "payloads" for the execute_reply message
const execute_payloads = Dict[]

function execute_request(socket, msg)
    code = msg.content["code"]
    @vprintln("EXECUTING ", code)
    global execute_msg = msg
    global n, In, Out, ans
    silent = msg.content["silent"]
    store_history = get(msg.content, "store_history", !silent)
    empty!(execute_payloads)

    if !silent
        n += 1
        send_ipython(publish[],
                     msg_pub(msg, "execute_input",
                             Dict("execution_count" => n,
                                          "code" => code)))
    end

    silent = silent || ismatch(r";\s*$", code)
    if store_history
        In[n] = code
    end

    # "; ..." cells are interpreted as shell commands for run
    code = replace(code, r"^\s*;.*$",
                   m -> string(replace(m, r"^\s*;", "Base.repl_cmd(`"),
                               "`, STDOUT)"), 0)

    # a cell beginning with "? ..." is interpreted as a help request
    hcode = replace(code, r"^\s*\?", "")
    if hcode != code
        code = helpcode(hcode)
    end

    try
        for hook in preexecute_hooks
            hook()
        end

        #run the code!
        ans = result = ismatch(magics_regex, code) ? magics_help(code) :
            include_string(code, "In[$n]")

        if silent
            result = nothing
        elseif result != nothing
            if store_history
                if result != Out # workaround for Julia #3066
                    Out[n] = result
                end
            end
        end

        user_expressions = Dict()
        for (v,ex) in msg.content["user_expressions"]
            user_expressions[v] = eval(Main,parse(ex))
        end

        for hook in postexecute_hooks
            hook()
        end

        # flush pending stdio
        flush_all()

        undisplay(result) # dequeue if needed, since we display result in pyout
        display() # flush pending display requests

        if result != nothing
            # Work around for Julia issue #265 (see # #7884 for context)
            # We have to explicitly invoke the correct metadata method.
            result_metadata = invoke(metadata, (typeof(result),), result)
            send_ipython(publish[],
                         msg_pub(msg, "execute_result",
                                 Dict("execution_count" => n,
                                              "metadata" => result_metadata,
                                              "data" => display_dict(result))))

        end
        send_ipython(requests[],
                     msg_reply(msg, "execute_reply",
                               Dict("status" => "ok",
                                            "payload" => execute_payloads,
                                            "execution_count" => n,
                                            "user_expressions" => user_expressions)))
        empty!(execute_payloads)
    catch e
        bt = catch_backtrace()
        try
            # flush pending stdio
            flush_all()
            for hook in posterror_hooks
                hook()
            end
        catch
        end
        empty!(displayqueue) # discard pending display requests on an error
        content = error_content(e,bt)
        send_ipython(publish[], msg_pub(msg, "error", content))
        content["status"] = "error"
        content["execution_count"] = n
        send_ipython(requests[], msg_reply(msg, "execute_reply", content))
    end
end


text:
using Test
using LinearAlgebra
using ProximalCore
using ProximalCore: prox, gradient, convex_conjugate
using ProximalCore: Zero, IndZero
import ProximalCore: prox!, is_convex, is_generalized_quadratic

@testset "Package sanity checks" begin
    @test isempty(detect_unbound_args(ProximalCore))
end

struct UnitInfNormBall end

is_convex(::Type{UnitInfNormBall}) = true

function (::UnitInfNormBall)(x)
    R = eltype(y)
    return norm(x, Inf) <= 1 ? R(0) : R(Inf)
end

function prox!(y, ::UnitInfNormBall, x, gamma)
    R = eltype(y)
    y .= max.(R(-1), min.(R(1), x))
    return R(0)
end

@testset "Basic types" begin

    @testset "Zero" begin

        @inferred (f -> Val(is_convex(f)))(Zero())
        @inferred (f -> Val(is_generalized_quadratic(f)))(Zero())

        @test is_convex(Zero())
        @test is_generalized_quadratic(Zero())

        for T in [Float32, Float64]
            @test let x = T[1.0, 2.0, 3.0]
                prox(Zero(), x, T(42)) == (x, T(0))
                gradient(Zero(), x) == (T[0, 0, 0], T(0))
            end
        end

    end

    @testset "IndZero" begin

        @inferred (f -> Val(is_convex(f)))(IndZero())
        @inferred (f -> Val(is_generalized_quadratic(f)))(IndZero())

        @test is_convex(IndZero())
        @test is_generalized_quadratic(IndZero())

        for T in [Float32, Float64]
            @test let x = T[1.0, 2.0, 3.0]
                prox(IndZero(), x, T(42)) == (T[0, 0, 0], T(0))
            end
        end

    end

    @testset "Conjugation" begin

        @test convex_conjugate(Zero()) isa IndZero
        @test convex_conjugate(IndZero()) isa Zero

    end

end

@testset "Custom types" begin

    @testset "UnitInfNormBall" begin
        f = UnitInfNormBall()
        x = Float64[-2, -1, 0, 1, 2]
        y, v = @inferred prox(f, x)

        @test is_convex(f)
        @test !is_generalized_quadratic(f)

        @test typeof(v) == real(eltype(x))
        @test y == Float64[-1, -1, 0, 1, 1]
        @test v == 0

        f_conjugate = convex_conjugate(f)
        y, v = @inferred prox(f_conjugate, x)

        @test is_convex(f_conjugate)
        @test !is_generalized_quadratic(f_conjugate)

        @test typeof(v) == real(eltype(x))
        @test y == Float64[-1, 0, 0, 0, 1]
        @test v == norm(y, 1)

        y, v = @inferred prox(f_conjugate, x, 0.5)

        @test typeof(v) == real(eltype(x))
        @test y == Float64[-1.5, -0.5, 0, 0.5, 1.5]
        @test v == norm(y, 1)

        f_biconjugate = convex_conjugate(f_conjugate)
        y, v = @inferred prox(f_biconjugate, x, 0.5)

        @test is_convex(f_biconjugate)
        @test !is_generalized_quadratic(f_biconjugate)

        @test typeof(v) == real(eltype(x))
        @test y == Float64[-1, -1, 0, 1, 1]
        @test v == 0
    end

end


text:
for (func, typ) in [(:clblasStrmm, cl.CL_float),
                    (:clblasDtrmm, cl.CL_double),
                    (:clblasCtrmm, CL_float2),
                    (:clblasZtrmm, CL_double2)]

    @eval @blasfun $func(order::clblasOrder, side::clblasSide,
                              uplo::clblasUplo, transA::clblasTranspose,
                              diag::clblasDiag,
                              M::Csize_t, N::Csize_t,
                              alpha::$typ,
                              A::cl.CL_mem, offA::Csize_t, lda::Csize_t,
                              B::cl.CL_mem, offB::Csize_t, ldb::Csize_t,
                              n_queues::cl.CL_uint,
                              queues::Ptr{cl.CL_command_queue},
                              n_events_in_wait_list::cl.CL_uint,
                              event_wait_list::Ptr{cl.CL_event},
                              events::Ptr{cl.CL_event})

    @eval @blasfun2 $func(order::clblasOrder, side::clblasSide,
                               uplo::clblasUplo, transA::clblasTranspose,
                               diag::clblasDiag,
                               M::Csize_t, N::Csize_t,
                               alpha::$typ,
                               A::cl.CL_mem, offA::Csize_t, lda::Csize_t,
                               B::cl.CL_mem, offB::Csize_t, ldb::Csize_t)

end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Glibmm_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Glibmm")
JLLWrappers.@generate_main_file("Glibmm", UUID("5d85a9da-21f7-5855-afec-cdc5039c46e8"))
end  # module Glibmm_jll


text:
using Distributions, Parameters, DynamicHMC, LogDensityProblems, TransformVariables
using Random, StatsFuns
import Distributions: pdf, logpdf, rand

export LBA,pdf,logpdf,rand

Base.@kwdef struct LBA{T1,T2,T3,T4} <: ContinuousUnivariateDistribution
    ν::T1
    A::T2
    k::T3
    τ::T4
    σ::Float64 = 1.0
end

Base.broadcastable(x::LBA) = Ref(x)

###
### simulation
###

function selectWinner(dt)
    if any(x -> x > 0,dt)
        mi, mv = 0, Inf
        for (i, t) in enumerate(dt)
            if (t > 0) && (t < mv)
                mi = i
                mv = t
            end
        end
    else
        return 1,-1.0
    end
    return mi,mv
end

function sampleDriftRates(ν,σ)
    noPositive=true
    v = similar(ν)
    while noPositive
        v = [rand(Normal(d,σ)) for d in ν]
        any(x->x>0, v) ? noPositive=false : nothing
    end
    return v
end

function rand(d::LBA)
    @unpack τ,A,k,ν,σ = d
    b=A+k
    N = length(ν)
    v = sampleDriftRates(ν, σ)
    a = rand(Uniform(0, A), N)
    dt = @. (b-a)/v
    choice,mn = selectWinner(dt)
    rt = τ .+ mn
    return choice,rt
end

function rand(d::LBA, N::Int)
    choice = fill(0, N)
    rt = fill(0.0, N)
    for i in 1:N
        choice[i], rt[i] = rand(d)
    end
    return (choice=choice, rt=rt)
end

function simulateLBA(;Nd, v=[1.0,1.5,2.0], A=.8, k=.2, tau=.4, kwargs...)
    return (rand(LBA(ν=v, A=A, k=k, τ=tau), Nd)..., N=Nd, Nc=length(v))
end

###
### log densities
###

function logpdf(d::LBA,data::T) where {T<:NamedTuple}
    return sum(logpdf.(d,data...))
end

logpdf(dist::LBA,data::Array{<:Tuple,1}) = sum(d -> logpdf(dist, d), data)

function logpdf(d::LBA, c, rt)
    @unpack τ,A,k,ν,σ = d
    b = A + k
    logden = 0.0
    rt < τ && return -Inf
    for (i,v) in enumerate(ν)
        if c == i
            logden += log_dens_win(d, v, rt)
        else
            logden += log1mexp(log_dens_lose(d, v, rt))
        end
    end
    return logden - log1mexp(logpnegative(d))
end

logpdf(d::LBA, data::Tuple) = logpdf(d, data...)

function log_dens_win(d::LBA, v, rt)
    @unpack τ,A,k,σ = d
    dt = rt-τ; b=A+k
    n1 = (b-A-dt*v)/(dt*σ)
    n2 = (b-dt*v)/(dt*σ)
    Δcdfs = cdf(Normal(0,1),n2) - cdf(Normal(0,1),n1)
    Δpdfs = pdf(Normal(0,1),n1) - pdf(Normal(0,1),n2)
    println(Δpdfs)
    return -log(A) + logaddexp(log(σ) + log(Δpdfs), log(v) + log(Δcdfs))
end

function log_dens_lose(d::LBA, v, rt)
    @unpack τ,A,k,σ = d
    dt = rt-τ; b=A+k
    n1 = (b-A-dt*v)/(dt*σ)
    n2 = (b-dt*v)/(dt*σ)
    cm = 1 + ((b-A-dt*v)/A)*cdf(Normal(0, 1), n1) -
        ((b-dt*v)/A)*cdf(Normal(0, 1), n2) + ((dt*σ)/A)*pdf(Normal(0, 1), n1) -
        ((dt*σ)/A)*pdf(Normal(0, 1), n2)
    cm = max(cm, 1e-10)
    return log(cm)
end

function logpnegative(d::LBA)
    @unpack ν,σ=d
    sum(v -> logcdf(Normal(0, 1), -v/σ), ν)
end

struct LBAProb{T}
    data::T
    N::Int
    Nc::Int
end

function (problem::LBAProb)(θ)
    @unpack data=problem
    @unpack v,A,k,tau=θ
    d = LBA(ν=v, A=A, k=k, τ=tau)
    minRT = minimum(last, data)
    logprior = (sum(logpdf.(TruncatedNormal(0, 3, 0, Inf), v)) +
                logpdf(TruncatedNormal(.8, .4, 0, Inf) ,A) +
                logpdf(TruncatedNormal(.2, .3, 0, Inf), k) +
                logpdf(TruncatedNormal(.4, .1, 0, minRT), tau))
    loglikelihood = logpdf(d, data)
end

function sampleDHMC(choice, rt, N, Nc, nsamples)
    data = [(c,r) for (c,r) in zip(choice,rt)]
    return sampleDHMC(data, N, Nc, nsamples)
end

Random.seed!(54548)
N = 10
data = simulateLBA(Nd = N)
p = LBAProb(collect(zip(data.choice, data.rt)), N, data.Nc)
p((v=fill(.5, data.Nc),A=.8, k=.2, tau=.4))
minRT = minimum(data.rt)
trans = as((v=as(Array,asℝ₊,data.Nc),A=asℝ₊,k=asℝ₊,tau=as(Real,0,minRT)))
P = TransformedLogDensity(trans, p)
∇P = ADgradient(:ForwardDiff, P)
results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇P, 2000)#; warmup_stages =
    #default_warmup_stages(; local_optimization = nothing, M = DynamicHMC.Symmetric))
posterior = trans.(results.chain)

max_stars_repo_name: itsdfish/Miscellaneous

max_issues_repo_name: itsdfish/Miscellaneous

max_forks_repo_name: itsdfish/Miscellaneous


text:
@testset "Plotting" begin

    fname = joinpath(dirname(@__FILE__), "../data", "test_Hz19.5-testing.bdf")
    s = read_SSR(fname)
    s = rereference(s, "Cz")
    s = filter_highpass(s)
    s = extract_epochs(s)
    s = create_sweeps(s, epochsPerSweep = 2)
    s = ftest(s)

    @testset "Spectrum" begin
        p = plot_spectrum(s, "20Hz_SWN_70dB_R", targetFreq = 3.0)
        #= display(p) =#

        p = plot_spectrum(vec(s.data[:, 1]), samplingrate(Int, s), dBPlot = false)
        #= display(p) =#

        p = plot_spectrum(s, 3, targetFreq = 40.0390625)
        #= display(p) =#
    end

    @testset "Filter reponse" begin

        p = plot_filter_response(s.processing["filter1"], samplingrate(Int, s))
        #= display(p) =#
    end
end


text:
include("Points.jl")
include("Vectors.jl")


text:
using NLopt

println("Use NLOpt to fit the parameter")

obj = build_loss_objective(prob1,Tsit5(),L2Loss(t,data),
                           maxiters=10000,verbose=false)

opt = Opt(:LN_COBYLA, 1)
min_objective!(opt, obj)
(minf,minx,ret) = NLopt.optimize(opt,[1.3])
@test minx[1] ≈ 1.5 atol=1e-3

opt = Opt(:GN_ESCH, 1)
min_objective!(opt, obj.cost_function2)
lower_bounds!(opt,[1.0])
upper_bounds!(opt,[3.0])
xtol_rel!(opt,1e-3)
maxeval!(opt, 10000)
(minf,minx,ret) = NLopt.optimize(opt,[1.3])
@test minx[1] ≈ 1.5 atol=1e-1

opt = Opt(:GN_ISRES, 1)
min_objective!(opt, obj.cost_function2)
lower_bounds!(opt,[1.0])
upper_bounds!(opt,[3.0])
xtol_rel!(opt,1e-3)
maxeval!(opt, 100000)
(minf,minx,ret) = NLopt.optimize(opt,[1.2])
@test minx[1] ≈ 1.5 atol=1e-1

# test differentiation

for autodiff in (false, true)
   global obj = build_loss_objective(prob1, Tsit5(), L2Loss(t,data);
                              mpg_autodiff = autodiff, maxiters = 10000)

   global opt = Opt(:LD_MMA, 1)
   min_objective!(opt, obj.cost_function2)
   xtol_rel!(opt,1e-3)
   maxeval!(opt, 10000)
   global minf,minx,ret
   (minf,minx,ret) = NLopt.optimize(opt, [1.3])
   @test minx[1] ≈ 1.5 atol=5e-1 #take a look at this, it behaves weirdly
end


text:
@testset "AgeCouplingParams" begin
  @testset "agegroup() assings to the right groups" begin

    age_thresholds = Int[0, 5, 12, 18, 30, 40, 50, 60, 70]
    reference_implementation(age) =
      age < 5 ? 1 :
      age < 12 ? 2 :
      age < 18 ? 3 :
      age < 30 ? 4 :
      age < 40 ? 5 :
      age < 50 ? 6 :
      age < 60 ? 7 :
      age < 70 ? 8 :
      9

      for age in 0:100
        @test MocosSim.agegroup(age_thresholds, age) == reference_implementation(age)
      end

  end
end


text:
using UnzipUser
using Documenter

DocMeta.setdocmeta!(UnzipUser, :DocTestSetup, :(using UnzipUser); recursive=true)

makedocs(;
    modules=[UnzipUser],
    authors="Lilith Orion Hafner <60898866+LilithHafner@users.noreply.github.com> and contributors",
    repo="https://github.com/LilithHafner/UnzipUser.jl/blob/{commit}{path}#{line}",
    sitename="UnzipUser.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://LilithHafner.github.io/UnzipUser.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/LilithHafner/UnzipUser.jl",
    devbranch="main",
)


text:
"""
    setup_problem(N, Npairs, pair_order)

Create the list of objects necessary for the PAPA reconstruction.
"""
function setup_problem(N::Int,Npairs::Int,pair_order::Array{Int,2})

    # Dimensions and parameter numbers
    dim = 2^N
    Choi_dim = dim^2
    dim2q = 16
    Nel = convert(Int,dim2q*(dim2q+1)/2)

    # Two-qubit bases
    Bases = zeros(Npairs,dim2q,dim,dim)
    B2 = zeros(1,16,4,4)

    bases_create!(Bases,N,Npairs,pair_order)
    Bases = convert(Array{ComplexF64,4},Bases)

    p_temp = zeros(Int,1,2)
    p_temp[1,:] = [1 2]
    # p_temp = convert(Array{Int,2},p_temp)
    bases_create!(B2,2,1,p_temp)
    B2 = convert(Array{ComplexF64,4},B2)

    # Maximally entangled states
    mrho = zeros(dim^2,dim^2)
    mrho2 = zeros(16,16)
    max_ent_state!(dim;mrho=mrho)
    max_ent_state!(4,mrho=mrho2)

    # Identities
    IN = iden(dim)
    I4 = iden(4)

    # Choi and chi matrices
    Choi = zeros(ComplexF64,Npairs,Choi_dim,Choi_dim)
    Choi2q = zeros(ComplexF64,dim2q,dim2q)
    chitemp = zeros(ComplexF64,dim2q,dim2q)
    rhotemp = zeros(ComplexF64,dim,dim)
    tempE = zeros(ComplexF64,16)
    tempF = zeros(ComplexF64,Nel)

    # Counters and qubit labels
    start = 0
    stop = 0
    qstemp = ones(Int,N)

    return (Bases,B2,N,Npairs,pair_order,mrho,mrho2,IN,I4,Choi_dim,dim2q,Nel,Choi,Choi2q,chitemp,rhotemp,tempE,tempF,start,stop,qstemp)
end

"""
Functionalized version of the full process solver
"""
function process_solver(chilist::Vector{Float64},sigmaVec::Vector{Float64},setup_tuple::Tuple,CP_penalty::Float64)
    return process_solver_full(chilist,sigmaVec,CP_penalty,setup_tuple...)
end

function process_solver!(F::Vector{Float64},chilist::Vector{Float64},sigmaVec::Vector{Float64},setup_tuple::Tuple,CP_penalty::Float64)
    F[:] = process_solver_full(chilist,sigmaVec,CP_penalty,setup_tuple...)
    nothing
end

"""
Full process solver that implements the PAPA_BEAR bootstrapping
"""
function process_solver_full(chilist::Vector{Float64},sigmaVec::Vector{Float64},CP_penalty::Float64,Bases::Array{ComplexF64,4},B2::Array{ComplexF64,4},N::Int,Npairs::Int,pair_order::Array{Int,2},mrho::Array{Float64,2},mrho2::Array{Float64,2},IN::Array{ComplexF64,2},I4::Array{ComplexF64,2},Choi_dim::Int,dim2q::Int,Nel::Int,Choi::Array{ComplexF64,3},Choi2q::Array{ComplexF64,2},chitemp::Array{ComplexF64,2},rhotemp::Array{ComplexF64,2},tempE::Array{ComplexF64,1},tempF::Array{ComplexF64,1},start::Int,stop::Int,qstemp::Vector{Int})::Vector{Float64}

    F = []

    fill!(Choi,0.)

    for nn = 1:1:Npairs
        fill!(chitemp,0.)
        start = 1 + (nn-1)*Nel
        stop = nn*Nel;
        chitemp[triu(trues(dim2q,dim2q))] = chilist[start:1:stop] + 1im*chilist[(start+Npairs*Nel):1:(stop+Npairs*Nel)]
        chitemp[:,:] = chitemp+chitemp'
        chitemp[:,:] = chitemp - diagm(0 => diag(chitemp))/2

        fill!(Choi2q,0.)
        for aa = 1:1:dim2q
            for bb = 1:1:dim2q
                # Choi state for the 2-qubit process on all N qubits
                Choi[nn,:,:] = Choi[nn,:,:] + chitemp[aa,bb]*kron(IN,Bases[nn,aa,:,:])*mrho*kron(IN,Bases[nn,bb,:,:]')

                # Choi state for the 2-qubit process on the 2-qubits only (for constraints use)
                Choi2q[:,:] = Choi2q[:,:] + chitemp[aa,bb]*kron(I4,B2[1,aa,:,:])*mrho2*kron(I4,B2[1,bb,:,:]')
            end
        end

        tempE = eigvals(Choi2q)

        F = vcat(F,[real(4*partial_trace(Choi2q,[0,1],[4,4],2)[:]-I4[:]);imag(4*partial_trace(Choi2q,[0,1],[4,4],2)[:]-I4[:]);CP_penalty*abs(sum(tempE[tempE.<0]))])

    end

    for pp = 1:1:Npairs
        fill!(chitemp,0.)  #ignore the name, used as a temp array

        fill!(qstemp,1)
        qstemp[pair_order[pp,1]] = 0
        qstemp[pair_order[pp,2]] = 0

        # Calculate the reduced process
        for dd = 1:1:dim2q
            fill!(rhotemp,0.) #ignore the name, used as a temp array
            full_process(rhotemp,Bases[pp,dd,:,:],Choi,IN)
            chitemp[:,:] = chitemp[:,:] + kron(B2[1,dd,:,:],partial_trace(rhotemp,qstemp,2*ones(Int,N),N));
        end
        chitemp[:,:] = chitemp[:,:]/2^N

        # Compare to input data
        start = 1 + (pp-1)*Nel
        stop = pp*Nel
        tempF[:] = chitemp[triu(trues(dim2q,dim2q))]

        F = vcat(F,[real(tempF)-sigmaVec[start:1:stop];imag(tempF)-sigmaVec[(start+Npairs*Nel):1:(stop+Npairs*Nel)]])
    end

    return abs.(F)

end


text:
const peri_str = Dict(
            Minute => "minute",
            Day => "day",
            Month => "month",
            Year => "year",
           )
const freq_str = Dict(
            Minute => "minute",
            Day => "daily",
            Week => "weekly",
            Month => "monthly"
           )

"""
ms_since_epoch(t::DateTime)

get the milisecond since UNIX EPOCH of a DateTime object.
"""
ms_since_epoch(t::DateTime) = Dates.value(t) - Dates.UNIXEPOCH

"""
    price_history(ticker; kwargs...)

Get price history of a given ticker, everything but `ticker` have
default value specified by TD api:
https://developer.tdameritrade.com/price-history/apis/get/marketdata/{symbol}/pricehistory

# Optional kwargs:
periodType
period
frequencyType
frequency
endDate
startDate
needExtendedHoursData

# Examples
```julia
using Dates

price_history("GE", Minute(1), Day(1)) #a tick = a day for the past day
price_history("GE", Day(1), Year(2)) #a tick = a day for the past 2 years
price_history("AAPL", Day(1), now()-Day(1), now())
```
"""
function price_history(ticker; kwargs...)
    kwargs = (kwargs...,
               apikey = AUTH_KEY.CONSUMER_KEY,
               Authorization = "Bearer "*AUTH_KEY.ACCESS_TOKEN
          )
    uri = construct_api("marketdata/$ticker/pricehistory", kwargs)
    return @pipe HTTP.get(uri; retries = 1).body |> JSON3.read |> _[:candles]
end

"""
price_history(ticker, freq, peri)

see TD API for default and compatible frequency vs. tickers
"""
function price_history(ticker, freq, peri)
    Dates.toms(freq) > Dates.toms(peri) && ArgumentError("Period must be larger than frequency")

    frequencyType = freq_str[typeof(freq)]
    frequency = freq.value

    periodType = peri_str[typeof(peri)]
    period = peri.value

    price_history(ticker; periodType=periodType,
                        period=period,
                        frequencyType=frequencyType,
                        frequency=frequency)
end

"""
price_history(ticker, freq, start::DateTime=now()-Day(1), stop::DateTime=now())

get the price history with `interval` ticks from `start` to `stop`

valid frequency:
minute: 1*, 5, 10, 15, 30
"""
function price_history(ticker, freq, start::DateTime=now()-Day(1), stop::DateTime=now())
    @assert stop > start
    startDate, endDate = ms_since_epoch.((start, stop))
    frequencyType = freq_str[typeof(freq)]
    frequency = freq.value
    price_history(ticker; startDate=startDate,
                endDate=endDate,
                frequencyType=frequencyType,
                frequency=frequency)
end


text:
@inline function two_inv(b::T) where {T<:FloatWithFMA}
     hi = inv(b)
     lo = fma(-hi, b, one(T))
     lo /= b
     return hi, lo
end

@inline function two_dvi(a::T, b::T) where {T<:FloatWithFMA}
     hi = a / b
     lo = fma(-hi, b, a)
     lo /= b
     return hi, lo
end

@inline function two_sqrt(a::T) where {T<:FloatWithFMA}
    hi = sqrt(a)
    lo = fma(-hi, hi, a)
    lo /= 2
    lo /= hi
    return hi, lo
end

function two_cbrt(a::T) where {T<:FloatWithFMA}
     return HILO(DoubleFloat{T}(a, zero(T))^(one(T)/3))
end


text:
# RK4 is a function which performs one step of the Runge-Kutta 4 ray tracing
#
# EndPoint = RK4(StartPoint, GradientVolume, StepSize);
#
# inputs :
#      StartPoint: 2D or 3D location in vectorfield
#      GradientVolume: Vectorfield
#      Stepsize : The stepsize
#
# outputs :
#      EndPoint : The new location (zero if outside image)
#
# Function is written by D.Kroon University of Twente (July 2008)
# Port from C to Julia by Lars Hellemo (May 2014)

#function mindex2(x, y, sizx)
#    return y*sizx+x
#end

function mindex3(x::Int, y::Int, z::Int, sizx::Int, sizy::Int)
    return z*sizy*sizx+y*sizx+x
end

function checkBounds2d( point::Array{Float64,1}, Isize::Array{Int,1})
    if (point[1]<1)||(point[2]<1)||(point[1]>(Isize[1]))||(point[2]>(Isize[2]))
        return false
    end
    return true
end

# TODO check dimensions of point, seems to be bogus
function checkBounds3d( point::Array{Float64,3}, Isize::Array{Int,1})
    if (point[1,1,1]<1)||(point[2,1,1]<1)||(point[3,1,1]<1)||(point[1,1,1]>(Isize[1]))||(point[2,1,1]>(Isize[2]))||(point[3,1,1]>(Isize[3]))
         return false
    end
    return true
end

# TODO: Replace with (wrapped) built in norm?
function norm2(a::Array{Float64,1})
    return sqrt(a[1]*a[1]+a[2]*a[2])
end

# TODO: Replace with (wrapped) built in norm?
function norm3(a::Array{Float64,1})
    return sqrt(a[1]*a[1]+a[2]*a[2]+a[3]*a[3])
end

#function interpgrad2d(k1, gradientArray, gradientArraySize, startPoint,tmp)
#    tmp = e1(copy(startPoint),gradientArray,0.5)
#    k1[1] = tmp[1]
#    k1[2] = tmp[2]
#end

function interpgrad2d(Ireturn::Array{Float64,1},I::Array{Float64,3}, Isize::Array{Int64,1}, point::Array{Float64,1})
    #  Linear interpolation variables

    perc = zeros(4,)
    #index = [1::Int64 for i = 1:4]
    index = fill(4,4,1)
    fTlocalx = floor(point[1])
    fTlocaly = floor(point[2])
    xBas0 = Int(fTlocalx)
    yBas0 = Int(fTlocaly)
    xBas1=xBas0+1
    yBas1=yBas0+1

    # Linear interpolation constants (percentages)
    xCom=point[1]-fTlocalx
    yCom=point[2]-fTlocaly
    xComi=(1-xCom)
    yComi=(1-yCom)
    perc[1]=xComi * yComi
    perc[2]=xComi * yCom
    perc[3]=xCom * yComi
    perc[4]=xCom * yCom

    # Stick to boundary
    if xBas0<1
        xBas0=1
        if xBas1<1
            xBas1=1
        end
    end

    if yBas0<1
        yBas0=1
        if yBas1<1
            yBas1=1
        end
    end
    if xBas1>Isize[1]
        xBas1=Isize[1]
        if xBas0>Isize[1]
            xBas0=Isize[1]
        end
    end

    if yBas1>Isize[2]
        yBas1=Isize[2]
        if yBas0>Isize[2]
            yBas0=Isize[2]
        end
    end

    # Get the neighbour intensities
    index[1] = LinearIndices(size(I))[xBas0,yBas0,1] #TODO Why is this dim 3,3,2?
    index[2] = LinearIndices(size(I))[xBas0,yBas1,1]
    index[3] = LinearIndices(size(I))[xBas1,yBas0,1]
    index[4] = LinearIndices(size(I))[xBas1,yBas1,1]

   f=Isize[1]*Isize[2]

    # the interpolated color
    Ireturn[1]=I[index[1]]*perc[1]+I[index[2]]*perc[2]+I[index[3]]*perc[3]+I[index[4]]*perc[4]
    Ireturn[2]=I[index[1]+f]*perc[1]+I[index[2]+f]*perc[2]+I[index[3]+f]*perc[3]+I[index[4]+f]*perc[4]

    return Ireturn
end
# TODO: Port 3D
# __inline void interpgrad3d(double *Ireturn, double *I, int *Isize, double *point) {
#     /*  Linear interpolation variables */
#     int xBas0, xBas1, yBas0, yBas1, zBas0, zBas1;
#     double perc[8];
#     double xCom, yCom, zCom;
#     double xComi, yComi, zComi;
#     double fTlocalx, fTlocaly, fTlocalz;
#     int f0, f1;
#     int index[8];
#     double temp;

#     fTlocalx = floor(point[0]); fTlocaly = floor(point[1]); fTlocalz = floor(point[2]);
#     xBas0=(int) fTlocalx; yBas0=(int) fTlocaly; zBas0=(int) fTlocalz;
#     xBas1=xBas0+1; yBas1=yBas0+1; zBas1=zBas0+1;

#     /* Linear interpolation constants (percentages) */
#     xCom=point[0]-fTlocalx;  yCom=point[1]-fTlocaly;   zCom=point[2]-fTlocalz;
#     xComi=(1-xCom); yComi=(1-yCom); zComi=(1-zCom);
#     perc[0]=xComi * yComi; perc[1]=perc[0] * zCom; perc[0]=perc[0] * zComi;
#     perc[2]=xComi * yCom;  perc[3]=perc[2] * zCom; perc[2]=perc[2] * zComi;
#     perc[4]=xCom * yComi;  perc[5]=perc[4] * zCom; perc[4]=perc[4] * zComi;
#     perc[6]=xCom * yCom;   perc[7]=perc[6] * zCom; perc[6]=perc[6] * zComi;

#     /* Stick to boundary */
#     if(xBas0<0) { xBas0=0; if(xBas1<0) { xBas1=0; }}
#     if(yBas0<0) { yBas0=0; if(yBas1<0) { yBas1=0; }}
#     if(zBas0<0) { zBas0=0; if(zBas1<0) { zBas1=0; }}
#     if(xBas1>(Isize[0]-1)) { xBas1=Isize[0]-1; if(xBas0>(Isize[0]-1)) { xBas0=Isize[0]-1; }}
#     if(yBas1>(Isize[1]-1)) { yBas1=Isize[1]-1; if(yBas0>(Isize[1]-1)) { yBas0=Isize[1]-1; }}
#     if(zBas1>(Isize[2]-1)) { zBas1=Isize[2]-1; if(zBas0>(Isize[2]-1)) { zBas0=Isize[2]-1; }}

#    /*Get the neighbour intensities */
#     index[0]=mindex3(xBas0, yBas0, zBas0, Isize[0], Isize[1]);
#     index[1]=mindex3(xBas0, yBas0, zBas1, Isize[0], Isize[1]);
#     index[2]=mindex3(xBas0, yBas1, zBas0, Isize[0], Isize[1]);
#     index[3]=mindex3(xBas0, yBas1, zBas1, Isize[0], Isize[1]);
#     index[4]=mindex3(xBas1, yBas0, zBas0, Isize[0], Isize[1]);
#     index[5]=mindex3(xBas1, yBas0, zBas1, Isize[0], Isize[1]);
#     index[6]=mindex3(xBas1, yBas1, zBas0, Isize[0], Isize[1]);
#     index[7]=mindex3(xBas1, yBas1, zBas1, Isize[0], Isize[1]);
#     f0=Isize[0]*Isize[1]*Isize[2];
#     f1=f0+f0;

#    /*the interpolated color */
#     temp=I[index[0]]*perc[0]+I[index[1]]*perc[1]+I[index[2]]*perc[2]+I[index[3]]*perc[3];
#     Ireturn[0]=temp+I[index[4]]*perc[4]+I[index[5]]*perc[5]+I[index[6]]*perc[6]+I[index[7]]*perc[7];
#     temp=I[index[0]+f0]*perc[0]+I[index[1]+f0]*perc[1]+I[index[2]+f0]*perc[2]+I[index[3]+f0]*perc[3];
#     Ireturn[1]=temp+I[index[4]+f0]*perc[4]+I[index[5]+f0]*perc[5]+I[index[6]+f0]*perc[6]+I[index[7]+f0]*perc[7];
#     temp=I[index[0]+f1]*perc[0]+I[index[1]+f1]*perc[1]+I[index[2]+f1]*perc[2]+I[index[3]+f1]*perc[3];
#     Ireturn[2]=temp+I[index[4]+f1]*perc[4]+I[index[5]+f1]*perc[5]+I[index[6]+f1]*perc[6]+I[index[7]+f1]*perc[7];
# }
function rk4(startPoint::Array{Float64,1}, GradientVolume::Array{Float64,3}, Stepsize::Float64)

    # Perform the RK4 raytracing step

    #gradientArraySize[0]=gradientArraySizeC[0];
    #    gradientArraySize[1]=gradientArraySizeC[1];

   #gradientArray = copy(GradientVolume)
   gradientArray = GradientVolume # TODO Test if this works (without copy)
    gradientArraySize = [size(gradientArray)...]
    nextPoint = startPoint

    startPoint1 = startPoint
#    startPoint1[1]=startPoint[1]-1.0
#    startPoint1[2]=startPoint[2]-1.0
#    if RK4STEP_2D(gradientArray, gradientArraySize, startPoint1, nextPoint, Stepsize)
    nextPoint =  RK4STEP_2D(gradientArray, startPoint1, ones(2,), Stepsize)
 #       nextPoint[1]=nextPoint[1]+1.0
 #       nextPoint[2]=nextPoint[2]+1.0
#    else
#        nextPoint[1]=1.
#        nextPoint[2]=1.
 #   end
    # TODO: 3D
    # else if(PointLength==3) {
    #     gradientArraySize[0]=gradientArraySizeC[0];
    #     gradientArraySize[1]=gradientArraySizeC[1];
    #     gradientArraySize[2]=gradientArraySizeC[2];
    #     startPoint1[0]=startPoint[0]-1.0;
    #     startPoint1[1]=startPoint[1]-1.0;
    #     startPoint1[2]=startPoint[2]-1.0;
    #     if(RK4STEP_3D(gradientArray, gradientArraySize, startPoint1, nextPoint, stepSize)) {
    #         nextPoint[0]=nextPoint[0]+1.0;
    #         nextPoint[1]=nextPoint[1]+1.0;
    #         nextPoint[2]=nextPoint[2]+1.0;
    #     }
    #     else {
    #         nextPoint[0]=0; nextPoint[1]=0; nextPoint[2]=0;
    #     }

    # }
    # else {
    #     mexErrMsgTxt("Starting Point must be 2D or 3D");
    # }

   return nextPoint
end

function RK4STEP_2D(gradientArray::Array{Float64,3}, startPoint::Array{Float64,1}, nextPoint::Array{Float64,1},stepSize::Float64)

    # Perform one step of the RK4 algorithm
    k1 = ones(2,)
    k2 = ones(2,)
    k3 = ones(2,)
    k4 = ones(2,)
    tempPoint = ones(2,)

    gradientArraySize = [size(gradientArray)...]

    # Calculate k1
    k1 = interpgrad2d(k1, gradientArray, gradientArraySize, startPoint)
    tempnorm=max(norm2(k1),1e-6)
    k1[1] = k1[1]*stepSize/tempnorm
    k1[2] = k1[2]*stepSize/tempnorm

    tempPoint[1]=startPoint[1] - k1[1]*0.5
    tempPoint[2]=startPoint[2] - k1[2]*0.5

    # Check if still inside the domain
    if !checkBounds2d([tempPoint...], gradientArraySize)
        return ones(2,)
    end

    # Calculate k2
    k2 = interpgrad2d(k2, gradientArray, gradientArraySize, tempPoint)
    tempnorm=max(norm2(k2),1e-6)
    k2[1] = k2[1]*stepSize/tempnorm
    k2[2] = k2[2]*stepSize/tempnorm

    tempPoint[1]=startPoint[1] - k2[1]*0.5
    tempPoint[2]=startPoint[2] - k2[2]*0.5

    # Check the if are still inside the domain
    if !checkBounds2d([tempPoint...], gradientArraySize)
        return ones(2,)
    end

    # Calculate k3
    k3 = interpgrad2d(k3, gradientArray, gradientArraySize, tempPoint)
    tempnorm=max(norm2(k3),1e-6)
    k3[1] = k3[1]*stepSize/tempnorm
    k3[2] = k3[2]*stepSize/tempnorm

    tempPoint[1]=startPoint[1] - k3[1]
    tempPoint[2]=startPoint[2] - k3[2]

    # Check the if are still inside the domain
    if !checkBounds2d([tempPoint...], gradientArraySize)
        return ones(2,)
    end

    # Calculate k4
    k4 = interpgrad2d(k4, gradientArray, gradientArraySize, tempPoint)
    tempnorm=max(norm2(k4),1e-6)
    k4[1] = k4[1]*stepSize/tempnorm
    k4[2] = k4[2]*stepSize/tempnorm

    # Calculate final point
    nextPoint[1] = startPoint[1] - (k1[1] + k2[1]*2.0 + k3[1]*2.0 + k4[1])/6.0
    nextPoint[2] = startPoint[2] - (k1[2] + k2[2]*2.0 + k3[2]*2.0 + k4[2])/6.0

    # Set step to step size
    # /*
    # D[0]=(nextPoint[0]-startPoint[0]);
    # D[1]=(nextPoint[1]-startPoint[1]);
    # dl=stepSize/(sqrt(D[0]*D[0]+D[1]*D[1])+1e-15);
    # D[0]*=dl; D[1]*=dl;
    # nextPoint[0]=startPoint[0]+D[0];
    # nextPoint[1]=startPoint[1]+D[1];
    #*/

    # Check the if are still inside the domain
    if !checkBounds2d([nextPoint...], gradientArraySize)
        return ones(2,)
    end

    return nextPoint
end


text:
"""
    MeanSquareError(NeuralNetwork)

Sum of square error loss function.

    Usage:
```julia-repl
MeanSquareError(NeuralNetwork)
```
Arguments:
- `NeuralNetwork` : The NeuralNetwork structure.
"""
function MeanSquareError(NeuralNetwork)
	ΣSE = 0.0

	for i in 1:size(NeuralNetwork.y, 1)
		ΣSE += (NeuralNetwork.y[i] - NeuralNetwork.output[i]) ^ 2.0
    end

    loss = 1.0 / size(NeuralNetwork.y) * ΣSE
    return loss
end


text:
module Parse

using RBNF
using RBNF: Token

using ..Types

struct QASMLang end

second((a, b)) = b
second(vec::V) where {V<:AbstractArray} = vec[2]

# roses are red
# violets are blue
# pirates are good
RBNF.crate(::Type{Symbol}) = gensym(:qasm)
RBNF.crate(::Type{VersionNumber}) = VersionNumber("0.0.0")

Base.convert(::Type{VersionNumber}, t::Token) = VersionNumber(t.str)
Base.convert(::Type{String}, t::Token) = t.str
Base.convert(::Type{Int}, t::Token{:int}) = Base.parse(Int, t.str)
Base.convert(::Type{Float64}, t::Token{:float64}) = Base.parse(Float64, t.str)
Base.convert(::Type{Symbol}, t::Token{:id}) = Symbol(t.str)
Base.convert(::Type{Symbol}, t::Token{:reserved}) = Symbol(t.str)
Base.convert(::Type{String}, t::Token{:str}) = String(t.str[2:end-1])

RBNF.@parser QASMLang begin
    # define ignorances
    ignore{space, comment}

    @grammar
    # define grammars
    mainprogram::MainProgram := ["OPENQASM", version = float64, ';', prog = program]
    program = statement{*}
    statement = (regdecl | gate | opaque | qop | ifstmt | barrier | inc)
    # stmts
    ifstmt::IfStmt := [:if, '(', left = id, :(==), right = int, ')', body = qop]
    opaque::Opaque := [:opaque, name = id, ['(', [cargs = idlist].?, ')'].?, qargs = idlist, ';']
    barrier::Barrier := [:barrier, qargs = bitlist, ';']
    regdecl::RegDecl := [type = :qreg | :creg, name = id, '[', size = int, ']', ';']
    inc::Include := [:include, file = str, ';']
    # gate
    gate::Gate := [decl = gatedecl, [body = goplist].?, '}']
    gatedecl::GateDecl := [:gate, name = id, ['(', [cargs = idlist].?, ')'].?, qargs = idlist, '{']

    goplist = (uop | barrier){*}

    # qop
    qop = (uop | measure | reset)
    reset::Reset := [:reset, qarg = bit, ';']
    measure::Measure := [:measure, qarg = bit, :(->), carg = bit, ';']

    uop = (inst | ugate | csemantic_gate)
    inst::Instruction := [name = id, ['(', [cargs = explist].?, ')'].?, qargs = bitlist, ';']
    ugate::UGate := [:U, '(', z1 = exp, ',', y = exp, ',', z2 = exp, ')', qarg = bit, ';']
    csemantic_gate::CXGate := [:CX, ctrl = bit, ',', qarg = bit, ';']

    idlist = @direct_recur begin
        init = [id]
        prefix = [recur..., (',', id) % second]
    end

    bit::Bit := [name = id, ['[', address = int, ']'].?]
    bitlist = @direct_recur begin
        init = [bit]
        prefix = [recur..., (',', bit) % second]
    end

    explist = @direct_recur begin
        init = [exp]
        prefix = [recur..., (',', exp) % second]
    end

    con = (float64 | int | :pi | id | call)
    num = (['(', exp, ')'] % second) | neg | con
    add = ( :+ | :- )
    mul = ( :* | :/ )
    exp = @direct_recur begin
        init = term
        prefix = (recur, add, term)
    end
    term = @direct_recur begin
        init = num
        prefix = (recur, mul, term)
    end
    # term = (add | sub | num)
    neg::Neg := [:-, val = num]
    call::Call := [name=fn, "(", args = exp, ")"]
    fn = (:sin | :cos | :tan | :exp | :ln | :sqrt)
    # binop = (:+ | :- | :* | :/)

    # define tokens
    @token
    id := r"\G[a-z]{1}[A-Za-z0-9_]*"
    float64 := r"\G([0-9]+\.[0-9]*|[0-9]*\.[0.9]+)([eE][-+]?[0-9]+)?"
    int := r"\G([1-9]+[0-9]*|0)"
    space := r"\G\s+"
    comment := r"\G//.*"
    str := @quote ("\"", "\\\"", "\"")
end

end # Parse


text:
using MultivariateStats
using Base.Test

srand(34568)

## LinearDiscriminant

@assert LinearDiscriminant <: Discriminant

w = [1., 2., 3., 4., 5.]
b = 2.5
x = [4., 5., 2., 3., 1.]

f = LinearDiscriminant(w, b)
@test length(f) == 5
@test evaluate(f, x) == 39.5
@test evaluate(f, -x) == -34.5

@test predict(f, x) == true
@test predict(f, -x) == false

X = rand(5, 8)
Y = evaluate(f, X)
@test size(Y) == (8,)
for i = 1:8
	@test_approx_eq Y[i] evaluate(f, X[:,i])
end
@test predict(f, X) == (Y .> 0)

## prepare data

t1 = deg2rad(30)
t2 = deg2rad(75)

R1 = [cos(t1) -sin(t1); sin(t1) cos(t1)]
R2 = [cos(t2) -sin(t2); sin(t2) cos(t2)]

n = 20
Xp = Diagonal([1.2, 3.6]) * randn(2, n) .+ [1.0, -3.0]
Xn = Diagonal([2.8, 1.8]) * randn(2, n) .+ [-5.0, 2.0]

up = vec(mean(Xp, 2))
un = vec(mean(Xn, 2))
if VERSION < v"0.5.0-dev+660"
    Cp = cov(Xp; vardim=2)
    Cn = cov(Xn; vardim=2)
else
    Cp = cov(Xp, 2)
    Cn = cov(Xn, 2)
end
C = 0.5 * (Cp + Cn)

w_gt = C \ (up - un)
w_gt .*= (2 / (dot(w_gt, up) - dot(w_gt, un)))
b_gt = 1.0 - dot(w_gt, up)

@test_approx_eq dot(w_gt, up) + b_gt 1.0
@test_approx_eq dot(w_gt, un) + b_gt -1.0

## LDA

f = ldacov(C, up, un)
@test_approx_eq f.w w_gt
@test_approx_eq f.b b_gt

f = ldacov(Cp, Cn, up, un)
@test_approx_eq f.w w_gt
@test_approx_eq f.b b_gt

f = fit(LinearDiscriminant, Xp, Xn)
@test_approx_eq f.w w_gt
@test_approx_eq f.b b_gt


text:
# println("loading packages... ")
#
# using OrthographicProjection
#
# println("packages OK")
#
# function main()
#
# 	txtpotreedirs = raw"C:\Users\marte\Documents\potreeDirectory\pointclouds\CAVA"
# 	project_name = "TEST_SLICES"
# 	output_folder = raw"C:\Users\marte\Documents\Julia_package\package_test\TEST\SCRIPT_SLICING"
# 	bbin = "458117.68 4493768.53 196.68 458452.43 4494171.78 237.49"
# 	steps = nothing
# 	step = 5.
# 	n = 20
# 	p1_ = "458145.180 4493834.030 224.250"
# 	p2_ = "458255.180 4493775.530 226.050"
# 	axis_ = "0 0 1"
# 	thickness = 1.0
#
# 	b = tryparse.(Float64,split(bbin, " "))
# 	if length(b) == 6
# 		bbin = OrthographicProjection.AABB(b[4],b[1],b[5],b[2],b[6],b[3])
# 	end
#
# 	p1 = tryparse.(Float64,split(p1_, " "))
# 	@assert length(p1) == 3 "a 3D point needed"
# 	p2 = tryparse.(Float64,split(p2_, " "))
# 	@assert length(p2) == 3 "a 3D point needed"
# 	axis_y = tryparse.(Float64,split(axis_, " "))
# 	@assert length(axis_y) == 3 "a 3D axis needed"
#
# 	if !isnothing(step) && !isnothing(n)
# 		steps = fill(step,n)
# 	elseif !isnothing(steps)
# 		steps = tryparse.(Float64,split(steps, " "))
# 	else
# 		steps = Float64[]
# 	end
#
# 	prepend!(steps,0.0)
#
#
# 	proj_folder, plane, model = OrthographicProjection.preprocess(project_name, output_folder, bbin, p1, p2, axis_y, thickness)
# 	OrthographicProjection.get_parallel_sections(txtpotreedirs, project_name, proj_folder, bbin, steps, plane, model)
#
#
# end
#
# @time main()


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct DataDisk <: SwaggerModel
    lun::Any # spec type: Union{ Nothing, Int32 } # spec name: lun
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    vhd::Any # spec type: Union{ Nothing, VirtualHardDisk } # spec name: vhd
    image::Any # spec type: Union{ Nothing, VirtualHardDisk } # spec name: image
    caching::Any # spec type: Union{ Nothing, Caching } # spec name: caching
    writeAcceleratorEnabled::Any # spec type: Union{ Nothing, Bool } # spec name: writeAcceleratorEnabled
    createOption::Any # spec type: Union{ Nothing, CreateOption } # spec name: createOption
    diskSizeGB::Any # spec type: Union{ Nothing, Int32 } # spec name: diskSizeGB
    managedDisk::Any # spec type: Union{ Nothing, ManagedDiskParameters } # spec name: managedDisk
    toBeDetached::Any # spec type: Union{ Nothing, Bool } # spec name: toBeDetached
    diskIOPSReadWrite::Any # spec type: Union{ Nothing, Int64 } # spec name: diskIOPSReadWrite
    diskMBpsReadWrite::Any # spec type: Union{ Nothing, Int64 } # spec name: diskMBpsReadWrite

    function DataDisk(;lun=nothing, name=nothing, vhd=nothing, image=nothing, caching=nothing, writeAcceleratorEnabled=nothing, createOption=nothing, diskSizeGB=nothing, managedDisk=nothing, toBeDetached=nothing, diskIOPSReadWrite=nothing, diskMBpsReadWrite=nothing)
        o = new()
        validate_property(DataDisk, Symbol("lun"), lun)
        setfield!(o, Symbol("lun"), lun)
        validate_property(DataDisk, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(DataDisk, Symbol("vhd"), vhd)
        setfield!(o, Symbol("vhd"), vhd)
        validate_property(DataDisk, Symbol("image"), image)
        setfield!(o, Symbol("image"), image)
        validate_property(DataDisk, Symbol("caching"), caching)
        setfield!(o, Symbol("caching"), caching)
        validate_property(DataDisk, Symbol("writeAcceleratorEnabled"), writeAcceleratorEnabled)
        setfield!(o, Symbol("writeAcceleratorEnabled"), writeAcceleratorEnabled)
        validate_property(DataDisk, Symbol("createOption"), createOption)
        setfield!(o, Symbol("createOption"), createOption)
        validate_property(DataDisk, Symbol("diskSizeGB"), diskSizeGB)
        setfield!(o, Symbol("diskSizeGB"), diskSizeGB)
        validate_property(DataDisk, Symbol("managedDisk"), managedDisk)
        setfield!(o, Symbol("managedDisk"), managedDisk)
        validate_property(DataDisk, Symbol("toBeDetached"), toBeDetached)
        setfield!(o, Symbol("toBeDetached"), toBeDetached)
        validate_property(DataDisk, Symbol("diskIOPSReadWrite"), diskIOPSReadWrite)
        setfield!(o, Symbol("diskIOPSReadWrite"), diskIOPSReadWrite)
        validate_property(DataDisk, Symbol("diskMBpsReadWrite"), diskMBpsReadWrite)
        setfield!(o, Symbol("diskMBpsReadWrite"), diskMBpsReadWrite)
        o
    end
end # type DataDisk

const _property_map_DataDisk = Dict{Symbol,Symbol}(Symbol("lun")=>Symbol("lun"), Symbol("name")=>Symbol("name"), Symbol("vhd")=>Symbol("vhd"), Symbol("image")=>Symbol("image"), Symbol("caching")=>Symbol("caching"), Symbol("writeAcceleratorEnabled")=>Symbol("writeAcceleratorEnabled"), Symbol("createOption")=>Symbol("createOption"), Symbol("diskSizeGB")=>Symbol("diskSizeGB"), Symbol("managedDisk")=>Symbol("managedDisk"), Symbol("toBeDetached")=>Symbol("toBeDetached"), Symbol("diskIOPSReadWrite")=>Symbol("diskIOPSReadWrite"), Symbol("diskMBpsReadWrite")=>Symbol("diskMBpsReadWrite"))
const _property_types_DataDisk = Dict{Symbol,String}(Symbol("lun")=>"Int32", Symbol("name")=>"String", Symbol("vhd")=>"VirtualHardDisk", Symbol("image")=>"VirtualHardDisk", Symbol("caching")=>"Caching", Symbol("writeAcceleratorEnabled")=>"Bool", Symbol("createOption")=>"CreateOption", Symbol("diskSizeGB")=>"Int32", Symbol("managedDisk")=>"ManagedDiskParameters", Symbol("toBeDetached")=>"Bool", Symbol("diskIOPSReadWrite")=>"Int64", Symbol("diskMBpsReadWrite")=>"Int64")
Base.propertynames(::Type{ DataDisk }) = collect(keys(_property_map_DataDisk))
Swagger.property_type(::Type{ DataDisk }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_DataDisk[name]))}
Swagger.field_name(::Type{ DataDisk }, property_name::Symbol) =  _property_map_DataDisk[property_name]

function check_required(o::DataDisk)
    (getproperty(o, Symbol("lun")) === nothing) && (return false)
    (getproperty(o, Symbol("createOption")) === nothing) && (return false)
    true
end

function validate_property(::Type{ DataDisk }, name::Symbol, val)
end


text:
using DataStructures

input = joinpath(@__DIR__, "input")
lines = readlines(input)

m = hcat(map(line -> map(n -> parse(Int, n), collect(line)), lines)...)

function expand(m)
    acc = [m]
    for i in 1:4
        push!(acc, map(x -> mod1(x, 9), m.+ i))
    end
    tmp = hcat(acc...)
    acc = [tmp]
    for i in 1:4
        n = tmp.+ 1
        push!(acc, map(x -> mod1(x, 9), tmp.+ i))
    end
    return vcat(acc...)
end

function neighbors(m, idx)
    x, y = Tuple(idx)
    return filter(c -> checkbounds(Bool, m, c), [
        CartesianIndex(x + 1, y),
        CartesianIndex(x - 1, y),
        CartesianIndex(x, y + 1),
        CartesianIndex(x, y - 1)
    ])
end

function path(m, origin, dest)
    unvisited = Set{}(collect(CartesianIndices(m)))
    candidates = PriorityQueue{CartesianIndex, Int}()
    node = origin
    # initial node is free
    cost = 0

    while dest in unvisited
        delete!(unvisited, node)
        for n in neighbors(m, node)
            if n in unvisited
                haskey(candidates, n) ?
                    candidates[n] = minimum([candidates[n], cost + m[n]]) :
                    candidates[n] = cost + m[n]
            end
        end
        node, cost = dequeue_pair!(candidates)
        node == dest && delete!(unvisited, node)
    end
    return cost
end

p1 = path(m, CartesianIndex(1, 1), CartesianIndex(size(m)))
m = expand(m)
p2 = path(m, CartesianIndex(1, 1), CartesianIndex(size(m)))

println("-----------------------------------------------------------------------")
println("chiton -- part one :: $p1")
println("chiton -- part two :: $p2")
println("-----------------------------------------------------------------------")

@assert(p1 == 790)
@assert(p2 == 2998)


text:
using Turing
using LinearAlgebra

using BenchmarkHelper

include("lr_helper.jl")

if !haskey(BenchmarkSuite, "nuts")
    BenchmarkSuite["nuts"] = BenchmarkGroup(["nuts"])
end

X, Y = readlrdata()

@model function lr_nuts(x, y, σ)

    N,D = size(x)

    α ~ Normal(0, σ)
    β ~ MvNormal(zeros(D), σ^2 * I)

    for n = 1:N
        y[n] ~ BinomialLogit(1, dot(x[n,:], β) + α)
    end
end

# Sampling parameter settings
n_samples = 1_000
n_adapts = 1_000

# Sampling
BenchmarkSuite["nuts"]["lr"] = @benchmarkable sample(lr_nuts(X, Y, 100),
                                                     NUTS(0.65), n_samples)


text:
module SimLynx

include("stub.jl")

export greet

end


text:
#####
##### Repository of Available Games
#####

const AVAILABLE_GAMES = ["tictactoe", "connect-four", "tictac-go" #=, "mancala" =#]

for game in AVAILABLE_GAMES
  file = "../games/$game/main.jl"
  @eval include($file)
end

# Convert a string identifier in upper snake case
function snake(str)
  ret = Char[]
  new_word = true
  for c in str
    if c ∈ ['-', '_']
      new_word = true
    else
      push!(ret, new_word ? uppercase(c) : c)
      new_word = false
    end
  end
  return String(ret)
end

const GAME_MODULE = Dict([g => eval(Symbol(snake(g))) for g in AVAILABLE_GAMES])


text:
int_rules_5_3_3 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 5.3.3*(d+e*x)^m*(a+b*arctan(c*x^n))^p =#
    @apply_utils Antiderivative((~(a') + ~(b') * atan(~(c') * ~x)) ^ ~(p') / (~d + ~(e') * ~x), ~x) => (-((~a + ~b * atan(~c * ~x)) ^ ~p) * log(2 / (1 + (~e * ~x) / ~d))) / ~e + ((~b * ~c * ~p) / ~e) * Antiderivative(((~a + ~b * atan(~c * ~x)) ^ (~p - 1) * log(2 / (1 + (~e * ~x) / ~d))) / (1 + (~c) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (IGtQ(~p, 0) && EqQ((~c) ^ 2 * (~d) ^ 2 + (~e) ^ 2, 0))
    @apply_utils Antiderivative((~(a') + ~(b') * acot(~(c') * ~x)) ^ ~(p') / (~d + ~(e') * ~x), ~x) => (-((~a + ~b * acot(~c * ~x)) ^ ~p) * log(2 / (1 + (~e * ~x) / ~d))) / ~e - ((~b * ~c * ~p) / ~e) * Antiderivative(((~a + ~b * acot(~c * ~x)) ^ (~p - 1) * log(2 / (1 + (~e * ~x) / ~d))) / (1 + (~c) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (IGtQ(~p, 0) && EqQ((~c) ^ 2 * (~d) ^ 2 + (~e) ^ 2, 0))
    @apply_utils Antiderivative((~(a') + ~(b') * atan(~(c') * ~x)) / (~d + ~(e') * ~x), ~x) => ((-((~a + ~b * atan(~c * ~x))) * log(2 / (1 - im * ~c * ~x))) / ~e + ((~b * ~c) / ~e) * Antiderivative(log(2 / (1 - im * ~c * ~x)) / (1 + (~c) ^ 2 * (~x) ^ 2), ~x) + ((~a + ~b * atan(~c * ~x)) * log((2 * ~c * (~d + ~e * ~x)) / ((~c * ~d + im * ~e) * (1 - im * ~c * ~x)))) / ~e) - ((~b * ~c) / ~e) * Antiderivative(log((2 * ~c * (~d + ~e * ~x)) / ((~c * ~d + im * ~e) * (1 - im * ~c * ~x))) / (1 + (~c) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && NeQ((~c) ^ 2 * (~d) ^ 2 + (~e) ^ 2, 0)
    @apply_utils Antiderivative((~(a') + ~(b') * acot(~(c') * ~x)) / (~d + ~(e') * ~x), ~x) => ((-((~a + ~b * acot(~c * ~x))) * log(2 / (1 - im * ~c * ~x))) / ~e - ((~b * ~c) / ~e) * Antiderivative(log(2 / (1 - im * ~c * ~x)) / (1 + (~c) ^ 2 * (~x) ^ 2), ~x)) + ((~a + ~b * acot(~c * ~x)) * log((2 * ~c * (~d + ~e * ~x)) / ((~c * ~d + im * ~e) * (1 - im * ~c * ~x)))) / ~e + ((~b * ~c) / ~e) * Antiderivative(log((2 * ~c * (~d + ~e * ~x)) / ((~c * ~d + im * ~e) * (1 - im * ~c * ~x))) / (1 + (~c) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && NeQ((~c) ^ 2 * (~d) ^ 2 + (~e) ^ 2, 0)
    @apply_utils Antiderivative((~(a') + ~(b') * atan(~(c') * ~x)) ^ 2 / (~d + ~(e') * ~x), ~x) => (((((-((~a + ~b * atan(~c * ~x)) ^ 2) * log(2 / (1 - im * ~c * ~x))) / ~e + (im * ~b * (~a + ~b * atan(~c * ~x)) * PolyLog(2, 1 - 2 / (1 - im * ~c * ~x))) / ~e) - ((~b) ^ 2 * PolyLog(3, 1 - 2 / (1 - im * ~c * ~x))) / (2 * ~e)) + ((~a + ~b * atan(~c * ~x)) ^ 2 * log((2 * ~c * (~d + ~e * ~x)) / ((~c * ~d + im * ~e) * (1 - im * ~c * ~x)))) / ~e) - (im * ~b * (~a + ~b * atan(~c * ~x)) * PolyLog(2, 1 - (2 * ~c * (~d + ~e * ~x)) / ((~c * ~d + im * ~e) * (1 - im * ~c * ~x)))) / ~e) + ((~b) ^ 2 * PolyLog(3, 1 - (2 * ~c * (~d + ~e * ~x)) / ((~c * ~d + im * ~e) * (1 - im * ~c * ~x)))) / (2 * ~e) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && NeQ((~c) ^ 2 * (~d) ^ 2 + (~e) ^ 2, 0)
    @apply_utils Antiderivative((~(a') + ~(b') * acot(~(c') * ~x)) ^ 2 / (~d + ~(e') * ~x), ~x) => (((-((~a + ~b * acot(~c * ~x)) ^ 2) * log(2 / (1 - im * ~c * ~x))) / ~e - (im * ~b * (~a + ~b * acot(~c * ~x)) * PolyLog(2, 1 - 2 / (1 - im * ~c * ~x))) / ~e) - ((~b) ^ 2 * PolyLog(3, 1 - 2 / (1 - im * ~c * ~x))) / (2 * ~e)) + ((~a + ~b * acot(~c * ~x)) ^ 2 * log((2 * ~c * (~d + ~e * ~x)) / ((~c * ~d + im * ~e) * (1 - im * ~c * ~x)))) / ~e + (im * ~b * (~a + ~b * acot(~c * ~x)) * PolyLog(2, 1 - (2 * ~c * (~d + ~e * ~x)) / ((~c * ~d + im * ~e) * (1 - im * ~c * ~x)))) / ~e + ((~b) ^ 2 * PolyLog(3, 1 - (2 * ~c * (~d + ~e * ~x)) / ((~c * ~d + im * ~e) * (1 - im * ~c * ~x)))) / (2 * ~e) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && NeQ((~c) ^ 2 * (~d) ^ 2 + (~e) ^ 2, 0)
    @apply_utils Antiderivative((~(a') + ~(b') * atan(~(c') * ~x)) ^ 3 / (~d + ~(e') * ~x), ~x) => ((((((-((~a + ~b * atan(~c * ~x)) ^ 3) * log(2 / (1 - im * ~c * ~x))) / ~e + (3 * im * ~b * (~a + ~b * atan(~c * ~x)) ^ 2 * PolyLog(2, 1 - 2 / (1 - im * ~c * ~x))) / (2 * ~e)) - (3 * (~b) ^ 2 * (~a + ~b * atan(~c * ~x)) * PolyLog(3, 1 - 2 / (1 - im * ~c * ~x))) / (2 * ~e)) - (3 * im * (~b) ^ 3 * PolyLog(4, 1 - 2 / (1 - im * ~c * ~x))) / (4 * ~e)) + ((~a + ~b * atan(~c * ~x)) ^ 3 * log((2 * ~c * (~d + ~e * ~x)) / ((~c * ~d + im * ~e) * (1 - im * ~c * ~x)))) / ~e) - (3 * im * ~b * (~a + ~b * atan(~c * ~x)) ^ 2 * PolyLog(2, 1 - (2 * ~c * (~d + ~e * ~x)) / ((~c * ~d + im * ~e) * (1 - im * ~c * ~x)))) / (2 * ~e)) + (3 * (~b) ^ 2 * (~a + ~b * atan(~c * ~x)) * PolyLog(3, 1 - (2 * ~c * (~d + ~e * ~x)) / ((~c * ~d + im * ~e) * (1 - im * ~c * ~x)))) / (2 * ~e) + (3 * im * (~b) ^ 3 * PolyLog(4, 1 - (2 * ~c * (~d + ~e * ~x)) / ((~c * ~d + im * ~e) * (1 - im * ~c * ~x)))) / (4 * ~e) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && NeQ((~c) ^ 2 * (~d) ^ 2 + (~e) ^ 2, 0)
    @apply_utils Antiderivative((~(a') + ~(b') * acot(~(c') * ~x)) ^ 3 / (~d + ~(e') * ~x), ~x) => ((((-((~a + ~b * acot(~c * ~x)) ^ 3) * log(2 / (1 - im * ~c * ~x))) / ~e - (3 * im * ~b * (~a + ~b * acot(~c * ~x)) ^ 2 * PolyLog(2, 1 - 2 / (1 - im * ~c * ~x))) / (2 * ~e)) - (3 * (~b) ^ 2 * (~a + ~b * acot(~c * ~x)) * PolyLog(3, 1 - 2 / (1 - im * ~c * ~x))) / (2 * ~e)) + (3 * im * (~b) ^ 3 * PolyLog(4, 1 - 2 / (1 - im * ~c * ~x))) / (4 * ~e) + ((~a + ~b * acot(~c * ~x)) ^ 3 * log((2 * ~c * (~d + ~e * ~x)) / ((~c * ~d + im * ~e) * (1 - im * ~c * ~x)))) / ~e + (3 * im * ~b * (~a + ~b * acot(~c * ~x)) ^ 2 * PolyLog(2, 1 - (2 * ~c * (~d + ~e * ~x)) / ((~c * ~d + im * ~e) * (1 - im * ~c * ~x)))) / (2 * ~e) + (3 * (~b) ^ 2 * (~a + ~b * acot(~c * ~x)) * PolyLog(3, 1 - (2 * ~c * (~d + ~e * ~x)) / ((~c * ~d + im * ~e) * (1 - im * ~c * ~x)))) / (2 * ~e)) - (3 * im * (~b) ^ 3 * PolyLog(4, 1 - (2 * ~c * (~d + ~e * ~x)) / ((~c * ~d + im * ~e) * (1 - im * ~c * ~x)))) / (4 * ~e) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && NeQ((~c) ^ 2 * (~d) ^ 2 + (~e) ^ 2, 0)
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~(q') * (~(a') + ~(b') * atan(~(c') * ~x)), ~x) => ((~d + ~e * ~x) ^ (~q + 1) * (~a + ~b * atan(~c * ~x))) / (~e * (~q + 1)) - ((~b * ~c) / (~e * (~q + 1))) * Antiderivative((~d + ~e * ~x) ^ (~q + 1) / (1 + (~c) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~q], ~x) && NeQ(~q, -1)
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~(q') * (~(a') + ~(b') * acot(~(c') * ~x)), ~x) => ((~d + ~e * ~x) ^ (~q + 1) * (~a + ~b * acot(~c * ~x))) / (~e * (~q + 1)) + ((~b * ~c) / (~e * (~q + 1))) * Antiderivative((~d + ~e * ~x) ^ (~q + 1) / (1 + (~c) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~q], ~x) && NeQ(~q, -1)
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~(q') * (~(a') + ~(b') * atan(~(c') * ~x)) ^ ~p, ~x) => ((~d + ~e * ~x) ^ (~q + 1) * (~a + ~b * atan(~c * ~x)) ^ ~p) / (~e * (~q + 1)) - ((~b * ~c * ~p) / (~e * (~q + 1))) * Antiderivative(ExpandIntegrand((~a + ~b * atan(~c * ~x)) ^ (~p - 1), (~d + ~e * ~x) ^ (~q + 1) / (1 + (~c) ^ 2 * (~x) ^ 2), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (IGtQ(~p, 1) && (IntegerQ(~q) && NeQ(~q, -1)))
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~(q') * (~(a') + ~(b') * acot(~(c') * ~x)) ^ ~p, ~x) => ((~d + ~e * ~x) ^ (~q + 1) * (~a + ~b * acot(~c * ~x)) ^ ~p) / (~e * (~q + 1)) + ((~b * ~c * ~p) / (~e * (~q + 1))) * Antiderivative(ExpandIntegrand((~a + ~b * acot(~c * ~x)) ^ (~p - 1), (~d + ~e * ~x) ^ (~q + 1) / (1 + (~c) ^ 2 * (~x) ^ 2), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (IGtQ(~p, 1) && (IntegerQ(~q) && NeQ(~q, -1)))
    @apply_utils Antiderivative((~(a') + ~(b') * atan(~(c') * (~x) ^ ~n)) / (~d + ~(e') * ~x), ~x) => (log(~d + ~e * ~x) * (~a + ~b * atan(~c * (~x) ^ ~n))) / ~e - ((~b * ~c * ~n) / ~e) * Antiderivative(((~x) ^ (~n - 1) * log(~d + ~e * ~x)) / (1 + (~c) ^ 2 * (~x) ^ (2 * ~n)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && IntegerQ(~n)
    @apply_utils Antiderivative((~(a') + ~(b') * acot(~(c') * (~x) ^ ~n)) / (~d + ~(e') * ~x), ~x) => (log(~d + ~e * ~x) * (~a + ~b * acot(~c * (~x) ^ ~n))) / ~e + ((~b * ~c * ~n) / ~e) * Antiderivative(((~x) ^ (~n - 1) * log(~d + ~e * ~x)) / (1 + (~c) ^ 2 * (~x) ^ (2 * ~n)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && IntegerQ(~n)
    @apply_utils Antiderivative((~(a') + ~(b') * atan(~(c') * (~x) ^ ~n)) / (~d + ~(e') * ~x), ~x) => With([k = Denominator(~n)], k * Subst(Antiderivative(((~x) ^ (k - 1) * (~a + ~b * atan(~c * (~x) ^ (k * ~n)))) / (~d + ~e * (~x) ^ k), ~x), ~x, (~x) ^ (1 / k))) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && FractionQ(~n)
    @apply_utils Antiderivative((~(a') + ~(b') * acot(~(c') * (~x) ^ ~n)) / (~d + ~(e') * ~x), ~x) => With([k = Denominator(~n)], k * Subst(Antiderivative(((~x) ^ (k - 1) * (~a + ~b * acot(~c * (~x) ^ (k * ~n)))) / (~d + ~e * (~x) ^ k), ~x), ~x, (~x) ^ (1 / k))) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && FractionQ(~n)
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * atan(~(c') * (~x) ^ ~n)), ~x) => ((~d + ~e * ~x) ^ (~m + 1) * (~a + ~b * atan(~c * (~x) ^ ~n))) / (~e * (~m + 1)) - ((~b * ~c * ~n) / (~e * (~m + 1))) * Antiderivative(((~x) ^ (~n - 1) * (~d + ~e * ~x) ^ (~m + 1)) / (1 + (~c) ^ 2 * (~x) ^ (2 * ~n)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x) && NeQ(~m, -1)
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * acot(~(c') * (~x) ^ ~n)), ~x) => ((~d + ~e * ~x) ^ (~m + 1) * (~a + ~b * acot(~c * (~x) ^ ~n))) / (~e * (~m + 1)) + ((~b * ~c * ~n) / (~e * (~m + 1))) * Antiderivative(((~x) ^ (~n - 1) * (~d + ~e * ~x) ^ (~m + 1)) / (1 + (~c) ^ 2 * (~x) ^ (2 * ~n)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x) && NeQ(~m, -1)
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * atan(~(c') * (~x) ^ ~n)) ^ ~p, ~x) => Antiderivative(ExpandIntegrand((~a + ~b * atan(~c * (~x) ^ ~n)) ^ ~p, (~d + ~e * ~x) ^ ~m, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (IGtQ(~p, 1) && IGtQ(~m, 0))
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * acot(~(c') * (~x) ^ ~n)) ^ ~p, ~x) => Antiderivative(ExpandIntegrand((~a + ~b * acot(~c * (~x) ^ ~n)) ^ ~p, (~d + ~e * ~x) ^ ~m, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (IGtQ(~p, 1) && IGtQ(~m, 0))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * atan(~(c') * (~x) ^ ~n)) ^ ~(p'), ~x) => Unintegrable((~d + ~e * ~x) ^ ~m * (~a + ~b * atan(~c * (~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n, ~p], ~x)
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * acot(~(c') * (~x) ^ ~n)) ^ ~(p'), ~x) => Unintegrable((~d + ~e * ~x) ^ ~m * (~a + ~b * acot(~c * (~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n, ~p], ~x)
end

max_stars_repo_path: src/rules/5 Inverse trig functions/5.3 Inverse tangent/5.3.3 (d+e x)^m (a+b arctan(c x^n))^p.jl


text:
using Distributed
using Random

workers = 16

if nprocs() <= workers
    addprocs(workers + 1 - nprocs())
end

@everywhere include("models/Dicke.jl")
@everywhere include("modules/ClassicalDynamics.jl")

""" Calculates Poincaré sections with Lyapunov exponents for various energies
    Parallel calculation, takes usually days to finish
"""

@everywhere function SolveItem(energy, parameters, dimension; file="Dicke.txt", path="", alreadySolved=[], sectionPlane=2)
    for (E, λ) in alreadySolved
        if isapprox(E, energy) && isapprox(parameters[1], λ)
            println("Skipped $parameters, E=$energy, dim=$dimension (already calculated)")
            return
        end
    end

    time = @elapsed averageLyapunov, maximumLyapunov, freg, trajectories = SolveEnergy(energy, parameters, dimension, savePath=path, sectionPlane=sectionPlane)

    chaos = 0
    total = 0
    error = 0
    meanLyapunov = 0
    meanLyapunovChaos = 0

    for x in averageLyapunov
        if x > 0.0
            total += 1
            meanLyapunov += x
        end
        if x == 0.0
            error += 1
        end
        if x > 0.01
            chaos += 1
            meanLyapunovChaos += x
        end
    end

    result = [energy, parameters[1], total, chaos, error, freg, total > 0 ? meanLyapunov / total : 0, chaos > 0 ? meanLyapunovChaos / chaos : 0, maximumLyapunov, myid(), time, trajectories]

    open(path * file, "a") do io
        println(io, result)
    end

    return
end

function RunMap(; δ=1.0, ω=1.0, ω₀=1.0, path="", dimension=101, step=0.1, sectionPlane=2)
    path *= "Dicke_"

    file = "Map_dim=$(dimension)_$([δ, ω, ω₀]).txt"
    alreadySolved = ReadMap(path * file)

    println("Already calculated $(length(alreadySolved)) points.")

    λᵪ = sqrt(ω * ω₀) / (1.0 + δ)
    println("Critical value λc=$λᵪ")

    input = shuffle([(energy, [λ * λᵪ, δ, ω, ω₀], dimension) for energy in 6:-step:-4, λ in 4:-step:0])
    pmap((args)->SolveItem(args...; file=file, path=path, alreadySolved=alreadySolved, sectionPlane=sectionPlane), input)

    return
end

""" Calculates freg for one given lambda """
function RunLambda(; λ=2.0, δ=1.0, ω=1.0, ω₀=1.0, path="", dimension=101, step=0.1, sectionPlane=2)
    path *= "Dicke_"

    file = "Energy_dim=$(dimension)_$([λ, δ, ω, ω₀]).txt"
    alreadySolved = ReadMap(path * file)

    println("Already calculated $(length(alreadySolved)) points.")

    λᵪ = sqrt(ω * ω₀) / (1.0 + δ)
    println("Critical value λc=$λᵪ")

    input = [(energy, [λ * λᵪ, δ, ω, ω₀], dimension) for energy in -4:step:30]

    pmap((args)->SolveItem(args...; file=file, path=path, alreadySolved=alreadySolved, sectionPlane=sectionPlane), input)

    return
end

function ReadMap(file="")
    Eλ = []

    try
        for line in eachline(file)
            line = replace(line, "[" => "")
            line = replace(line, "]" => "")
            line = replace(line, "," => "")
            elements = split(line)

            append!(Eλ, [(parse(Float64, elements[1]), parse(Float64, elements[2]))])
        end
    catch x
    end

    return Eλ
end

max_stars_repo_name: PavelStransky/ClassicalDynamics

max_issues_repo_name: PavelStransky/ClassicalDynamics

max_forks_repo_name: PavelStransky/ClassicalDynamics


text:
using HoJBot:
    DiscourseData, DiscoursePost,
    discourse_run_query, discourse_load, discourse_save,
    current, message, next, previous

import JSON3

@testset "discourse" begin

    json = discourse_run_query("search.json", "pluto")

    let
        @test json isa JSON3.Object
        @test hasproperty(json, :topics)
        @test length(json.topics) > 0
        @test hasproperty(json.topics[1], :id)
        @test hasproperty(json.topics[1], :slug)
    end

    let data = DiscourseData(json)
        @test length(data) == length(json.topics)
        @test 1 <= data.index <= length(data)
        @test current(data) isa DiscoursePost

        let index = data.index
            # cycle once using `next`
            for i in 1:length(data)
                data = next(data)
            end
            @test data.index === index

            # cycle once using `pevious`
            for i in 1:length(data)
                data = previous(data)
            end
            @test data.index === index
        end

        let msg = message(data)
            @test length(msg) > 0
            @test occursin("http", msg)
        end

        # load/save
        let
            id = UInt64(0)
            @test_nowarn discourse_save(id, data)
            loaded = @test_nowarn discourse_load(id)
            @test length(data) == length(loaded)
            @test data.index == loaded.index
            @test data.posts == loaded.posts
        end
    end
end


text:
to_scale(x::Function) = x
to_scale(x::AbstractArray) = x
to_scale(x::Any) = nothing
to_scale(x::Nothing) = nothing


text:
# UnsafeStridedView
struct UnsafeStridedView{T,N,PT,F<:Union{FN,FC,FA,FT}} <: AbstractStridedView{T,N,F}
    ptr::Ptr{PT}
    size::NTuple{N,Int}
    strides::NTuple{N,Int}
    offset::Int
    op::F
end
function UnsafeStridedView(a::Ptr{PT}, size::NTuple{N,Int}, strides::NTuple{N,Int},
    offset::Int, op::F) where {PT,N,F}

    @assert isbitstype(PT)
    T = Base.promote_op(op, PT)
    UnsafeStridedView{T,N,PT,F}(a, size, _normalizestrides(size, strides), offset, op)
end

UnsafeStridedView(a::Ptr{T}, size::NTuple{N,Int}, strides::NTuple{N,Int},
    offset::Int = 0) where {T,N} = UnsafeStridedView(a, size, strides, offset, identity)
UnsafeStridedView(a::DenseArray) = UnsafeStridedView(pointer(a), size(a), strides(a))
UnsafeStridedView(a::StridedView) = UnsafeStridedView(pointer(a), size(a), strides(a))

UnsafeStridedView(a::Adjoint{<:Any, <:Any}) = UnsafeStridedView(a')'
UnsafeStridedView(a::Transpose{<:Any, <:Any}) =
    transpose(UnsafeStridedView(transpose(a)))
UnsafeStridedView(a::Base.SubArray) = sview(UnsafeStridedView(a.parent), a.indices...)
UnsafeStridedView(a::Base.ReshapedArray) = sreshape(UnsafeStridedView(a.parent), a.dims)

# Methods for UnsafeStridedView
Base.size(a::UnsafeStridedView) = a.size
Base.strides(a::UnsafeStridedView) = a.strides
Base.stride(a::UnsafeStridedView{<:Any, 0}, n::Int) = 1
Base.stride(a::UnsafeStridedView{<:Any, N}, n::Int) where {N} =
    (n <= N) ? a.strides[n] : a.strides[N]*a.size[N]
offset(a::UnsafeStridedView) = a.offset

function blasstrides(a::UnsafeStridedView{<:Any,2})
    # canonialize strides to make compatible with gemm
    if size(a, 2) == 1 && stride(a, 1) == 1
        return UnsafeStridedView(a.ptr, a.size, (1, size(a,1)), a.offset, a.op)
    else
        return a
    end
end

Base.similar(a::UnsafeStridedView, ::Type{T}, dims::NTuple{N,Int}) where {N,T} =
    StridedView(Array{T}(undef, dims))
Base.copy(a::UnsafeStridedView) = copyto!(similar(a), a)

Base.unsafe_convert(::Type{Ptr{T}}, a::UnsafeStridedView{T}) where {T} = a.ptr + Base.elsize(a)*a.offset

Base.dataids(a::UnsafeStridedView) = (UInt(a.ptr),)

Base.IndexStyle(::Type{<:UnsafeStridedView}) = Base.IndexCartesian()

# Indexing with N integer arguments
@inline function Base.getindex(a::UnsafeStridedView{<:Any,N}, I::Vararg{Int,N}) where N
    @boundscheck checkbounds(a, I...)
    r = a.op(unsafe_load(a.ptr, a.offset+_computeind(I, a.strides)))
    return r
end
@inline function Base.setindex!(a::UnsafeStridedView{<:Any,N}, v, I::Vararg{Int,N}) where N
    @boundscheck checkbounds(a, I...)
    unsafe_store!(a.ptr, a.op(v), a.offset+_computeind(I, a.strides))
    return a
end

function Base.getindex(a::UnsafeStridedView{<:Any,N}, I::Vararg{SliceIndex,N}) where {N}
    UnsafeStridedView(a.ptr, _computeviewsize(a.size, I), _computeviewstrides(a.strides, I), a.offset + _computeviewoffset(a.strides, I), a.op)
end

@propagate_inbounds Base.getindex(a::UnsafeStridedView, I::ParentIndex) =
    a.op(unsafe_load(a.ptr, I.i))
@propagate_inbounds Base.setindex!(a::UnsafeStridedView, v, I::ParentIndex) =
    (unsafe_store!(a.ptr, a.op(v), I.i); return a)

# Specialized methods for `UnsafeStridedView` which produce views/share data
Base.conj(a::UnsafeStridedView{<:Real}) = a
Base.conj(a::UnsafeStridedView) =
    UnsafeStridedView(a.ptr, a.size, a.strides, a.offset, _conj(a.op))

function Base.permutedims(a::UnsafeStridedView{<:Any,N}, p) where {N}
    _isperm(N, p) || throw(ArgumentError("Invalid permutation of length $N: $p"))
    newsize = TupleTools._permute(a.size, p)
    newstrides = TupleTools._permute(a.strides, p)
    return UnsafeStridedView(a.ptr, newsize, newstrides, a.offset, a.op)
end

LinearAlgebra.transpose(a::UnsafeStridedView{<:Number,2}) = permutedims(a, (2,1))
LinearAlgebra.adjoint(a::UnsafeStridedView{<:Number,2}) = permutedims(conj(a), (2,1))
LinearAlgebra.adjoint(a::UnsafeStridedView{<:Any,2}) = # act recursively, like Base
    permutedims(UnsafeStridedView(a.ptr,a.size,a.strides,a.offset, _adjoint(a.op)), (2,1))
LinearAlgebra.transpose(a::UnsafeStridedView{<:Any,2}) = # act recursively, like Base
    permutedims(UnsafeStridedView(a.ptr,a.size,a.strides,a.offset, _transpose(a.op)), (2,1))

Base.map(::FC, a::UnsafeStridedView{<:Real}) = a
Base.map(::FT, a::UnsafeStridedView{<:Number}) = a
Base.map(::FA, a::UnsafeStridedView{<:Number}) = conj(a)
Base.map(::FC, a::UnsafeStridedView) =
    UnsafeStridedView(a.ptr, a.size, a.strides, a.offset, _conj(a.op))
Base.map(::FT, a::UnsafeStridedView) =
    UnsafeStridedView(a.ptr, a.size, a.strides, a.offset, _transpose(a.op))
Base.map(::FA, a::UnsafeStridedView) =
    UnsafeStridedView(a.ptr, a.size, a.strides, a.offset, _adjoint(a.op))

@inline function sreshape(a::UnsafeStridedView, newsize::Dims)
    if any(isequal(0), newsize)
        any(isequal(0), size(a)) || throw(DimensionMismatch())
        newstrides = one.(newsize)
    else
        newstrides = _computereshapestrides(newsize, _simplify(size(a), strides(a))...)
    end
    UnsafeStridedView(a.ptr, newsize, newstrides, a.offset, a.op)
end


text:
using Test
using MeshArrays

@testset "MeshArrays tests:" begin
for nTopo=1:3
    if nTopo==1; grTopo="cs"; nFaces=6; N=200;
    elseif nTopo==2; grTopo="llc"; nFaces=5; N=200;
    elseif nTopo==3; grTopo="ll"; nFaces=1; N=500;
    end;
    Npt=nFaces*N*N
    GridVariables=GridOfOnes(grTopo,nFaces,N)
    mygrid=GridVariables["XC"].grid
    @test mygrid.class == grTopo
    Rini= 0.; Rend= 0.;
    (Rini,Rend,DXCsm,DYCsm)=MeshArrays.demo2(GridVariables);
    @test isa(Rend,MeshArray)
    @test sum(isfinite.(Rend)) == Npt
    Sini=sqrt(sum(Rini*Rini)/(Npt-1.0))
    Send=sqrt(sum(Rend*Rend)/(Npt-1.0))
    #println([Sini Send])
    @test isapprox(Sini,1.000; atol=1e-2)
    @test isapprox(Send,0.093; atol=1e-2)
end;
end;


text:
using StaticArrays
using BenchmarkTools

combination(x) = combination_(x...)
combination_(x) = x
combination_(x,y) = (x,y,(x,y))
combination_(x,y...) = (x,combination(y)...,map(z->(x,z...),combination(y))...)

combination((1))
combination((2,3))
combination((2,3,4))
combination((1,2,3,4))

for i in 1:15
    print(count_ones(i)," ")
end

function one_positions_array(x::Int64, bits=64)
    out = Vector{Int}([])
    sizehint!(out,count_ones(x))
    for position in 1:bits-leading_zeros(x)
        if x&1 == 1
            push!(out,position)
        end
        x>>>=1
    end
    out
end

function one_positions_mvector(x::Int64, bits=64)
    out = zeros(MVector{count_ones(x),Int64})
    i = 1
    for position in 1:bits-leading_zeros(x)
        if x&1 == 1
            out[i]=position
            i+=1
        end
        x>>>=1
    end
    out
end

function one_positions_array2(x::Int64, bits=64)
    out = zeros(Int64,count_ones(x))
    i = 1
    for position in 1:bits-leading_zeros(x)
        if x&1 == 1
            out[i]=position
            i+=1
        end
        x>>>=1
    end
    out
end

@benchmark one_positions_array(x) setup=(x=rand(Int))
@benchmark one_positions_mvector(x) setup=(x=rand(Int))
@benchmark one_positions_array2(x) setup=(x=rand(Int))

for i in 1:20
    println(i," => ",one_positions(i))
end
for i in 1:20
    print(i," => ")
    for p in OnePositions(i)
        print(p," ")
    end
    println()
end

@benchmark sum(one_positions_array2(x)) setup=(x=rand(Int))
@benchmark sum(OnePositions(x)) setup=(x=rand(Int))

struct Strategy{T}
    groups :: Vector{T}
end

mutable struct Solution{T<:Integer}
    s :: T
    n :: Int
    a :: Array{T,1}
end
solution(s,n) = Solution{Int}(s,n,[])
function canonical!(solution::Solution)
    sort!(solution.a)
end
function isvalid(solution::Solution)
    s = solution.s
    n = solution.n
    a = solution.a
    s<3 && return false
    n<16 && return false
    n>40 && return false
    length(a)<2 && return false
    length(a)>(s-1) && return false
    for x in a
        x<1 && return false
        x>(s-1) && return false
    end
    true
end

function Base.string(sol::Solution)
    out = string(sol.s,"^",sol.n,"=>{",sol.a[1])
    for value in sol.a[2:end]
        out = out*string(",",value)
    end
    out = out*"}"
    out
end

function err(sol::Solution, acc=BigInt(0))
    acc = acc + sol.s^sol.n
    for a in sol.a
        acc = acc - a^sol.n
    end
    acc
end

function random_search(s, n, iterations=1e5,
                             best_e=s^(BigInt(n)),
                             powers=powers_tuple_old(s,n))
    s_to_n = s^(BigInt(n))
    best_a = Vector{Int}()
    all_a = Set(1:s-1)
    gotone = false
    for i in 1:iterations
        x = copy(all_a)
        a = Vector{Int}()
        e = s_to_n
        while true
            length(x)==0 && break
            new_a = rand(x)
            pop!(x,new_a)
            new_e = e - powers[s-new_a]
            abs(new_e)>abs(e) && break
            push!(a,new_a)
            e = new_e
        end
        if abs(e)<abs(best_e)
            gotone = true
            best_e = e
            best_a = a
        end
        if best_e == 0
            break # can't get better
        end
    end
    if gotone
        return Solution(s,n,sort(best_a,rev=true)),best_e
    else
        return nothing
    end
end

function improve!(solution::Solution, strategy::Strategy;
                  powers=powers_tuple(solution.s,solution.n),
                  initial_error=err(solution),
                  s_to_n=solution.s^(BigInt(solution.n))
                  )
    current_error = BigInt(initial_error)
    new_error = BigInt(initial_error)
    current_solution = solution
    while true
        for group in strategy.groups
            new_error = best_group!(current_solution, group,
                                    powers=powers,
                                    initial_error=current_error,
                                    s_to_n = s_to_n)
        end
        new_error>=current_error && break
        current_error = new_error
    end
    solution = current_solution
    current_error
end

function best_group!(solution::Solution, group;
                     powers=powers_tuple(solution.s,solution.n),
                     initial_error=err(solution),
                     s_to_n=solution.s^(BigInt(solution.n))
                    )
    lhs = BigInt(s_to_n)
    a_not_in_group = setdiff(solution.a,group)
    for i in a_not_in_group
        lhs -= powers[i]
    end
    best_e = BigInt(initial_error)
    best_c = Array{Int,1}[]
    for c in combinations(group)
        e = BigInt(lhs)
        for i in c
            e=e-powers[i]
            e<0 && abs(e)>abs(best_e) && break
        end
        if abs(e)<=abs(best_e)
            best_e = e
            best_c = c
        end
        best_e==0 && break
    end
    best_a = sort!(union(best_c,a_not_in_group))
    solution.a = best_a
    best_e
end

powers_tuple(s,n) = ntuple(x->x^BigInt(n),s-1)
powers_tuple_old(s,n) = ntuple(x->(s-x)^BigInt(n),s-1)
function best(s,n,powers=powers_tuple_old(s,n))
    s_to_n = s^(BigInt(n))
    best_e = s_to_n
    all = collect(s-1:-1:1)
    e = s_to_n
    for i in all
        e = e - powers[s-i]
    end
    if e>=0
        return Solution(s,n,all),e
    end
    best_c = Array{Int,1}[]
    for c = combinations(s-1:-1:1)
        e = s_to_n
        for i in c
            e = e - powers[s-i]
            e<0 && abs(e)>abs(best_e) && break # teminate if hopeless
        end
        if abs(e)<abs(best_e)
            best_e = e
            best_c = c
        end
        if best_e == 0
            break # can't get better
        end
    end
    return Solution(s,n,best_c),best_e
end

function max_all_a(n)
    s=2
    while true
        powers=powers_tuple_old(s,n)
        all = collect(s-1:-1:1)
        e = s^(BigInt(n))
        for i in all
            e = e - powers[s-i]
        end
        if e<0
            return s-1
        end
        s+=1
    end
end

function print_best(io::IO,s,n,powers=powers_tuple_old(s,n))
    sol,e = best(s,n,powers)
    println(io,string(sol),",e=",string(e))
end

function print_best_to_file(s_start,s_stop,n,file="data/n$(n)best.txt")
    open(file,"a") do io
        for s in s_start:s_stop
            print_best(io,s,n)
            flush(io)
        end
    end
end

function ns()
    n = 2:17
    s = [4,5,7,8,10,11,12,14,15,17,18,20,21,23,24,26]
    n,s
end

end # module


text:
using Test
using StaticArrays
using ClimateMachine.VariableTemplates
using CUDA
using KernelAbstractions
using KernelAbstractions.Extras: @unroll

include("complex_models.jl")

run_gpu = CUDA.has_cuda_gpu()
if run_gpu
    CUDA.allowscalar(false)
else
    CUDA.allowscalar(true)
end

get_device() = run_gpu ? CPU() : CUDADevice()
device_array(a, ::CUDADevice) = CuArray(a)
device_array(a, ::CPU) = Array(a)
device_rand(::CUDADevice, args...) = CUDA.rand(args...)
device_rand(::CPU, args...) = rand(args...)
number_states(m) = varsize(state(m, Int))

@kernel function mem_copy_kernel!(
    m::AbstractModel,
    dst::AbstractArray{FT, N},
    src::AbstractArray{FT, N},
) where {FT, N}
    @uniform begin
        ns = number_states(m)
        vs = state(m, FT)
        local_src = MArray{Tuple{ns}, FT}(undef)
        local_dst = MArray{Tuple{ns}, FT}(undef)
    end
    i = @index(Group, Linear)
    @inbounds begin
        @unroll for s in 1:ns
            local_src[s] = src[s, i]
        end
        mem_copy!(m, Vars{vs}(local_dst), Vars{vs}(local_src))
        @unroll for s in 1:ns
            dst[s, i] = local_dst[s]
        end
    end
end

function mem_copy!(m::ScalarModel, dst::Vars, src::Vars)
    dst.x = src.x
end

function mem_copy!(m::NTupleContainingModel{N}, dst::Vars, src::Vars) where {N}
    dst.vector_model.x = src.vector_model.x
    dst.scalar_model.x = src.scalar_model.x

    up = vuntuple(i -> src.ntuple_model[i].scalar_model.x, N)
    up_v = vuntuple(i -> src.ntuple_model[i].vector_model.x, N)
    up_sv = SVector(up...)

    @unroll_map(N) do i
        dst.ntuple_model[i].scalar_model.x = up[i]    # index into tuple
        dst.ntuple_model[i].scalar_model.x = up_sv[i] # index into SArray
        dst.ntuple_model[i].vector_model.x = up_v[i]
    end
end

@testset "ScalarModel" begin
    FT = Float32
    device = get_device()
    n_elem = 10
    m = ScalarModel()
    ns = number_states(m)
    a_src = Array{FT}(undef, ns, n_elem)
    a_dst = Array{FT}(undef, ns, n_elem)
    d_src = device_array(a_src, device)
    d_dst = device_array(a_dst, device)
    d_src .= device_rand(device, FT, ns, n_elem)
    fill!(d_dst, 0)

    work_groups = (1,)
    ndrange = (n_elem,)
    kernel! = mem_copy_kernel!(device, work_groups)
    event = kernel!(m, d_dst, d_src, ndrange = ndrange)
    wait(device, event)

    a_src = Array(d_src)
    a_dst = Array(d_dst)
    @test a_src == a_dst
end

@testset "NTupleContainingModel" begin
    FT = Float32
    device = get_device()
    n_elem = 10
    Nv = 4
    N = 3
    m = NTupleContainingModel(N, Nv)
    ns = number_states(m)
    a_src = Array{FT}(undef, ns, n_elem)
    a_dst = Array{FT}(undef, ns, n_elem)
    d_src = device_array(a_src, device)
    d_dst = device_array(a_dst, device)
    d_src .= device_rand(device, FT, ns, n_elem)
    fill!(d_dst, 0)

    work_groups = (1,)
    ndrange = (n_elem,)
    kernel! = mem_copy_kernel!(device, work_groups)
    event = kernel!(m, d_dst, d_src, ndrange = ndrange)
    wait(device, event)

    a_src = Array(d_src)
    a_dst = Array(d_dst)
    @test a_src == a_dst
end


text:
module DCDC
    using LinearAlgebra, DataFrames,Optim, ForwardDiff, BenchmarkTools,Distributions,Expectations, QuantEcon, Statistics, GLM
    using Distributions: invsqrt2π, log2π, sqrt2, invsqrt2
    # using ScikitLearn,JLD,PyCall
    # @sk_import kernel_ridge: KernelRidge
    # @sk_import svm: SVR
    # @sk_import neighbors: KNeighborsRegressor
    # using ScikitLearn: fit!, predict
    using Distributed, Suppressor
    using Flux #Get derivatives
    using Flux.Tracker
    using Flux.Tracker: update!

    # export Utility,ApproxFn,UpdateVal,UpdateData, DynamicDecisionProcess, Transition
    export RealVector,ekernel4, ekernel2, Kernel,estimate,
    forecast,bw_constant,compute_w,ApproxFn,State,
    DynamicDecisionProcess,dynamic_decision_process,simulate_ddc,
    Utility, Transition, State, find_optim, UpdateVal!,computeEquilibrium,check_ee,
    convert_data, EulerEquation, simulate_ee, _moment
    # Parameter, State, profit, ProfitFn,

    include("kernel.jl")
    # include("param.jl")
    include("HouseholdProblem.jl")
    include("misc.jl")
end # module]

a = zeros(10)
Threads.@threads for i = 1:10
   a[i] = Threads.threadid()
end


text:
import Base: ==

mutable struct Point{T<:Number}
    x::T
    y::T
end

==(p::Point, q::Point) = p.x == q.x && p.y == q.y

@enum Heading NORTH=1 EAST=2 SOUTH=3 WEST=4

mutable struct Robot
    position::Point
    heading::Heading
end

function Robot(position::Tuple{T, T}, heading::Heading) where T<:Number
    Robot(Point(position...), heading)
end

heading(r::Robot) = r.heading
position(r::Robot) = r.position

function turn_right!(r::Robot)
    r.heading = Heading(mod1(Int(r.heading) + 1, 4))
    r
end

function turn_left!(r::Robot)
    r.heading = Heading(mod1(Int(r.heading) - 1, 4))
    r
end

function advance!(r::Robot)
    if heading(r) == NORTH
        r.position.y += 1
    elseif heading(r) == SOUTH
        r.position.y -= 1
    elseif heading(r) == EAST
        r.position.x += 1
    else
        r.position.x -= 1
    end
    r
end

function move!(r::Robot, instructions::AbstractString)
    moves = Dict('A' => advance!, 'R' => turn_right!, 'L' => turn_left!)
    for move in instructions
        move in keys(moves) && moves[move](r)
    end
    r
end


text:
function convert{N, T, PointType}(::Type{Array{T, 2}}, points::Vector{Point{N, PointType}})
    A = Array{T}(N, length(points))
    for i = 1:N
        for j = 1:length(points)
            A[i,j] = points[j][i]
        end
    end
    A
end

function convert{N, T, FaceType, Offset}(::Type{Array{T, 2}}, faces::Vector{GeometryTypes.Face{N, FaceType, Offset}})
    A = Array{T}(N, length(faces))
    for i = 1:N
        for j = 1:length(faces)
            A[i,j] = faces[j][i]
        end
    end
    A
end

function geom_data(mesh::AbstractMesh; color=[1;0;0;0.5])
    msg = lcmdrake[:lcmt_viewer_geometry_data]()
    msg[:type] = msg[:MESH]
    msg[:position] = [0;0;0]
    msg[:quaternion] = [1;0;0;0]
    msg[:color] = color
    msg[:string_data] = ""
    msg[:float_data] = Float64[length(vertices(mesh));
        length(faces(mesh));
        convert(Array{Float64,2}, vertices(mesh))[:];
        convert(Array{Float64,2}, map(f -> convert(Face{3, Int, -1}, f), faces(mesh)))[:];
    ]
    msg[:num_float_data] = length(msg[:float_data])
    msg
end

function link_data(mesh::AbstractMesh, name::AbstractString, robot_num::Integer=1; color=[1;0;0;0.5])
    msg = lcmdrake[:lcmt_viewer_link_data]()
    msg[:name] = name
    msg[:robot_num] = robot_num
    msg[:num_geom] = 1
    push!(msg["geom"], geom_data(mesh, color=color))
    msg
end

function link_data_list(robot::FixedObject, state::FixedObjectState, robot_num::Integer=1)
    lb = Vec(-1., -1, -1)
    ub = Vec(1.5, 1.5, 1.5)
    widths = ub - lb
    bounds = HyperRectangle(lb, widths)
    mesh = barrier_mesh(bounds, state.barrier)
    [link_data(mesh, "barrier", robot_num, color=[0;1;0;0.2])]
end

function link_data_list(robot::SoftRobot, state::SoftRobotState, robot_num::Integer=1)
    [link_data(HomogenousMesh(state.positions, robot.faces), "body", robot_num, color=[1;1;0;0.5])
     link_data(barrier_mesh(state), "barrier", robot_num, color=[1;0;0;0.2])]
 end

function draw(lcm::PyLCM.LCM, world, world_state)
    msg = lcmdrake[:lcmt_viewer_load_robot]()
    for (i, robot) in enumerate(world.objects)
        append!(msg["link"], link_data_list(robot, world_state[robot], i))
    end
    msg[:num_links] = length(msg[:link])
    publish(lcm, "DRAKE_VIEWER_LOAD_ROBOT", msg)
end

function draw(world, world_state)
    lcm = LCM()
    draw(lcm, world, world_state)
    finalize(lcm)
end

function draw{T}(robot, state::SoftRobotState{T})
    world = World3D()
    update_barrier!(robot, state)
    push!(world.objects, robot)
    world_state = Dict(robot => state)
    draw(world, world_state)
end