text:
module Consecutive

export consecutive

import Base: iterate, IteratorSize, IteratorEltype

macro ifsomething(ex)
    return quote
        result = $(esc(ex))
        result === nothing && return nothing
        result
    end
end

struct ConsecIter{I}
    xs::I
end

"""
    consecutive(xs)

Return an iterator that returns consecutive elements of `xs`.

```jldoctest
julia> a = 1:2:11
1:2:11

julia> collect(consecutive(a))
5-element Array{Tuple{Int64,Int64},1}:
 (1, 3)
 (3, 5)
 (5, 7)
 (7, 9)
 (9, 11)
```
"""
consecutive(i) = ConsecIter(i)

function iterate(it::ConsecIter)
    first, state = @ifsomething iterate(it.xs)
    second, state = @ifsomething iterate(it.xs, state)

    return ((first, second), (second, state))
end

function iterate(it::ConsecIter, s::Tuple)
    first, state = s
    second, state = @ifsomething iterate(it.xs, state)

    return ((first, second), (second, state))
end

restrict(x) = x
restrict(::Base.HasShape) = Base.HasLength()

# Forwards
Base.length(it::ConsecIter) = length(it.xs) - 1
Base.eltype(it::ConsecIter) = Tuple{eltype(it.xs),eltype(it.xs)}

# Shape gets destroyed when doing consecutive iteration. Since "HasLength()" is "weaker"
# than HasShape, restrict HasShape to HasLength. Pass everhthing else through.
IteratorSize(::ConsecIter{I}) where {I} = restrict(IteratorSize(I))
IteratorEltype(::ConsecIter{I}) where {I} = IteratorEltype(I)

end # module


text:
# Autogenerated wrapper script for cmocka_jll for x86_64-linux-gnu
export libcmocka

JLLWrappers.@generate_wrapper_header("cmocka")
JLLWrappers.@declare_library_product(libcmocka, "libcmocka.so.0")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libcmocka,
        "lib/libcmocka.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
using Bokeh
using Compat
using Base.Dates

start = now()
x = map(d -> start + Dates.Hour(d), 1:(30 * 24))
println(typeof(x))
# x: Array{DateTime,1}
day = 1000 * 3600 * 24
y = 0.5sin.(Dates.value.(x) * pi / day) .+ cos.(Dates.value.(x) * pi / (day * 7))
# note x needs to be an array of type DateTime or Date
plot(x, y)
showplot("dates_xy.html")

start = Date(1980, 1, 1)
x = map(d -> start + Dates.Day(d), 1:(365 * 30))
println(typeof(x))
# x: Array{Date,1}
y = 0.5sin.(Dates.value.(x) * pi / 365) .+ cos.(Dates.value.(x) * pi / (365 * 7))

plot(y, x)
showplot("dates_yx.html")


text:
# lets do some sort of exact gradient test

prob = StateTransferProblem(
    B = [Sx, Sy],
    A = Sz,
    X_init = ρinit,
    X_target = ρfin,
    duration = 1.0,
    n_timeslices = 10,
    n_controls = 2,
    initial_guess = rand(2, 10),
)
dt = prob.duration / prob.n_timeslices

pw_ham_save_bad(prob.A, prob.B, prob.initial_guess, prob.n_controls, dt, prob.n_timeslices)

function pw_ham_save_bad(
    H₀::T,
    Hₓ_array::Array{T,1},
    x_arr::Array{Float64},
    n_pulses,
    timeslices,
) where {T}
    out = Vector{T}(undef, timeslices)
    for i = 1:timeslices
        Htot = H₀
        for j = 1:n_pulses
            @views Htot = Htot + Hₓ_array[j] * x_arr[j, i]
        end
        out[i] = Htot
    end
    out
end


text:
module ThetisSupport

Base.@ccallable function bar(x :: Cint) :: Cint
    println("bar:", x)
    return x+3
end

precompile(bar, (Cint,))

end

max_stars_repo_name: AkillesAILimited/Thetis

max_issues_repo_name: AkillesAILimited/Thetis

max_forks_repo_name: AkillesAILimited/Thetis


text:
using Widgets, Observables, OrderedCollections, InteractBase, WebIO
using Colors
using Widgets: Widget, @layout, widgettype
import Widgets: observe, @auto
using Test

@testset "utils" begin
    d = Widget{:test}(Dict(:a => 1, :b => Observable(2), :c => Widget{:test}(; output = Observable(5))))
    m = d[:a] + d[:b][] + d[:c][]
    n = Observables.@map d[:a] + &d[:b] + d[:c][]
    @test m == 8
    @test n[] == 8
    d[:b][] = 3
    sleep(0.1)
    @test m == 8
    @test n[] == 9
    @test isa(Observables.@map(&d[:c]), Observable)

    v = [1]
    t = Widget{:test}(Dict(:a => Observable(2), :b => Observable(50)));
    Observables.@on v[1] += &t[:b]
    observe(t, :b)[] = 15
    sleep(0.1)
    @test v[1] == 16
    observe(t, :b)[] = 30
    sleep(0.1)
    @test v[1] == 46

    d = Widget{:test}(Dict(:a => 1, :b => Observable(2), :c => Widget{:test}(; output = Observable(5))))
    m = d |> @layout :a + :b[]
    n = d |> @layout Observables.@map :a + &(:b)
    @test m == 3
    @test n[] == 3
    d[:b][] = 3
    sleep(0.1)
    @test m == 3
    @test n[] == 4
    @test isa(@layout(d, :c), Widget)
end

function myui(x)
    a = x + 1
    b = Observable(10)
    output = Observables.@map &b + a
    @auto x = "aa"
    Widget{:myui}(
        ["a" => a, "b" => b, "x" => x],
        output = output,
        layout = t -> Observables.@map "The sum is "*string(&t)
    )
end

Widgets.widget(::Val{:myui}, args...; kwargs...) = myui(args...; kwargs...)

@testset "widget" begin
    ui = myui(5)
    @test ui[:a] == 6
    @test ui[:b][] == 10
    @test ui.output[] == 16
    @test ui.layout(ui)[] == "The sum is 16"
    @test widgettype(ui[:x]) == :textbox
    @test observe(ui[:x])[] == "aa"

    ui = Widgets.widget(Val(:myui), 5)
    @test ui[:a] == 6
    @test ui[:b][] == 10
    @test ui.output[] == 16
    @test ui.layout(ui)[] == "The sum is 16"

    ui = Widgets.@nodeps myui(5)
    @test ui[:a] == 6
    @test ui[:b][] == 10
    @test ui.output[] == 16
    @test ui.layout(ui)[] == "The sum is 16"

    ui[:b][] = 11
    sleep(0.1)
    @test ui.output[] == 17
    @test ui.layout(ui)[] == "The sum is 17"
end

@testset "auto" begin
    @auto x = 10
    @test observe(x)[] == 10
    @test widgettype(x) == :spinbox
end

@testset "layout" begin
    wdg = slider(1:100)
    wdg2 = wdg(style = Dict("color" => "red"))
    n = WebIO.render(wdg2)
    @test props(n)[:style]["color"] == "red"
    @test Widgets.node(wdg) isa Node
end

@testset "observable" begin
    t = Widget{:test}(Dict(:a => Observable(2), :b => Observable(50)), output = Observable(12));
    s = map(x->x-1, t)
    @test t[] == 12
    @test s[] == 11
