text:
import Base: @kwdef

# structs defining terms of a composable model for the independent-particle
# susceptibility χ0. The struct define a call operator, which does some setup
# and returns an `apply!(δρ, δV, α=1)` function to from `δρ .+= α * χ0 * δV`.
# `δV[1]` holds the spin-up (or total) potential change, `δV[2]` the
# spin-down potential as a `RealFourierArray`. `δρ[:, :, :, 1]` is the spin-up
# or the total density change and `δρ[:, :, :, 2]` the spin-down density change.

abstract type χ0Model end

@doc raw"""
Represents the LDOS-based ``χ_0`` model
```math
χ_0(r, r') = (-D_\text{loc}(r) δ(r, r') + D_\text{loc}(r) D_\text{loc}(r') / D)
```
where ``D_\text{loc}`` is the local density of states and ``D`` the density of states.
For details see Herbst, Levitt 2020 arXiv:2009.01665
"""
struct LdosModel <: χ0Model end
function (::LdosModel)(basis; eigenvalues, ψ, εF, kwargs...)
    n_spin = basis.model.n_spin_components
    dVol   = basis.model.unit_cell_volume / prod(basis.fft_size)

    # Catch cases that will yield no contribution
    iszero(basis.model.temperature) && return nothing
    ldos = [compute_ldos(εF, basis, eigenvalues, ψ, spins=[σ]) for σ in 1:n_spin]
    if maximum(maximum(abs, ldos[σ]) for σ in 1:n_spin) < eps(eltype(basis))
        return nothing
    end

    dos = sum(sum, ldos) * dVol  # Integrate LDOS to form total DOS
    function apply!(δρ, δV, α=1)
        dotldosδV = sum(dot(ldos[σ], δV[σ].real) for σ = 1:n_spin)
        for σ in 1:n_spin
            δρ[:, :, :, σ] .+= α .* (-ldos[σ] .* δV[σ].real
                                     .+ ldos[σ] .* dotldosδV .* dVol ./ dos)
        end
        δρ
    end
end

@doc raw"""
A localised dielectric model for ``χ_0``:
```math
\sqrt{L(x)} \text{IFFT} \frac{C_0 G^2}{4π (1 - C_0 G^2 / k_{TF}^2)} \text{FFT} \sqrt{L(x)}
```
where ``C_0 = 1 - ε_r``, `L(r)` is a real-space localisation function
and otherwise the same conventions are used as in [`DielectricMixing`](@ref).
"""
@kwdef struct DielectricModel <: χ0Model
    εr::Real  = 10
    kTF::Real = 0.8
    localisation::Function = identity
end
function (χ0::DielectricModel)(basis; kwargs...)
    T   = eltype(basis)
    εr  = T(χ0.εr)
    kTF = T(χ0.kTF)
    C0  = 1 - εr
    n_spin = basis.model.n_spin_components
    iszero(C0) && return nothing  # Will yield no contribution

    Gsq = [sum(abs2, G) for G in G_vectors_cart(basis)]
    apply_sqrtL = identity
    if χ0.localisation != identity
        sqrtL = sqrt.(χ0.localisation.(r_vectors(basis)))
        apply_sqrtL = x -> from_real(basis, sqrtL .* x.real)
    end

    function apply!(δρ, δV, α=1)
        for σ in 1:n_spin
            loc_δV = apply_sqrtL(δV[σ]).fourier
            dielectric_loc_δV = @. C0 * kTF^2 * Gsq / 4T(π) / (kTF^2 - C0 * Gsq) * loc_δV
            δρ[:, :, :, σ] .+= α * apply_sqrtL(from_fourier(basis, dielectric_loc_δV)).real
        end
        δρ
    end
end

"""
Full χ0 application, optionally dropping terms or disabling Sternheimer.
All keyword arguments passed to [`apply_χ0`](@ref).
"""
struct Applyχ0Model <: χ0Model
    kwargs_apply_χ0
end
Applyχ0Model(; kwargs_apply_χ0...) = Applyχ0Model(kwargs_apply_χ0)
function (χ0::Applyχ0Model)(basis; ham, eigenvalues, ψ, εF, n_ep_extra, kwargs...)
    # self_consistent_field uses a few extra bands, which are not converged by the eigensolver
    # For the χ0 application, bands need to be perfectly converged, so we discard them here
    ψ_cvg = [@view ψk[:, 1:end-n_ep_extra]  for ψk in ψ]
    eigenvalues_cvg = [εk[1:end-n_ep_extra] for εk in eigenvalues]

    function apply!(δρ, δV, α=1)
        # χ0δV[1] is total, χ0δV[2] is spin
        χ0δV = apply_χ0(ham, ψ_cvg, εF, eigenvalues_cvg, δV...; χ0.kwargs_apply_χ0...)
        if basis.model.n_spin_components == 1
            δρ[:, :, :, 1] .+= α .* χ0δV[1].real
        else
            δρ[:, :, :, 1] .+= α .* (χ0δV[1].real .+ χ0δV[2].real) ./ 2
            δρ[:, :, :, 2] .+= α .* (χ0δV[1].real .- χ0δV[2].real) ./ 2
        end
        δρ
    end
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export grid

"""
    grid(;kwargs...)
    grid(children::Any;kwargs...)
    grid(children_maker::Function;kwargs...)

A Grid component.
Carbon Grid
Keyword arguments:
- `children` (Array of a list of or a singular dash component, string or numbers | a list of or a singular dash component, string or number; optional): The children of the element
- `style` (Dict; optional): The inline styles
- `id` (String; optional): The ID of this component, used to identify dash components
in callbacks. The ID needs to be unique across all of the
components in an app.
- `className` (String; optional): The class of the element
"""
function grid(; kwargs...)
        available_props = Symbol[:children, :style, :id, :className]
        wild_props = Symbol[]
        return Component("grid", "Grid", "dash_carbon_components", available_props, wild_props; kwargs...)
end

grid(children::Any; kwargs...) = grid(;kwargs..., children = children)
grid(children_maker::Function; kwargs...) = grid(children_maker(); kwargs...)

ext: jl

max_stars_repo_name: Matheus-Rangel/dash-carbon-components

max_issues_repo_name: Matheus-Rangel/dash-carbon-components

max_forks_repo_name: Matheus-Rangel/dash-carbon-components


text:
import Serialization

mutable struct RunLogger <: Logger
    trial_id
    logged_organisms_alive
    logged_organisms_dead

    RunLogger(trial_id) = new(string(trial_id), Dict(), Dict())
end

function log_step(logger::RunLogger, model)
    time = get_time(model)
    if time % LOG_FREQUENCY != 0 return end

    snapshot_id = string(time_ns())

    current_organisms = get_organisms(model)

    for organism in current_organisms
        id = get_id(model, organism)

        if haskey(logger.logged_organisms_alive, id)
            logger.logged_organisms_alive[id].snapshot_id = snapshot_id
        end
    end

    save_snapshot(model, snapshot_id)
    save_log(logger, snapshot_id)
end

function log_birth(logger::RunLogger, model, child, parents = nothing)
    time = get_time(model)

    if parents === nothing return end

    parent_ids = [get_id(model, parent) for parent in parents]
    parent_genotypes = [get_genotype(model, parent) for parent in parents]

    abstracted_organism = get_abstracted_organism(model, child, parent_ids, parent_genotypes)
    abstracted_organism.time_birth = time

    logger.logged_organisms_alive[abstracted_organism.id] = abstracted_organism
end

function log_death(logger::RunLogger, model, organism)
    time = get_time(model)
    organism_id = get_id(model, organism)

    if haskey(logger.logged_organisms_alive, organism_id) == false return end

    logged_organism = logger.logged_organisms_alive[organism_id]
    logged_organism.time_death = time

    delete!(logger.logged_organisms_alive, organism_id)
    logger.logged_organisms_dead[organism_id] = logged_organism
end

function save_snapshot(model, snapshot_id)
    logger = get_logger(model)
    path = SNAPSHOTS_FOLDER * logger.trial_id * "_" * string(snapshot_id)

    set_logger!(model, DoNothingLogger())
    Serialization.serialize(path, model)
    set_logger!(model, logger)
end

function save_log(logger::RunLogger, snapshot_id)
    path = LOGGER_FOLDER * logger.trial_id * "_" * snapshot_id
    Serialization.serialize(path, logger.logged_organisms_dead)
    logger.logged_organisms_dead = Dict()
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
module Plot2D

using GLMakie
export plot2d, buttons

LOG = nothing
const P1 = [Node(Vector{Point2f0}(undef, 6000))]
const P2 = [Node(Vector{Point2f0}(undef, 6000))]

function autoscale(ax, x, y)
    xlims!(ax, x[1], x[end])
    range = maximum(y) - minimum(y)
    if range > 0.0f0
        ylims!(ax, minimum(y)-0.05*range, maximum(y)+0.05*range)
    end
end

function plot2d(se, ax, label1, log, points, field, lower=false)
    global LOG, P1, P2
    LOG=log
    if lower
        P2[1]=points
    else
        P1[1]=points
    end
    unit = ""
    y = [1f0]
    if field == :height
        unit = "[m]"
        y    = log.extlog.z ./ se().zoom
    elseif field == :elevation
        unit = "[°]"
        y    = log.syslog.elevation/pi*180.0
    elseif field == :azimuth
        unit = "[°]"
        y    = log.syslog.azimuth/pi*180.0
    elseif field == :v_reelout
        unit = "[m/s]"
        y    = log.syslog.v_reelout
    elseif field == :force
        unit = "[N]"
        y    = log.syslog.force
    elseif field == :depower
        unit = "[%]"
        y    = log.syslog.depower * 100.0f0
    elseif field == :v_app
        unit = "[m/s]"
        y    = log.syslog.v_app
    elseif field == :l_tether
        unit = "[m]"
        y    = log.syslog.l_tether
    elseif field == :power
        unit = "[kW]"
        y    = log.syslog.v_reelout .* log.syslog.force * 0.001f0
    end
    x       = log.syslog.time
    if field == :power
        label1[] = "mechanical " * string(field) * " " * unit
    else
        label1[] = string(field) * " " * unit
    end
    points[]    =  Point2f0.(x, y)
    autoscale(ax, x, y)
end

function buttons(fig, bg, se, ax1, ax2, label1, label2, reset)
    textsize=14
    btn_height         = Button(fig, label = "height", textsize=textsize)
    btn_elevation      = Button(fig, label = "elevation", textsize=textsize)
    btn_azimuth        = Button(fig, label = "azimuth", textsize=textsize)
    btn_v_reelout      = Button(fig, label = "v_reelout", textsize=textsize)
    btn_force          = Button(fig, label = "force", textsize=textsize)
    btn_depower        = Button(fig, label = "depower", textsize=textsize)
    btn_v_app          = Button(fig, label = "v_app", textsize=textsize)
    btn_l_tether       = Button(fig, label = "l_tether", textsize=textsize)
    btn_power          = Button(fig, label = "power", textsize=textsize)
    sw = Toggle(fig, active = false) # active means plot in the upper area
    bg[1, 1:10] = [btn_height, btn_elevation, btn_azimuth, btn_v_reelout, btn_force, btn_depower, btn_v_app, btn_l_tether, btn_power, sw]

    on(btn_height.clicks) do c
        if sw.active[]
            plot2d(se, ax1, label1, LOG, P1[1], :height)
        else
            plot2d(se, ax2, label2, LOG, P2[1], :height, true)
        end
        reset()
    end
    on(btn_elevation.clicks) do c
        if sw.active[]
            plot2d(se, ax1, label1, LOG, P1[1], :elevation)
        else
            plot2d(se, ax2, label2, LOG, P2[1], :elevation, true)
        end
        reset()
    end
    on(btn_azimuth.clicks) do c
        if sw.active[]
            plot2d(se, ax1, label1, LOG, P1[1], :azimuth)
        else
            plot2d(se, ax2, label2, LOG, P2[1], :azimuth, true)
        end
        sleep(0.05)
        reset()
    end
    on(btn_v_reelout.clicks) do c
        if sw.active[]
            plot2d(se, ax1, label1, LOG, P1[1], :v_reelout)
        else
            plot2d(se, ax2, label2, LOG, P2[1], :v_reelout, true)
        end
        reset()
    end
    on(btn_force.clicks) do c
        if sw.active[]
            plot2d(se, ax1, label1, LOG, P1[1], :force)
        else
            plot2d(se, ax2, label2, LOG, P2[1], :force, true)
        end
        reset()
    end
    on(btn_depower.clicks) do c
        if sw.active[]
            plot2d(se, ax1, label1, LOG, P1[1], :depower)
        else
            plot2d(se, ax2, label2, LOG, P2[1], :depower, true)
        end
        reset()
    end
    on(btn_v_app.clicks) do c
        if sw.active[]
            plot2d(se, ax1, label1, LOG, P1[1], :v_app)
        else
            plot2d(se, ax2, label2, LOG, P2[1], :v_app, true)
        end
        reset()
    end
    on(btn_l_tether.clicks) do c
        if sw.active[]
            plot2d(se, ax1, label1, LOG, P1[1], :l_tether)
        else
            plot2d(se, ax2, label2, LOG, P2[1], :l_tether, true)
        end
        reset()
    end
    on(btn_power.clicks) do c
        if sw.active[]
            plot2d(se, ax1, label1, LOG, P1[1], :power)
        else
            plot2d(se, ax2, label2, LOG, P2[1], :power)
        end
        reset()
    end
end

end


text:
# Iterated extended Kalman smoothing
mutable struct IEKS <: AbstractEK
    prior::Symbol
    order::Int
    diffusionmodel::Symbol
    smooth::Bool
    linearize_at
end

"""
    IEKS(; prior=:ibm, order=1, diffusionmodel=:dynamic, linearize_at=nothing)

Gaussian ODE filtering with iterated extended Kalman smoothing. To use it, use
`solve_ieks(prob, IEKS(), args...)`
instead of
`solve(prob, IEKS(), args...)`,
since it is implemented as an outer loop around the solver.

Currently, only the integrated Brownian motion prior `:ibm` is supported.
For the diffusionmodel, chose one of
`[:dynamic, :dynamicMV, :fixed, :fixedMV, :fixedMAP]`.
Just like the [`EK1`](@ref) it requires that the Jacobian of the rhs function is available.

See also: [`EK0`](@ref), [`EK1`](@ref), [`solve_ieks`](@ref)

# References:
- F. Tronarp, S. Särkkä, and P. Hennig: **Bayesian ODE Solvers: The Maximum A Posteriori Estimate**
"""
function IEKS(; prior=:ibm, order=1, diffusionmodel=:dynamic, linearize_at=nothing)
    if !isnothing(linearize_at)
        @assert linearize_at isa ProbODESolution
        @assert linearize_at.alg.prior == prior
        @assert linearize_at.alg.order == order
        @assert linearize_at.alg.diffusionmodel == diffusionmodel
        @assert linearize_at.alg.smooth == true
    end
    return IEKS(prior, order, diffusionmodel, true, linearize_at)