end

@testset "manipulate" begin
    ui = @manipulate for r = 0:.05:1, g = 0:.05:1, b = 0:.05:1
        RGB(r,g,b)
    end
    @test observe(ui)[] == RGB(0.5, 0.5, 0.5)
    observe(ui, :r)[] = 0.1
    sleep(0.1)
    @test observe(ui)[] == RGB(0.1, 0.5, 0.5)
    ui = @manipulate throttle = 1 for r = 0:.05:1, g = 0:.05:1, b = 0:.05:1
        RGB(r,g,b)
    end
    observe(ui, :r)[] = 0.1
    sleep(0.1)
    observe(ui, :r)[] = 0.3
    sleep(0.1)
    observe(ui, :g)[] = 0.1
    sleep(0.1)
    observe(ui, :g)[] = 0.3
    sleep(0.1)
    observe(ui, :b)[] = 0.1
    sleep(0.1)
    observe(ui, :b)[] = 0.3
    sleep(0.1)
    @test observe(ui)[] != RGB(0.3, 0.3, 0.3)
    sleep(1.5)
    @test observe(ui)[] == RGB(0.3, 0.3, 0.3)
end


text:
#mutable clock structure to store time info
mutable struct Clock{T <: Real, N <: Integer}
    n_iter::N
    time::T
end

#clock constructor
function Clock(;
                n_iter = 0,
                time = 0)
    return Clock(n_iter, time)
end


text:
# legacy TU SpecialEuclidean

# TODO -- Change to type and overload the operators
# TODO -- uncomment
# function vee!(rv::Vector{Float64,1},T::SE2)
#   rv[1] = T.t[1]
#   rv[2] = T.t[2]
#   rv[3] = wrapRad(atan(-T.R.R[1,2], T.R.R[1,1]))
#   nothing
# end
# function *(a::SE2, b::SE2)
#   return SE2(R.R*b.R, vec(a.R.R*b.t + a.t))
# end

function SE2(X::AbstractArray{P,1}) where P <: Real
    T = Matrix{P}(LinearAlgebra.I, 3,3)
    T[1:2,1:2] = R(X[3])
    T[1,3] = X[1]
    T[2,3] = X[2]
    return T
end

function se2vee!(retval::AbstractArray{<:Real,1}, T::AbstractArray{<:Real,2})
    retval[1] = T[1,3]
    retval[2] = T[2,3]
    retval[3] = wrapRad(atan(-T[1,2], T[1,1]))
    nothing
end

function se2vee(T::AbstractArray{<:Real,2})
    # retval = zeros(3)
    # se2vee!(retval, T)
    # return retval
    return [T[1,3], T[2,3], wrapRad(atan(-T[1,2], T[1,1]))]
end

function veePose3(s::SE3)
  TransformUtils.veeEuler(s)
end
function veePose(s::SE3)
  TransformUtils.veeEuler(s)
end


text:
module MaxHelpingHandMovingFlagTouchSwitch

using ..Ahorn, Maple

@mapdef Entity "MaxHelpingHand/MovingFlagTouchSwitch" MovingFlagTouchSwitch(x::Integer, y::Integer,
    flag::String="moving_flag_touch_switch", icon::String="vanilla", persistent::Bool=false,
    inactiveColor::String="5FCDE4", movingColor::String="FF8080", activeColor::String="FFFFFF", finishColor::String="F141DF")

const placements = Ahorn.PlacementDict(
    "Flag Touch Switch (Moving)\n(max480's Helping Hand + Outback Helper)" => Ahorn.EntityPlacement(
        MovingFlagTouchSwitch
    )
)

const bundledIcons = String["vanilla", "tall", "triangle", "circle", "diamond", "double", "heart", "square", "wide", "winged"]

Ahorn.editingOptions(entity::MovingFlagTouchSwitch) = Dict{String,Any}(
    "icon" => bundledIcons
)

Ahorn.editingOrder(entity::MovingFlagTouchSwitch) = String["x", "y", "width", "height", "inactiveColor", "movingColor", "activeColor", "finishColor"]

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::MovingFlagTouchSwitch)
    px, py = Ahorn.position(entity)

    sprite = "collectables/outback/movingtouchswitch/container.png"

    for node in get(entity.data, "nodes", ())
        nx, ny = Int.(node)

        theta = atan(py - ny, px - nx)
        Ahorn.drawArrow(ctx, px, py, nx + cos(theta) * 8, ny + sin(theta) * 8, Ahorn.colors.selection_selected_fc, headLength=6)
        Ahorn.drawSprite(ctx, sprite, nx, ny)

        px, py = nx, ny
    end
end

function Ahorn.selection(entity::MovingFlagTouchSwitch)
    nodes = get(entity.data, "nodes", ())
    x, y = Ahorn.position(entity)

    sprite = "collectables/outback/movingtouchswitch/container.png"

    res = Ahorn.Rectangle[Ahorn.getSpriteRectangle(sprite, x, y)]

    for node in nodes
        nx, ny = Int.(node)

        push!(res, Ahorn.getSpriteRectangle(sprite, nx, ny))
    end

    return res
end

Ahorn.nodeLimits(entity::MovingFlagTouchSwitch) = 0, -1

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::MovingFlagTouchSwitch)
    icon = get(entity.data, "icon", "vanilla")

    iconPath = "objects/touchswitch/icon00.png"
    if icon != "vanilla"
        iconPath = "objects/MaxHelpingHand/flagTouchSwitch/$(icon)/icon00.png"
    end

    Ahorn.drawSprite(ctx, "collectables/outback/movingtouchswitch/container.png", 0, 0)
    Ahorn.drawSprite(ctx, iconPath, 0, 0)
end

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %% [markdown] slideshow={"slide_type": "slide"}
# # Van der Pol 振動子
#
# * Author: 黒木玄
# * Date: 2019-04-05～2019-04-15, 2021-07-14
# * Repository: https://github.com/genkuroki/DifferentialEquations
# $
# \newcommand\Z{{\mathbb Z}}
# \newcommand\R{{\mathbb R}}
# \newcommand\C{{\mathbb C}}
# \newcommand\eps{\varepsilon}
# \newcommand\QED{\text{□}}
# $
#
# このファイルは [nbviewer](https://nbviewer.jupyter.org/github/genkuroki/DifferentialEquations/blob/master/01-1%20Van%20der%20Pol%20oscillator.ipynb) でも閲覧できる.
#
# [Julia言語](https://julialang.org/) と [Jupyter環境](https://jupyter.org/) の簡単な解説については次を参照せよ:
#
# * [JuliaとJupyterのすすめ](https://nbviewer.jupyter.org/github/genkuroki/msfd28/blob/master/msfd28genkuroki.ipynb?flush_cached=true)
#
# [Julia言語](https://julialang.org/) 環境の整備の仕方については次を参照せよ:
#
# * [Julia v1.1.0 の Windows 8.1 へのインストール](https://nbviewer.jupyter.org/github/genkuroki/msfd28/blob/master/install.ipynb)

# %% [markdown] slideshow={"slide_type": "-"} toc=true
#
目次

#

1  Van der Pol 振動子

2  Van der Pol 方程式のベクトル場と流れの図