end

"""
    solve_ieks(prob::DiffEqBase.AbstractODEProblem, alg::IEKS, args...; iterations=10, kwargs...)

Solve method to be used with the [`IEKS`](@ref). The IEKS works essentially by solving the
ODE multiple times. `solve_ieks` therefore wraps a call to the standard `solve` method,
passing `args...` and `kwargs...`.
"""
function solve_ieks(prob::DiffEqBase.AbstractODEProblem, alg::IEKS, args...;
                    iterations=10, kwargs...)
    sol = nothing
    for i in 1:iterations
        alg.linearize_at = sol
        sol = solve(prob, alg, args...; kwargs...)
    end
    return sol
end


text:
########################################
## File Name: state_transition_test.jl
## Author: Haruki Nishimura (hnishimura@stanford.edu)
## Date Created: 2020/05/12
## Description: Test code for transition functions in
## src/state_transition_observation_models.jl
########################################

using LinearAlgebra
using Random

@testset "State Transition Test" begin

# TransModel_Pos
transPos = TransModel_Pos();
xPos = PhysPos(0.,zeros(2));
uPos = PosControl(0.,[2.,3.2]);
p = trans(transPos,xPos,uPos,0.01);
@test p.t == 0.01
@test p.pos == [0.02, 0.032]
@test p.dim == 2

p2 = trans(transPos, xPos, uPos, 0.01, Matrix(1.0I, 2, 2), MersenneTwister(1234));

J = trans_jacobi(transPos,xPos,uPos);
@test J == zeros(2, 2);

J2 = trans_jacobi_euler(transPos,xPos,uPos,0.01)
@test J2 == Matrix(1.0I, 2, 2)

C = trans_u_coeff(transPos, xPos)
@test C == Matrix(1.0I, 2, 2)

# TransModel_Manipulate2D
transMP = TransModel_Manipulate2D();
xMP = PhysManipulate2D(0.,[0.1,0.2],0.5,[0.3,0.4],0.6,3.,5.,[10.,11.],0.7);
uMP = MControl2D(0.,6.,7.,8.);
m = trans(transMP,xMP,uMP,0.01)
@test m.t == 0.01
@test all(isapprox.(m.pos, [0.103, 0.204]))
@test isapprox(m.θ, 0.506)
@test all(isapprox.(m.vel, [0.3193, 0.4224]))
@test isapprox(m.ω, 0.4916580629)
@test m.m == 3.0
@test m.J == 5.0
@test m.r == [10., 11.]
@test m.μ == 0.7
@test m.dim == 2

m2 = trans(transMP, xMP, uMP, 0.01, Matrix(1.0I, 11, 11), MersenneTwister(1234))

J3 = trans_jacobi(transMP, xMP, uMP)
@test all(isapprox.(trans_jacobi_auto(transMP, xMP, uMP), J3))

J4 = trans_jacobi_euler(transMP, xMP, uMP, 0.01)
@test all(isapprox.(J4, J3.*0.01 + Matrix(1.0I, 11, 11)))

C2 = trans_u_coeff(transMP, xMP)
@test all(isapprox.(C2, [0.0       0.0       0.0
                         0.0       0.0       0.0
                         0.0       0.0       0.0
                         0.333333  0.0       0.0
                         0.0       0.333333  0.0
                        -2.88953   0.700429  0.2
                         0.0       0.0       0.0
                         0.0       0.0       0.0
                         0.0       0.0       0.0
                         0.0       0.0       0.0
                         0.0       0.0       0.0], atol=1e-5))

end


text:
using ArrayAllez
using Test

@static if Sys.isapple()
    using AppleAccelerate
    @info "testing with AppleAccelerate (always done on Apple machines)"
elseif Sys.isunix()
    using IntelVectorMath
    @info "testing with IntelVectorMath (always done on Linux machines)"
else
    @info "testing without AppleAccelerate nor IntelVectorMath, to check fallbacks"
end

@testset "exp/log/inv/scale" begin
    @testset "small" begin

        m = rand(3,7)
        v = randn(3)
        r = randn(7)'

        @test exp!(copy(m)) ≈ exp_(m) ≈ exp0(m) ≈ exp_(:test, m)
        @test exp!(copy(v)) ≈ exp_(v) ≈ exp0(v) ≈ exp_(:test, v)
        @test exp!(copy(r)) ≈ exp_(r) ≈ exp0(r) ≈ exp_(:test, r)

        @test exp!(copy(m)) ≈ exp_(m) ≈ exp0(m) ≈ exp_(:test, m)
        @test exp!(copy(v)) ≈ exp_(v) ≈ exp0(v) ≈ exp_(:test, v)
        @test exp!(copy(r)) ≈ exp_(r) ≈ exp0(r) ≈ exp_(:test, r)

        @test inv!(copy(m)) ≈ inv_(m) ≈ inv0(m) ≈ inv_(:test, m)
        @test inv!(copy(v)) ≈ inv_(v) ≈ inv0(v) ≈ inv_(:test, v)
        @test inv!(copy(r)) ≈ inv_(r) ≈ inv0(r) ≈ inv_(:test, r)

        @test scale!(copy(m),π) ≈ scale_(m,π) ≈ scale0(m,π)
        @test scale!(copy(v),π) ≈ scale_(v,π) ≈ scale0(v,π)
        @test scale!(copy(r),π) ≈ scale_(r,π) ≈ scale0(r,π)

        @test scale!(copy(m),v) ≈ scale_(m,v) ≈ scale0(m,v)
        @test scale!(copy(v),v) ≈ scale_(v,v) ≈ scale0(v,v)

        @test scale!(copy(m),r) ≈ scale_(m,r) ≈ scale0(m,r)
        @test scale!(copy(r),r) ≈ scale_(r,r) ≈ scale0(r,r)

        @test iscale!(copy(m),π) ≈ iscale_(m,π) ≈ iscale0(m,π)
        @test iscale!(copy(v),π) ≈ iscale_(v,π) ≈ iscale0(v,π)
        @test iscale!(copy(r),π) ≈ iscale_(r,π) ≈ iscale0(r,π)

        @test iscale!(copy(m),v) ≈ iscale_(m,v) ≈ iscale0(m,v)
        @test iscale!(copy(v),v) ≈ iscale_(v,v) ≈ iscale0(v,v)

        @test iscale!(copy(m),r) ≈ iscale_(m,r) ≈ iscale0(m,r)
        @test iscale!(copy(r),r) ≈ iscale_(r,r) ≈ iscale0(r,r)

    end
    @testset "large" begin # needed because some functions switch on threading

        m = rand(300,700);
        v = randn(300);
        r = randn(700)';

        @test exp!(copy(m)) ≈ exp_(m) ≈ exp0(m) ≈ exp_(:test, m)
        @test exp!(copy(v)) ≈ exp_(v) ≈ exp0(v) ≈ exp_(:test, v)
        @test exp!(copy(r)) ≈ exp_(r) ≈ exp0(r) ≈ exp_(:test, r)

        @test exp!(copy(m)) ≈ exp_(m) ≈ exp0(m) ≈ exp_(:test, m)
        @test exp!(copy(v)) ≈ exp_(v) ≈ exp0(v) ≈ exp_(:test, v)
        @test exp!(copy(r)) ≈ exp_(r) ≈ exp0(r) ≈ exp_(:test, r)

        @test inv!(copy(m)) ≈ inv_(m) ≈ inv0(m) ≈ inv_(:test, m)
        @test inv!(copy(v)) ≈ inv_(v) ≈ inv0(v) ≈ inv_(:test, v)
        @test inv!(copy(r)) ≈ inv_(r) ≈ inv0(r) ≈ inv_(:test, r)

    end
end

@testset "dropdims" begin

    @dropdims begin
        a = sum(ones(3,7), dims=2)
        b = sum(10 .* randn(2,10); dims=2) do x
            trunc(Int, x)
        end
    end
    @test a isa Vector
    @test b isa Vector

end

@info "loading Tracker"
using Tracker
using Tracker: TrackedArray, gradcheck, back!, data, grad

gradtest(f, dims...) = gradtest(f, rand.(Float64, dims)...) ## from Flux tests
gradtest(f, xs::AbstractArray...) = gradcheck((xs...) -> sum(sin.(f(xs...))), xs...)

using ForwardDiff

mycheck(f, x) = ForwardDiff.gradient(z -> sum(sin,f(z)), x) ≈ Tracker.gradient(z -> sum(sin,f(z)), x)[1]

@testset "Tracker gradients" begin
    @testset "exp + log" begin

        @test gradtest(exp0, (2,3))
        @test gradtest(sum∘exp_, (2,3))
        @test gradtest(sum∘exp!∘copy, (2,3))

        p = param(randn(2,3));
        back!(sum(exp.(p)))
        pg = p.grad
        p.grad[:] .= 0;
        back!(sum(exp!(p)))
        @test p.grad ≈ pg

        p.grad[:] .= 0;
        back!(sum(exp!!(p)))
        @test p.grad ≈ pg

        @test gradtest(log0, rand(2,3))
        @test gradtest(log_, rand(2,3))
        @test gradtest(log!∘copy, rand(2,3))

        p = param(rand(2,3));
        back!(sum(log.(p)))
        pg = p.grad
        p.grad[:] .= 0;
        back!(sum(log!(p)))
        @test p.grad ≈ pg

        # p.grad[:] .= 0;
        # back!(sum(log!!(p)))
        # @test p.grad ≈ pg

        @test gradcheck(A -> scale0(A,4) |> sum, rand(2,3))
        @test gradcheck(A -> scale_(A,4) |> sum, rand(2,3))

    end
    @testset "exp + log II" begin # using ForwardDiff, no problem to test exp! etc.

        m = rand(3,7)

        mycheck(z -> log0(z), m)
        mycheck(z -> log_(z), m)
        mycheck(z -> log!(z), m)
        mycheck(z -> log!!(z), m)

        mycheck(z -> exp0(z), m)
        mycheck(z -> exp_(z), m)
        mycheck(z -> exp!(z), m)
        mycheck(z -> exp!!(z), m)

    end
    @testset "scale + inv" begin

        m = rand(3,7)
        v = randn(3)
        r = randn(7)'

        @test gradtest(z -> scale_(z,9), m)
        @test gradtest(z -> scale_(z,v), m)
        @test gradtest(z -> scale_(z,r), m)

        # @test gradtest(z -> iscale_(z,9), m)
        # @test gradtest(z -> iscale_(z,v), m)
        # @test gradtest(z -> iscale_(z,r), m)

        # @test gradcheck(z -> sum(inv_(z)), m)
        # @test gradcheck(z -> sum(inv_(z,9)), m)

        # @test gradcheck(z -> sum(scale_(m,z)), v) # crash?
        # @test gradcheck(z -> sum(scale_(m,z)), r)
        #
        # @test gradcheck(z -> sum(iscale_(m,z)), v)
        # @test gradcheck(z -> sum(iscale_(m,z)), r)

    end
    @testset "scale + inv II" begin

        m = rand(3,7)
        v = randn(3)
        r = randn(7)'

        @test mycheck(z -> scale0(z,9), m)
        @test mycheck(z -> scale_(z,9), m)
        @test mycheck(z -> scale!(z,9), m)

        @test mycheck(z -> scale0(z,v), m)
        @test mycheck(z -> scale_(z,v), m)
        @test mycheck(z -> scale!(z,v), m)

        @test mycheck(z -> scale0(z,r), m)
        @test mycheck(z -> scale_(z,r), m)
        # @test mycheck(z -> scale!(z,r), m) # ambiguous

        @test mycheck(z -> scale0(z,m), m)
        # @test mycheck(z -> scale_(z,m), m) # no method
        # @test mycheck(z -> scale!(z,m), m)

    end
    @testset "prod + cumprod" begin # https://github.com/FluxML/Flux.jl/pull/524

        @test gradtest(x -> prod(x, dims=(2, 3)), (3,4,5))
        @test gradtest(x -> prod(x, dims=1), (3,4,5))
        @test gradtest(x -> prod(x, dims=1), (3,))
        @test gradtest(x -> prod(x), (3,4,5))
        @test gradtest(x -> prod(x), (3,))

        rzero(dims...) = (x = rand(dims...); x[2]=0; x)
        @test gradtest(x -> prod(x, dims=(2, 3)), rzero(3,4,5))
        @test gradtest(x -> prod(x, dims=1), rzero(3,4,5))
        @test gradtest(x -> prod(x, dims=1), rzero(3,))
        @test gradtest(x -> prod(x), rzero(3,4,5))
        @test gradtest(x -> prod(x), rzero(3,))

        @test gradtest(x -> cumsum(x, dims=2), (3,4,5))
        @test gradtest(x -> cumsum(x, dims=1), (3,))
        @test gradtest(x -> cumsum(x), (3,))

        @test gradtest(x -> cumprod(x, dims=2), (3,4,5))
        @test gradtest(x -> cumprod(x, dims=1), (3,))
        @test gradtest(x -> cumprod(x), (3,))
        @test gradtest(x -> cumprod(x, dims=2), rzero(3,4,5))
        @test gradtest(x -> cumprod(x, dims=1), rzero(3,))
        @test gradtest(x -> cumprod(x), rzero(3,))

    end
end

#=
@info "loading Zygote"
using Zygote: Zygote

@testset "Zygote gradients" begin
    @testset "* left & right" begin

        @test Zygote.gradient(*ˡ, 2,3) == (3, nothing)
        @test Zygote.gradient(sum∘*ˡ, rand(2,2), ones(2,2)) == ([2 2; 2 2], nothing)

        @test Zygote.gradient(*ʳ, 2,3) == (nothing, 2)
        @test Zygote.gradient(sum∘*ʳ, ones(2,2), rand(2,2)) == (nothing, [2 2; 2 2])

    end
    @testset "odot" begin

        @test Zygote.gradient(⊙ˡ, 2,3) == (3, nothing)
        @test Zygote.gradient(sum∘⊙ˡ, rand(2,2), ones(2,2,2,2)) == ([8 8; 8 8], nothing)

        @test Zygote.gradient(⊙ʳ, 2,3) == (nothing, 2)
        @test Zygote.gradient(sum∘⊙ʳ, ones(2,2,2), rand(2,2)) == (nothing, [4 4; 4 4])

    end
end
=#


text:
@testset "Witness Sets" begin

    @testset "affine" begin
        @var x y

        F = System([x^2 + y^2 - 5], [x, y])

        W = witness_set(F)

        @test dim(W) == 1
        @test codim(W) == 1
        @test degree(W) == 2
        @test solutions(W) isa Vector{Vector{ComplexF64}}
        @test results(W) isa Vector{PathResult}

        L = LinearSubspace([1 1], [-1])

        W_L = witness_set(W, L)
        @test degree(W_L) == 2
        @test sort(real.(solutions(W_L))) ≈ [[-2, 1], [1, -2]]
        @test linear_subspace(W_L) == L

        @test trace_test(W) < 1e-8
    end

    @testset "projective" begin
        @var x y z

        F = System([x^2 + y^2 - 5z^2], [x, y, z])

        W = witness_set(F; compile = false)

        @test dim(W) == 1
        @test codim(W) == 2
        @test degree(W) == 2
        @test solutions(W) isa Vector{Vector{ComplexF64}}
        @test results(W) isa Vector{PathResult}
        @test trace_test(W) < 1e-8

        L = LinearSubspace([1 1 1])
        W_L = witness_set(W, L; compile = false)
        @test degree(W_L) == 2

        L = rand_subspace(3; codim = 1, affine = false)
        W_L = witness_set(W, L; compile = false)
        @test degree(W_L) == 2

        L = rand_subspace([x, y, z]; codim = 1, affine = false)
        W_L = witness_set(W, L; compile = false)
        @test degree(W_L) == 2

        L = rand_subspace(3; codim = 1, affine = true)
        @test_throws ErrorException witness_set(W, L; compile = false)
    end

    @testset "dim / codim" begin
        @var x[1:6]
        homogeneous = true

        f = System([
            rand_poly(x, 2; homogeneous = homogeneous),
            rand_poly(x, 2; homogeneous = homogeneous),
            rand_poly(x, 2; homogeneous = homogeneous),
            rand_poly(x, 2; homogeneous = homogeneous),
        ])

        # quadratic_surface = rand_poly([x,y,z,w], 2; homogeneous = true)
        # sextic_curve = System([cubic_surface, quadratic_surface])
        @test degree(witness_set(f; dim = 1, compile = false)) == 16
        @test degree(witness_set(f; codim = 4, compile = false)) == 16
        @test degree(witness_set(f; compile = false)) == 16

        homogeneous = false
        f = System([
            rand_poly(x, 2; homogeneous = homogeneous),
            rand_poly(x, 2; homogeneous = homogeneous),
            rand_poly(x, 2; homogeneous = homogeneous),
            rand_poly(x, 2; homogeneous = homogeneous),
        ])
        # quadratic_surface = rand_poly([x,y,z,w], 2; homogeneous = true)
        # sextic_curve = System([cubic_surface, quadratic_surface])
        @test degree(witness_set(f; dim = 2, compile = false)) == 16
        @test degree(witness_set(f; codim = 4, compile = false)) == 16
        @test degree(witness_set(f; compile = false)) == 16
    end
end


text:
function printanalysis0(paramvals,logposteriors)
    println("Parameter Values: $(round(paramvals,4))")
    println("Mean (STD / %) of LogPosteriors: $(mean(logposteriors)) ($(round(std(logposteriors),2)) / $(abs(round(std(logposteriors)/mean(logposteriors)*100,2)))%)")
end

function printanalysis1(functionname,paramvals,logposteriors)
    println("====================================================")
    println("Results of $functionname for $(length(logposteriors)) iterations")
    println("====================================================")
    printanalysis0(paramvals,logposteriors)
end

function printanalysis2(functionname,paramvals,logposteriors)
    println("====================================================")
    println("Results of $functionname for $(size(logposteriors,1)) iterations and $(size(logposteriors,2)) parameter sets")
    println("====================================================")
    for i=1:size(paramvals,2)
        printanalysis0(paramvals[:,i],logposteriors[:,i])
    end
end

function currentplot(model::GMHPhotoReceptor.PhotoReceptorModel,current::Array)
    PyPlot.plot(dataindex(model),current)
    PyPlot.plot(dataindex(model),vec(mean(current,2));linewidth=1,color="black",label="mean")
    PyPlot.xlim(dataindex(model)[1],dataindex(model)[end])
    PyPlot.legend(loc="upper right",fancybox="true")
    PyPlot.grid("on")
    PyPlot.xlabel("Time (s)")
    PyPlot.ylabel("Current (nA)")
    PyPlot.title("PhotoReceptorModel Light-Induced Current")
end

function meancurrentplot(model::GMHPhotoReceptor.PhotoReceptorModel,meancurrent::Vector)
    PyPlot.plot(dataindex(model),measurements(model);linewidth=2,color="yellow",label="measured")
    PyPlot.plot(dataindex(model),meancurrent;linewidth=1,color="black",label="mean")
    PyPlot.legend(loc="upper right",fancybox="true")
    PyPlot.grid("on")
    PyPlot.xlabel("Time (s)")
    PyPlot.ylabel("Current (nA)")
    PyPlot.title("Measured vs. Mean PhotoReceptorModel Light-Induced Current")
end

function meancurrentsplot(model::GMHPhotoReceptor.PhotoReceptorModel,meancurrents::Array)
    PyPlot.plot(dataindex(model),meancurrents)
    PyPlot.xlim(dataindex(model)[1],dataindex(model)[end])
    PyPlot.grid("on")
    PyPlot.xlabel("Time (s)")
    PyPlot.ylabel("Current (nA)")
    PyPlot.title("Mean PhotoReceptorModel Light-Induced Currents")
end

function stdcurrentplot(model::GMHPhotoReceptor.PhotoReceptorModel,stdcurrent::Array)
    PyPlot.plot(dataindex(model),stdcurrent)
    PyPlot.xlim(dataindex(model)[1],dataindex(model)[end])
    PyPlot.grid("on")
    PyPlot.xlabel("Time (s)")
    PyPlot.ylabel("Current (nA)")
    PyPlot.title("STD of PhotoReceptorModel Light-Induced Current")
end

function logposteriorplot(model::GMHPhotoReceptor.PhotoReceptorModel,logposteriors::Vector)
    nbins = round(Int,sqrt(length(logposteriors)))
    e,c = hist(logposteriors,nbins)
    f = Distributions.fit(Distributions.Normal,logposteriors)
    PyPlot.plt[:hist](logposteriors,nbins)
    PyPlot.grid("on")
    PyPlot.plot(e,Distributions.pdf(f,e)*length(logposteriors)*(maximum(e)-minimum(e))/nbins)
    PyPlot.title("Histogram of Log-Posterior Values")
    f
end