2.1  $\mu = 0$ の場合

2.2  $\mu$ が正の場合

2.3  $\mu$ が負の場合

3  Van der Pol 方程式の数値解

3.1  $\mu=0$

3.2  $\mu = 0.1$

3.3  $\mu=0.5$

3.4  $\mu=1$

4  正弦波による強制力付きの場合

# %% slideshow={"slide_type": "-"}
using PyPlot: PyPlot, plt

# %% slideshow={"slide_type": "-"}
function plot_stream(f, g;
        x = range(-4, 4, length=201),
        y = range(-4, 4, length=201),
        density = 1.2, figtitle="", sign=1.0)

    # meshgrid
    xx, yy = reim(complex.(x', y))
    XX = f.(xx, sign*yy)
    YY = g.(xx, sign*yy) * sign

    plt.streamplot(xx, yy, XX, YY, linewidth=0.5, density=density, color="blue")
    plt.xlim(extrema(x)...)
    plt.ylim(extrema(y)...)
    plt.grid(ls=":")
    figtitle == "" || plt.title(figtitle)
end

function plot_vector_field(f, g;
        x = range(-4, 4, length=21),
        y = range(-4, 4, length=21),
        figtitle="", scale=2.0, sign=1.0)

    # meshgrid
    xx, yy = reim(complex.(x', y))
    XX = f.(xx, sign*yy)
    YY = g.(xx, sign*yy) * sign
    MM = maximum(@.(√(XX^2+YY^2)))
    #println(MM) # for debug

    plt.quiver(xx, yy, XX, YY, scale=scale*MM, color="red")
    plt.xlim(extrema(x)...)
    plt.ylim(extrema(y)...)
    plt.grid(ls=":")
    figtitle == "" || plt.title(figtitle)
end

# %% [markdown]
# ## Van der Pol 振動子
#
# [**Van der Pol (ファン・デル・ポル) 振動子**](https://en.wikipedia.org/wiki/Van_der_Pol_oscillator) とは以下の非線形常微分方程式で記述される系のことである:
#
# $$
# \ddot{x} - \mu(1-x^2)\dot{x} + x = 0.
# $$
#
# ここで $x=x(t)$ の時刻 $t$ に関する1階と2階の導函数を $\dot{x}$, $\ddot{x}$ と表した.  この微分方程式を **Van der Pol 方程式** と呼ぶ.

# %% [markdown] slideshow={"slide_type": "slide"}
# ## Van der Pol 方程式のベクトル場と流れの図
#
# $v = \dot{x}$ とおくと, Van der Pol 方程式は次の連立常微分方程式に書き直される:
#
# $$
# \begin{cases}
# \dot{x} = v, \\
# \dot{v} = \mu(1-x^2)v - x.
# \end{cases}
# $$
#
# この微分方程式は位置 $x$ と速度 $v$ を平面上の点 $(x,v)$ で表すとき, その点の動きの速度ベクトルが $(v, \mu(1-x^2)v-x)$ となることを意味している.  平面上の各点にその点における速度ベクトルを対応させる函数を平面上の **ベクトル場** と呼ぶ. 上の連立常微分方程式はベクトル場に沿って平面上の点が流れて行く様子を表している.

# %% slideshow={"slide_type": "-"}
function plot_Van_der_Pol(μ; scale=1/√abs(μ), sign=1.0)
    # Van der Pol equation
    f(x,v) = v
    g(x,v) = μ*(1-x^2)*v - x

    plt.figure(figsize=(7,3.5))
    plt.subplot("121")
    plot_stream(f, g, figtitle="Van der Pol: μ = $μ", sign=sign)
    plt.subplot("122")
    plot_vector_field(f, g, figtitle="Van der Pol: μ = $μ", scale=scale, sign=sign)
    plt.tight_layout()
end

# %% [markdown] slideshow={"slide_type": "slide"}
# ### $\mu = 0$ の場合
#
# $\mu=0$ のとき, Van der Pol 方程式は
#
# $$
# \begin{cases}
# \dot{x} = v, \\
# \dot{v} = - x.
# \end{cases}
# $$
#
# の形になる.  これが定める流れは時計回りで回転する流れになる. その一般解は, $r\geqq 0$ と実数 $a$ を使って次のように書ける:
#
# $$
# \begin{cases}
# \dot{x} = r\sin(t-a), \\
# \dot{v} = r\cos(t-a).
# \end{cases}
# $$

# %% slideshow={"slide_type": "subslide"} tags=[]
plot_Van_der_Pol(0.0, scale=10.0)

# %% [markdown] slideshow={"slide_type": "slide"} tags=[]
# ### $\mu$ が正の場合

# %% slideshow={"slide_type": "-"}
plot_Van_der_Pol(0.1)

# %% slideshow={"slide_type": "subslide"}
plot_Van_der_Pol(0.5)

# %% slideshow={"slide_type": "subslide"}
plot_Van_der_Pol(1.0)

# %% [markdown] slideshow={"slide_type": "slide"}
# ### $\mu$ が負の場合
#
# Van der Pol 方程式
#
# $$
# \ddot{x} - \mu(1-x^2)\dot{x} + x = 0.
# $$
#
# で時刻 $t$ が進む向きを逆転させると, すなわち, $t$ に $-t$ を代入すると,
#
# $$
# \ddot{x} + \mu(1-x^2)\dot{x} + x = 0.
# $$
#
# の形になる. すなわち, 時間の向きを逆転させることと, パラメーター $\mu$ をその $-1$ 倍で置き換えることは等しい.
#
# 以下の $\mu=0.2$ と $\mu=-0.2$ の場合のプロットを比較すると, 矢印の向きが逆になっていることがわかる.

# %% slideshow={"slide_type": "-"}
plot_Van_der_Pol(0.2)

# %% slideshow={"slide_type": "-"}
plot_Van_der_Pol(-0.2, sign=-1.0)

# %% [markdown] slideshow={"slide_type": "-"}
# **注意:** 上のセルの `sign=-1.0` はプロットするときに縦軸($v$ 軸)の向きを反転してプロットすることを意味している.

# %% [markdown] slideshow={"slide_type": "slide"}
# ## Van der Pol 方程式の数値解

# %%
using DifferentialEquations
using Plots
default(:bglegend, plot_color(default(:bg), 0.7))
default(:fglegend, plot_color(ifelse(isdark(plot_color(default(:bg))), :white, :black), 0.6))

gr()
ENV["PLOTS_TEST"] = "true";

# %%
### Van der Pol 方程式の記述
function vanderpol!(du, u, p, t)
    # u[1] = x
    # u[2] = v
    # du[1] = dx/dt
    # du[2] = dv/dt
    # p[1] = μ

    # dx/dt = v
    du[1] = u[2]
    # dv/dt = μ (1 - x^2) v - x
    du[2] = p[1]*(1 - u[1]^2)*u[2] - u[1]
end

function solve_vanderpol(μ, x0, v0, tmax)
    u0 = [x0, v0] # initial values
    p = [μ] # parameters
    tspan = (0.0, tmax) # time span
    prob = ODEProblem(vanderpol!, u0, tspan, p) # ordinary differential equation problem
    solve(prob)
end

function plot_vanderpol_sol(sol, μ; l1=:best, l2=:best)
    plot(size=(600, 250))
    plot!(sol, label=["x" "v"], legend=l1, lw=1)
    title!("Van der Pol equation (mu = $μ)", titlefontsize=12) |> display

    plot(size=(300, 300))
    plot!(sol, vars=(1,2), label="(x,v)", legend=l2, lw=1)
    title!("Van der Pol equation (mu = $μ)", titlefontsize=10) |> display
end

function solve_and_plot_vanderpol(; μ=0.5, x0=0.1, v0=0.0, tmax=50.0, l1=:topleft, l2=:best)
    sol =  solve_vanderpol(μ, x0, v0, tmax)
    plot_vanderpol_sol(sol, μ; l1=l1, l2=l2)
end

# %% [markdown]
# ### $\mu=0$
#
# $\mu = 0$ の場合に van der Pol 方程式の解は $(x,v)$ 平面上の円運動になる.

# %%
solve_and_plot_vanderpol(; μ=0.0, x0=1.0, v0=0.0, tmax=10.0, l1=:topright)

# %% [markdown]
# ### $\mu = 0.1$
#
# $\mu = 0.1$ の場合に, van der Pol 方程式の解は $(x,v)$ 平面上の歪んだ円に内側と外側から巻き付くような運動になる. 巻き付く先の歪んだ円は **limit cycle** と呼ばれる.

# %%
solve_and_plot_vanderpol(; μ=0.1, x0=0.1, v0=0.0, tmax=100.0)

# %%
solve_and_plot_vanderpol(; μ=0.1, x0=-8.0, v0=0.0, tmax=40.0, l1=:bottomright, l2=:bottomleft)

# %% [markdown]
# ### $\mu=0.5$
#
# $\mu$ を大きくすると limit cycle の形が円から離れて行く.

# %%
solve_and_plot_vanderpol(; μ=0.5, x0=0.01, v0=0.0, tmax=60.0, l2=:topleft)

# %%
solve_and_plot_vanderpol(; μ=0.5, x0=-8.0, v0=0.0, tmax=60.0, l1=:bottomright, l2=:bottomleft)

# %% [markdown]
# ### $\mu=1$
#
# $\mu = 1$ での limit cycle は斜めに傾いた平行四辺形に近い形になる.

# %%
solve_and_plot_vanderpol(; μ=1.0, x0=0.1, v0=0.0, tmax=30.0, l2=:topleft)

# %%
solve_and_plot_vanderpol(; μ=1.0, x0=-8.0, v0=0.0, tmax=60.0, l1=:bottomright, l2=:topleft)

# %% [markdown]
# ## 正弦波による強制力付きの場合
#
# $$
# \begin{cases}
# \dot{x} = v, \\
# \dot{v} = \mu(1-x^2)v - x + A\sin\omega t.
# \end{cases}
# $$

# %%
using DifferentialEquations
using Plots
default(:bglegend, plot_color(default(:bg), 0.7))
default(:fglegend, plot_color(ifelse(isdark(plot_color(default(:bg))), :white, :black), 0.6))

gr()
ENV["PLOTS_TEST"] = "true";

# %%
### sin による強制力付きの Van der Pol 方程式の記述
function vanderpolsin!(du, u, p, t)
    # x = u[1]
    # v = u[2]
    # dx/dt = du[1]
    # dv/dt = du[2]
    # μ = p[1]
    # A = p[2]
    # ω = p[3]

    # dx/dt = v
    du[1] = u[2]
    # dv/dt = μ (1 - x^2) v - x + A sin ωt
    du[2] = p[1]*(1 - u[1]^2)*u[2] - u[1] + p[2]*sin(p[3]*t)
end

function solve_vanderpolsin(; μ=1.0, A=1.2, ω=2π/10, x0=0.1, v0=0.0, tmax=100)
    u0 = [x0, v0] # initial values
    p = [μ, A, ω] # parameters
    tspan = (0.0, tmax) # time span
    prob = ODEProblem(vanderpolsin!, u0, tspan, p) # ordinary differential equation problem
    solve(prob)
end

function plot_vanderpol_sol(sol, μ, A, ω; l1=:best, l2=:best)
    plot(size=(600, 200))
    plot!(sol, vars=(0,1), label="x", legend=l1, lw=1)
    title!("Forced Van der Pol (mu = $μ, A = $A, omega = $(round(ω, digits=2)))", titlefontsize=10) |> display

    sleep(0.1)
    plot(size=(300, 300))
    plot!(sol, vars=(1,2), label="(x,v)", legend=l2, lw=1)
    title!("mu = $μ, A = $A, omega = $(round(ω, digits=2))", titlefontsize=10) |> display
end

function solve_and_plot_vanderpolsin(; μ=1.0, A=1.2, ω=2π/5, x0=0.1, v0=0.0, tmax=50.0, l1=false, l2=false)
    sol =  solve_vanderpolsin(μ=μ, A=A, ω=ω, x0=x0, v0=v0, tmax=tmax)
    plot_vanderpol_sol(sol, μ, A, ω; l1=l1, l2=l2)
end

# %%
solve_and_plot_vanderpolsin(; μ=8.53, A=1.2, ω=2π/10, x0=0.1, v0=0.0, tmax=300.0)

# %% [markdown]
# カオス的な動き方をしている.

max_stars_repo_path: 2021-07-14 Van der Pol.jl

max_stars_repo_name: genkuroki/DIfferentialEquations

max_issues_repo_path: 2021-07-14 Van der Pol.jl

max_issues_repo_name: genkuroki/DIfferentialEquations

max_forks_repo_path: 2021-07-14 Van der Pol.jl

max_forks_repo_name: genkuroki/DIfferentialEquations


text:
using Images, ShapeFromShading
using Test

@testset "ShapeFromShading" begin
    include("deterministic.jl")
    include("deterministic_custom.jl")
    include("discreteshape.jl")
    include("discreteshapebound.jl")
    include("multiresolutionhybrid.jl")
    include("pentland.jl")
    include("photometric.jl")
    include("shah.jl")
    # include("SimulatedAnnealing.jl")
end

@testset "Anisotropicdiffusion" begin
    include("anisotropicdiffusion.jl")
end
@testset "Durou" begin
    include("durou.jl")
end
@testset "Frankot" begin
    include("frankot.jl")
end
@testset "Horn" begin
    include("horn.jl")
end
@testset "Mumfordshah" begin
    include("mumfordshah.jl")
end
@testset "Nonconvex1" begin
    include("nonconvex1.jl")
end
@testset "Path" begin
    include("path.jl")
end
@testset "Splitpath" begin
    include("splitpath.jl")
end
@testset "Totalvariation" begin
    include("totalvariation.jl")
end
@testset "Quadratic" begin
    include("quadratic.jl")
end

@testset "SyntheticSurface" begin
    include("cake.jl")
    include("cake2.jl")
    include("dem.jl")
    include("prism.jl")
    include("ripple.jl")
    include("ripple2.jl")
    include("supergaussian.jl")
    include("synthsphere.jl")
    include("synthvase.jl")
    include("synthgaussian.jl")
    include("tent.jl")
end

@testset "ImageGeneration" begin
    include("estimatealbedo.jl")
    include("generate_surface.jl")
end

@testset "OtherFunctions" begin
    include("mask.jl")
end

@testset "Benchmarks" begin
    include("benchmark_shape.jl")
    include("benchmark_integration.jl")
    include("benchmark_normals.jl")
    include("benchmark_image.jl")
end


text:
using Documenter, squirrel

makedocs(
    sitename    = "squirrel.jl",
    format      = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    authors     = "Justin C. Feng and Filip Hejda and Sante Carloni",
    pages = [
                "Home"               => "index.md",
                "Metrics"            => "metrics.md",
                "Geodesic solver"   => "geodesics.md",
                "Squirrel algorithm" => "squirrel.md",
                "Evaluation"         => "evaluation.md",
            ]
        )


text:
@enum Surface begin
    Shiny        = 1
    CheckerBoard = 2
end

struct Vector
    x::Float64
    y::Float64
    z::Float64
end

struct Color
    r::Float64
    g::Float64
    b::Float64
end

struct RGBColor
    b::UInt8
    g::UInt8
    r::UInt8
    a::UInt8

    function RGBColor(b::UInt8, g::UInt8, r::UInt8, a::UInt8)
        new(b, g, r, a)
    end
end

struct Image
    width::Int
    height::Int
    data::Array{RGBColor,1}

    function Image(w::Int, h::Int)
        data = Array{RGBColor,1}(undef, w * h);
        new(w, h, data)
    end
end

struct Camera
    pos::Vector
    forward::Vector
    right::Vector
    up::Vector

    function Camera(pos::Vector, lookAt::Vector)
        _down = Vector(0.0, -1.0, 0.0);
        _forward = lookAt - pos
        pos = pos
        forward = Norm(_forward)
        right = Norm(Cross(forward, _down)) * 1.5
        up = Norm(Cross(forward, right)) * 1.5
        new(pos, forward, right, up)
    end
end

struct Ray
    start::Vector
    dir::Vector
end

abstract type Thing end

struct Plane <: Thing
    normal::Vector
    offset::Float64
    surface::Surface
end

struct Sphere <: Thing
    center::Vector
    radius2::Float64
    surface::Surface
end

struct Intersection
    thing::Thing
    ray::Ray
    dist::Float64
end

struct Light
    pos::Vector
    color::Color
end

struct SurfaceProperties
    diffuse::Color
    specular::Color
    reflect::Float64
    roughness::Float64
end

const WORD  = UInt16
const DWORD = UInt32
const LONG  = Int32

struct BITMAPINFOHEADER
    biSize::DWORD
    biWidth::LONG
    biHeight::LONG
    biPlanes::WORD
    biBitCount::WORD
    biCompression::DWORD
    biSizeImage::DWORD
    biXPelsPerMeter::LONG
    biYPelsPerMeter::LONG
    biClrUsed::DWORD
    biClrImportant::DWORD

    function BITMAPINFOHEADER(w::LONG, h::LONG)
        size = sizeof(BITMAPINFOHEADER)
        sizeImage = w * h * 4
        new(size, w, -h, 1, 32, 0, sizeImage, 0, 0, 0, 0)
    end
end

struct BITMAPFILEHEADER
    bfType::WORD
    bfSize::DWORD
    bfReserved::DWORD
    bfOffBits::DWORD

    function BITMAPFILEHEADER(imageSize::DWORD)
        bfType =  0x4D42 # 'B' + ('M' << 8)
        offBits = 54
        size = convert(DWORD, offBits + imageSize);
        new(bfType, size, 0, offBits)
    end
end

struct Scene
    things::Array{Thing,1}
    lights::Array{Light,1}
    camera::Camera

    function Scene()
        things = [
            Plane(Vector(0.0, 1.0, 0.0), 0.0,  CheckerBoard)
            Sphere(Vector(0.0, 1.0, -0.25), 1.0, Shiny)
            Sphere(Vector(-1.0, 0.5, 1.5), 0.5*0.5, Shiny)
        ]
        lights = [
            Light(Vector(-2.0, 2.5, 0.0), Color(0.49, 0.07, 0.07))
            Light(Vector(1.5, 2.5, 1.5), Color(0.07, 0.07, 0.49))
            Light(Vector(1.5, 2.5, -1.5), Color(0.07, 0.49, 0.071))
            Light(Vector(0.0, 3.5, 0.0), Color(0.21, 0.21, 0.35))
        ]
        camera = Camera(Vector(3.0, 2.0, 4.0), Vector(-1.0, 0.5, 0.0));

        new(things, lights, camera)
    end
end

const maxDepth = 5
const White = Color(1.0, 1.0, 1.0)
const Grey  = Color(0.5, 0.5, 0.5)
const Black = Color(0.0, 0.0, 0.0)
const Background = Black
const Defaultcolor = Black

Length(v::Vector) = sqrt(v.x^2 + v.y^2 + v.z^2)
Scale(v::Vector, k::Float64) =  Vector(k * v.x, k * v.y, k * v.z)
Dot(a::Vector, b::Vector) =  a.x * b.x + a.y * b.y + a.z * b.z

Base.:*(v::Vector, k::Float64) = Scale(v, k)
Base.:*(a::Vector, b::Vector) = Dot(a, b)
Base.:+(a::Vector, b::Vector) = Vector(a.x + b.x, a.y + b.y, a.z + b.z)
Base.:-(a::Vector, b::Vector) = Vector(a.x - b.x, a.y - b.y, a.z - b.z)

function Norm(v::Vector)
    len = Length(v)
    div = (len == 0) ? floatmax(Float32) : 1.0 / len
    Scale(v, div)
end

function Cross(a::Vector, b::Vector)
    return Vector(
        a.y * b.z - a.z * b.y,
        a.z * b.x - a.x * b.z,
        a.x * b.y - a.y * b.x
    )
end

Base.:+(a::Color, b::Color) = Color(a.r + b.r, a.g + b.g, a.b + b.b)
Base.:-(a::Color, b::Color) = Color(a.r - b.r, a.g - b.g, a.b - b.b)
Base.:*(a::Color, b::Color) = Color(a.r * b.r, a.g * b.g, a.b * b.b)

Scale(v::Color, k::Float64) = Color(k * v.r, k * v.g, k * v.b)
GetNormal(p::Plane, pos::Vector) = p.normal
GetNormal(s::Sphere, pos::Vector) = Norm(pos - s.center)

function ToRgbColor(c::Color)::RGBColor
    b = floor(UInt8, clamp(c.b, 0.0, 1.0) * 255.0)
    g = floor(UInt8, clamp(c.g, 0.0, 1.0) * 255.0)
    r = floor(UInt8, clamp(c.r, 0.0, 1.0) * 255.0)
    a = convert(UInt8, 255)
    return RGBColor(b, g, r, a)
end

function GetPoint(camera::Camera, x::Int, y::Int, screenWidth::Int, screenHeight::Int)
    recenterX = (x - (screenWidth / 2.0)) / 2.0 / screenWidth
    recenterY = -(y - (screenHeight / 2.0)) / 2.0 / screenHeight
    return Norm(camera.forward + ((camera.right * recenterX) + (camera.up * recenterY)))
end

function GetInteresection(p::Plane, ray::Ray)::Union{Intersection,Nothing}
    denom = Dot(p.normal, ray.dir)
    if (denom > 0.0)
        return nothing
    end
    dist = (Dot(p.normal, ray.start) + p.offset) / (-denom)
    Intersection(p, ray, dist)
end

function GetInteresection(s::Sphere, ray::Ray)::Union{Intersection,Nothing}
    eo = (s.center - ray.start)
    v = Dot(eo, ray.dir)
    dist = 0.0
    if (v >= 0.0)
        disc = s.radius2 - (Dot(eo, eo) - v * v)
        dist = (disc >= 0.0) ? v - sqrt(disc) : dist
    end
    (dist == 0.0) ? nothing : Intersection(s, ray, dist)
end

function GetSurfaceProperties(s::Surface, pos::Vector)::SurfaceProperties
    if s == Shiny
        return SurfaceProperties(White, Grey, 0.7, 250.0)
    else
        condition = Int32(floor(pos.x) + floor(pos.z)) & 1 != 0
        color = condition ? White : Black
        reflect = condition ? 0.1 : 0.7
        SurfaceProperties(color, White, reflect, 250.0)
    end
end

function Save(image::Image, fileName::String)

    w = convert(LONG, image.width)
    h = convert(LONG, image.height)

    infoHeader = BITMAPINFOHEADER(w, h)

    offBits = convert(DWORD, 54)
    size  = convert(DWORD, offBits + infoHeader.biSizeImage)

    fileHeader = BITMAPFILEHEADER(infoHeader.biSizeImage)

    f = open(fileName, "w")
    write(f, fileHeader.bfType)
    write(f, fileHeader.bfSize)
    write(f, fileHeader.bfReserved)
    write(f, fileHeader.bfOffBits)

    write(f, infoHeader.biSize)
    write(f, infoHeader.biWidth)
    write(f, infoHeader.biHeight)
    write(f, infoHeader.biPlanes)
    write(f, infoHeader.biBitCount)
    write(f, infoHeader.biCompression)
    write(f, infoHeader.biSizeImage)
    write(f, infoHeader.biXPelsPerMeter)
    write(f, infoHeader.biYPelsPerMeter)
    write(f, infoHeader.biClrUsed)
    write(f, infoHeader.biClrImportant)

    for c in image.data
        write(f, c.b)
        write(f, c.g)
        write(f, c.r)
        write(f, c.a)
    end
    close(f)
end

function GetClosestIntersection(scene::Scene, ray::Ray)::Union{Intersection,Nothing}
    closest = floatmax(Float64)
    closestInter = nothing

    for thing in scene.things
        inter = GetInteresection(thing, ray)
        if (inter != nothing && inter.dist < closest)
            closestInter = inter
            closest = inter.dist
        end
    end
    return closestInter;
end

function TraceRay(scene::Scene, ray::Ray, depth::Int)
    isect = GetClosestIntersection(scene, ray);
    isect != nothing ? Shade(scene, isect, depth) : Background
end

function Shade(scene::Scene, isect::Intersection, depth::Int)
    d = isect.ray.dir;
    pos = (d * isect.dist) + isect.ray.start;
    normal = GetNormal(isect.thing, pos);
    reflectDir = Norm(d - ((normal * (normal * d)) * 2.0))

    surface =  GetSurfaceProperties(isect.thing.surface, pos);

    naturalColor = Background + GetNaturalColor(scene, surface, pos, normal, reflectDir);
    reflectedColor = (depth >= maxDepth) ? Grey : GetReflectionColor(scene, surface, pos, reflectDir, depth);

    return naturalColor + reflectedColor;
end

function GetReflectionColor(scene::Scene, surface::SurfaceProperties, pos::Vector, reflectDir::Vector, depth::Int)
    ray = Ray(pos, reflectDir)
    color = TraceRay(scene, ray, depth + 1)
    return Scale(color, surface.reflect)
end

function GetNaturalColor(scene::Scene, surface::SurfaceProperties, pos::Vector, norm::Vector, reflectDir::Vector)
    result = Black

    for light in scene.lights
        ldist = light.pos - pos
        livec = Norm(ldist)
        ray = Ray(pos, livec)

        neatIsect = GetClosestIntersection(scene, ray)
        isInShadow = neatIsect != nothing && neatIsect.dist < Length(ldist)

        if (!isInShadow)
            illum    = livec * norm
            specular = livec * reflectDir

            lcolor = (illum > 0) ? Scale(light.color, illum) : Defaultcolor
            scolor = (specular > 0) ? Scale(light.color, specular^surface.roughness) : Defaultcolor
            result = result + (lcolor * surface.diffuse) + (scolor * surface.specular)
        end
    end
    return result
end

function Render(scene::Scene)
    camera = scene.camera
    image = Image(500, 500)

    w = image.width - 1
    h = image.height - 1

    for y in 0:h, x in 0:w
        pt    = GetPoint(camera, x, y, w, h)
        ray   = Ray(scene.camera.pos, pt)
        color = TraceRay(scene, ray, 0)
        image.data[y * image.width + x + 1] = ToRgbColor(color)
    end
    return image
end

scene = Scene()
image = Render(scene)
@time begin
    scene = Scene()
    image = Render(scene)
end
Save(image, "julia-ray.bmp")

max_stars_repo_name: alkalinin/raytracer

max_issues_repo_name: alkalinin/raytracer

max_forks_repo_name: alkalinin/raytracer


text:
# ω0 => fator inercial inicial
# ωf => fator inercial final
# ϕp => coeficiente de aceleração partícula
# ϕg => coeficiente de aceleração grupo
# η => passo incremental para criação de nova solução
# λ => decaimento da temperatura

function hpso(func::Function, particles::Particles, constraints, args, kwargs,
    ω0, ωf , ϕp, ϕg, η, λ , maxiter, maxiterSA, minstep, minfunc, verbose)

    obj = x -> func(x, args...; kwargs...)
    cons = make_constraints(constraints, args, kwargs, verbose)
    is_feasible = x -> all(cons(x) .>= 0)

    S = particles.swarmsize
    D = particles.particlelength
    verbose && println("Initialize current particle function values.")
    particles.fx = [obj(particles.x[i]) for i = 1:S]  # current particle function values
    particles.fs = [is_feasible(particles.x[i]) for i = 1:S]  # feasibility of each particle

    particles.history = []
    particles.fp = copy(particles.fx)  # best particle function values
    particles.order = sortperm(particles.fx, rev=true) #particle neighborhood order

    particles.g = copy(particles.x)  # best neighborhood swarm position
    particles.fg = ones(S) * Inf  # best neighborhood swarm position starting value

    # procurar alguma solução feasible e depois
    # pegar a com o menor valor dentre estas
    if .|(particles.fs...)
        aux1, aux2 = findmin(particles.fx[particles.fs])
        particles.best_position = copy(particles.x[particles.fs][aux2]) # best position
        particles.best_value = aux1 # best position value
    else
        aux1, aux2 = findmin(particles.fx)
        particles.best_position = copy(particles.x[aux2]) # best position
        particles.best_value = aux1 # best position value
    end

    # Store particle's best position (if constraints are satisfied)
    verbose && println("Store particle's best position.")
    update_position!(particles)

    # Update swarm's best position
    verbose && println("Update swarm's best position.")
    findBestGlobal!(particles)

    #temperature
    temperature = -(maximum(particles.fp) - minimum(particles.fp) )/ log(0.1)

    it = 1
    it_best = 1
    verbose && println("Start iteration.")
    while it <= maxiter

        # Update the particles' velocities and positions using equations 1 and 2
        # Eq(1): v = ω*v + ϕp . rp . (p -x) + ϕg . rg . (g - x)
        # Eq(2): x = x +v
        ω = ωf + (ω0 - ωf) * ( (maxiter - it) / (maxiter - 1) )
        rp = rand(particles.swarmsize) # random value
        rg = rand(particles.swarmsize) # random value
        particles.v = (ω .* particles.v) .+ (ϕp .* rp .* (particles.p .- particles.x)) .+
        (ϕg .* rg .* (particles.g .- particles.x))
        particles.x = particles.x .+ particles.v

        # Correct for bound violations
        maskl = [particles.x[i] .< particles.lb for i = 1:S]
        masku = [particles.x[i] .> particles.ub for i = 1:S]
        mask = [ maskl[i] .| masku[i] for i = 1:S]
        particles.x = map((x,y) -> x.*(.~(y)), particles.x, mask)
        particles.x .+= [particles.lb.*maskl[i] for i = 1:S]
        particles.x .+= [particles.ub.*masku[i] for i=1:S]

        # Update objectives and constraints
        particles.fx = map(x-> obj(x), particles.x)
        particles.fs = map(x-> is_feasible(x), particles.x)

        # Store particle's best position (if constraints are satisfied)
        update_position!(particles)

        # Busca local pelo algoritm SA conforme artigo
        verbose && println("Local Search.")
        # TODO: lembrar de no IHPSO fazer essa parte aqui ser um loop com as melhores posições
        indextemp = argmin(particles.fp)
        gtemp = copy(particles.p[indextemp])

        p_feasible = [is_feasible(particles.p[i]) for i = 1:S]
        if .|(p_feasible...)
            indextemp = argmin(particles.fp[p_feasible])
            gtemp = copy(particles.p[p_feasible][indextemp])
            indextemp = findfirst(x -> x == gtemp, particles.p)
        end

        patemp = 1.

        itSA = 1.
        # TODO: remover esse print aqui depois de apresentar os resultados
        verbose && println("Local Search best at start $(gtemp) $(obj(gtemp)).")
        while itSA <= maxiterSA
            xtemp = gtemp + η * (particles.ub - particles.lb) .* randn(D)

            # Correct for bound violations
            maskltemp = xtemp .< particles.lb
            maskutemp = xtemp .> particles.ub
            mask = maskltemp .| maskutemp
            xtemp = xtemp.*(.~(mask))
            xtemp .+= particles.lb.*maskltemp
            xtemp .+= particles.ub.*maskutemp

            fxtemp = obj(xtemp)
            if (is_feasible(xtemp) && is_feasible(gtemp))
                patemp = min(1., (exp((obj(gtemp) - fxtemp)/temperature)) )
            elseif (!is_feasible(xtemp) && is_feasible(gtemp))
                patemp = 0.0
            elseif (is_feasible(xtemp) && !is_feasible(gtemp))
                patemp = 1.
            end

            if patemp > rand()
                #verbose && println("$(itSA) $(patemp) e $(obj(gtemp)) e $(obj(xtemp)) $(temperature)")
                gtemp = xtemp
            end

            itSA += 1
        end
        particles.p[indextemp] = gtemp
        particles.fp[indextemp] = obj(gtemp)

        verbose && println("Local Search best at end $(gtemp) $(obj(gtemp)).")

        # Update swarm's best position
        findBestGlobal!(particles)

        # Compare swarm's best position with global best position
        i_min = findmin(particles.fp)
        position_temp = copy(particles.p[i_min[2]])

        if .|(p_feasible...)
            i_min = findmin(particles.fp[p_feasible])
            position_temp = particles.p[p_feasible][i_min[2]]
        end

        if i_min[1] < particles.best_value
            it_best = it
            stepsize = √(sum((particles.best_position .- position_temp).^2))
            if abs.(particles.best_value .- i_min[1]) <= minfunc
                verbose && println("\n Stopping search: Swarm best objective change
                less than $(minfunc)")
                break
            end
            if stepsize <= minstep
                verbose && println("Stopping search: Swarm best position change
                less than $(minstep)")
                break
            end
            # if feasible change value
            if is_feasible(position_temp)
                particles.best_position = copy(position_temp) # best position
                particles.best_value = i_min[1] # best position value
            elseif !is_feasible(position_temp) && !is_feasible(particles.best_position)
                particles.best_position = copy(position_temp) # best position
                particles.best_value = i_min[1] # best position value
            end
        end
        push!(particles.history, particles.best_value)
        verbose && (print("\r Iteration : $(it)/$(maxiter) | best-iteration $(it_best)
        | Best value: $(particles.best_value) "))
        it += 1
    end
    verbose && println("\n Stopping search: maximum iterations reached --> $(maxiter)")
    is_feasible(particles.best_position) || print("However, the optimization couldn't find a
     feasible design. Sorry")
    return nothing
end

function hpso(func::Function, lb::Vector, ub::Vector; constraints=nothing, args=(), kwargs=Dict(),
    swarmsize=100, omegastart=0.9, omegaend=0.4, phip=0.5, phig=0.5, eta = 0.001, lambda = 0.94,
     maxiter=100, maxiterSA = 30, minstep=1e-8, minfunc=1e-8, verbose=false)

@assert length(ub) == length(lb) "ub and lb must have same dimension"
@assert all(ub .> lb) "Each value of 'ub' must be greater than 'lb'"
particles = Particles(swarmsize, lb, ub) #create particles

hpso(func, particles, constraints, args, kwargs, omegastart, omegaend, phip,
 phig, eta, lambda, maxiter, maxiterSA, minstep, minfunc, verbose)
return particles
end

function hpso(particles:: Particles, func::Function; constraints=nothing, args=(), kwargs=Dict(),
    omegastart=0.9, omegaend=0.4, phip=0.5, phig=0.5, eta = 0.001, lambda = 0.94,
    maxiter=100, maxiterSA = 30, minstep=1e-8, minfunc=1e-8, verbose=false)

hpso(func, particles, constraints, args, kwargs, omegastart, omegaend, phip,
 phig, eta, lambda, maxiter, maxiterSA, minstep, minfunc, verbose)

end


text:
struct DefaultLogger{IO,Fallback<:AbstractLogger} <: AbstractLogger
    io::IO
    fallback::Fallback
end

DefaultLogger() = DefaultLogger(stderr, current_logger())

function Logging.handle_message(
    logger::DefaultLogger,
    level::Logging.LogLevel,
    message,
    _module,
    group,
    id,
    filepath,
    line;
    exception = nothing,
    kwargs...,
)
    if exception !== nothing
        handle_message(
            logger.fallback,
            level,
            message,
            _module,
            group,
            id,
            filepath,
            line;
            exception,
            kwargs...,
        )
    end
    printing_in_oneline(logger.io) do io, should_stop
        if level != Logging.Info
            printstyled(io, level; color = :blue)
            print(io, ": ")
        end
        print(io, message)
        should_stop() && return
        for (k, v) in kwargs
            print(io, ' ', k, '=')
            print(io, v)
            should_stop() && return
        end
    end
end

# These are unused but just for making it a proper logger:
Logging.shouldlog(logger::DefaultLogger, level, _module, group, id) = true
Logging.min_enabled_level(::DefaultLogger) = Logging.Info


text:
import FrankWolfe
import LinearAlgebra
import Random

using Test

n = Int(1e3)
k = 10000

s = rand(1:100)
@info "Seed: $s"

Random.seed!(s)

xpi = rand(n);
total = sum(xpi);
const xp = xpi ./ total;

f(x) = LinearAlgebra.norm(x - xp)^2
function grad!(storage, x)
    @. storage = 2 * (x - xp)
    return nothing
end

lmo = FrankWolfe.UnitSimplexOracle(1.0)
x0 = FrankWolfe.compute_extreme_point(lmo, rand(n))

@time xblas, _ = FrankWolfe.frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Shortstep(),
    L=2,
    emphasis=FrankWolfe.blas,
    verbose=false,
    trajectory=false,
    momentum=0.9,
)

@time xmem, _ = FrankWolfe.frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Shortstep(),
    L=2,
    emphasis=FrankWolfe.memory,
    verbose=true,
    trajectory=true,
    momentum=0.9,
)

@test f(xblas) ≈ f(xmem)


text:
"""

General Framework for the data augmented Gaussian Processes

"""
module AugmentedGaussianProcesses

const AGP = AugmentedGaussianProcesses
export AGP
export AbstractGP, GP, VGP, SVGP, VStP, MCGP, MOVGP, MOSVGP, MOARGP, OnlineSVGP # All models
export AbstractLikelihood,
    RegressionLikelihood, ClassificationLikelihood, MultiClassLikelihood, EventLikelihood # All categories of likelihoods
export GaussianLikelihood, StudentTLikelihood, LaplaceLikelihood, HeteroscedasticLikelihood # Regression Likelihoods
export LogisticLikelihood, BayesianSVM # Classification Likelihoods
export SoftMaxLikelihood, LogisticSoftMaxLikelihood # Multiclass Classification Likelihoods
export PoissonLikelihood, NegBinomialLikelihood # Event Likelihoods
export AbstractInference, Analytic, AnalyticVI, AnalyticSVI # Inference objects
export GibbsSampling, HMCSampling # Sampling inference
export NumericalVI,
    NumericalSVI, MCIntegrationVI, MCIntegrationSVI, QuadratureVI, QuadratureSVI # Numerical inference
export PriorMean, ZeroMean, ConstantMean, EmpiricalMean, AffineMean # Prior means
#Useful functions
export train!, sample
export predict_f, predict_y, proba_y
export fstar
export ELBO
export covariance, diag_covariance, prior_mean
export @augmodel

#General modules
using Reexport
@reexport using KernelFunctions
@reexport using Flux.Optimise
@reexport using InducingPoints

using AbstractMCMC
using AdvancedHMC
using ChainRulesCore: ChainRulesCore, NoTangent
using Distributions:
    Distributions,
    Distribution,
    dim,
    cov,
    mean,
    var,
    pdf,
    logpdf,
    loglikelihood,
    Normal,
    Poisson,
    NegativeBinomial,
    InverseGamma,
    Laplace,
    MvNormal,
    Gamma
using FastGaussQuadrature: gausshermite
using Flux: params, destructure
using ForwardDiff
using KernelFunctions: ColVecs, RowVecs
using LinearAlgebra
using ProgressMeter
using Random
using StatsBase
using SimpleTraits
using StatsFuns
using SpecialFunctions
using Zygote

#Include custom module for additional distributions
include(joinpath("ComplementaryDistributions", "ComplementaryDistributions.jl"))
using .ComplementaryDistributions

# Main classes
abstract type AbstractInference{T<:Real} end
abstract type VariationalInference{T} <: AbstractInference{T} end
abstract type SamplingInference{T} <: AbstractInference{T} end
abstract type AbstractLikelihood{T<:Real} end
abstract type AbstractLatent{T<:Real,Tpr,Tpo} end

include(joinpath("mean", "priormean.jl"))
include(joinpath("data", "datacontainer.jl"))
include(joinpath("functions", "utils.jl"))

# Models
include(joinpath("models", "AbstractGP.jl"))
include(joinpath("gpblocks", "latentgp.jl"))
include(joinpath("models", "GP.jl"))
include(joinpath("models", "VGP.jl"))
include(joinpath("models", "MCGP.jl"))
include(joinpath("models", "SVGP.jl"))
include(joinpath("models", "VStP.jl"))
include(joinpath("models", "MOSVGP.jl"))
include(joinpath("models", "MOVGP.jl"))
include(joinpath("models", "OnlineSVGP.jl"))
include(joinpath("models", "single_output_utils.jl"))
include(joinpath("models", "multi_output_utils.jl"))

include(joinpath("inference", "inference.jl"))
include(joinpath("likelihood", "likelihood.jl"))
include(joinpath("likelihood", "generic_likelihood.jl"))

include(joinpath("functions", "KLdivergences.jl"))
include(joinpath("functions", "ELBO.jl"))
include(joinpath("data", "utils.jl"))
include(joinpath("functions", "plotting.jl"))

# Training and prediction functions
include(joinpath("training", "training.jl"))
include(joinpath("training", "sampling.jl"))
include(joinpath("training", "onlinetraining.jl"))
include(joinpath("hyperparameter", "autotuning.jl"))
include(joinpath("training", "predictions.jl"))
include("ar_predict.jl")

end #End Module


text:
#=
  Package: FishEBM.jl
  File: summary.jl
  Devin Rose
  Functions used for summarizing the simulation results upon completion
  Created: September 2016
=#

"""
  Description: Function used to generate summary files from an agent datavase for completed
    simulations for archiving and analysis. Saves all available data after a completed simulation.

  Returns: N/A

  Last update: September 2016
"""
function simSummary(adultAssumpt::AdultAssumptions, agentAssumpt::AgentAssumptions,
  agentDB::Vector{EnviroAgent}, effort::Vector{Float64}, finalWeek::Int64, initStock::Vector{Int64},
  carryingCap::Vector{Float64}, stageDataFrame::DataFrame, adultDataFrame::DataFrame,
  harvestDataFrame::DataFrame, harvestZoneData::DataFrame, spawnDataFrame::DataFrame,
  killedDataFrame::DataFrame, userInput::String)

  # Find directory for the results
  simDir()
  path = runDir(dateDir(resultsDir(setProjPath())[1])[1])[2]

  # Save all available data from the completed simulation

  # Weekly population (by stage)
  aliveData(stageDataFrame, path)

  # Weekly adult age (by year) distribution
  ageData(adultDataFrame, path)

  # Weekly harvest total by location (harvest zone number) and age
  harvestData(harvestDataFrame, harvestZoneData, path)

  # Weekly spawning total from each age
  spawnData(spawnDataFrame, path)

  # Wekly mortality totals from each layer of mortality
  killedData(killedDataFrame, path)

  # Create a readme file to save the simulation paramaters and description
  simReadme(adultAssumpt, agentAssumpt, effort, initStock, carryingCap, path, userInput)

  # Saving an empty file inside of the run for quickly checking the simulation check
  cd(path)
  if contains(userInput, "\n") && startswith(userInput, "\n") == false
    touch(split(userInput, "\n")[1])
  else
    touch(userInput)
  end

end