function meanlogposteriorplot(model::GMHPhotoReceptor.PhotoReceptorModel,logposteriors::Array)
    meanlogpost,varlogpost = streamingstats(logposteriors')
    PyPlot.plot(1:size(meanlogpost,2),meanlogpost';linewidth=1)
    PyPlot.grid("on")
    PyPlot.xlabel("Iteration")
    PyPlot.ylabel("Log-Posterior")
    PyPlot.title("Mean of Log-Posterior Distribution")
end

function analysis1(model::GMHPhotoReceptor.PhotoReceptorModel,niterations::Integer;plotmode=:repl)
    #check the value of plotmode by converting it into a trait
    plotmodetrait = GeneralizedMetropolisHastings.trait(:plotmode,plotmode)
    #perform the required analysis
    paramvals,macrocurrents,logposteriors = evaluate1(model,niterations)
    meancurrent = vec(mean(macrocurrents,2))
    stdcurrent = vec(std(macrocurrents,2))
    #print the results
    printanalysis1("analysis1",paramvals,logposteriors)
    #plot in suitable manner for repl or notebook
    if traitvalue(plotmodetrait) == :repl
        PyPlot.figure()
        PyPlot.subplot(221)
    elseif traitvalue(plotmodetrait) == :notebook
        PyPlot.figure()
        PyPlot.subplot(211)
    end
    currentplot(model,macrocurrents)
    if traitvalue(plotmodetrait) == :repl
        PyPlot.subplot(223)
    elseif traitvalue(plotmodetrait) == :notebook
        PyPlot.figure()
        PyPlot.subplot(211)
    end
    meancurrentplot(model,meancurrent)
    if traitvalue(plotmodetrait) == :repl
        PyPlot.subplot(222)
    elseif traitvalue(plotmodetrait) == :notebook
        PyPlot.figure()
        PyPlot.subplot(211)
    end
    stdcurrentplot(model,stdcurrent)
    if traitvalue(plotmodetrait) == :repl
        PyPlot.subplot(224)
    elseif traitvalue(plotmodetrait) == :notebook
        PyPlot.figure()
    end
    f = logposteriorplot(model,logposteriors)
    #return the results
    paramvals,macrocurrents,logposteriors,f
end

function analysis2(model::GMHPhotoReceptor.PhotoReceptorModel,niterations::Integer;plotmode=:repl)
    #check the value of plotmode by converting it into a trait
    plotmodetrait = GeneralizedMetropolisHastings.trait(:plotmode,plotmode)
    #perform the required analysis
    paramvals,macrocurrentstats,logposteriors = evaluate2(model,niterations)
    meancurrent = mean(macrocurrentstats)
    stdcurrent = std(macrocurrentstats)
    #print the results
    printanalysis1("analysis2",paramvals,logposteriors)
    #plot in suitable manner for repl or notebook
    if traitvalue(plotmodetrait) == :repl
        PyPlot.figure()
        PyPlot.subplot(221)
    elseif traitvalue(plotmodetrait) == :notebook
        PyPlot.figure()
        PyPlot.subplot(211)
    end
    meancurrentplot(model,meancurrent)
    if traitvalue(plotmodetrait) == :repl
        PyPlot.subplot(222)
    elseif traitvalue(plotmodetrait) == :notebook
        PyPlot.figure()
        PyPlot.subplot(211)
    end
    stdcurrentplot(model,stdcurrent)
    if traitvalue(plotmodetrait) == :repl
        PyPlot.subplot(223)
    elseif traitvalue(plotmodetrait) == :notebook
        PyPlot.figure()
        PyPlot.subplot(121)
    end
    f = logposteriorplot(model,logposteriors)
    if traitvalue(plotmodetrait) == :repl
        PyPlot.subplot(224)
    elseif traitvalue(plotmodetrait) == :notebook
        PyPlot.subplot(122)
    end
    meanlogposteriorplot(model,logposteriors)
    #return the results
    paramvals,macrocurrentstats,logposteriors,f
end

function analysis3(model::GMHPhotoReceptor.PhotoReceptorModel,niterations::Integer;plotmode=:repl)
    #check the value of plotmode by converting it into a trait
    plotmodetrait = GeneralizedMetropolisHastings.trait(:plotmode,plotmode)
    #perform the required analysis
    nparamsets = 4
    paramvals,meancurrents,stdcurrents,logposteriors = evaluate3(model,niterations,nparamsets)
    #print the results
    printanalysis2("analysis3",paramvals,logposteriors)
    #plot the results
    if traitvalue(plotmodetrait) == :repl
        PyPlot.figure()
        PyPlot.subplot(221)
    elseif traitvalue(plotmodetrait) == :notebook
        PyPlot.figure()
        PyPlot.subplot(211)
    end
    meancurrentsplot(model,meancurrents)
    if traitvalue(plotmodetrait) == :repl
        PyPlot.subplot(222)
    elseif traitvalue(plotmodetrait) == :notebook
        PyPlot.figure()
        PyPlot.subplot(211)
    end
    stdcurrentplot(model,stdcurrents)
    if traitvalue(plotmodetrait) == :repl
        for i=1:4
            PyPlot.subplot(4,4,8+i)
            logposteriorplot(model,logposteriors[:,i])
            PyPlot.subplot(4,4,12+i)
            meanlogposteriorplot(model,logposteriors[:,i])
        end
    elseif traitvalue(plotmodetrait) == :notebook
        for s in Any[[1,2],[3,4]]
            PyPlot.figure()
            for (i,p) in enumerate(s)
                PyPlot.subplot(220+i)
                logposteriorplot(model,logposteriors[:,p])
                PyPlot.subplot(222+i)
                meanlogposteriorplot(model,logposteriors[:,p])
            end
        end
    end

    paramvals,meancurrents,stdcurrents,logposteriors
end

function analysis4(model::GMHPhotoReceptor.PhotoReceptorModel,niterations::Integer,nparamsets::Integer;plotmode=:repl)
    #check the value of plotmode by converting it into a trait
    plotmodetrait = GeneralizedMetropolisHastings.trait(:plotmode,plotmode)
    #perform the required analysis
    paramvals,meancurrents,stdcurrents,logposteriors = evaluate3(model,niterations,nparamsets)
    meanlogposteriors = mean(logposteriors,1)
    stdlogposteriors = std(logposteriors,1)
    #print the results
    printanalysis2("analysis4",paramvals,logposteriors)
    #plot the results
    if traitvalue(plotmodetrait) == :repl
        PyPlot.figure()
        PyPlot.subplot(221)
    elseif traitvalue(plotmodetrait) == :notebook
        PyPlot.figure()
        PyPlot.subplot(211)
    end
    meancurrentsplot(model,meancurrents)
    PyPlot.plot(dataindex(model),measurements(model);linewidth=2,color="yellow",label="measured")
    PyPlot.plot(dataindex(model),mean(meancurrents,2);linewidth=2,color="white",label="mean")
    if traitvalue(plotmodetrait) == :repl
        PyPlot.subplot(222)
    elseif traitvalue(plotmodetrait) == :notebook
        PyPlot.figure()
        PyPlot.subplot(211)
    end
    stdcurrentplot(model,stdcurrents)
    PyPlot.plot(dataindex(model),mean(stdcurrents,2);linewidth=2,color="white",label="mean")
    if traitvalue(plotmodetrait) == :repl
        PyPlot.subplot(223)
    elseif traitvalue(plotmodetrait) == :notebook
        PyPlot.figure(figsize=(10.0,10.0))
    end
    PyPlot.scatter(meanlogposteriors,stdlogposteriors)
    PyPlot.xlabel("Mean Log-Posterior")
    PyPlot.ylabel("STD Log-Posterior")
    if traitvalue(plotmodetrait) == :repl
        PyPlot.subplot(224)
    elseif traitvalue(plotmodetrait) == :notebook
        PyPlot.figure(figsize=(10.0,10.0))
    end
    PyPlot.scatter(meanlogposteriors,100*abs(stdlogposteriors./meanlogposteriors))
    PyPlot.ylim([0.0,5.0])
    PyPlot.xlabel("Mean Log-Posterior")
    PyPlot.ylabel("STD/Mean Log-Posterior")
    paramvals,meancurrents,stdcurrents,logposteriors
end

function analysis5(model::GMHPhotoReceptor.PhotoReceptorModel,niterations::Integer,gradienteps=1e-4;plotmode=:repl)
    #check the value of plotmode by converting it into a trait
    plotmodetrait = GeneralizedMetropolisHastings.trait(:plotmode,plotmode)
    #perform the required analysis
    paramvals,meancurrents,stdcurrents,logposteriors = gradient1(model,niterations,gradienteps)
    #print the results
    printanalysis2("analysis5",paramvals,logposteriors)
    printanalysis2("analysis5",paramvals,logposteriors-mean(logposteriors[:,1]))
    #plot the results
    if traitvalue(plotmodetrait) == :repl
        PyPlot.figure()
        PyPlot.subplot(221)
    elseif traitvalue(plotmodetrait) == :notebook
        PyPlot.figure()
        PyPlot.subplot(211)
    end
    meancurrentsplot(model,meancurrents)
    PyPlot.plot(dataindex(model),measurements(model);linewidth=2,color="yellow",label="measured")
    if traitvalue(plotmodetrait) == :repl
        PyPlot.subplot(222)
    elseif traitvalue(plotmodetrait) == :notebook
        PyPlot.figure()
        PyPlot.subplot(211)
    end
    stdcurrentplot(model,stdcurrents)
    if traitvalue(plotmodetrait) == :repl
        PyPlot.subplot(212)
    elseif traitvalue(plotmodetrait) == :notebook
        PyPlot.figure()
        PyPlot.subplot(211)
    end
    meanlogposteriorplot(model,logposteriors)
    paramvals,meancurrents,stdcurrents,logposteriors
end


text:
i = 1
hostname = "localhost"
port = 5000

max_stars_repo_name: keithseahus/msgpack-rpc-julia

max_issues_repo_name: keithseahus/msgpack-rpc-julia

max_forks_repo_name: keithseahus/msgpack-rpc-julia


text:
include("filtering.jl");
################################ Poke and Travel analysis ##################################
list = ["PreVehicle",
    "Altanserin",
    "SB242084",
    "Way_100135",
    "Methysergide",
    "Citalopram",
    "SB242084_opt",
    "Saline"]
filter!(r->r.Treatment in list &&
    r.Trial < 51 &&
    r.MouseID != "pc7",
    streaks)
streaks[streaks.Treatment .== "PreVehicle",:Treatment] .= "Control"
streaks[streaks.Treatment .== "Saline",:Treatment] = [o ? "Optogenetic" : "Control" for o in streaks[streaks.Treatment .== "Saline",:Stim]]
streaks[streaks.Treatment .== "SB242084_opt",:Phase] .=  "SB242084_opt"
streaks[streaks.Treatment .== "SB242084_opt",:Treatment] = [o ? "SB242084_opt" : "Control" for o in streaks[streaks.Treatment .== "SB242084_opt",:Stim]]
streaks[!,:ROILeavingTime] = streaks.Stop_trial .- streaks.Stop_poking
##
s = filter(r -> r.Phase in ["Optogenetic", "Citalopram"], streaks)
##
df1 = combine(groupby(s,:Phase)) do dd
    subdf = unstack(dd,:Treatment,:Num_pokes)
    current_drug = Symbol(subdf[1,:Phase])
    rename!(subdf, current_drug => :Drug)
end
open_html_table(df1)
df2 = combine(groupby(df1,[:Phase])) do dd
    wilcoxon(dd,:Control, :Drug; f = x -> mean(skipmissing(x)))
end
Drug_colors!(df2)
plot_wilcoxon(df2; sorting = Plotting_position)
##
savefig(joinpath(figs_loc,"FENS/Plot5_NumPokes_50trials.pdf"))
##
df1 = combine(groupby(s,:Phase)) do dd
    subdf = unstack(dd,:Treatment,:AfterLast)
    current_drug = Symbol(subdf[1,:Phase])
    rename!(subdf, current_drug => :Drug)
end
df2 = combine(groupby(df1,[:Phase])) do dd
    wilcoxon(dd,:Control, :Drug; f = x -> mean(skipmissing(x)))
end
Drug_colors!(df2)
plot_wilcoxon(df2; sorting = Plotting_position)
##
savefig(joinpath(figs_loc,"FENS/Plot9_Afterlast_50trials.pdf"))
##
df1 = combine(groupby(s,:Phase)) do dd
    subdf = unstack(dd,:Treatment,:Poking_Travel_to)
    current_drug = Symbol(subdf[1,:Phase])
    rename!(subdf, current_drug => :Drug)
end
df2 = combine(groupby(df1,[:Phase])) do dd
    wilcoxon(dd,:Control, :Drug; f = x -> mean(skipmissing(x)))
end
Drug_colors!(df2)
plot_wilcoxon(df2; sorting = Plotting_position)
##
savefig(joinpath(figs_loc,"FENS/Plot7_PokingTravel_50trials.pdf"))
##
df1 = combine(groupby(s,:Phase)) do dd
    subdf = unstack(dd,:Treatment,:Trial_Travel_to)
    current_drug = Symbol(subdf[1,:Phase])
    rename!(subdf, current_drug => :Drug)
end
df2 = combine(groupby(df1,[:Phase])) do dd
    wilcoxon(dd,:Control, :Drug; f = x -> mean(skipmissing(x)))
end
Drug_colors!(df2)
plot_wilcoxon(df2; sorting = Plotting_position)
##
savefig(joinpath(figs_loc,"FENS/Plot8_ROITravel_50trials.pdf"))
##
df1 = combine(groupby(s,:Phase)) do dd
    subdf = unstack(dd,:Treatment,:ROILeavingTime)
    current_drug = Symbol(subdf[1,:Phase])
    rename!(subdf, current_drug => :Drug)
end
df2 = combine(groupby(df1,[:Phase])) do dd
    wilcoxon(dd,:Control, :Drug; f = x -> mean(skipmissing(x)))
end
Drug_colors!(df2)
plot_wilcoxon(df2; sorting = Plotting_position)
##
savefig(joinpath(figs_loc,"FENS/Plot14_ROILeaving_50trials.pdf"))
################################ Leaving analysis ##################################
odc = ODC(pokes)
@df filter(r -> r.PokeDuration < 4,odc) density(:PokeDuration,group = :Reward)
plot!(repeat([0.43],2), [0,8])
plot!(repeat([0.7],2), [0,8])
odc[!,:ClearDuration] = [r.Reward ? r.PokeDuration > 0.7 : r.PokeDuration < 0.4 for r in eachrow(odc)]
list = ["PreVehicle",
    "Altanserin",
    "SB242084",
    "Way_100135",
    "Methysergide",
    "Citalopram",
    "SB242084_opt",
    "Saline"]

filter!(r -> r.Treatment in list &&
    !ismissing(r.Pre_Interpoke) &&
    r.Pre_Interpoke < 0.5 &&
    r.ClearDuration &&
    r.Trial < 51 &&
    r.MouseID != "d5" &&
    r.MouseID != "pc7",# ["pc1","pc2","pc3","pc4","pc5","pc6","pc8","pc9","pc10"],
    odc)

odc[odc.Treatment .== "PreVehicle",:Treatment] .= "Control"
odc[odc.Treatment .== "Saline",:Treatment] = [o ? "Optogenetic" : "Control" for o in odc[odc.Treatment .== "Saline",:Stim]]
odc[odc.Treatment .== "SB242084_opt",:Phase] .=  "SB242084_opt"
odc[odc.Treatment .== "SB242084_opt",:Treatment] = [o ? "SB242084_opt" : "Control" for o in odc[odc.Treatment .== "SB242084_opt",:Stim]]

df1 = filter(r -> !r.Reward &&
    r.ExpDay > 3 &&
    !ismissing(r.PokeFromLastRew) &&
    0 < r.PokeFromLastRew &&
    !ismissing(r.ODC) &&
    r.Treatment in ["Control","Citalopram", "Optogenetic", "Altanserin","SB242084"] &&
    !in(r.Phase,["Methysergide", "Way_100135","SB242084_opt"]) ,odc)

categorical!(df1,[:MouseID,:Treatment,:Phase])

levels!(df1.Treatment, ["Control",
    "Citalopram",
    "Optogenetic",
    "Altanserin",
    "SB242084"])
# Normalise omission duty cicle with animal session mean
transform!(groupby(df1,[:MouseID,:Phase,:Treatment]), :ODC => mean => :ODC_mean)
df1[!,:NODC] = df1.ODC ./ df1.ODC_mean
# Convert Normalise ODC in 15 quantiles
transform!(groupby(df1,[:MouseID,:Phase,:Treatment]), :NODC => binquantile => :QODC)
df1[!,:Leave] = [r == 0 for r in df1.PokeFromLeaving]
##
filtered = filter(r -> r.PokeFromLeaving < 7 && r.Phase in ["Optogenetic"], df1)
ps = plot_QODC(filtered,:PokeFromLeaving; xflip = true, ylims = :auto)
##
savefig(joinpath(figs_loc,"Report_sept_2020/QODCfromLeaving_opto_50trials.pdf"))
##
filtered = filter(r -> r.PokeFromLeaving == 0 &&
    r.Phase in ["Optogenetic","Citalopram","SB242084","Altanserin"], df1)
gd2 = groupby(filtered,[:MouseID,:Phase,:Treatment])
df2 = combine(gd2, :QODC => mean => :QODC)
gd3 = groupby(df2,[:Phase])
df3 = combine(gd3) do dd
    dd[:,:Treatment] = [t == "Control" ? t : "Manipulation" for t in dd.Treatment]
    unstack(dd,:Treatment, :QODC)
end
dropmissing!(df3)
gd4 = groupby(df3,:Phase)
df4 = combine(gd4) do dd
    wilcoxon(dd,:Control,:Manipulation)
end
df4[!,:Phase] = String.(df4.Phase)
Drug_colors!(df4)
plot_wilcoxon(df4; sorting = Plotting_position)
##
savefig(joinpath(figs_loc,"FENS/Plot13_QODCatLeaving_50trials.pdf"))


text:
"""
# script bench_dim

- Julia version: 1.4
- Author: niallcullinane
- Date: 2020-08-20

Script to benchmark model with different input sizes for CPU and GPU.
"""

using DrWatson
@quickactivate "MEngProject"
using MEngProject,
    CUDA,
    DifferentialEquations,
    PyPlot,
    NNlib,
    ImageFiltering,
    Images,
    MEngProject.LaminartKernels,
    MEngProject.LaminartInitFunc,
    MEngProject.Utils,
    BenchmarkTools,
    Test

using OrdinaryDiffEq,
    ParameterizedFunctions, LSODA, Sundials, DiffEqDevTools, Noise

batch = 1000

# files = readdir(datadir("res_test"))[2:end]
files = [
    "ver_lines_025_gs.png",
    "ver_lines_050_gs.png",
    "ver_lines_075_gs.png",
    "ver_lines_100_gs.png",
    "ver_lines_200_gs.png",
    "ver_lines_300_gs.png",
    "ver_lines_400_gs.png",
]

global benchm_gpu
global benchm_cpu
global y1Res_gpu
global y1Res_cpu
global prob_d

benchm_gpu = []
benchm_cpu = []
y1Res_gpu = []
y1Res_cpu = []

tspan = (0.0f0, 10f0)

batch_ = string(batch, "_", rand(1000:9999))
mkdir(plotsdir(string("bench_dim", batch_)))

test_name = ["025", "050", "075", "100", "200", "300", "400"]
test_name_plt = [
    "\$25×25\$",
    "\$50×50\$",
    "\$75×75\$",
    "\$100×100\$",
    "\$200×200\$",
    "\$300×300\$",
    "\$400×400\$",
]

for file in enumerate(files)
    # try
    p = LaminartInitFunc.parameterInit_conv_gpu(
        datadir("res_test", file[2]),
        Parameters.parameters_f32,
    )

    u0 = cu(reshape(
        zeros(Float32, p.dim_i, p.dim_j * (5 * p.K + 2)),
        p.dim_i,
        p.dim_j,
        5 * p.K + 2,
        1,
    ))

    arr1 = similar(@view u0[:, :, 1:2, :])
    arr2 = similar(@view u0[:, :, 1:1, :])

    f = LaminartFunc.LamFunction(
        arr1, #x
        similar(arr1), #m
        similar(arr1), #s
        arr2, #x_lgn,
        similar(arr1), #C,
        similar(arr1), #H_z,
        similar(arr1), # dy_temp,
        similar(arr1), # dm_temp,
        similar(arr1), # dz_temp,
        similar(arr1), # ds_temp,
        similar(arr2), # dv_temp,
        similar(arr1), # H_z_temp,
        similar(arr2), #  V_temp_1,
        similar(arr2), #  V_temp_2,
        similar(arr1), #  Q_temp,
        similar(arr1), #   P_temp
    )

    prob_d = ODEProblem(f, u0, tspan, p)
    push!(benchm_gpu, @benchmark solve(prob_d))
    sol = solve(prob_d)

    # @inbounds begin
    #     t = 10
    #     v0 = @view sol(t)[:, :, :, 1]
    #     axMax = findmax(v0)[1]
    #
    #
    #     k = 7
    #
    #     fig, ax = plt.subplots()
    #
    #     v1 = @view sol(t)[:, :, k, 1]
    #     v2 = @view sol(t)[:, :, k+1, 1]
    #     im = ax.imshow(
    #         v1,
    #         cmap = matplotlib.cm.PRGn,
    #         vmax = axMax,
    #         vmin = -axMax,
    #     )
    #     im2 = ax.imshow(
    #         v2,
    #         cmap = matplotlib.cm.RdBu_r,
    #         vmax = axMax,
    #         vmin = -axMax,
    #         alpha = 0.5,
    #     )
    #
    #     cbar = fig.colorbar(im2, shrink = 0.9, ax = ax)
    #     cbar.ax.set_xlabel("\$k=2\$")
    #     cbar = fig.colorbar(im, shrink = 0.9, ax = ax)
    #     cbar.set_alpha(0.5)
    #     cbar.draw_all()
    #     cbar.ax.set_xlabel("\$k=1\$")
    #     layer = Utils.layers[k]
    #     plt.title(string(
    #         "Layer: $layer, \$t=$t\$, resolution=",
    #         test_name_plt[file[1]],
    #     ))
    #     plt.axis("off")
    #     fig.tight_layout()
    #     plt.savefig(plotsdir(
    #         string("bench_dim", batch_),
    #         string(
    #             file[2],
    #             "_res_",
    #             test_name[file[1]],
    #             "_t",
    #             t,
    #             "_",
    #             Utils.la[k],
    #             ".png",
    #         ),
    #     ))
    #     close("all")
    #
    #     v3 = @view sol[:, :, 7:7, :, :]
    #     push!(y1Res, Array(v3))

    #     end
    # 	catch err
    # finally
    # 	u0 = nothing
    #         p = nothing
    #         arr1 = nothing
    #         arr2 = nothing
    #         f = nothing
    #         prob_d = nothing
    #         sol = nothing
    # end

    try

        p = LaminartInitFunc.parameterInit_conv_cpu(
            datadir("res_test", file[2]),
            Parameters.parameters_f32,
        )

        u0 = reshape(
            zeros(Float32, p.dim_i, p.dim_j * (5 * p.K + 2)),
            p.dim_i,
            p.dim_j,
            5 * p.K + 2,
            1,
        )

        arr1 = similar(@view u0[:, :, 1:2, :])
        arr2 = similar(@view u0[:, :, 1:1, :])

        f = LaminartFunc.LamFunction(
            arr1, #x
            similar(arr1), #m
            similar(arr1), #s
            arr2, #x_lgn,
            similar(arr1), #C,
            similar(arr1), #H_z,
            similar(arr1), # dy_temp,
            similar(arr1), # dm_temp,
            similar(arr1), # dz_temp,
            similar(arr1), # ds_temp,
            similar(arr2), # dv_temp,
            similar(arr1), # H_z_temp,
            similar(arr2), #  V_temp_1,
            similar(arr2), #  V_temp_2,
            similar(arr1), #  Q_temp,
            similar(arr1), #   P_temp
        )

        prob_d = ODEProblem(f, u0, tspan, p)
        push!(benchm_cpu, @benchmark solve(prob_d))
        sol = solve(prob_d)

        push!(y1Res_cpu, sol[:, :, 7:7, :, :])
    catch
    finally
        u0 = nothing
        p = nothing
        arr1 = nothing
        arr2 = nothing
        f = nothing
        prob_d = nothing
        sol = nothing
    end
end

# # time plot
# fig, axs = plt.subplots()
#
# for result ∈ enumerate(y1Res_gpu)
#     lab = "$test_name[result[1]]"
#     axs.plot(
#         result[2][
#             findmax(result[2][:, :, 1, 1, end])[2][1],
#             findmax(result[2][:, :, 1, 1, end])[2][2],
#             k,
#             1,
#             :,
#         ],
#         c = Utils.Colour[result[1]],
#         "--",
#         label = "$lab GPU",
#         alpha=0.8
#     )
# end
#
# for result ∈ enumerate(y1Res_cpu)
#     lab = "$test_name[result[1]]"
#     axs.plot(
#         result[2][
#             findmax(result[2][:, :, 1, 1, end])[2][1],
#             findmax(result[2][:, :, 1, 1, end])[2][2],
#             k,
#             1,
#             :,
#         ],
#         c = Utils.Colour[result[1]],
#         ":",
#         label = "$lab CPU",
#         alpha=0.8
#     )
# end
# axs.set_xlabel("Steps")
# axs.set_ylabel("Activation")
# plt.title("L2/3, \$k=1\$")
# plt.legend()
# fig.tight_layout()
# plt.savefig(plotsdir(
#     string("bench_dim", batch_),
#     string(file[2], "_para_", test_name[file[1]], "_time.png"),
# ))
# close("all")

# benchmark plot

fig, ax = plt.subplots()
for bm_d ∈ enumerate(benchm_gpu)
    ax.scatter(
        test_name_plt[bm_d[1]],
        median(bm_d[2].times) * 1e-9,
        label = "GPU",
        color = Utils.colours[1],
        edgecolors = "none",
    )
end

for bm_d ∈ enumerate(benchm_cpu)
    ax.scatter(
        test_name_plt[bm_d[1]],
        median(bm_d[2].times) * 1e-9,
        label = "CPU",
        color = Utils.colours[2],
        edgecolors = "none",
    )
end

ax.legend()
axs.set_xlabel("Resolution (\$px\$)")
axs.set_ylabel("Time (\$s\$)")
ax.grid(true)
ax.set_ylim(ymin = 0)
fig.tight_layout()
plt.savefig(plotsdir(
    string("bench_dim", batch_),
    string(file[2], "_para_", test_name[file[1]], "_time.png"),
))
close("all")

# memory
fig, ax = plt.subplots()
for bm_d ∈ enumerate(benchm_gpu)
    ax.scatter(
        test_name_plt[bm_d[1]],
        bm_d[2].memory * 1e-6,
        color = Utils.colours[1],
        label = "GPU",
        edgecolors = "none",
    )
end

ax.legend()
axs.set_xlabel("Resolution (\$px\$)")
axs.set_ylabel("Memory")
ax.grid(true)
ax.set_ylim(ymin = 0)
fig.tight_layout()
plt.savefig(plotsdir(
    string("bench_dim", batch_),
    string(file[2], "_para_", test_name[file[1]], "_time.png"),
))
close("all")

benchm_gpu = nothing
benchm_cpu = nothing
y1Res_gpu = nothing
y1Res_cpu = nothing
prob_d = nothing
CUDA.reclaim()


text:
function lagrangian_multiplier(model::MOI.ModelLike, p, certificate, index, domain, T::Type)
    monos = Certificate.get(certificate, Certificate.MultiplierBasis(), index, domain, p)
    MCT = SOS.matrix_cone_type(typeof(certificate))
    return SOS.add_gram_matrix(model, MCT, monos)..., monos
end

struct SOSPolynomialInSemialgebraicSetBridge{
    T, F <: MOI.AbstractVectorFunction, DT <: SemialgebraicSets.AbstractSemialgebraicSet,
    CT <: Certificate.AbstractIdealCertificate,
    UMCT <: MOI.ConstraintIndex{MOI.VectorOfVariables},
    UMST,
    MT <: MP.AbstractMonomial, MVT <: AbstractVector{MT}} <: MOIB.Constraint.AbstractBridge
    lagrangian_monomials::Vector{MVT}
    lagrangian_variables::Vector{Vector{MOI.VariableIndex}}
    lagrangian_constraints::Vector{UMCT}
    constraint::MOI.ConstraintIndex{F, SOS.SOSPolynomialSet{DT, MT, MVT, CT}}
    monomials::MVT
end

function MOI.Bridges.Constraint.bridge_constraint(
    ::Type{SOSPolynomialInSemialgebraicSetBridge{T, F, DT, CT, UMCT, UMST, MT, MVT}},
    model::MOI.ModelLike, f::MOI.AbstractVectorFunction,
    set::SOS.SOSPolynomialSet{<:SemialgebraicSets.BasicSemialgebraicSet}) where {T, F, DT, CT, UMCT, UMST, MT, MVT}

    @assert MOI.output_dimension(f) == length(set.monomials)
    p = MP.polynomial(collect(MOIU.eachscalar(f)), set.monomials)
    n = length(set.domain.p)
    λ_monos     = MVT[]
    λ_variables = Vector{MOI.VariableIndex}[]
    λ_constraints = UMCT[]
    for index in Certificate.get(set.certificate, Certificate.PreorderIndices(), set.domain)
        λ, λ_variable, λ_constraint, λ_mono = lagrangian_multiplier(
            model, p, set.certificate, index, set.domain, T)
        push!(λ_variables, λ_variable)
        push!(λ_constraints, λ_constraint)
        push!(λ_monos, λ_mono)
        # As `*(::MOI.ScalarAffineFunction{T}, ::S)` is only defined if `S == T`, we
        # need to call `changecoefficienttype`. This is critical since `T` is
        # `Float64` when used with JuMP and the coefficient type is often `Int` if
        # `set.domain.V` is `FullSpace` or `FixedPolynomialsSet`.
        g = Certificate.get(set.certificate, Certificate.Generator(), index, set.domain)
        p -= λ * MP.changecoefficienttype(g, T)
    end
    new_set = SOS.SOSPolynomialSet(
        set.domain.V, MP.monomials(p),
        Certificate.get(set.certificate, Certificate.IdealCertificate()))
    constraint = MOI.add_constraint(model, MOIU.vectorize(MP.coefficients(p)),
                                    new_set)

    return SOSPolynomialInSemialgebraicSetBridge{
        T, F, DT, CT, UMCT, UMST, MT, MVT}(λ_monos, λ_variables, λ_constraints, constraint, set.monomials)
end

function MOI.supports_constraint(::Type{SOSPolynomialInSemialgebraicSetBridge{T}},
                                 ::Type{<:MOI.AbstractVectorFunction},
                                 ::Type{<:SOS.SOSPolynomialSet{<:SemialgebraicSets.BasicSemialgebraicSet}}) where T
    return true
end
function MOIB.added_constrained_variable_types(::Type{<:SOSPolynomialInSemialgebraicSetBridge{T, F, DT, CT}}) where {T, F, DT, CT}
    return constrained_variable_types(SOS.matrix_cone_type(CT))
end
function MOIB.added_constraint_types(
    ::Type{SOSPolynomialInSemialgebraicSetBridge{T, F, DT, CT, UMCT, UMST, MT, MVT}}) where {T, F, DT, CT, UMCT, UMST, MT, MVT}
    return [(F, SOS.SOSPolynomialSet{DT, MT, MVT, CT})]
end
function MOIB.Constraint.concrete_bridge_type(
    ::Type{<:SOSPolynomialInSemialgebraicSetBridge{T}},
    F::Type{<:MOI.AbstractVectorFunction},
    ::Type{<:SOS.SOSPolynomialSet{SemialgebraicSets.BasicSemialgebraicSet{S, PS, AT}, MT, MVT, CT}}) where {T, S, PS, AT, CT, MT, MVT}

    # promotes VectorOfVariables into VectorAffineFunction, it should be enough
    # for most use cases
    G = MOIU.promote_operation(-, T, F, MOI.VectorOfVariables)
    MCT = SOS.matrix_cone_type(CT)
    UMCT = union_constraint_types(MCT)
    UMST = union_set_types(MCT)
    IC = Certificate.get(CT, Certificate.IdealCertificate())
    return SOSPolynomialInSemialgebraicSetBridge{T, G, AT, IC, UMCT, UMST, MT, MVT}
end

# Attributes, Bridge acting as an model
function MOI.get(bridge::SOSPolynomialInSemialgebraicSetBridge,
                 ::MOI.NumberOfVariables)
    return reduce((n, vis) -> n + length(vis), bridge.lagrangian_variables, init=0)
end
function MOI.get(bridge::SOSPolynomialInSemialgebraicSetBridge, ::MOI.ListOfVariableIndices)
    return Iterators.flatten(bridge.lagrangian_variables)
end
function MOI.get(bridge::SOSPolynomialInSemialgebraicSetBridge{T, F, DT, CT, UMCT, UMST},
                 ::MOI.NumberOfConstraints{MOI.VectorOfVariables, S}) where {T, F, DT, CT, UMCT, UMST, S<:UMST}
    return reduce((n, ci) -> n + (ci isa MOI.ConstraintIndex{MOI.VectorOfVariables, S}),
                  bridge.lagrangian_constraints, init=0)
end
function MOI.get(b::SOSPolynomialInSemialgebraicSetBridge{T, F, DT, CT, UMCT, UMST},
                 ::MOI.ListOfConstraintIndices{MOI.VectorOfVariables, S}) where {T, F, DT, CT, UMCT, UMST, S<:UMST}
    return filter(ci -> ci isa MOI.ConstraintIndex{MOI.VectorOfVariables, S}, bridge.lagrangian_constraints)
end
function MOI.get(::SOSPolynomialInSemialgebraicSetBridge{T, F, DT, CT, UMCT, UMST, MT, MVT},
                 ::MOI.NumberOfConstraints{F, SOS.SOSPolynomialSet{DT, MT, MVT, CT}}) where {T, F, DT, CT, UMCT, UMST, MT, MVT}
    return 1
end
function MOI.get(b::SOSPolynomialInSemialgebraicSetBridge{T, F, DT, CT, UMCT, UMST, MT, MVT},
                 ::MOI.ListOfConstraintIndices{F, SOS.SOSPolynomialSet{DT, MT, MVT, CT}}) where {T, F, DT, CT, UMCT, UMST, MT, MVT}
    return [b.constraint]
end

# Indices
function MOI.delete(model::MOI.ModelLike,
                    bridge::SOSPolynomialInSemialgebraicSetBridge)
    MOI.delete(model, bridge.constraint)
    for variables in bridge.lagrangian_variables
        MOI.delete(model, variables)
    end
end

# Attributes, Bridge acting as a constraint

# The monomials might be different from the ones of the original polynomial
# because of the ∑ λ_i s_i(x) so we don't define ConstraintPrimal and
# ConstraintDual, as the caller won't know how to reshape it
function MOI.get(::MOI.ModelLike,
                 ::MOI.ConstraintPrimal,
                 ::SOSPolynomialInSemialgebraicSetBridge)
    throw(SOS.ValueNotSupported())
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintDual,
                 bridge::SOSPolynomialInSemialgebraicSetBridge)
    dual = MOI.get(model, attr, bridge.constraint)
    set = MOI.get(model, MOI.ConstraintSet(), bridge.constraint)
    μ = MultivariateMoments.measure(dual, set.monomials)
    return [dot(mono, μ) for mono in bridge.monomials]
end
function MOI.get(model::MOI.ModelLike, attr::PolyJuMP.MomentsAttribute,
                 bridge::SOSPolynomialInSemialgebraicSetBridge)
    return MOI.get(model, attr, bridge.constraint)
end

function MOI.get(model::MOI.ModelLike,
                 attr::Union{SOS.CertificateMonomials, SOS.GramMatrixAttribute,
                             SOS.MomentMatrixAttribute},
                 bridge::SOSPolynomialInSemialgebraicSetBridge)
    return MOI.get(model, attr, bridge.constraint)
end
function MOI.get(model::MOI.ModelLike, attr::SOS.LagrangianMultipliers,
                 bridge::SOSPolynomialInSemialgebraicSetBridge)
    @assert eachindex(bridge.lagrangian_variables) == eachindex(bridge.lagrangian_monomials)
    map(i -> SOS.build_gram_matrix(MOI.get(model, MOI.VariablePrimal(attr.N),
                                           bridge.lagrangian_variables[i]),
                                   bridge.lagrangian_monomials[i]),
        eachindex(bridge.lagrangian_variables))
end


text:
using Base.Test

include("hamming.jl")

@testset "identical strands" begin
    @test distance("A", "A") == 0
    @test distance("GGACTGA", "GGACTGA") == 0
end

@testset "complete distance" begin
    @test distance("A", "G") == 1
    @test distance("AG", "CT") == 2
end

@testset "small distance" begin
    @test distance("AT", "CT") == 1
    @test distance("GGACG", "GGTCG") == 1
    @test distance("ACCAGGG", "ACTATGG") == 2
end

@testset "non-unique characters" begin
    @test distance("AGA", "AGG") == 1
    @test distance("AGG", "AGA") == 1
end

@testset "same nucleotides in different positions" begin
    @test distance("TAG", "GAT") == 2
end

@testset "large distance" begin
    @test distance("GATACA", "GCATAA") == 4
    @test distance("GGACGGATTCTG", "AGGACGGATTCT") == 9
end

@testset "empty strands" begin
    @test distance("", "") == 0
end

@testset "different strand lengths" begin
    @test_throws ArgumentError distance("AATG", "AAA")
    @test_throws ArgumentError distance("ATA", "AGTG")
end

max_stars_repo_name: Devlious/LearningJulia

max_issues_repo_name: Devlious/LearningJulia

max_forks_repo_name: Devlious/LearningJulia


text:
@testset "NumField/NonSimpleNumField" begin
  include("NonSimpleNumField/Conjugates.jl")
  include("NonSimpleNumField/ComplexEmbeddings.jl")
end


text:
"""
$(TYPEDEF)

Dense Cholesky-based implementation for the sparse positive semidefinite
cone [`PosSemidefTriSparse`](@ref).
"""
struct PSDSparseDense <: PSDSparseImpl end

mutable struct PSDSparseDenseCache{T <: Real, R <: RealOrComplex{T}} <:
    PSDSparseCache{T, R}
    mat::Matrix{R}
    mat2::Matrix{R}
    inv_mat::Matrix{R}
    fact_mat
    PSDSparseDenseCache{T, R}() where {T <: Real, R <: RealOrComplex{T}} =
        new{T, R}()
end

function setup_extra_data!(
    cone::PosSemidefTriSparse{PSDSparseDense, T, R},
    ) where {T, R}
    cone.cache = cache = PSDSparseDenseCache{T, R}()
    cache.mat = zeros(R, cone.side, cone.side)
    cache.mat2 = zero(cache.mat)
    cache.inv_mat = zero(cache.mat)
    return
end

function update_feas(cone::PosSemidefTriSparse{PSDSparseDense})
    @assert !cone.feas_updated
    cache = cone.cache
    Λ = cache.mat

    svec_to_smat_sparse!(Λ, cone.point, cone)
    cache.fact_mat = cholesky!(Hermitian(Λ, :L), check = false)
    cone.is_feas = isposdef(cache.fact_mat)

    cone.feas_updated = true
    return cone.is_feas
end

function update_grad(cone::PosSemidefTriSparse{PSDSparseDense})
    @assert !cone.grad_updated && cone.is_feas
    cache = cone.cache
    inv_mat = cache.inv_mat

    inv_fact!(inv_mat, cache.fact_mat)
    copytri!(inv_mat, 'L', true)
    smat_to_svec_sparse!(cone.grad, inv_mat, cone)
    cone.grad .*= -1

    cone.grad_updated = true
    return cone.grad
end

function update_hess(
    cone::PosSemidefTriSparse{PSDSparseDense, T, T},
    ) where {T <: Real}
    @assert cone.grad_updated
    isdefined(cone, :hess) || alloc_hess!(cone)
    rt2 = cone.rt2
    H = cone.hess.data
    fill!(H, 0)
    Λi = cone.cache.inv_mat

    @inbounds for (idx2, (i2, j2)) in enumerate(zip(cone.row_idxs, cone.col_idxs))
        for idx1 in 1:idx2
            (i1, j1) = (cone.row_idxs[idx1], cone.col_idxs[idx1])
            H[idx1, idx2] = begin
                if (i1 == j1) && (i2 == j2)
                    abs2(Λi[i1, i2])
                elseif xor(i1 == j1, i2 == j2)
                    rt2 * Λi[i1, i2] * Λi[j1, j2]
                else
                    Λi[i1, i2] * Λi[j1, j2] + Λi[i1, j2] * Λi[j1, i2]
                end
            end
        end
    end

    cone.hess_updated = true
    return cone.hess
end

function update_hess(
    cone::PosSemidefTriSparse{PSDSparseDense, T, Complex{T}},
    ) where {T <: Real}
    @assert cone.grad_updated
    isdefined(cone, :hess) || alloc_hess!(cone)
    rt2 = cone.rt2
    H = cone.hess.data
    fill!(H, 0)
    Λi = cone.cache.inv_mat

    idx2 = 1
    @inbounds for (i2, j2) in zip(cone.row_idxs, cone.col_idxs)
        idx1 = 1
        if i2 == j2
            for (i1, j1) in zip(cone.row_idxs, cone.col_idxs)
                if i1 == j1
                    H[idx1, idx2] = abs2(Λi[i1, i2])
                    idx1 += 1
                else
                    c = rt2 * Λi[i2, i1] * Λi[j1, j2]
                    H[idx1, idx2] = real(c)
                    idx1 += 1
                    H[idx1, idx2] = -imag(c)
                    idx1 += 1
                end
                (idx1 > idx2) && break
            end
            idx2 += 1
        else
            for (i1, j1) in zip(cone.row_idxs, cone.col_idxs)
                if i1 == j1
                    c = rt2 * Λi[i1, i2] * Λi[j2, j1]
                    H[idx1, idx2] = real(c)
                    H[idx1, idx2 + 1] = -imag(c)
                    idx1 += 1
                else
                    b1 = Λi[i1, i2] * Λi[j2, j1]
                    b2 = Λi[j1, i2] * Λi[j2, i1]
                    c1 = b1 + b2
                    H[idx1, idx2] = real(c1)
                    H[idx1, idx2 + 1] = -imag(c1)
                    idx1 += 1
                    c2 = b1 - b2
                    H[idx1, idx2] = imag(c2)
                    H[idx1, idx2 + 1] = real(c2)
                    idx1 += 1
                end
                (idx1 > idx2) && break
            end
            idx2 += 2
        end
    end

    cone.hess_updated = true
    return cone.hess
end

function hess_prod_slow!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::PosSemidefTriSparse{PSDSparseDense},
    )
    cone.use_hess_prod_slow_updated || update_use_hess_prod_slow(cone)
    @assert cone.hess_updated
    cone.use_hess_prod_slow || return hess_prod!(prod, arr, cone)
    @assert is_feas(cone)
    cache = cone.cache

    @inbounds @views for j in 1:size(arr, 2)
        Λ = svec_to_smat_sparse!(cache.mat2, arr[:, j], cone)
        copytri!(Λ, 'L', true)
        ldiv!(cache.fact_mat, Λ)
        rdiv!(Λ, cache.fact_mat)
        smat_to_svec_sparse!(prod[:, j], Λ, cone)
    end

    return prod
end

function dder3(
    cone::PosSemidefTriSparse{PSDSparseDense},
    dir::AbstractVector,
    )
    @assert is_feas(cone)
    cache = cone.cache

    Λ = svec_to_smat_sparse!(cache.mat2, dir, cone)
    copytri!(Λ, 'L', true)
    ldiv!(cache.fact_mat.L, Λ)
    rdiv!(Λ, cache.fact_mat)
    outer_prod_vec_sparse!(cone.dder3, Λ, cone)

    return cone.dder3
end

function outer_prod_vec_sparse!(
    vec::AbstractVector{T},
    mat::AbstractMatrix{T},
    cone::PosSemidefTriSparse{PSDSparseDense, T, T},
    ) where {T <: Real}
    @assert length(vec) == length(cone.row_idxs)
    @inbounds for (idx, (i, j)) in enumerate(zip(cone.row_idxs, cone.col_idxs))
        @views x = dot(mat[:, i], mat[:, j])
        if i != j
            x *= cone.rt2
        end
        vec[idx] = x
    end
    return vec
end

function outer_prod_vec_sparse!(
    vec::AbstractVector{T},
    mat::AbstractMatrix{Complex{T}},
    cone::PosSemidefTriSparse{PSDSparseDense, T, Complex{T}},
    ) where {T <: Real}
    idx = 1
    @inbounds for (i, j) in zip(cone.row_idxs, cone.col_idxs)
        @views x = dot(mat[:, i], mat[:, j])
        if i == j
            vec[idx] = real(x)
            idx += 1
        else
            x *= cone.rt2
            vec[idx] = real(x)
            vec[idx + 1] = imag(x)
            idx += 2
        end
    end
    @assert idx == length(vec) + 1
    return vec
end

function svec_to_smat_sparse!(
    mat::AbstractMatrix{T},
    vec::AbstractVector{T},
    cone::PosSemidefTriSparse{PSDSparseDense, T, T},
    ) where {T <: Real}
    @assert length(vec) == length(cone.row_idxs)
    fill!(mat, 0)
    @inbounds for (idx, (i, j)) in enumerate(zip(cone.row_idxs, cone.col_idxs))
        x = vec[idx]
        if i != j
            x /= cone.rt2
        end
        mat[i, j] = x
    end
    return mat
end

function svec_to_smat_sparse!(
    mat::AbstractMatrix{Complex{T}},
    vec::AbstractVector{T},
    cone::PosSemidefTriSparse{PSDSparseDense, T, Complex{T}},
    ) where {T <: Real}
    fill!(mat, 0)
    idx = 1
    @inbounds for (i, j) in zip(cone.row_idxs, cone.col_idxs)
        if i == j
            mat[i, j] = vec[idx]
            idx += 1
        else
            mat[i, j] = Complex(vec[idx], vec[idx + 1]) / cone.rt2
            idx += 2
        end
    end
    @assert idx == length(vec) + 1
    return mat
end

function smat_to_svec_sparse!(
    vec::AbstractVector{T},
    mat::AbstractMatrix{T},
    cone::PosSemidefTriSparse{PSDSparseDense, T, T},
    ) where {T <: Real}
    @assert length(vec) == length(cone.row_idxs)
    fill!(vec, 0)
    @inbounds for (idx, (i, j)) in enumerate(zip(cone.row_idxs, cone.col_idxs))
        x = mat[i, j]
        if i != j
            x *= cone.rt2
        end
        vec[idx] = x
    end
    return vec
end

function smat_to_svec_sparse!(
    vec::AbstractVector{T},
    mat::AbstractMatrix{Complex{T}},
    cone::PosSemidefTriSparse{PSDSparseDense, T, Complex{T}},
    ) where {T <: Real}
    fill!(vec, 0)
    idx = 1
    @inbounds for (i, j) in zip(cone.row_idxs, cone.col_idxs)
        x = mat[i, j]
        if i == j
            vec[idx] = real(x)
            idx += 1
        else
            x *= cone.rt2
            vec[idx] = real(x)
            vec[idx + 1] = imag(x)
            idx += 2
        end
    end
    @assert idx == length(vec) + 1
    return vec
end


text:
export GLGM06

struct GLGM06 <: ContinuousPost
    options::TwoLayerOptions

    function GLGM06(𝑂::Options)
        𝑂new = validate_and_wrap_options(𝑂, options_GLGM06())
        return new(𝑂new)
    end
end

# convenience constructor from pairs of symbols
GLGM06(𝑂::Pair{Symbol,<:Any}...) = GLGM06(Options(Dict{Symbol,Any}(𝑂)))

# default options (they are added in the function validate_and_wrap_options)
GLGM06() = GLGM06(Options())

include("init.jl")
include("post.jl")
include("reach.jl")


text:
""""
Add LaTeX, Emoji, HTML, and Unicode names to Str string literals

Copyright 2018-2021 Gandalf Software, Inc., Scott P. Jones
Licensed under MIT License, see LICENSE.md
"""
module StrEntities

using ModuleInterfaceTools, StrTables
using HTML_Entities, LaTeX_Entities, Emoji_Entities, Unicode_Entities

@api extend! StrLiterals

# Export from here as well so that people can do just `using StrEntities`
export @f_str, @pr_str, @F_str, @PR_str, @sym_str

function _parse_entity(io, str, pos, begseq, fin, tab, nam)
    beg = pos # start location
    chr, pos = iterate(str, pos)
    while chr != fin
        check_done(str, pos, "\\$begseq missing ending $fin")
        chr, pos = iterate(str, pos)
    end
    seq = lookupname(tab, str[beg:pos-2])
    seq == "" && throw_arg_err("Invalid $nam name in ", str)
    print(io, seq)
    pos
end

_parse_html(io, str, pos, chr)  =
    _parse_entity(io, str, pos, "&", ';', HTML_Entities.default, "HTML")
_parse_latex(io, str, pos, chr) =
    _parse_entity(io, str, pos, "<", '>', LaTeX_Entities.default, "LaTeX")
_parse_emoji(io, str, pos, chr) =
    _parse_entity(io, str, pos, ":", ':', Emoji_Entities.default, "Emoji")

function _parse_unicode(io, str, pos, chr)
    check_done(str, pos, "\\N incomplete")
    chr, pos = iterate(str, pos)
    chr == '{' || throw_arg_err("\\N missing initial { in ", str)
    check_done(str, pos, "\\N{ incomplete")
    _parse_entity(io, str, pos, "N{", '}', Unicode_Entities.default, "Unicode")
end

function __init__()
    StrLiterals.parse_chr[':'] = _parse_emoji
    StrLiterals.parse_chr['&'] = _parse_html
    StrLiterals.parse_chr['<'] = _parse_latex
    StrLiterals.parse_chr['N'] = _parse_unicode
end

end # module StrEntities


text:
type FastCodingData{T}
  Dtx::Array{T, 2};
  DtD::Array{T, 2};
end
type FastCodingOpts{T}
  z::Array{T, 2};
  t::Array{T, 2};
  lambda_t::Array{T, 2};
  beta::T;       # 1
  mu_t::T;       # 1
  max_iter::Int; # 1000
  tol::T;        # 1e-10
  DEBUG::Bool;   # false
  debugData::Dict;
end
FastCodingOpts(T) = FastCodingOpts{T}(zeros(0, 0), zeros(0, 0), zeros(0, 0), 1, 1, 1000, 1e-10, false, Dict());

type LearningData{T}
  Ztx::Array{T, 2};
  ZtZ::Array{T, 2};
end
type LearningOpts{T}
  Dinit::Array{T, 2};
  s::Array{T, 2};
  lambda_s::Array{T, 2};
  mu_s::T;       # 1
  max_iter::Int; # 1000
  tol::T;        # 1e-10
  DEBUG::Bool;   # false
  debugData::Dict;
end
LearningOpts(T) = LearningOpts{T}(zeros(0, 0), zeros(0, 0), zeros(0, 0), 1., 1000, 1e-10, false, Dict());

type FSCOpts{T}
  Dinit::Array{T, 2};
  beta::T;       # 1
  mu_s::T;       # 1
  mu_t::T;       # 1
  max_iter::Int; # 1000
  tol::T;        # 1e-10
  FIXED_D::Bool; # false
  DEBUG::Bool;   # false
  debugData::Dict;
end
FSCOpts(T) = FSCOpts{T}(zeros(0, 0), 1., 1., 1., 1000, 1e-10, false, false, Dict());

max_stars_repo_name: bkong/FSC

max_issues_repo_name: bkong/FSC

max_forks_repo_name: bkong/FSC


text:
"""
    VectorCacheModule

Module to manage vector caches.
"""
module VectorCacheModule

using ..FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict

"""
    VectorCache{T<:Number, F<:Function}

Type for caching vectors.

`T` = type of the entries of the vector,
`F` = type of the function to update the entries of the vector.

Signature of the function to fill the cache with the value of the vector
at any given point `XYZ`, using the columns of the Jacobian matrix of
the element, `tangents`, and, if convenient, also the finite element
label, `fe_label`. Finally, the value of the vector may also depend on
the `time` (or the load factor):
```
fillcache!(cacheout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt; time::FFlt = 0.0)
```
The cache `cacheout` is filled with the value  of the vector.
"""
struct VectorCache{T<:Number, F<:Function}
	# Function to update and retrieve the vector
	fillcache!::F
    # Cache where the current value of the vector can be retrieved
    cache::Vector{T}
	# Current time (or current load factor). Do not set directly. Use `settime!`
	_time::Vector{FFlt}
end

"""
    VectorCache(::Type{T}, nentries::FInt, fillcache!::F) where {T<:Number, F<:Function}

Construct vector cache. The function to fill the vector cache is given.

This constructor is intended for *time-independent* vector caches.
This function needs to have a signature of
```
fillcache!(cacheout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
    Calculate the vector and copy it into the cache....
    return forceout
end
```
and it needs to  fill in the cache `cacheout` with the current vector at
the location `XYZ`, using the information supplied in the Jacobian
matrix `tangents`, and the label of the finite element, `fe_label`, if
appropriate.
"""
function VectorCache(::Type{T}, nentries::FInt, fillcache!::F) where {T<:Number, F<:Function}
	function fillcachenotime!(cacheout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt; time::FFlt = 0.0)
        return fillcache!(cacheout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
    end
    # Allocate the cache to be ready for the first call
    return VectorCache(fillcachenotime!, zeros(T, nentries), [0.0]);
end

"""
    VectorCache(::Type{T}, nentries::FInt, fillcache!::F, time::FFlt) where {T<:Number, F<:Function}

Construct vector cache. The function to fill the vector cache is given.

This constructor is intended for *time-dependent* vector caches.
This function needs to have a signature of
```
fillcache!(cacheout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt; time::FFlt = 0.0)
    Calculate the vector and copy it into the cache....
    return forceout
end
```
and it needs to  fill in the cache `cacheout` with the current vector at
the location `XYZ`, using the information supplied in the Jacobian
matrix `tangents`, and the label of the finite element, `fe_label`, if
appropriate. The time can also be supplied (keyword argument `time`).
"""
function VectorCache(::Type{T}, nentries::FInt, fillcache!::F, time::FFlt) where {T<:Number, F<:Function}
	# Allocate the cache to be ready for the first call
    return VectorCache(fillcache!, zeros(T, nentries), [time]);
end

"""
    VectorCache(vector::FVec{T}) where {T<:Number}

Construct vector cache. The *constant* vector is given.
"""
function VectorCache(vector::FVec{T}) where {T<:Number}
    function fillcache!(cacheout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt; time::FFlt = 0.0)
        # do nothing:  the vector is already in the cache
        return cacheout
    end
    return VectorCache(fillcache!, deepcopy(vector), [0.0]);
end

"""
    updateretrieve!(self::VectorCache, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)

Update the cache and retrieve the vector.

After the return from this function the updated vector can be read from
the cache as `self.cache` (also returned). If the vector depends on
time, the vector cache time first needs to be set as
```
settime!(c, t)
```
"""
function updateretrieve!(self::VectorCache, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
    self.fillcache!(self.cache, XYZ, tangents, fe_label; time = self._time[1])
    return self.cache
end

"""
    settime!(self::VectorCache, time::FFlt)

Set the current time for the vector cache.
"""
function settime!(self::VectorCache, time::FFlt)
	self._time[1] = time
	return self
end

end


text:
facts("asset constructors work") do

    context("Stocks") do

        s  = Stock(Ticker("AAPL"), 100.00)
        ls = LongStock(Ticker("AAPL"), 100.00)
        ss = ShortStock(Ticker("AAPL"), 100.00)

        @fact typeof(s)  => Stock
        @fact typeof(ls) => LongStock
        @fact typeof(ss) => ShortStock

        @fact s.ticker.id  => "AAPL"
        @fact ls.ticker.id => "AAPL"
        @fact ss.ticker.id => "AAPL"

        @fact s.basis  => 100.00
        @fact ls.basis => 100.00
        @fact ss.basis => 100.00

        @fact s.shares  => 1
        @fact ls.shares => 1
        @fact ss.shares => 1

        @fact s.currency  => USD
        @fact ls.currency => USD
        @fact ss.currency => USD

        @fact s.tick  => .01
        @fact ls.tick => .01
        @fact ss.tick => .01

        @fact s.multiplier  => 1.0
        @fact ls.multiplier => 1.0
        @fact ss.multiplier => 1.0

        @fact typeof(s.id)  => Nullable{FinancialID}
        @fact typeof(ls.id) => Nullable{FinancialID}
        @fact typeof(ss.id) => Nullable{FinancialID}
    end

    context("Futures") do

        fut  = Future(Ticker("ES"), 2000.00)
        lfut = LongFuture(Ticker("ES"), 2000.00)
        sfut = ShortFuture(Ticker("ES"), 2000.00)

        @fact typeof(fut)  => Future
        @fact typeof(lfut) => LongFuture
        @fact typeof(sfut) => ShortFuture

        @fact fut.ticker.id  => "ES"
        @fact lfut.ticker.id => "ES"
        @fact sfut.ticker.id => "ES"

        @fact fut.basis  => 2000.00
        @fact lfut.basis => 2000.00
        @fact sfut.basis => 2000.00

        @fact fut.suffix  => ""
        @fact lfut.suffix => ""
        @fact sfut.suffix => ""

        @fact fut.contracts  => 1
        @fact lfut.contracts => 1
        @fact sfut.contracts => 1

        @fact fut.currency  => USD
        @fact lfut.currency => USD
        @fact sfut.currency => USD

        @fact typeof(fut.tick)  => Nullable{Float64}
        @fact typeof(lfut.tick) => Nullable{Float64}
        @fact typeof(sfut.tick) => Nullable{Float64}

        @fact typeof(fut.multiplier)  => Nullable{Float64}
        @fact typeof(lfut.multiplier) => Nullable{Float64}
        @fact typeof(sfut.multiplier) => Nullable{Float64}

        @fact typeof(fut.expiry)  => Nullable{Date}
        @fact typeof(lfut.expiry) => Nullable{Date}
        @fact typeof(sfut.expiry) => Nullable{Date}

        @fact typeof(fut.id)  => Nullable{FinancialID}
        @fact typeof(lfut.id) => Nullable{FinancialID}
        @fact typeof(sfut.id) => Nullable{FinancialID}
    end

    context("Options") do

        lp = LongPut(Ticker("SPX"), 1.45, 2000.00)
        lc = LongCall(Ticker("SPX"), 1.45, 2000.00)
        sp = ShortPut(Ticker("SPX"), 1.45, 2000.00)
        sc = ShortCall(Ticker("SPX"), 1.45, 2000.00)

        @fact typeof(lp) => LongPut
        @fact typeof(lc) => LongCall
        @fact typeof(sp) => ShortPut
        @fact typeof(sc) => ShortCall

        @fact lp.ticker.id => "SPX"
        @fact lc.ticker.id => "SPX"
        @fact sp.ticker.id => "SPX"
        @fact sc.ticker.id => "SPX"

        @fact lp.basis => 1.45
        @fact lc.basis => 1.45
        @fact sp.basis => 1.45
        @fact sc.basis => 1.45

        @fact lp.strike => 2000.00
        @fact lc.strike => 2000.00
        @fact sp.strike => 2000.00
        @fact sc.strike => 2000.00

        @fact typeof(lp.expiry) => Nullable{Date}
        @fact typeof(lc.expiry) => Nullable{Date}
        @fact typeof(sp.expiry) => Nullable{Date}
        @fact typeof(sc.expiry) => Nullable{Date}

        @fact lp.contracts => 1
        @fact lc.contracts => 1
        @fact sp.contracts => 1
        @fact sc.contracts => 1

        @fact lp.currency => USD
        @fact lc.currency => USD
        @fact sp.currency => USD
        @fact sc.currency => USD

        @fact typeof(lp.style) => Nullable{OptionExercise}
        @fact typeof(lc.style) => Nullable{OptionExercise}
        @fact typeof(sp.style) => Nullable{OptionExercise}
        @fact typeof(sc.style) => Nullable{OptionExercise}

        @fact typeof(lp.tick) => Nullable{Float64}
        @fact typeof(lc.tick) => Nullable{Float64}
        @fact typeof(sp.tick) => Nullable{Float64}
        @fact typeof(sc.tick) => Nullable{Float64}

        @fact typeof(lp.multiplier) => Nullable{Float64}
        @fact typeof(lc.multiplier) => Nullable{Float64}
        @fact typeof(sp.multiplier) => Nullable{Float64}
        @fact typeof(sc.multiplier) => Nullable{Float64}

        @fact typeof(lp.id) => Nullable{FinancialID}
        @fact typeof(lc.id) => Nullable{FinancialID}
        @fact typeof(sp.id) => Nullable{FinancialID}
        @fact typeof(sc.id) => Nullable{FinancialID}
    end

    context("Bond") do
    end

    context("CurrencyPair") do

        pair = CurrencyPair(EUR, USD)

        @fact pair.baseside          => EUR
        @fact pair.quoteside         => USD
        @fact typeof(pair.baseside)  => Currency
        @fact typeof(pair.quoteside) => Currency

    end
end


text:
# # Tutorial 2
#
# ## Transformations, extruded geometries, volumes
#
import gmsh
#
# If `ARGS` is passed to `gmsh.initialize()`, Gmsh will parse the command line
# in the same way as the standalone Gmsh app:
gmsh.initialize(append!(["gmsh"], ARGS))
#
gmsh.model.add("t2")
#
# Copied from Tutorial 1:
lc = 1e-2
gmsh.model.geo.addPoint(0, 0, 0, lc, 1)
gmsh.model.geo.addPoint(.1, 0,  0, lc, 2)
gmsh.model.geo.addPoint(.1, .3, 0, lc, 3)
gmsh.model.geo.addPoint(0, .3, 0, lc, 4)
gmsh.model.geo.addLine(1, 2, 1)
gmsh.model.geo.addLine(3, 2, 2)
gmsh.model.geo.addLine(3, 4, 3)
gmsh.model.geo.addLine(4, 1, 4)
gmsh.model.geo.addCurveLoop([4, 1, -2, 3], 1)
gmsh.model.geo.addPlaneSurface([1], 1)
gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(0, [1, 2], 1)
gmsh.model.addPhysicalGroup(1, [1, 2], 2)
gmsh.model.addPhysicalGroup(2, [1], 6)
gmsh.model.setPhysicalName(2, 6, "My surface")
#
# We can then add new points and curves in the same way as we did in `t1.jl`:
gmsh.model.geo.addPoint(0, .4, 0, lc, 5)
gmsh.model.geo.addLine(4, 5, 5)
#
#
# But Gmsh also provides tools to transform (translate, rotate, etc.)
# elementary entities or copies of elementary entities.  Geometrical
# transformations take a vector of pairs of integers as first argument, which
# contains the list of entities, represented by (dimension, tag) pairs.  For
# example, the point 5 (dimension=0, tag=5) can be moved by 0.02 to the left
# (dx=-0.02, dy=0, dz=0) with
gmsh.model.geo.translate([(0, 5)], -0.02, 0, 0)
# And it can be further rotated by -Pi/4 around (0, 0.3, 0) (with the rotation
# along the z axis) with:
gmsh.model.geo.rotate([(0, 5)], 0,0.3,0, 0,0,1, -pi/4)
#
#
# Note that there are no units in Gmsh: coordinates are just numbers - it's
# up to the user to associate a meaning to them.
#
# Point 3 can be duplicated and translated by 0.05 along the y axis by using the
# copy() function, which takes a vector of (dim, tag) pairs as input, and
# returns another vector of (dim, tag) pairs:
ov = gmsh.model.geo.copy([(0, 3)])
gmsh.model.geo.translate(ov, 0, 0.05, 0)
#
# The new point tag is available in `ov[0][1]`, and can be used to create new
# lines:
gmsh.model.geo.addLine(3, ov[1][2], 7)
gmsh.model.geo.addLine(ov[1][2], 5, 8)
gmsh.model.geo.addCurveLoop([5,-8,-7,3], 10)
gmsh.model.geo.addPlaneSurface([10], 11)
#
# In the same way, we can translate copies of the two surfaces 1 and 11 to the
# right with the following command:
ov = gmsh.model.geo.copy([(2, 1), (2, 11)])
gmsh.model.geo.translate(ov, 0.12, 0, 0)
#
println("New surfaces ", ov[1][2], " and ", ov[2][2])
#
# Volumes are the fourth type of elementary entities in Gmsh. In the same way
# one defines curve loops to build surfaces, one has to define surface loops
# (i.e. *shells*) to build volumes. The following volume does not have holes and
# thus consists of a single surface loop:
gmsh.model.geo.addPoint(0., 0.3, 0.12, lc, 100)
gmsh.model.geo.addPoint(0.1, 0.3, 0.12, lc, 101)
gmsh.model.geo.addPoint(0.1, 0.35, 0.12, lc, 102)
#
# We would like to retrieve the coordinates of point 5 to create point 103, so
# we synchronize the model, and use `getValue()`
gmsh.model.geo.synchronize()
xyz = gmsh.model.getValue(0, 5, [])
gmsh.model.geo.addPoint(xyz[1], xyz[2], 0.12, lc, 103)
#
gmsh.model.geo.addLine(4, 100, 110)
gmsh.model.geo.addLine(3, 101, 111)
gmsh.model.geo.addLine(6, 102, 112)
gmsh.model.geo.addLine(5, 103, 113)
gmsh.model.geo.addLine(103, 100, 114)
gmsh.model.geo.addLine(100, 101, 115)
gmsh.model.geo.addLine(101, 102, 116)
gmsh.model.geo.addLine(102, 103, 117)
#
gmsh.model.geo.addCurveLoop([115, -111, 3, 110], 118)
gmsh.model.geo.addPlaneSurface([118], 119)
gmsh.model.geo.addCurveLoop([111, 116, -112, -7], 120)
gmsh.model.geo.addPlaneSurface([120], 121)
gmsh.model.geo.addCurveLoop([112, 117, -113, -8], 122)
gmsh.model.geo.addPlaneSurface([122], 123)
gmsh.model.geo.addCurveLoop([114, -110, 5, 113], 124)
gmsh.model.geo.addPlaneSurface([124], 125)
gmsh.model.geo.addCurveLoop([115, 116, 117, 114], 126)
gmsh.model.geo.addPlaneSurface([126], 127)
#
gmsh.model.geo.addSurfaceLoop([127, 119, 121, 123, 125, 11], 128)
gmsh.model.geo.addVolume([128], 129)
#
# When a volume can be extruded from a surface, it is usually easier to use the
# `extrude()` function directly instead of creating all the points, curves and
# surfaces by hand. For example, the following command extrudes the surface 11
# along the z axis and automatically creates a new volume (as well as all the
# needed points, curves and surfaces). As expected, the function takes a vector
# of (dim, tag) pairs as input as well as the translation vector, and returns a
# vector of (dim, tag) pairs as output:
ov2 = gmsh.model.geo.extrude([ov[2]], 0, 0, 0.12)
#
#
# Mesh sizes associated to geometrical points can be set by passing a vector of
# (dim, tag) pairs for the corresponding points:
gmsh.model.geo.mesh.setSize([(0,103), (0,105), (0,109), (0,102), (0,28),
                             (0, 24), (0,6), (0,5)], lc * 3)
#
# We finish by synchronizing the data from the built-in CAD kernel with the Gmsh
# model:
gmsh.model.geo.synchronize()
#
# We group volumes 129 and 130 in a single physical group with tag `1` and name
# "The volume":
gmsh.model.addPhysicalGroup(3, [129,130], 1)
gmsh.model.setPhysicalName(3, 1, "The volume")
#
gmsh.model.mesh.generate(3)
#
# We finally generate and save the mesh:
gmsh.write("t2.msh")
#
# Note that, if the transformation tools are handy to create complex geometries,
# it is also sometimes useful to generate the *flat* geometry, with an explicit
# representation of all the elementary entities.
#
# With the built-in CAD kernel, this can be achieved by saving the model in the
# *Gmsh Unrolled GEO* format:
#
# `gmsh.write("t2.geo_unrolled");`
#
# With the OpenCASCADE CAD kernel, unrolling the geometry can be achieved by
# exporting in the `OpenCASCADE BRep' format:
#
# `gmsh.write("t2.brep");`
#
# (OpenCASCADE geometries can also be exported as STEP files.)
#
# It is important to note that Gmsh never translates geometry data into a common
# representation: all the operations on a geometrical entity are performed
# natively with the associated CAD kernel. Consequently, one cannot export a
# geometry constructed with the built-in kernel as an OpenCASCADE BRep file; or
# export an OpenCASCADE model as an Unrolled GEO file.
gmsh.finalize()


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMSparse.jl/blob/master/LICENSE

using Test

@testset "test FEMSparse.jl" begin
    @testset "test SparseMatrixCOO" begin
        include("test_sparsematrixcoo.jl")
    end
    @testset "test SparseMatrixCOO" begin
        include("test_sparsematrixcsc.jl")
    end
end


text:
@static if VERSION < v"0.7.0-DEV.481"
    const _ldiv! = isdefined(:ldiv!) ? ldiv! : Base.A_ldiv_B!
else
    const _ldiv! = (@isdefined ldiv!) ? ldiv! : Base.A_ldiv_B!
end
const _lu = @static VERSION < v"0.7.0-DEV.3449" ? Base.lufact : lu

let
    d  = 2 .+ rand(5)
    dl = rand(4)
    du = rand(4)
    M = Tridiagonal(dl, d, du)
    F = _lu(M)
    src = rand(5,5,5)
    for dim = 1:3
        dest1 = mapslices(x->_ldiv!(F, x), copy(src), dims=dim)
        dest2 = similar(src)
        AxisAlgorithms.A_ldiv_B_md!(dest2, F, src, dim)
        @test dest1 ≈ dest2
    end
end


text:
abstract type AbstractTransform end

struct BoxCox <: AbstractTransform end

transform(::BoxCox, λ, y) = λ[1] == 0 ? log(y) : expm1(log(y) * λ[1]) / λ[1]
log∂λ(::BoxCox, λ, y) = (λ[1] - 1) * log(y)

struct IdentityTransform <: AbstractTransform end

transform(::IdentityTransform, _, y) = y
log∂(::IdentityTransform, _, _) = 0

struct YeoJohnson <: AbstractTransform end

function transform(::YeoJohnson, λ, y)
    λ′ = y >= 0 ? λ[1] : 2 - λ[1]
    return sign(y) * (λ′ == 0 ? log1p(abs(y)) : expm1(log1p(abs(y)) * λ′) / λ′)
end
log∂(::YeoJohnson, λ, y) = ((y >= 0 ? λ[1] : 2 - λ[1]) - 1) *  log(abs(y))

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
module OneDimensionalUtils

using PyPlot

using Oceananigans.Grids

defaultcolors = plt.rcParams["axes.prop_cycle"].by_key()["color"]
removespine(side) = gca().spines[side].set_visible(false)
removespines(sides...) = [removespine(side) for side in sides]

""" Unpack a vector of results associated with a convergence test. """
function unpack_solutions(results)
    c_ana = map(r -> r.cx.analytical[:], results)
    c_sim = map(r -> r.cx.simulation[:], results)
    return c_ana, c_sim
end

function unpack_errors(results)
     u_L₁ = map(r -> r.u.L₁, results)
     v_L₁ = map(r -> r.v.L₁, results)
    cx_L₁ = map(r -> r.cx.L₁, results)
    cy_L₁ = map(r -> r.cy.L₁, results)

     u_L∞ = map(r ->  r.u.L∞, results)
     v_L∞ = map(r ->  r.v.L∞, results)
    cx_L∞ = map(r -> r.cx.L∞, results)
    cy_L∞ = map(r -> r.cy.L∞, results)

    return (u_L₁, v_L₁, cx_L₁, cy_L₁,
            u_L∞, v_L∞, cx_L∞, cy_L∞)
end

""" Unpack a vector of grids associated with a convergence test. """
unpack_grids(results) = map(r -> r.grid, results)

function plot_solutions!(axs, all_results, names, linestyles, specialcolors)

    for j = 1:length(all_results)

        results = all_results[j]
        name = names[j]
        linestyle = linestyles[j]

        c_ana, c_sim = unpack_solutions(results)

        grids = unpack_grids(results)

        x = xnodes(Cell, grids[end])[:]

        sca(axs[1])

        plot(x, c_ana[end], "-", color=specialcolors[j], alpha=0.2, linewidth=3, label="$name analytical")

        for i in 1:length(c_sim)
            x = xnodes(Cell, grids[i])[:]

            plot(x, c_sim[i], linestyle, color=defaultcolors[i], alpha=0.8, linewidth=1,
                 label="$name simulated, \$ N_x \$ = $(grids[i].Nx)")

            xlim(minimum(x), maximum(x))
        end

        sca(axs[2])

        for i in 1:length(c_sim)
            x = xnodes(Cell, grids[i])[:]

            semilogy(x, abs.(c_sim[i] .- c_ana[i]), linestyle, color=defaultcolors[i], alpha=0.8,
                     label="$name, \$ N_x \$ = $(grids[i].Nx)")

            xlim(minimum(x), maximum(x))
        end

    end

    sca(axs[1])
    ylabel(L"c")
    removespines("top", "right", "bottom")
    axs[1].tick_params(bottom=false, labelbottom=false)
    lgd1 = legend(loc="upper left", prop=Dict(:size=>7), bbox_to_anchor=(-0.35, 0, 1.5, 1.0))

    sca(axs[2])
    removespines("top", "right")
    xlabel(L"x")
    ylabel("Absolute error \$ | c_\\mathrm{sim} - c_\\mathrm{analytical} | \$")
    lgd2 = legend(loc="upper right", prop=Dict(:size=>10), bbox_to_anchor=(-0.2, 0.4, 1.5, 1.0))

    return (lgd1, lgd2)
end

function plot_error_convergence!(axs, Nx, all_results, names)

    for j = 1:length(all_results)
        results = all_results[j]
        name = names[j]
        u_L₁, v_L₁, cx_L₁, cy_L₁, u_L∞, v_L∞, cx_L∞, cy_L∞  = unpack_errors(results)

        common_kwargs = (linestyle="None", color=defaultcolors[j], mfc="None", alpha=0.8)
        loglog(Nx,  u_L₁; basex=2, marker="o", label="\$L_1\$-norm, \$u\$ $name", common_kwargs...)
        loglog(Nx,  v_L₁; basex=2, marker="2", label="\$L_1\$-norm, \$v\$ $name", common_kwargs...)
        loglog(Nx, cx_L₁; basex=2, marker="*", label="\$L_1\$-norm, \$x\$ tracer $name", common_kwargs...)
        loglog(Nx, cy_L₁; basex=2, marker="+", label="\$L_1\$-norm, \$y\$ tracer $name", common_kwargs...)

        loglog(Nx,  u_L∞; basex=2, marker="1", label="\$L_\\infty\$-norm, \$u\$ $name", common_kwargs...)
        loglog(Nx,  v_L∞; basex=2, marker="_", label="\$L_\\infty\$-norm, \$v\$ $name", common_kwargs...)
        loglog(Nx, cx_L∞; basex=2, marker="^", label="\$L_\\infty\$-norm, \$x\$ tracer $name", common_kwargs...)
        loglog(Nx, cy_L∞; basex=2, marker="s", label="\$L_\\infty\$-norm, \$y\$ tracer $name", common_kwargs...)
    end

    # Guide line to confirm second-order scaling
    u_L₁, v_L₁, cx_L₁, cy_L₁, u_L∞, v_L∞, cx_L∞, cy_L∞  = unpack_errors(all_results[1])
    loglog(Nx, cx_L₁[1] .* (Nx[1] ./ Nx).^2, "k-", basex=2, alpha=0.8, label=L"\sim N_x^{-2}")

    xlabel(L"N_x")
    ylabel("\$L\$-norms of \$ | c_\\mathrm{sim} - c_\\mathrm{analytical} |\$")
    removespines("top", "right")
    lgd = legend(loc="upper right", bbox_to_anchor=(1.4, 1.0), prop=Dict(:size=>6))

    return lgd
end

end # module


text:
"""
      parsimony(tree::N, char::Dict{String,String}, gap::String="-")::Float64 where N<: GeneralNode

Do parsimony reconstruction for a tree and a set of characters. The characters are
supplied as the char dictionary with the key being the name of a leave and the value
string representation of the character.
"""
function parsimony(tree::N, char::Dict{String,String}, gap::String="-")::Float64 where N<: GeneralNode
      po = post_order(tree)
      states = filter(x -> x != "-", unique(collect(values(char))))
      nStates = length(states)
      if nStates == 0
            return 0.0
      end
      nNodes = length(po)
      costMatrix = zeros(nNodes, nStates) .- 1
      for node in po
            if node.nchild == 0
                  # i am a leave node
                  s = char[node.name]
                  costMatrix[node.num, :] .= s == gap ? 0 : abs.(log.(states .== s))
            else
                  # i am an internal node
                  daughters = node.children
                  localCosts = zeros(nStates, nStates, length(daughters))
                  for (j, s1) in enumerate(states),
                        (k, s2) in enumerate(states),
                        (l, d) in enumerate(daughters)

                        localCosts[j, k, l] = costMatrix[d.num, k] + Int(s1 != s2)
                  end
                  costMatrix[i, :] = vec(
                        mapslices(
                              sum,
                              mapslices(minimum, localCosts, dims = 2),
                              dims = 3,
                        ),
                  )
            end
      end
      minimum(costMatrix[tree.num, :])
end


text:
#####################################################################
# simulate.jl
#
# TODO: eliminate steps and make simulate with a simunit (maybe)
######################################################################

# sample an observation
#  1. receives an observation
#  2. updates belief (contained in filter)
#  3. updates vehicle position according to some policy
# For right now, assume a random policy
function step!(m::SearchDomain, x::Vehicle, f::AbstractFilter, p::Policy)
	# observe, update belief, select action, and act
	o = observe(m,x)
	update!(f, x, o)
	a = action(m, x, o, f, p)
	act!(m,x,a)
end

function simulate(m::SearchDomain,
                  x::Vehicle,
                  f::AbstractFilter,
                  p::Policy,
                  n_steps::Int=10
                 )

    tc = StepThreshold(n_steps)

	# Then go through steps
	step_count = 0
	while !is_complete(f, tc, step_count)
		step!(m,x,f,p)
		step_count += 1
	end
end

function simulate(m::SearchDomain, su::SimUnit)

    # Don't reset the filter, vehicle, and policy
    # I think I assume the SimUnit comes in clean and ready to go

    # What was the cost to getting this first observation?
    cost_sum = get_cost(su, m)

    # before doing anything else, we observe
    #  and update filter once
    o = observe(m, su.x)
    update!(su, o)

    # This was our first step; steps count number of observations
    step_count = 1

    while !is_complete(su.f, su.tc, step_count)
        # act
        a = action(m, su, o)
        act!(m, su.x, a)

        move_target!(m)

        # get cost and update step count
        cost_sum += get_cost(su, m, a)
        step_count += 1

        # observe and update
        o = observe(m, su.x)
        update!(su, o)
    end

    return cost_sum
end

# for running batches of simulations on a single sim unit
function simulate(m::SearchDomain, su::SimUnit, ns::Int)
    return vec( simulate(m, [su], ns) )
end

# for running batches of simulations on a vector of sim units
function simulate(m::SearchDomain, suv::Vector{SimUnit}, n_sims::Int;
                  random_theta::Bool=true
                )

    theta_start = m.theta

    costs = zeros(n_sims, length(suv))
    for sim_ind = 1:n_sims

        # reset jammer to random location
        m.theta = theta_start
        if random_theta
            theta!(m)
        end

        print("Simulation ", sim_ind, ": ")
        for (su_ind, su) in enumerate(suv)
            print(su_ind, ",")

            # reset the filter, vehicle, and policy
            reset!(m, su)

            costs[sim_ind, su_ind] = simulate(m, su)
        end
        println("complete.")
    end
    return costs
end


text:
using Reduce
using ReduceLinAlg
using Test

# write your own tests here
@test mat_jacobian((:x,),(:x,)) == mat_jacobian([:x],[:x])
@test hessian(:x,(:x,)) == hessian(:x,[:x])
@test jordan_block(:x,1) == jordan_block(Reduce.RExpr("x"),Reduce.RExpr("1")) |> parse