text:
using Test, QBase

@testset "./src/constructors/brakets.jl" begin

@testset "bloch_qubit_ket()" begin
    @testset "complex bloch qubit" begin
        ψ1 = bloch_qubit_ket(π/2,π/2)

        @test ψ1 isa Ket{Complex{Float64}}
        @test ψ1 ≈ [1,im]/sqrt(2)

        ψ2 = bloch_qubit_ket(π/2,3π/2)

        @test ψ2 isa Ket{Complex{Float64}}
        @test ψ2 ≈ [1,-im]/sqrt(2)

        @test_throws DomainError bloch_qubit_ket(-0.1,π/2)
        @test_throws DomainError bloch_qubit_ket(15π/7,π/2)
        @test_throws DomainError bloch_qubit_ket(π/2,-0.1)
        @test_throws DomainError bloch_qubit_ket(π/2,15π/7)
    end

    @testset "real bloch qubit" begin
        ψ1 = bloch_qubit_ket(π/2)

        @test ψ1 isa Ket{Float64}
        @test ψ1 ≈ [1,1]/sqrt(2)

        ψ2 = bloch_qubit_ket(-π/2)

        @test ψ2 isa Ket{Float64}
        @test ψ2 ≈ [1,-1]/sqrt(2)

        @test_throws DomainError bloch_qubit_ket(-2π - 0.1)
        @test_throws DomainError bloch_qubit_ket(2π + 0.1)
    end

    @testset "equivalence between real and complex versions" for θ in 0:0.2:2π
        @test bloch_qubit_ket(θ) ≈ bloch_qubit_ket(θ,0)
        @test isapprox(bloch_qubit_ket(-θ), bloch_qubit_ket(θ,π), atol=1e-7)
    end
end

@testset "computational_basis_kets()" begin
    comp_basis_kets = computational_basis_kets(3)
    @test comp_basis_kets isa Vector{Ket{Int64}}
    @test comp_basis_kets == [[1,0,0],[0,1,0],[0,0,1]]
    @test computational_basis_kets(8) == [
        [1,0,0,0,0,0,0,0],[0,1,0,0,0,0,0,0],[0,0,1,0,0,0,0,0],[0,0,0,1,0,0,0,0],
        [0,0,0,0,1,0,0,0],[0,0,0,0,0,1,0,0],[0,0,0,0,0,0,1,0],[0,0,0,0,0,0,0,1]
    ]
end

@testset "mirror_symmetric_qubit_kets()" begin
    @test mirror_symmetric_qubit_kets(0) isa Vector{Ket{Float64}}

    @test isapprox(mirror_symmetric_qubit_kets(0), [[1,0],[1,0],[1,0]])
    @test isapprox(mirror_symmetric_qubit_kets(π/4), [[1,0],[1,1]/sqrt(2),[1,-1]/sqrt(2)])
    @test isapprox(mirror_symmetric_qubit_kets(π/3),[[1,0],[0.5,(sqrt(3)/2)],[0.5,-sqrt(3)/2]])
    @test isapprox(mirror_symmetric_qubit_kets(π/2), [[1,0],[0,1],[0,-1]], atol=1e-16)

    @test_throws DomainError mirror_symmetric_qubit_kets(-0.1)
    @test_throws DomainError mirror_symmetric_qubit_kets(3π/5)
end

@testset "planar_symmetric_qubit_kets()" begin
    @testset "orthogonal qubits" begin
        kets = planar_symmetric_qubit_kets(2)
        @test kets isa Vector{Ket{Float64}}
        @test isapprox(kets, [[1,0],[0,1]], atol=1e-7)
    end

    @testset "scanning over cases verifying right number are computed" begin
        for n in 2:100
            kets = planar_symmetric_qubit_kets(n)

            @test length(kets) == n
            @test kets isa Vector{Ket{Float64}}
        end
    end

    @test_throws(DomainError, planar_symmetric_qubit_kets(1))
end

@testset "trine_qubit_kets()" begin
    trine_kets = trine_qubit_kets()
    @test trine_kets isa Vector{Ket{Float64}}
    @test isapprox(trine_kets, mirror_symmetric_qubit_kets(π/3), atol=1e-7)
    @test isapprox(trine_kets, planar_symmetric_qubit_kets(3), atol=1e-7)
end

@testset "sic_qubit_kets()" begin
    sic_kets = sic_qubit_kets()
    @test sic_kets isa Vector{Ket{Complex{Float64}}}

    θ_sic = 2*acos(1/sqrt(3))

    @test sic_kets[1] ≈ bloch_qubit_ket(0, 0)
    @test sic_kets[2] ≈ bloch_qubit_ket(θ_sic, 0)
    @test sic_kets[3] ≈ bloch_qubit_ket(θ_sic, 2π/3)
    @test sic_kets[4] ≈ bloch_qubit_ket(θ_sic, 4π/3)
end

@testset "bb84_qubit_kets" begin
    bb84_kets = bb84_qubit_kets()
    @test bb84_kets isa Vector{Ket{Float64}}
    @test isapprox(bb84_kets, planar_symmetric_qubit_kets(4), atol=1e-7)
end

@testset "bell_kets()" begin
    kets = bell_kets()
    k0 = [1,0]
    k1 = [0,1]

    @test kets isa Vector{Ket{Float64}}
    @test length(kets) == 4
    @test kets[1] == 1/sqrt(2)*(kron(k0,k0)+kron(k1,k1))
    @test kets[2] == 1/sqrt(2)*(kron(k0,k0)-kron(k1,k1))
    @test kets[3] == 1/sqrt(2)*(kron(k0,k1)+kron(k1,k0))
    @test kets[4] == 1/sqrt(2)*(kron(k0,k1)-kron(k1,k0))
end

@testset "generalized_bell_kets()" begin
    @testset "qubit bell states" begin
        kets = generalized_bell_kets(2)

        @test kets isa Vector{Ket{Complex{Float64}}}
        @test isapprox(kets, bell_kets(), atol=1e-7)
    end

    @testset "qutrit bell states" begin
        kets = generalized_bell_kets(3)

        basis = computational_basis_kets(3)
        ψ1 = (kron(basis[1],basis[1]) + kron(basis[2],basis[2]) + kron(basis[3],basis[3]))/sqrt(3)
        ψ2 = (kron(basis[1],basis[1]) + exp(im*2*π/3)*kron(basis[2],basis[2]) + exp(im*4*π/3)*kron(basis[3],basis[3]))/sqrt(3)
        ψ3 = (kron(basis[1],basis[1]) + exp(im*4*π/3)*kron(basis[2],basis[2]) + exp(im*2*π/3)*kron(basis[3],basis[3]))/sqrt(3)

        ψ4 = (kron(basis[1],basis[2]) + kron(basis[2],basis[3]) + kron(basis[3],basis[1]))/sqrt(3)
        ψ5 = (kron(basis[1],basis[2]) + exp(im*2*π/3)*kron(basis[2],basis[3]) + exp(im*4*π/3)*kron(basis[3],basis[1]))/sqrt(3)
        ψ6 = (kron(basis[1],basis[2]) + exp(im*4*π/3)*kron(basis[2],basis[3]) + exp(im*2*π/3)*kron(basis[3],basis[1]))/sqrt(3)

        ψ7 = (kron(basis[1],basis[3]) + kron(basis[2],basis[1]) + kron(basis[3],basis[2]))/sqrt(3)
        ψ8 = (kron(basis[1],basis[3]) + exp(im*2*π/3)*kron(basis[2],basis[1]) + exp(im*4*π/3)*kron(basis[3],basis[2]))/sqrt(3)
        ψ9 = (kron(basis[1],basis[3]) + exp(im*4*π/3)*kron(basis[2],basis[1]) + exp(im*2*π/3)*kron(basis[3],basis[2]))/sqrt(3)

        ψ_match = [ψ1,ψ2,ψ3,ψ4,ψ5,ψ6,ψ7,ψ8,ψ9]

        @test kets isa Vector{Ket{Complex{Float64}}}
        @test kets ≈ ψ_match
    end
end

end


text:
"""
    isavailable(filter_or_id)

Given a subtype of `Filters.Filter` or the filter ID number as an integer,
return `true` if the filter is available and `false` otherwise.
"""
isavailable(filter_or_id) = API.h5z_filter_avail(filter_or_id)

"""
    isencoderenabled(filter_or_id)

Given a subtype of `Filters.Filter` or the filter ID number as an integer,
return `true` if the filter can encode or compress data.
"""
function isencoderenabled(filter_or_id)
    info = API.h5z_get_filter_info(filter_or_id)
    return info & API.H5Z_FILTER_CONFIG_ENCODE_ENABLED != 0
end

"""
    isdecoderenabled(filter_or_id)

Given a subtype of `Filters.Filter` or the filter ID number as an integer,
return `true` if the filter can decode or decompress data.
"""
function isdecoderenabled(filter_or_id)
    info = API.h5z_get_filter_info(filter_or_id)
    return info & API.H5Z_FILTER_CONFIG_DECODE_ENABLED != 0
end

ext: jl


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

mutable struct NextHopResult <: SwaggerModel
    nextHopType::Any # spec type: Union{ Nothing, String } # spec name: nextHopType
    nextHopIpAddress::Any # spec type: Union{ Nothing, String } # spec name: nextHopIpAddress
    routeTableId::Any # spec type: Union{ Nothing, String } # spec name: routeTableId

    function NextHopResult(;nextHopType=nothing, nextHopIpAddress=nothing, routeTableId=nothing)
        o = new()
        validate_property(NextHopResult, Symbol("nextHopType"), nextHopType)
        setfield!(o, Symbol("nextHopType"), nextHopType)
        validate_property(NextHopResult, Symbol("nextHopIpAddress"), nextHopIpAddress)
        setfield!(o, Symbol("nextHopIpAddress"), nextHopIpAddress)
        validate_property(NextHopResult, Symbol("routeTableId"), routeTableId)
        setfield!(o, Symbol("routeTableId"), routeTableId)
        o
    end
end # type NextHopResult

const _property_map_NextHopResult = Dict{Symbol,Symbol}(Symbol("nextHopType")=>Symbol("nextHopType"), Symbol("nextHopIpAddress")=>Symbol("nextHopIpAddress"), Symbol("routeTableId")=>Symbol("routeTableId"))
const _property_types_NextHopResult = Dict{Symbol,String}(Symbol("nextHopType")=>"String", Symbol("nextHopIpAddress")=>"String", Symbol("routeTableId")=>"String")
Base.propertynames(::Type{ NextHopResult }) = collect(keys(_property_map_NextHopResult))
Swagger.property_type(::Type{ NextHopResult }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_NextHopResult[name]))}
Swagger.field_name(::Type{ NextHopResult }, property_name::Symbol) =  _property_map_NextHopResult[property_name]

const _allowed_NextHopResult_nextHopType = ["Internet", "VirtualAppliance", "VirtualNetworkGateway", "VnetLocal", "HyperNetGateway", "None"]

function check_required(o::NextHopResult)
    true
end

function validate_property(::Type{ NextHopResult }, name::Symbol, val)
    if name === Symbol("nextHopType")
        Swagger.validate_param(name, "NextHopResult", :enum, val, _allowed_NextHopResult_nextHopType)
    end
end


text:
# In this file we define special-cases to prevent Cassette related inference issues

const BLACK_LIST = (
    Base.promote_op, Base.to_shape,
    Core.getfield,
    Core.:(===),
    Base.iterate,
    Broadcast.broadcasted,
    Broadcast.preprocess,
    Broadcast.combine_axes,
    Base.not_int,
    Base.size,
    Base.haskey,
    Base.reduced_indices,
    LinearAlgebra.gemv!,
    Tuple,
)

for F in BLACK_LIST
    @eval @inline Cassette.overdub(ctx::RecordingCtx, f::typeof($F), xs...) = f(xs...)
    @eval @inline Cassette.overdub(ctx::ReplayCtx, f::typeof($F), xs...) = f(xs...)
end

@inline Cassette.overdub(ctx::RecordingCtx, ::Type{Val}, x) = Val(x)
@inline Cassette.overdub(ctx::ReplayCtx, ::Type{Val}, x) = Val(x)

@inline Cassette.overdub(ctx::RecordingCtx, ::typeof(getindex), x::IdDict, key) = getindex(x, key)
@inline Cassette.overdub(ctx::ReplayCtx, ::typeof(getindex), x::IdDict, key) = getindex(x, key)


text:
# Automatically generated using Clang.jl

const Clp_Simplex = Cvoid
const clp_callback = Ptr{Cvoid}
const Sbb_Model = Cvoid
const Cbc_Model = Cvoid
const sbb_callback = Ptr{Cvoid}
const cbc_callback = Ptr{Cvoid}
const cbc_cut_callback = Ptr{Cvoid}
const CoinBigIndex = Cint


text:
"""
    GPA(score::Int)::Float64
北大 GPA 计算公式
"""
function GPA(score::Int)::Float64
    score < 0 || score > 100 && error("不存在的分数")
    score < 60 && return 0
    return 4 - 3 * (100 - score) * (100 - score) / 1600
end

const id_weight = [7, 9, 10, 5, 8, 4, 2, 1, 6, 3, 7, 9, 10, 5, 8, 4, 2]
const id_check_code = "10X98765432"

"""
    get_id_check_code(id17::AbstractString)::Char
根据前面17位身份证数字计算校验码
"""
function get_id_check_code(id17::AbstractString)::Char
    length(id17) < 17 && error("至少要给身份证前17位")
    sum = 0
    for i = 1:17
        sum += (id17[i] - '0') * id_weight[i]
    end
    return id_check_code[mod(sum, 11) + 1]
end

"""
    check_idcard(id::AbstractString)::Bool
检验身份证校验码是否正确
"""
function check_idcard(id::AbstractString)::Bool
    length(id) != 18 && error("请使用18位身份证")
    return get_id_check_code(id) == id[18]
end


text:
using StaticArrays

# Frequently-used array types
const Mat3{T} = SMatrix{3, 3, T, 9} where T
const Vec3{T} = SVector{3, T} where T
const AbstractArray3{T} = AbstractArray{T, 3}


text:
# Convenience functions
using PyCall
using DocStringExtensions
using Pipe
using CSV
py"""
import numpy as np
"""

add_dim(x::Array) = reshape(x, (size(x)...,1))

"""
drop missing values in pipe
$(SIGNATURES)
"""
function drop_na(df::DataFrame)
    return df[completecases(df),:]
end

function drop_na(ar::Array{Union{Missing,Float64},1})
    mask1 = .!ismissing.(ar)
    mask2 = isfinite.(ar)
    return ar[mask1.&mask2]
end

"""
Enumerate unique values
"""
function enique(x)
    return enumerate(unique(x))
end

"""
Read observed data
$(SIGNATURES)
"""
function read_observed(path::String; timevar=:tday, yvar=:count)
    observed = @pipe CSV.read(path, DataFrame) |> leftjoin(_, combine(groupby(_[_.treatment.=="Co",:], timevar), yvar => mean), on=timevar)
    observed.response = observed[:,yvar] ./ observed[:,Symbol(string(yvar)*"_mean")]
    observed.tday = convert.(Float64, observed.tday)
    return observed
end

"""
Check if timepoint in timepoints
"""
function t_in_timepoints(x)
    return x in timepoints
end

"""
Plot population trajectories in different treatments, given predicted and observed data.
$(SIGNATURES)
"""
function plot_all(predicted::DataFrame, observed::DataFrame; alpha=.2, co_mean=nothing)
    trts = unique(predicted.treatment)
    fig, ax = plt.subplots(ncols=size(trts)[1], figsize=(20, 4), sharey=true)
    for (i, trt) in enumerate(trts)
        m = predicted[predicted.treatment.==trt,:]
        d = dropmissing(observed[observed.treatment.==trt,:])
        sns.lineplot(m.tday, m.N, alpha=alpha, units=m.replicate, estimator=nothing, ax=ax[i])
        sns.lineplot(d.tday, d.count, ax=ax[i])
        if co_mean != nothing
            ax[i].axhline(co_mean)
        end
        ax[i].set(title=trt, xlabel="time (days)")
        if i==1
            ax[i].set_ylabel("# Daphnids")
        end
    end
    plt.tight_layout()
    sns.despine()
end

function plot_pdf(dist, range)
    p = [pdf(dist, x) for x in range]
    sns.lineplot(range, p)
    ax.set(xlabel="value", ylabel="probability density")
    ax.get_yaxis().set_visible(false)
    sns.despine()
end

"""
Plot population trajectories in different treatments, given predicted and observed data.
$(SIGNATURES)
"""
function nplot(predicted::DataFrame, observed::DataFrame; ci=95)

    trts = unique(predicted.treatment)
    fig, ax = plt.subplots(ncols=size(trts)[1], figsize=(20, 4), sharey=true)

    for (i, trt) in enumerate(trts)
        m = predicted[predicted.treatment.==trt,:]

        pest = by(m, :tday) do df
            est, lower, upper = get_estimate(df.N)
            DataFrame(
                est = est,
                lower = lower,
                upper = upper
            )
        end

        d = dropmissing(observed[observed.treatment.==trt,:])
        sns.scatterplot(d.tday, d.ind_per_L, ax=ax[i])

        sns.lineplot(pest.tday, pest.est, color="gray", estimator=nothing, ax=ax[i])
        sns.lineplot(pest.tday, pest.lower, color="lightgray", alpha=.5, ax=ax[i])
        sns.lineplot(pest.tday, pest.upper, color="lightgray", alpha=.5, ax=ax[i])

        ax[i].set(title=trt, xlabel="time (days)")

        if i==1
            ax[i].set_ylabel("# Daphnids")
        end

    end

    plt.tight_layout()
    sns.despine()
end

# functions for log(x+1)-transformed axes
logx = (py"lambda x:np.log(x+1)", py"lambda x:np.exp(x)-1")


text:
module BinaryPlatforms

using Base.BinaryPlatforms
export platform_key_abi, platform_dlext, valid_dl_path, arch, libc,
       libgfortran_version, libstdcxx_version, cxxstring_abi, parse_dl_name_version,
       detect_libgfortran_version, detect_libstdcxx_version, detect_cxxstring_abi,
       call_abi, wordsize, triplet, select_platform, platforms_match,
       CompilerABI, Platform, UnknownPlatform, Linux, MacOS, Windows, FreeBSD

import Base.BinaryPlatforms: libgfortran_version, libstdcxx_version, platform_name,
                             wordsize, platform_dlext, tags, arch, libc, call_abi,
                             cxxstring_abi

struct UnknownPlatform <: AbstractPlatform
    UnknownPlatform(args...; kwargs...) = new()
end
tags(::UnknownPlatform) = Dict{String,String}("os"=>"unknown")

struct CompilerABI
    libgfortran_version::Union{Nothing,VersionNumber}
    libstdcxx_version::Union{Nothing,VersionNumber}
    cxxstring_abi::Union{Nothing,Symbol}

    function CompilerABI(;libgfortran_version::Union{Nothing, VersionNumber} = nothing,
                         libstdcxx_version::Union{Nothing, VersionNumber} = nothing,
                         cxxstring_abi::Union{Nothing, Symbol} = nothing)
        return new(libgfortran_version, libstdcxx_version, cxxstring_abi)
    end
end

# Easy replacement constructor
function CompilerABI(cabi::CompilerABI; libgfortran_version=nothing,
                                        libstdcxx_version=nothing,
                                        cxxstring_abi=nothing)
    return CompilerABI(;
        libgfortran_version=something(libgfortran_version, Some(cabi.libgfortran_version)),
        libstdcxx_version=something(libstdcxx_version, Some(cabi.libstdcxx_version)),
        cxxstring_abi=something(cxxstring_abi, Some(cabi.cxxstring_abi)),
    )
end

libgfortran_version(cabi::CompilerABI) = cabi.libgfortran_version
libstdcxx_version(cabi::CompilerABI) = cabi.libstdcxx_version
cxxstring_abi(cabi::CompilerABI) = cabi.cxxstring_abi

for T in (:Linux, :Windows, :MacOS, :FreeBSD)
    @eval begin
        struct $(T) <: AbstractPlatform
            p::Platform
            function $(T)(arch::Symbol; compiler_abi=nothing, kwargs...)
                if compiler_abi !== nothing
                    kwargs = (; kwargs...,
                        :libgfortran_version => libgfortran_version(compiler_abi),
                        :libstdcxx_version => libstdcxx_version(compiler_abi),
                        :cxxstring_abi => cxxstring_abi(compiler_abi)
                    )
                end
                return new(Platform(string(arch), $(string(T)); kwargs..., validate_strict=true))
            end
        end
    end

    # First, methods we need to coerce to Symbol for backwards-compatibility
    for f in (:arch, :libc, :call_abi, :cxxstring_abi)
        @eval begin
            function $(f)(p::$(T))
                str = $(f)(p.p)
                if str === nothing
                    return nothing
                end
                return Symbol(str)
            end
        end
    end

    # Next, things we don't need to coerce
    for f in (:libgfortran_version, :libstdcxx_version, :platform_name, :wordsize, :platform_dlext, :tags, :triplet)
        @eval begin
            $(f)(p::$(T)) = $(f)(p.p)
        end
    end

    # Finally, add equality testing between these wrapper types and other AbstractPlatforms
    @eval begin
        Base.:(==)(a::$(T), b::AbstractPlatform) = b == a.p
    end
end

# Add one-off functions
MacOS(; kwargs...) = MacOS(:x86_64; kwargs...)
FreeBSD(; kwargs...) = FreeBSD(:x86_64; kwargs...)

function triplet(p::AbstractPlatform)
    # We are going to sub off to `Base.BinaryPlatforms.triplet()` here,
    # with the important exception that we override `os_version` to better
    # mimic the old behavior of `triplet()`
    if Sys.isfreebsd(p)
        p = deepcopy(p)
        p["os_version"] = "11.1.0"
    elseif Sys.isapple(p)
        p = deepcopy(p)
        p["os_version"] = "14.0.0"
    end
    return Base.BinaryPlatforms.triplet(p)
end

"""
    platform_key_abi(machine::AbstractString)

Returns the platform key for the current platform, or any other though the
the use of the `machine` parameter.

This method is deprecated, import `Base.BinaryPlatforms` and use either `HostPlatform()`
to get the current host platform, or `parse(Base.BinaryPlatforms.Platform, triplet)`
to parse the triplet for some other platform instead.
"""
platform_key_abi() = HostPlatform()
platform_key_abi(triplet::AbstractString) = parse(Platform, triplet)

"""
    valid_dl_path(path::AbstractString, platform::Platform)

Return `true` if the given `path` ends in a valid dynamic library filename.
E.g. returns `true` for a path like `"usr/lib/libfoo.so.3.5"`, but returns
`false` for a path like `"libbar.so.f.a"`.

This method is deprecated and will be removed in Julia 2.0.
"""
function valid_dl_path(path::AbstractString, platform::AbstractPlatform)
    try
        parse_dl_name_version(path, string(os(platform))::String)
        return true
    catch e
        if isa(e, ArgumentError)
            return false
        end
        rethrow(e)
    end
end

end # module BinaryPlatforms


text:
# Proportional selection using the stohcastic acceptance algorithm of Lipowski and Lipowska
# Fitnesses are defined in a variant table.
# This version came from evotech/spatial and was the default for continous variable evolution.
#using DataStructures
export propsel, propsel!, reverse_propsel

@doc """ function propsel()
Apply proportional selection to Population pop to generate a new population of size n.
Assumes that all elements (integers) in pop have their fitnesses defined in variant_table.
"""
function propsel( pop::Population, n::Int64, variant_table::Dict{Int64,variant_type} )
  new_pop = zeros(Int64,n)
  fmax = 0.0
  for v in pop
    fitv = variant_table[v].fitness
    if fitv < 0.0
      error("negative fitness in propsel")
    end
    if fitv > fmax
      fmax = fitv
    end
  end
  if fmax == 0.0
    # all elements have fitness zero
    println("Warning: all elements have fitness zero")
    return pop
  end

  N = length(pop)
  k = 0
  while k < n
    i = rand(1:N)
    w = variant_table[pop[i]].fitness/fmax
    if rand() < w
      new_pop[k + 1] = pop[i]
      k += 1
    end
  end
  new_pop
end

@doc """function reverse_propsel()
Apply reverse proportional selection to pop to return the population without the selected elements.
Reverse proportional selection uses rescaled netative fitness to favor lower fitness individuals.
The rescaling is  fit -> -fit + fmax + fmin.
Note that under this rescaling, fmax -> fmin,  fmin -> fmax.
Assumes that all elements (integers) in pop have their fitnesses defined in variant_table.
n  must be less than N/2.
"""
function reverse_propsel( pop::Population, n::Int64, variant_table::Dict{Int64,variant_type} )
  N = length(pop)
  #println("reverse_propsel N: ",N,"  n: ",n)
  if n > Int(ceil(N/2))
    error("n  (number selected) is greater than N/2 in reverse_propsel!")
  end
  fmax = 0.0
  fmin = typemax(Float64)
  for v in pop
    fitv = variant_table[v].fitness
    if fitv < 0.0
      error("negative fitness in reverse_propsel!")
    end
    if fitv > fmax
      fmax = fitv
    end
    if fitv < fmin
      fmin = fitv
    end
  end
  if fmax == 0.0
    # all elements have fitness zero
    println("Warning: all elements have fitness zero")
    return pop
  end

  selected = trues(N)
  k = 0
  while k < n
    i = rand(1:N)
    #println("i: ",i," s[i]: ",selected[i],"  fit: ",variant_table[pop[i]].fitness)
    while !selected[i]
      i = rand(1:N)
      #println("  i: ",i," s[i]: ",selected[i])
    end
    w = (-variant_table[pop[i]].fitness+fmax+fmin)/fmax
    #println("w: ",w)
    if rand() < w
      selected[i] = false
      k += 1
      #println("selected")
    end
  end
  pop[selected]
end

@doc """ function propsel_alt()
  Roulette wheel version of proportial selection
  Less efficient in almost all circumstances than the above version.
  Written as check of correctness for the above version
"""
function propsel_alt( pop::Population, dfe::Function)
  fitness_table = Dict{Int64,Float64}()
  N = length(pop)
  new_pop = zeros(Int64,N)
  fitdict = Dict{Int64,Float64}()
  popctr = pop_counter( pop )
  for p in pop
    fit = dfe_fitness(p, dfe, fitness_table )  # set fitness of p to be dfe(p).
  end
  sum_fitness = 0.0
  for p in keys(popctr)
    dfit = dfe_fitness( p, dfe, fitness_table )*popctr[p]
    sum_fitness += dfit
  end
  for p in keys(popctr)
    fitdict[p] = dfe_fitness( p, dfe, fitness_table )*popctr[p]/sum_fitness
  end
  for i = 1:N
    r = rand()
    sumfit = 0.0
    for p in keys(fitdict)
      sumfit += fitdict[p]
      if r < sumfit
        new_pop[i] = p
        break
      end
    end
  end
  new_pop
end

@doc """ function negative_propsel()
Apply proportional selection to Population pop to generate the population of size n to be deleted
Use rescaled netative fitness to favor lower fitness individuals.
The rescaling is  fit -> -fit + fmax + fmin.
Under this rescaling, fmax -> fmin,  fmin -> fmax.
Assumes that all elements (integers) in pop have their fitnesses defined in variant_table.
Not as useful as reverse_propsel() defined above.
"""
function negative_propsel( pop::Population, n::Int64, variant_table::Dict{Int64,variant_type} )
  new_pop = zeros(Int64,n)
  fmax = 0.0
  fmin = typemax(Float64)
  for v in pop
    fitv = variant_table[v].fitness
    if fitv < 0.0
      error("negative fitness in negative_propsel")
    end
    if fitv > fmax
      fmax = fitv
    end
    if fitv < fmin
      fmin = fitv
    end
  end
  if fmax == 0.0
    # all elements have fitness zero
    println("Warning: all elements have fitness zero")
    return pop
  end

  N = length(pop)
  k = 0
  while k < n
    i = rand(1:N)
    w = (-variant_table[pop[i]].fitness+fmax+fmin)/fmax
    if rand() < w
      new_pop[k + 1] = pop[i]
      k += 1
    end
  end
  new_pop
end


text:
# media(SyncPlot, Media.Plot)
# media(Plot, Media.Plot)
#
# Media.render(::Juno.PlotPane, plot::SyncPlot) = display_blink(plot)
#
# @render Juno.PlotPane plot::Plot SyncPlot(plot)


text:
module SkullSegment
using NIfTI, Compat
include(joinpath(Pkg.dir("NIfTI"), "examples", "register.jl"))

immutable XYZ
    x::Int
    y::Int
    z::Int
end

Base.setindex!{T<:Real}(X::Array{T}, val::T, x::XYZ) = (X[x.x, x.y, x.z] = val)
Base.getindex{T<:Real}(X::Array{T}, x::XYZ) = X[x.x, x.y, x.z]

function process!(processed, to_process, v::XYZ)
    if !processed[v]
        processed[v] = true
        push!(to_process, v)
    end
end

function floodfill(arr, seed::Tuple{Int, Int, Int}, lb, ub)
    lb <= arr[seed[1], seed[2], seed[3]] <= ub || warn("seed voxel not in range")
    processed = zeros(Bool, size(arr))
    out = zeros(Bool, size(arr))
    to_process = [XYZ(seed[1], seed[2], seed[3])]
    processed[seed[1], seed[2], seed[3]] = 1
    while !isempty(to_process)
        cur = pop!(to_process)
        if lb <= arr[cur] <= ub
            out[cur] = true
            cur.x > 1 && process!(processed, to_process, XYZ(cur.x-1, cur.y, cur.z))
            cur.x < size(arr, 1) && process!(processed, to_process, XYZ(cur.x+1, cur.y, cur.z))
            cur.y > 1 && process!(processed, to_process, XYZ(cur.x, cur.y-1, cur.z))
            cur.y < size(arr, 2) && process!(processed, to_process, XYZ(cur.x, cur.y+1, cur.z))
            cur.z > 1 && process!(processed, to_process, XYZ(cur.x, cur.y, cur.z-1))
            cur.z < size(arr, 3) && process!(processed, to_process, XYZ(cur.x, cur.y, cur.z+1))
        end
    end
    out
end

# Add a suffix to the part of a filename before the .
function addsuffix(path::String, suffix::String)
    dir = dirname(path)
    base = basename(path)
    dotindex = rsearch(base, '.')
    joinpath(dir, base[1:dotindex-1]*"_"*suffix*base[dotindex:end])
end

# Average several anatomical volumes
function avg{T<:String}(paths::Vector{T}, out::String)
	!isempty(paths) || error("no volume specified")
	ni = niread(paths[1], mmap=true)
	affine = getaffine(ni.header)
	niraw = isa(ni.raw, Array{Float32}) ? copy(ni.raw) : float32(ni.raw)
	for i = 2:length(paths)
		ni2 = niread(paths[i], mmap=true)
		size(ni2) == size(ni) || error("volume sizes do not match")
		affine == getaffine(ni2.header) || warn("volumes have different affine matrices")
		broadcast!(+, niraw, ni2.raw)
	end
	scale!(niraw, 1/length(paths))
	niwrite(out, NIVolume(ni.header, niraw))
end

function conform_header(origheader)
    header = deepcopy(origheader)
    header.pixdim = (header.pixdim[1], 1f0, 1f0, 1f0, header.pixdim[5:end]...)
    header.xyzt_units = 0x02 # NIfTI_UNITS_MM
    header.srow_x = map(x->x/origheader.pixdim[2], header.srow_x)
    header.srow_y = map(x->x/origheader.pixdim[3], header.srow_y)
    header.srow_z = map(x->x/origheader.pixdim[4], header.srow_z)
    header
end

# Normalize T1 scan using N3
function n3normalize(path::String)
    nu = addsuffix(path, "nu")
    wm110 = addsuffix(path, "wm110")

    # Pretend the dimensions are 1 mm in each direction because
    # FreeSurfer is stupid.
    # XXX This probably doesn't work if the slices aren't iso
    nutemp = tempname()*".nii"
    ni = niread(path, mmap=true)
    origheader = ni.header
    ni.header = conform_header(ni.header)
    niwrite(nutemp, ni)

    # Correct for inhomogeneity due to coils
    # Optimized parameters from Zheng, W., Chee, M. W. L., &
    # Zagorodnov, V. (2009). Improvement of brain segmentation accuracy
    # by optimizing non-uniformity correction using N3. NeuroImage,
    # 48(1), 73–83. doi:10.1016/j.neuroimage.2009.06.039
    try
        println(`mri_nu_correct.mni --i $nutemp --o $nu --proto-iters 1000 --distance 30`)
        run(`mri_nu_correct.mni --i $nutemp --o $nu --proto-iters 1000 --distance 30`)
    finally
        rm(nutemp)
    end

    # Normalize white matter to 110
    run(`mri_normalize -monkey $nu $wm110 -v -mprage -nosnr`)

    # Fix the dimensions
    for vol in (nu, wm110)
        ni = niread(vol, mmap=true)
        ni.header.pixdim = (ni.header.pixdim[1], origheader.pixdim[2:4]..., ni.header.pixdim[5:end]...)
        ni.header.xyzt_units = origheader.xyzt_units
        ni.header.srow_x = map(x->x*origheader.pixdim[2], ni.header.srow_x)
        ni.header.srow_y = map(x->x*origheader.pixdim[3], ni.header.srow_y)
        ni.header.srow_z = map(x->x*origheader.pixdim[4], ni.header.srow_z)
        f = open(vol, "r+")
        write(f, ni.header)
        close(f)
    end

    (nu, wm110)
end

# Generate brain mask for anatomical
function extractbrain(path::String; threshold=0.4, flood=20)
    norm = addsuffix(path, "wm110")
    isfile(norm) || n3normalize(path)

    masked = addsuffix(path, "masked")
    run(`mri_gcut $norm $masked -110 -T $threshold`)

    mask = addsuffix(path, "mask")
    ni = niread(masked)
    if flood != 0
        ni.raw[floodfill(ni.raw, (1, 1, 1), 0, flood)] = 0
        niwrite(masked, ni)
    end

    # Register back to the original nu
    maskvol = NIVolume(ni.header, 0xff*(ni.raw .!= 0))
    nu = niread(addsuffix(path, "nu"), mmap=true)
    maskvol = register(nu, maskvol)
    niwrite(mask, maskvol)
    mask
end

function postprocess_skull!(skull, nu, threshold)
    # Clean up cases where we get too much skin
    for k = 1:size(skull, 3), j = 1:size(skull, 2)
        for i = 1:size(skull, 1)
            if skull[i, j, k] != 0
                for i2 = i:size(skull, 1)
                    if nu[i2, j, k] >= threshold
                        skull[i2, j, k] = 0
                    else
                        break
                    end
                end
                break
            end
        end
        for i = size(skull, 1):-1:1
            if skull[i, j, k] != 0
                for i2 = i:-1:1
                    if nu[i2, j, k] >= threshold
                        skull[i2, j, k] = 0
                    else
                        break
                    end
                end
                break
            end
        end
    end
    skull
end

# Extract skull
function extractskull(path::String; lower_threshold=0, upper_threshold=0, erode_threshold=90)
    norm = addsuffix(path, "nu")
    isfile(norm) || n3normalize(path)
    mask = addsuffix(path, "mask")
    isfile(mask) || extractbrain(path)

    maskvol = niread(mask)
    center = size(maskvol, 2) >> 1
    slice = maskvol.raw[:, center, :]
    proj = vec(sum(slice, 1))
    bottom = findfirst(proj)-20
    top = findlast(proj)+20
    proj = vec(sum(slice, 3))
    left = findfirst(proj)-30
    right = findlast(proj)+30
    slice = maskvol.raw[left+div(right-left+1, 2), :, bottom+div(top-bottom+1, 2)]
    proj = vec(slice)
    front = findfirst(proj)-30
    back = findlast(proj)+50
    center = div(back-front+1, 2)

    vol = niread(norm)
    volsubset = vol.raw[left:right, front:back, bottom:top]
    fill!(vol.raw, 0)
    vol.raw[left:right, front:back, bottom:top] = volsubset
    cropped = addsuffix(path, "cropped")
    # origheader = vol.header
    # vol.header = conform_header(origheader)
    niwrite(cropped, vol)

    masksubset = maskvol.raw[left:right, front:back, bottom:top]
    fill!(maskvol.raw, 0)
    maskvol.raw[left:right, front:back, bottom:top] = masksubset
    cropped_mask = addsuffix(path, "cropped_mask")
    # maskvol.header = conform_header(maskvol.header)
    niwrite(cropped_mask, maskvol)

    skull = addsuffix(path, "skull")
    run(`skullfinder -i $cropped -o $skull -m $cropped_mask -v 2 -l $lower_threshold -u $upper_threshold --scalplabel 0 --skulllabel 255 --spacelabel 255 --brainlabel 0`)

    # skullfinder screws up the registration, so fix it
    skullvol = niread(skull)
    skullvol.header = vol.header
    postprocess_skull!(skullvol.raw, vol.raw, erode_threshold)

    niwrite(skull, skullvol)
end

# View brain
# function viewbrain(vol::NIVolume)
#     f = figure()
#     @manipulate for z = 1:size(vol, 3); withfig(f) do
#         imshow(vol.raw[:, :, z].', cmap="gray", origin="upper")
#     end end
# end
# viewbrain(vol::String) = viewbrain(niread(vol))
end


text:
using MosekTools, DifferentialEquations, PyPlot, LaTeXStrings, MarkovBounds
DE = DifferentialEquations
u = []
@polyvar(x[1:2])
@polyvar(t)

μ = 1/2
f = [x[2]; - x[1] + μ*(1 - x[1]^2) * x[2]]
g = [0; 0.5]
σ = g*g'

Tf = 2.0
x0 = [0.5; 0.5]
U = FullSpace()
X = @set(x[1] >= 0)

## data
order = 8
nt = 20
trange = range(Tf/nt, stop = Tf, length=nt)

## solve
MP = DiffusionProcess(x, f, σ, FullSpace())
CP = ControlProcess(MP, Tf, u, t, U, nothing, nothing, ExitProbability(X), 0)

bnds, V = optimal_control(CP, x0, order, trange, Mosek.Optimizer; value_function = true)

N = 100
drift(x,p,t) = [f[1](x[2]); f[2](x[1],x[2])]
diffusion(x,p,t) = [0; g[2]]
sde = DE.SDEProblem(drift, diffusion, x0, (0.0,Tf))
sde_ensemble = DE.EnsembleProblem(sde)

trajectories = DE.solve(sde_ensemble, trajectories=N, EM(), dt = 0.01)

fig, ax = subplots(figsize=(5,5))
stayed = 0
for sol in trajectories
    max_idx = findfirst(m -> m[1] < 0, sol.u)
    if max_idx == nothing
        global stayed += 1
    end
    max_idx = (max_idx == nothing ? length(sol.u) : max_idx)
    global ln = ax.plot([sol.u[i][1] for i in 1:max_idx], [sol.u[i][2] for i in 1:max_idx], color="blue", linewidth=0.5)
end
ax.scatter([x0[1]], [x0[2]], marker="o", color="black", s=5)
ax.set_xlabel(L"x_1")
ax.set_ylabel(L"x_2")
display(fig)

println("Upper Bound:  P(exit) ≦ ", bnds.bounds)
println("Sampled: P(exit) = ", stayed/N)


text:
export TracerProvider

"""
    TracerProvider(;kw...)

# Keyword Arguments

  - `sampler`::[`AbstractSampler`](@ref)=[`DEFAULT_ON`](@ref)
  - `resource`=[`Resource`](@ref)()
  - `span_processor`::[`AbstractSpanProcessor`](@ref)=[`CompositSpanProcessor`](@ref)()
  - `id_generator`::[`AbstractIdGenerator`](@ref)=[`RandomIdGenerator`](@ref)()
  - `limit_info`=[`LimitInfo`](@ref)()

The following extra methods are provided beyond those defined in [`AbstractTracerProvider`](@ref):

  - [`force_flush!(p::TracerProvider)`](@ref)
  - [`shut_down!(p::TracerProvider)`](@ref)
  - [`Base.push!(p::TracerProvider, sp::AbstractSpanProcessor)`](@ref)
"""
Base.@kwdef struct TracerProvider{
    S<:AbstractSampler,
    R<:Resource,
    SP<:AbstractSpanProcessor,
    RNG,
} <: AbstractTracerProvider
    sampler::S = DEFAULT_ON
    resource::R = Resource()
    span_processor::SP = CompositSpanProcessor()
    id_generator::RNG = RandomIdGenerator()
    limit_info::LimitInfo = LimitInfo()
    is_shut_down::Ref{Bool} = Ref(false)
end

OpenTelemetryAPI.resource(p::TracerProvider) = p.resource

"""
    force_flush!(p::TracerProvider)

Shorthand to force flush inner span processors
"""
force_flush!(p::TracerProvider) = force_flush!(p.span_processor)

"""
    shut_down!(p::TracerProvider)

Shut down inner span processors and then mark itself as shut down.
"""
function shut_down!(p::TracerProvider)
    shut_down!(p.span_processor)
    p.is_shut_down[] = true
end

"""
    Base.push!([p::TracerProvider], sp::AbstractSpanProcessor)

Add an extra span processor `sp` into the [`TracerProvider`](@ref) p.
"""
Base.push!(sp::AbstractSpanProcessor) = push!(global_tracer_provider(), sp)
Base.push!(p::TracerProvider, sp::AbstractSpanProcessor) = push!(p.span_processor, sp)


text:
using ImageCore, ImageView, TestImages

z = ones(10,50);
y = 8; x = 2;
z[y,x] = 0
guidict = imshow_now(z)
idx = annotate!(guidict, AnnotationText(x, y, "x", color=RGB(0,0,1), fontsize=3))
idx2 = annotate!(guidict, AnnotationPoint(x+10, y, shape='.', size=4, color=RGB(1,0,0)))
idx3 = annotate!(guidict, AnnotationPoint(x+20, y-6, shape='.', size=1, color=RGB(1,0,0), linecolor=RGB(0,0,0), scale=true))
idx4 = annotate!(guidict, AnnotationLine(x+10, y, x+20, y-6, linewidth=2, color=RGB(0,1,0)))
idx5 = annotate!(guidict, AnnotationBox(x+10, y, x+20, y-6, linewidth=2, color=RGB(0,0,1)))
idx6 = annotate!(guidict, AnnotationPoints([(x+25, y-3), (x+24, y-3), (x+23, y-3)],
                                           shape='o', size=3, color=RGB(1,0,1),
                                           linecolor=RGB(0,0,0), scale=true))

sleep(0.01)
delete!(guidict, idx)
sleep(0.01)

img = testimage("lighthouse")
guidict = imshow_now(img)
scalebar(guidict, 30; x = 0.1, y = 0.05)

# Grid of images (issue #202)
gd = imshow_gui((300, 300), (2, 2))
canvases = gd["canvas"]
anns = [annotations() annotations();
        annotations() annotations()]
makeimage(color) = fill(color, 100, 100)
imgs = makeimage.([colorant"red" colorant"green";
                   colorant"blue" colorant"purple"])
roidict = [imshow(canvases[1,1], imgs[1,1], anns[1,1]) imshow(canvases[1,2], imgs[1,2], anns[1,2]);
           imshow(canvases[2,1], imgs[2,1], anns[2,1]) imshow(canvases[2,2], imgs[2,2], anns[2,2])]
idx = annotate!(anns[2,2], canvases[2,2], roidict[2,2], AnnotationBox(5, 5, 30, 80, linewidth=3, color=RGB(1,1,0)))


text:
# ODEInterface.jl Algorithms

abstract type ODEInterfaceAlgorithm <: DiffEqBase.AbstractODEAlgorithm end
abstract type ODEInterfaceImplicitAlgorithm <: ODEInterfaceAlgorithm end
abstract type ODEInterfaceExplicitAlgroithm <: ODEInterfaceAlgorithm end
struct dopri5 <: ODEInterfaceExplicitAlgroithm end
struct dop853 <: ODEInterfaceExplicitAlgroithm end
struct odex <: ODEInterfaceExplicitAlgroithm end
struct seulex{T} <: ODEInterfaceImplicitAlgorithm
    jac_lower::T
    jac_upper::T
end
struct radau{T} <: ODEInterfaceImplicitAlgorithm
    jac_lower::T
    jac_upper::T
end
struct radau5{T} <: ODEInterfaceImplicitAlgorithm
    jac_lower::T
    jac_upper::T
end
struct rodas{T} <: ODEInterfaceImplicitAlgorithm
    jac_lower::T
    jac_upper::T
end
struct ddeabm <: ODEInterfaceExplicitAlgroithm end
struct ddebdf{T} <: ODEInterfaceImplicitAlgorithm
    jac_lower::T
    jac_upper::T
end

seulex(;jac_lower=nothing,jac_upper=nothing) = seulex(jac_lower,jac_upper)
radau(;jac_lower=nothing,jac_upper=nothing) = radau(jac_lower,jac_upper)
radau5(;jac_lower=nothing,jac_upper=nothing) = radau5(jac_lower,jac_upper)
rodas(;jac_lower=nothing,jac_upper=nothing) = rodas(jac_lower,jac_upper)
ddebdf(;jac_lower=nothing,jac_upper=nothing) = ddebdf(jac_lower,jac_upper)


text:
using LinearAlgebra
using PyPlot
function linesegment_intersect(p1, p2, q1, q2, δ)
    A = [p2-p1 q1-q2]
    b = q1-p1
    if rank(A) == 2
        st = A \ b
        if all(-δ .≤ st .≤ 1 + δ)
            pintersect = p1 + st[1] * (p2 - p1)
            return true, pintersect
        else
            return false, nothing
        end
    else
        # Warning: parallel lines, should check if they actually overlap
        return false, nothing
    end
end

# Example: n random lines, find and plot all intersections
n = 10
lines = [ [rand(2), rand(2)] for i = 1:n ]

# Plot all lines
for line in lines
    plot([line[1][1], line[2][1]], [line[1][2], line[2][2]])
end

# Find all intersections - note that we only consider j>i to avoid checking
# the same 2 line-segments twice
for i = 1:n
    for j = i+1:n
        isect, pq = linesegment_intersect(lines[i][1], lines[i][2], lines[j][1], lines[j][2], 0)
        if isect
            plot(pq[1], pq[2], "o")
        end
    end
end

axis("equal"); grid(true);


text:
using DDD
using Test
using Plots
gr()
cd(@__DIR__)

@testset "Plot nodes" begin
    # Load and create.
    fileSlipSystem = "./testData/BCC.json"
    fileDislocationLoop = "./testData/sampleDislocation.json"
    dictSlipSystem = loadJSON(fileSlipSystem)
    slipSystems = loadSlipSystem(dictSlipSystem)
    # There can be multiple dislocations per simulation parameters.
    dictDislocationLoop = loadJSON(fileDislocationLoop)
    loops = zeros(DislocationLoop, length(dictDislocationLoop))
    if typeof(dictDislocationLoop) <: AbstractArray
        loops = zeros(DislocationLoop, length(dictDislocationLoop))
        for i in eachindex(loops)
            loops[i] = loadDislocationLoop(dictDislocationLoop[i], slipSystems)
        end
    else
        loops = loadDislocationLoop(dictDislocationLoop, slipSystems)
    end
    network = DislocationNetwork(loops)

    function sumNodes(loops)
        totalNodes = 0
        for i in eachindex(loops)
            totalNodes += loops[i].numSides * loops[i].nodeSide * loops[i].numLoops
        end
        return totalNodes
    end
    totalNodes = sumNodes(loops)
    # Modifying current plot.
    fig = plot()
    plotNodes!(fig, network, m = 1, l = 3, legend = false)
    @test fig.n == totalNodes
    # Creating new plot.
    fig2 = plotNodes(network, m = 1, l = 3, legend = false)
    @test fig2.n == totalNodes

    function plotLoops(fig, loops)
        totalNodes = 0
        return for i in eachindex(loops)
            plotNodes!(fig, loops[i], m = 1, l = 3, legend = false)
            totalNodes += loops[i].numSides * loops[i].nodeSide * loops[i].numLoops
            @test fig.n == totalNodes
        end
    end
    function plotLoops(loops)
        return for i in eachindex(loops)
            fig = plotNodes(loops[i], m = 1, l = 3, legend = false)
            @test fig.n == loops[i].numSides * loops[i].nodeSide * loops[i].numLoops
        end
    end
    fig = plot()
    plotLoops(fig, loops)
    plotLoops(loops)

    DictMaterialParameters = loadJSON("./testData/sampleMaterialParameters.json")
    DictFEMParameters = loadJSON("./testData/sampleFEMParameters.json")
    matParams = loadMaterialParameters(DictMaterialParameters)
    femParams = loadFEMParameters(DictFEMParameters)
    regularCuboidMesh = buildMesh(matParams, femParams)

    fig = plot()
    plotNodes!(fig, regularCuboidMesh, network, m = 1, l = 3, legend = false)
    @test fig.n == totalNodes + 6
    # Creating new plot.
    fig2 = plotNodes(regularCuboidMesh, network, m = 1, l = 3, legend = false)
    @test fig2.n == totalNodes + 6

    fig3 = plotFEDomain(regularCuboidMesh)
    @test fig3.n == 19
end

max_stars_repo_name: dcelisgarza/DDD

max_forks_repo_name: dcelisgarza/DDD


text:
module FileIO

export read_csv, read_simple, read_big

import ..Advent2020

using DelimitedFiles: readdlm

"""
    read_csv(f_name)

Read a CSV file.
"""
read_csv(f_name) = open(f->readdlm(f, ',', Int), f_name)

"""
    read_simple(f_name)

Simple data reading function.
"""
read_simple(f_name, T=Int) = parse.(T, open(readlines, f_name))

"""
    read_big(f_name)

Read in single big integer
"""
read_big(f_name) = open(f->readlines(f)[1], f_name)

end


text:
#=
  From srl-pp-tutorial-wasp-stockholm.pdf
  "Statistical Relational learning and Probabilistic Programming"
  by Luc De Raedt, Anton Dries, and Angelika Kimmig
  https://dtai.cs.kuleuven.be/problog/wasp17-tutorial.html

  slide 394f:
  """
  - Flipping a coin with unknown weight
  - Prior: uniform distribution on [0,1]
  - Observation: 5x heads in a row
  """

  The ProbLog model return the following which corresponds to the
  density of the probabilist
  """
   weight(c1,0.1): 3.3994697e-13
   weight(c1,0.3): 2.1679411e-06
   weight(c1,0.5): 0.0041497433
   weight(c1,0.7): 0.1317485
   weight(c1,0.9): 0.86409959         <----

   weight(c2,0.1): 3.2276726e-06
   weight(c2,0.3): 0.024109997
   weight(c2,0.5): 0.66724754         <----
   weight(c2,0.7): 0.30628626
   weight(c2,0.9): 0.0023529733
  """

  Here we observe 13 tosses and detecting the probability of throwing a head.

  In this version we simulate the ProbLog model more closely
  than in coins_learning.jl

  Result:
  Distributions of variable param (num:0)
  0.90000 =>    8646  (0.864600)
  0.70000 =>    1313  (0.131300)
  0.50000 =>      41  (0.004100)

  Distributions of variable ix (num:0)
  5.00000 =>    8646  (0.864600)
  4.00000 =>    1313  (0.131300)
  3.00000 =>      41  (0.004100)

  Distributions of variable ix1 (num:0)
  0.00000 =>   10000  (1.000000)

  Distributions of variable ix2 (num:0)
  0.00000 =>   10000  (1.000000)

  Distributions of variable ix3 (num:0)
  0.00000 =>    9959  (0.995900)
  1.00000 =>      41  (0.004100)

  Distributions of variable ix4 (num:0)
  0.00000 =>    8687  (0.868700)
  1.00000 =>    1313  (0.131300)

  Distributions of variable ix5 (num:0)
  1.00000 =>    8646  (0.864600)
  0.00000 =>    1354  (0.135400)

  Cf ~/webppl/coins_learning2.wppl

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function coins_learning(obs)
    params = [0.1,0.3,0.5,0.7,0.9];
    params_prior = [0.05,0.2,0.5,0.2,0.05];
    len = length(params)

    #  Generate the index for param/weight interval
    #  according to the priors in params_prior
    ix ~ Categorical(simplex(params_prior)) # [0,1,2,3,4]

    #  generate the param value to use for this data point
    obs ~ filldist(flip(params[ix]),length(obs))

    param ~ Dirac(params[ix])
    ix1 ~ Dirac(ix == 1)
    ix2 ~ Dirac(ix == 2)
    ix3 ~ Dirac(ix == 3)
    ix4 ~ Dirac(ix == 4)
    ix5 ~ Dirac(ix == 5)
end

# p: 0.9344
obs = [true for _ in 1:13]
# p: 0.5984
# obs = [true,false,true,true,true,true,true,false,false,true,false,false,true]
model = coins_learning(obs)
num_chains = 4

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
# chns = sample(model, PG(5), 10_000)
chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

show_var_dist_pct(chns,:param)
show_var_dist_pct(chns,:ix)
show_var_dist_pct(chns,:ix1)
show_var_dist_pct(chns,:ix2)
show_var_dist_pct(chns,:ix3)
show_var_dist_pct(chns,:ix4)
show_var_dist_pct(chns,:ix5)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using Test
using LinearAlgebra
using SparseArrays
using ExtendableSparse
using Printf
using BenchmarkTools

using Pardiso
using AlgebraicMultigrid
using IncompleteLU

##############################################################
@testset "Constructors" begin
    function test_constructors()
        m1=ExtendableSparseMatrix(10,10)
        m2=ExtendableSparseMatrix(Float16,10,10)
        m3=ExtendableSparseMatrix(Float32,Int16,10,10)

        csc=sprand(10,10,0.1)
        m4=ExtendableSparseMatrix(csc)
        sprand!(m4,6)
        flush!(m4)
        true
    end
    @test test_constructors()
end

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

@testset "Updates" begin
    A=ExtendableSparseMatrix(10,10)
    @test nnz(A)==0
    A[1,3]=5
    updateindex!(A,+,6.0,4,5)
    updateindex!(A,+,0.0,2,3)
    @test nnz(A)==2
    rawupdateindex!(A,+,0.0,2,3)
    @test nnz(A)==3
    dropzeros!(A)
    @test nnz(A)==2
    rawupdateindex!(A,+,0.1,2,3)
    @test nnz(A)==3
    dropzeros!(A)
    @test nnz(A)==3
end
#################################################################
function test_timing(k,l,m)
    t1=@belapsed fdrand($k,$l,$m,matrixtype=$SparseMatrixCSC) seconds=0.1
    t2=@belapsed fdrand($k,$l,$m,matrixtype=$ExtendableSparseMatrix)  seconds=0.1
    t3=@belapsed fdrand($k,$l,$m,matrixtype=$SparseMatrixLNK)  seconds=0.1
    @printf("CSC: %.4f EXT: %.4f LNK: %.4f\n",t1*1000,t2*1000,t3*1000 )
    t3<t2<t1
end

@testset "Timing" begin
    @test test_timing(1000,1,1)
    @test test_timing(100,100,1)
    @test test_timing(20,20,20)
end

##################################################################################
function test_assembly(;m=1000,n=1000,xnnz=5000,nsplice=1)
    A=ExtendableSparseMatrix{Float64,Int64}(m,n)
    m,n=size(A)
    S=spzeros(m,n)
    for isplice=1:nsplice
        for inz=1:xnnz
            i=rand((1:m))
            j=rand((1:n))
            a=1.0+rand(Float64)
            S[i,j]+=a
            A[i,j]+=a
        end
        flush!(A)
        for j=1:n
            @assert(issorted(A.cscmatrix.rowval[A.cscmatrix.colptr[j]:A.cscmatrix.colptr[j+1]-1]))
        end
        @assert(nnz(S)==nnz(A))

        (I,J,V)=findnz(S)
        for inz=1:nnz(S)
            @assert(A[I[inz],J[inz]]==V[inz])
        end

        (I,J,V)=findnz(A)
        for inz=1:nnz(A)
            @assert(S[I[inz],J[inz]]==V[inz])
        end

    end
    return true
end

@testset "Assembly" begin
    @test test_assembly(m=10,n=10,xnnz=5)
    @test test_assembly(m=100,n=100,xnnz=500,nsplice=2)
    @test test_assembly(m=1000,n=1000,xnnz=5000,nsplice=3)

    @test test_assembly(m=20,n=10,xnnz=5)
    @test test_assembly(m=200,n=100,xnnz=500,nsplice=2)
    @test test_assembly(m=2000,n=1000,xnnz=5000,nsplice=3)

    @test test_assembly(m=10,n=20,xnnz=5)
    @test test_assembly(m=100,n=200,xnnz=500,nsplice=2)
    @test test_assembly(m=1000,n=2000,xnnz=5000,nsplice=3)

    for irun=1:10
        m=rand((1:10000))
        n=rand((1:10000))
        nnz=rand((1:10000))
        nsplice=rand((1:5))
        @test test_assembly(m=m,n=n,xnnz=nnz,nsplice=nsplice)
    end
end

#####################################################################
function test_transient_construction(;m=10,n=10,d=0.1) where {Tv,Ti}
    csc=sprand(m,n,d)
    lnk=SparseMatrixLNK(csc)
    csc2=SparseMatrixCSC(lnk)
    return csc2==csc
end

@testset "Transient Construction" begin
    for irun=1:10
        m=rand((1:1000))
        n=rand((1:1000))
        d=0.3*rand()
        @test test_transient_construction(m=m,n=n,d=d)
    end
end

#####################################################################
function test_addition(;m=10,n=10,d=0.1) where {Tv,Ti}
    csc=sprand(m,n,d)
    lnk=SparseMatrixLNK(csc)
    csc2=csc+lnk
    csc2==2*csc
end

function test_operations(n)
    A=ExtendableSparseMatrix(n,n)
    sprand_sdd!(A)
    b=rand(n)
    x=A\b
    Ax=A*x
    b≈ Ax
end

@testset "Operations" begin
    for irun=1:10
        m=rand((1:1000))
        n=rand((1:1000))
        d=0.3*rand()
        @test test_addition(m=m,n=n,d=d)
    end

    @test test_operations(10)
    @test test_operations(100)
    @test test_operations(1000)
end

##############################################
function test_fdrand(k,l,m)
    A=fdrand(k,l,m,matrixtype=ExtendableSparseMatrix)
    jacobi_iteration_matrix=I-inv(Diagonal(A))*A
    ext=extrema(real(eigvals(Matrix(jacobi_iteration_matrix))))
    mininv=minimum(inv(Matrix(A)))
    abs(ext[1])<1 &&  abs(ext[2])<1 && mininv>0
end

@testset "fdrand" begin
  @test  test_fdrand(1000,1,1)
  @test  test_fdrand(20,20,1)
  @test  test_fdrand(10,10,10)
end

##############################################
function test_fdrand_coo(k,l,m)
    Acsc=fdrand(k,l,m,rand=()->1,matrixtype=SparseMatrixCSC)
    Acoo=fdrand(k,l,m,rand=()->1,matrixtype=:COO)
    Acsc≈Acoo
end

@testset "fdrand_coo" begin
  @test  test_fdrand_coo(1000,1,1)
  @test  test_fdrand_coo(20,20,1)
  @test  test_fdrand_coo(10,10,10)
end

##############################################
function test_fdrand_update(k,l,m)
    A1=fdrand(k,l,m,rand=()->1,matrixtype=ExtendableSparseMatrix,update = (A,v,i,j)-> A[i,j]+=v)
    A2=fdrand(k,l,m,rand=()->1,matrixtype=ExtendableSparseMatrix,update = (A,v,i,j)-> rawupdateindex!(A,+,v,i,j))
    A3=fdrand(k,l,m,rand=()->1,matrixtype=ExtendableSparseMatrix,update = (A,v,i,j)-> updateindex!(A,+,v,i,j))

    A1≈A2 && A1 ≈ A3
end

@testset "fdrand_update" begin
  @test  test_fdrand_update(1000,1,1)
  @test  test_fdrand_update(20,20,1)
  @test  test_fdrand_update(10,10,10)
end

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

function test_precon(Precon,k,l,m;maxiter=10000)
    A=fdrand(k,l,m,matrixtype=ExtendableSparseMatrix, rand= ()-> 1.0)
    b=ones(size(A,2))
    exact=A\b
    Pl=Precon(A)
    it,hist=simple(A,b,Pl=Pl,maxiter=maxiter,reltol=1.0e-10,log=true)
    r=hist[:resnorm]
    nr=length(r)
    tail=min(100,length(r)÷2)
    all(x-> x<1,r[end-tail:end]./r[end-tail-1:end-1]),norm(it-exact)
end

function test_precon2(precon,k,l,m;maxiter=10000)
    A=fdrand(k,l,m,matrixtype=ExtendableSparseMatrix, rand= ()-> 1.0)
    b=ones(size(A,2))
    exact=A\b
    factorize!(precon,A)
    it,hist=simple(A,b,Pl=precon,maxiter=maxiter,reltol=1.0e-10,log=true)
    r=hist[:resnorm]
    nr=length(r)
    tail=min(100,length(r)÷2)
    all(x-> x<1,r[end-tail:end]./r[end-tail-1:end-1]),norm(it-exact)
end

@testset "preconditioners" begin
    @test   all(isapprox.(test_precon(ILU0Preconditioner,20,20,20),           (true, 1.3535160424212675e-5), rtol=1.0e-5))
    @test   all(isapprox.(test_precon(JacobiPreconditioner,20,20,20),         (true, 2.0406032775945658e-5), rtol=1.0e-5))
    @test   all(isapprox.(test_precon(ParallelJacobiPreconditioner,20,20,20), (true, 2.0406032775945658e-5), rtol=1.0e-5))
    @test   all(isapprox.(test_precon(ILUTPreconditioner,20,20,20),           (true, 1.2719511868322086e-5), rtol=1.0e-5))
    @test   all(isapprox.(test_precon(AMGPreconditioner,20,20,20),            (true, 6.863753664354144e-7), rtol=1.0e-2))
end

@testset "preconditioners 2" begin
    @test   all(isapprox.(test_precon2(ILU0Preconditioner(),20,20,20),           (true, 1.3535160424212675e-5), rtol=1.0e-5))
    @test   all(isapprox.(test_precon2(JacobiPreconditioner(),20,20,20),         (true, 2.0406032775945658e-5), rtol=1.0e-5))
    @test   all(isapprox.(test_precon2(ParallelJacobiPreconditioner(),20,20,20), (true, 2.0406032775945658e-5), rtol=1.0e-5))
    @test   all(isapprox.(test_precon2(ILUTPreconditioner(),20,20,20),           (true, 1.2719511868322086e-5), rtol=1.0e-5))
    @test   all(isapprox.(test_precon2(AMGPreconditioner(),20,20,20),            (true, 6.863753664354144e-7), rtol=1.0e-2))
end

##############################################
function test_symmetric(n,uplo)
    A=ExtendableSparseMatrix(n,n)
    sprand_sdd!(A)
    b=rand(n)
    flush!(A)
    SA=Symmetric(A,uplo)
    Scsc=Symmetric(A.cscmatrix,uplo)
    SA\b≈ Scsc\b
end

@testset "symmetric" begin
    @test test_symmetric(3,:U)
    @test test_symmetric(3,:L)
    @test test_symmetric(30,:U)
    @test test_symmetric(30,:L)
    @test test_symmetric(300,:U)
    @test test_symmetric(300,:L)
end

##############################################
function test_hermitian(n,uplo)
    A=ExtendableSparseMatrix{ComplexF64,Int64}(n,n)
    sprand_sdd!(A)
    flush!(A)
    A.cscmatrix.nzval.=(1.0+0.01im)*A.cscmatrix.nzval
    b=rand(n)
    HA=Hermitian(A,uplo)
    Hcsc=Hermitian(A.cscmatrix,uplo)
    HA\b≈ Hcsc\b
end

@testset "hermitian" begin
    @test test_hermitian(3,:U)
    @test test_hermitian(3,:L)
    @test test_hermitian(30,:U)
    @test test_hermitian(30,:L)
    @test test_hermitian(300,:U)
    @test test_hermitian(300,:L)
end

function test_lu1(k,l,m; lufac=LUFactorization())
    Acsc=fdrand(k,l,m,rand=()->1,matrixtype=SparseMatrixCSC)
    b=rand(k*l*m)
    LUcsc=lu(Acsc)
    x1csc=LUcsc\b
    for i=1:k*l*m
        Acsc[i,i]+=1.0
    end
    LUcsc=lu!(LUcsc,Acsc)
    x2csc=LUcsc\b

    Aext=fdrand(k,l,m,rand=()->1,matrixtype=ExtendableSparseMatrix)
    lu!(lufac,Aext)
    x1ext=lufac\b
    for i=1:k*l*m
        Aext[i,i]+=1.0
    end
    update!(lufac)
    x2ext=lufac\b
    x1csc≈x1ext && x2csc ≈ x2ext
end

function test_lu2(k,l,m;lufac=LUFactorization())
    Aext=fdrand(k,l,m,rand=()->1,matrixtype=ExtendableSparseMatrix)
    b=rand(k*l*m)
    lu!(lufac,Aext)
    x1ext=lufac\b
    for i=4:k*l*m-3
        Aext[i,i+3]-=1.0e-4
        Aext[i-3,i]-=1.0e-4
    end
    lufac=lu!(lufac,Aext)
    x2ext=lufac\b
    all(x1ext.< x2ext)
end

@testset "LUFactorization" begin
    @test test_lu1(10,10,10)
    @test test_lu1(25,40,1)
    @test test_lu1(1000,1,1)

    @test test_lu2(10,10,10)
    @test test_lu2(25,40,1)
    @test test_lu2(1000,1,1)
end

@testset "Cholesky" begin
    @test test_lu1(10,10,10,lufac=CholeskyFactorization())
    @test test_lu1(25,40,1,lufac=CholeskyFactorization())
    @test test_lu1(1000,1,1,lufac=CholeskyFactorization())

    @test test_lu2(10,10,10,lufac=CholeskyFactorization())
    @test test_lu2(25,40,1,lufac=CholeskyFactorization())
    @test test_lu2(1000,1,1,lufac=CholeskyFactorization())
end

@testset "mkl-pardiso" begin
    @test test_lu1(10,10,10,lufac=MKLPardisoLU())
    @test test_lu1(25,40,1,lufac=MKLPardisoLU())
    @test test_lu1(1000,1,1,lufac=MKLPardisoLU())

    @test test_lu2(10,10,10,lufac=MKLPardisoLU())
    @test test_lu2(25,40,1,lufac=MKLPardisoLU())
    @test test_lu2(1000,1,1,lufac=MKLPardisoLU())

    @test test_lu2(10,10,10,lufac=MKLPardisoLU(mtype=2))
    @test test_lu2(25,40,1,lufac=MKLPardisoLU(mtype=2))
    @test test_lu2(1000,1,1,lufac=MKLPardisoLU(mtype=2))
end

if Pardiso.PARDISO_LOADED[]
    @testset "pardiso" begin
        @test test_lu1(10,10,10,lufac=PardisoLU())
        @test test_lu1(25,40,1,lufac=PardisoLU())
        @test test_lu1(1000,1,1,lufac=PardisoLU())

        @test test_lu2(10,10,10,lufac=PardisoLU())
        @test test_lu2(25,40,1,lufac=PardisoLU())
        @test test_lu2(1000,1,1,lufac=PardisoLU())

        @test test_lu2(10,10,10,lufac=PardisoLU(mtype=2))
        @test test_lu2(25,40,1,lufac=PardisoLU(mtype=2))
        @test test_lu2(1000,1,1,lufac=PardisoLU(mtype=2))

    end
end


text:
### A Pluto.jl notebook ###
# v0.14.7

using Markdown
using InteractiveUtils

# ╔═╡ 86ddd2e2-5eda-11ec-22b8-e7be38483773
begin
	using HTTP
	res = HTTP.get("https://xkcd.com/552/info.0.json");
end

# ╔═╡ f0acd270-dbf3-4fde-b52e-d8443a203e02
begin
	using JSON3
	 JSON3.read(res.body)

end

# ╔═╡ 77cc4849-8f6c-400e-bf9d-e754c09eb416
begin
	struct Comic
    json::JSON3.Object
	end

    Comic() = Comic(JSON3.read(HTTP.get("https://xkcd.com/info.0.json").body))

    Comic(i::Int) =         Comic(JSON3.read(HTTP.get("https://xkcd.com/$i/info.0.json").body))
Comic(552)

end

# ╔═╡ Cell order:
# ╠═86ddd2e2-5eda-11ec-22b8-e7be38483773
# ╠═f0acd270-dbf3-4fde-b52e-d8443a203e02
# ╠═77cc4849-8f6c-400e-bf9d-e754c09eb416

max_stars_repo_name: BioMedAI/Pluto-Note

max_issues_repo_name: BioMedAI/Pluto-Note

max_forks_repo_name: BioMedAI/Pluto-Note


text:
using JuliaPetra
using Test
using TypeStability

include("TestUtil.jl")

# function based tests
include("DenseMultiVectorTests.jl")
include("BasicDirectoryTests.jl")
include("LocalCommTests.jl")

#need access to MPI comm to ensure that tests only run on the root process
#use distinct types
GID = UInt64
PID = UInt16
LID = UInt32
comm = MPIComm(GID, PID, LID)

pid = myPid(comm)
nProcs = numProc(comm)

#only print errors from one process
if pid != 1
    #redirect_stdout()
    #redirect_stderr()
end

#tries are to allow barriers to work correctly, even under erronious situtations
try
    @testset "MPI Tests" begin
        try
            @testset "Comm MPI Tests" begin
                include("MPICommTests.jl")
                include("MPIBlockMapTests.jl")
                include("MPIimport-export Tests.jl")

                include("LocalCommTests.jl")
                runLocalCommTests(comm)

                include("BasicDirectoryTests.jl")
                basicDirectoryTests(comm)
            end

            @testset "Data MPI Tests" begin
                denseMultiVectorTests(comm)

                include("CSRMatrixMPITests.jl")
            end

        finally
            #print results sequentially
            for i in 1:pid
                barrier(comm)
            end
        end
        @info "process $pid test results:"
    end

finally
    #print results sequentially
    for i in pid:nProcs
        barrier(comm)
    end
end


text:
import KitML
using KitML.PyCall

# load tensorflow model
cd(@__DIR__)
model = KitML.load_model("best_model.h5"; mode = :tf)
model.summary()

# prepare data
X = randn(Float32, 10, 4)
Y = randn(Float32, 10, 1)

# reinforcement training
KitML.sci_train!(model, (X, Y); epoch = 2)
# it's equivalent as
model.fit(X, Y, epochs=2)

model.predict(X)

KitML.neural_closure(model, X)


text:
using Base.Test
import COBRAUtils

@testset "manipulations" begin

@testset "remove reaction" begin
    model = COBRAUtils.loadEColiTestModel()
    n = length(model.rxns)
    model2 = COBRAUtils.remove_reaction(model, 1)
    @test length(model2.rxns) == length(model.rxns) - 1
    @test length(model2.lb) == length(model.lb) - 1
    @test length(model2.ub) == length(model.ub) - 1
    @test length(model2.c) == length(model.c) - 1
    @test size(model2.S) == (size(model.S, 1), n - 1)
end

@testset "remove reaction list" begin
    model = COBRAUtils.loadEColiTestModel()
    n = length(model.rxns)
    model2 = COBRAUtils.remove_reactions(model, [1,2])
    @test length(model2.rxns) == length(model.rxns) - 2
    @test length(model2.lb) == length(model.lb) - 2
    @test length(model2.ub) == length(model.ub) - 2
    @test length(model2.c) == length(model.c) - 2
    @test size(model2.S) == (size(model.S, 1), n - 2)
end

end


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

JLLWrappers.@generate_main_file_header("Qt6Base")
JLLWrappers.@generate_main_file("Qt6Base", UUID("c0090381-4147-56d7-9ebc-da0b1113ec56"))
end  # module Qt6Base_jll


text:
using DynamicGrids, Dispersal, Test, Unitful, Dates
using DynamicGrids: applyrule, SimData, extent
using Unitful: d

@testset "discrete growth" begin
    init =  [1.0 4.0 7.0;
             2.0 5.0 8.0;
             3.0 6.0 9.0]
    test1 = init
    test2 = init.*2.0
    test3 = init.*2.0.^2

    output = ArrayOutput(init; tspan=1:3)
    rule = Ruleset(DiscreteGrowth(intrinsicrate=2.0))
    sim!(output, rule)

    @test output[1] == test1
    @test output[2] == test2
    @test output[3] == test3
end

@testset "discrete growth map" begin
    init =  [1.0 4.0 7.0;
             2.0 5.0 8.0;
             3.0 6.0 9.0]
    suit =  [1.0 1.0 2.0;
             2.0 1.0 0.5;
             1.0 1.0 0.5]

    output = ArrayOutput(init; tspan=1:3, aux=(suit=suit,))
    rule = Ruleset(DiscreteGrowthMap(layerkey=Val(:suit)))
    sim!(output, rule)

    @test output[1] == init
    @test output[2] == init.*suit
    @test output[3] == init.*suit.^2
end

@testset "exponential growth" begin
    init =  [1.0 4.0 7.0;
             2.0 5.0 8.0;
             3.0 6.0 9.0]
    test1 = init.*exp(log(2.0)*0)
    # test1 = [ 1.0  4.0  7.0;
    #           2.0  5.0  8.0;
    #           3.0  6.0  9.0]
    test2 = init.*exp(log(2.0)*1)
    # test2 = [ 2.0  8.0 14.0;
    #           4.0 10.0 16.0;
    #           6.0 12.0 18.0]
    test3 = init.*exp(log(2.0)*2)
    # test3 = [ 4.0 16.0 28.0;
    #           8.0 20.0 32.0;
    #          12.0 24.0 36.0]

    output = ArrayOutput(init; tspan=1:3)
    rule = Ruleset(ExactExponentialGrowth(intrinsicrate=log(2.0)))
    sim!(output, rule)

    @test output[1] == test1
    @test output[2] == test2
    @test output[3] == test3

    output = ArrayOutput(init; tspan=1d:1d:3d)
    rule = Ruleset(ExactExponentialGrowth(intrinsicrate=log(2.0), timestep=1d); timestep=1d)
    sim!(output, rule)
    typeof(rule.rules[1].timestep)

    @test output[1] == test1
    @test output[2] == test2
    @test output[3] == test3

    output = ArrayOutput(init; tspan=DateTime(2001,1,1):Day(5):DateTime(2001,1,15))
    rule = Ruleset(ExactExponentialGrowth(intrinsicrate=log(2.0)/5, timestep=Day(1)); timestep=Day(5))
    sim!(output, rule)

    @test output[1] == test1
    @test output[2] == test2
    @test output[3] == test3
end

@testset "exponential growth with rate from growth map" begin
    init = [1.0 1.0 1.0;
            1.0 1.0 1.0;
            1.0 1.0 1.0]

    suit =  log.([1.0 1.0 2.0;
                  2.0 1.0 0.5;
                  1.0 1.0 0.5])

    output = ArrayOutput(init; tspan=1:3, aux=(suit=suit,))
    output.extent
    rule = Ruleset(ExactExponentialGrowthMap(layerkey=Val(:suit)))
    sim!(output, rule)

    @test output[1] == [1.0 1.0 1.0;
                        1.0 1.0 1.0;
                        1.0 1.0 1.0]
    @test output[2] ≈  [1.0 1.0 2.0;
                        2.0 1.0 0.5;
                        1.0 1.0 0.5]
    @test output[3] ≈  [1.0 1.0 4.0;
                        4.0 1.0 0.25;
                        1.0 1.0 0.25]

    @test DynamicGrids.normalise.(output[3], 0.25, 4) == [0.2  0.2  1.0;
                                                          1.0  0.2  0.0;
                                                          0.2  0.2  0.0]

end

@testset "logistic growth" begin

    init =  [1.0 4.0 7.0;
             2.0 5.0 8.0;
             3.0 6.0 9.0]

    test1 = [1.0  4.0  7.0;
             2.0  5.0  8.0;
             3.0  6.0  9.0]

    test2 = [1.81818  5.71429  8.23529;
             3.33333  6.66667  8.88889;
             4.61538  7.5      9.47368]

    test3 = [3.07692  7.27273  9.03226;
             5.0      8.0      9.41176;
             6.31579  8.57143  9.72973]

    rule = Ruleset(ExactLogisticGrowth(intrinsicrate=log(2.0), carrycap=10))
    output = ArrayOutput(init; tspan=1:3)

    sim!(output, rule)

    @test output[1] == test1
    @test output[2] ≈ test2 atol=1e-4
    @test output[3] ≈ test3 atol=1e-4
end

@testset "logistic growth with rate from a growth map" begin

    init = [1.0 1.0 1.0;
            1.0 1.0 1.0;
            1.0 1.0 1.0]

    test1 = [1.0 1.0 1.0;
             1.0 1.0 1.0;
             1.0 1.0 1.0]

    test2 = [1.0     1.0  1.81818;
             1.81818 1.0  0.5    ;
             1.0     1.0  0.5    ]

    test3 = [1.0      1.0  3.07692;
             3.07692  1.0  0.25   ;
             1.0      1.0  0.25   ]

    suit =  log.([1.0 1.0 2.0;
                  2.0 1.0 0.5;
                  1.0 1.0 0.5])

    output = ArrayOutput(init; tspan=1:3, aux=(suit=suit,))
    rule = Ruleset(ExactLogisticGrowthMap(layerkey=Val(:suit), carrycap=10))
    sim!(output, rule)

    @test output[1] == test1
    @test output[2] ≈ test2 atol=1e-4
    @test output[3] ≈ test3 atol=1e-4

end

@testset "growth map masking" begin
    init = [1.0 1.0 1.0;
            1.0 1.0 1.0;
            1.0 1.0 1.0]

    suit =  [1.0 1.0 2.0;
             2.0 2.0 0.5;
             1.0 1.0 0.5]

    test1 = [1.0 1.0 1.0;
             1.0 1.0 1.0;
             1.0 1.0 1.0]

    test2 = [0.0 0.0 1.0;
             1.0 1.0 0.0;
             0.0 0.0 0.0]

    output = ArrayOutput(init; tspan=1:3, aux=(suit=suit,))
    maskrule = MaskGrowthMap(layerkey=Val(:suit), threshold=1.1)
    ruleset = Ruleset(maskrule)
    data = SimData(extent(output), ruleset)

    @test applyrule(data, maskrule, 5, (1, 1)) === 0
    @test applyrule(data, maskrule, 5, (2, 2)) === 5

    @test applyrule(data, maskrule, 5.0, (1, 1)) === 0.0
    @test applyrule(data, maskrule, 5.0, (2, 2)) === 5.0

    sim!(output, Ruleset(maskrule))

    @test output[1] == test1
    @test output[2] == test2
end

@testset "Test double layers ExactLogisticGrowthMap2" begin
    popSizeInit = [ 1.0 4.0 7.0;
                    2.0 5.0 8.0;
                    3.0 6.0 9.0]

    intrinsicRate = cat([ 1.0 1.0 1.0;
                        1.0 1.0 1.0;
                        1.0 1.0 1.0],
                        [ 2.0 2.0 2.0;
                        2.0 2.0 2.0;
                        2.0 2.0 2.0],
                        [ 1.0 1.0 1.0;
                        1.0 1.0 1.0;
                        1.0 1.0 1.0]; dims=3)

    carryingCapacity = cat([ 10.0 10.0 10.0;
                            10.0 10.0 10.0;
                            10.0 10.0 10.0],
                            [ 10.0 10.0 10.0;
                            10.0 10.0 10.0;
                            10.0 10.0 10.0],
                            [ 10.0 10.0 10.0;
                            10.0 10.0 10.0;
                            10.0 10.0 10.0]; dims=3)

    popParameter = cat(intrinsicRate, carryingCapacity; dims = 4)

    popSizeGrids = ArrayOutput(popSizeInit; tspan=1:6, aux=(popParameter=popParameter,));
    growthRule = Ruleset(ExactLogisticGrowthMap2(layerkey=:popParameter));
    sim!(popSizeGrids, growthRule);
end

# @testset "Test double layers ExactLogisticGrowthMap3" begin
#     popSizeInit = [ 1.0 4.0 7.0;
#                     2.0 5.0 8.0;
#                     3.0 6.0 9.0]

#     intrinsicRate = cat([ 1.0 1.0 1.0;
#                         1.0 1.0 1.0;
#                         1.0 1.0 1.0],
#                         [ 2.0 2.0 2.0;
#                         2.0 2.0 2.0;
#                         2.0 2.0 2.0],
#                         [ 1.0 1.0 1.0;
#                         1.0 1.0 1.0;
#                         1.0 1.0 1.0]; dims=3)

#     carryingCapacity = cat([ 10.0 10.0 10.0;
#                             10.0 10.0 10.0;
#                             10.0 10.0 10.0],
#                             [ 10.0 10.0 10.0;
#                             10.0 10.0 10.0;
#                             10.0 10.0 10.0],
#                             [ 10.0 10.0 10.0;
#                             10.0 10.0 10.0;
#                             10.0 10.0 10.0]; dims=3)

#     popSizeGrids = ArrayOutput(popSizeInit; tspan=1:6);
#     growthRule = ExactLogisticGrowthMap3(ratekey=intrinsicRate, carrycapkey=carryingCapacity);
#     sim!(popSizeGrids, growthRule);
# end

@testset "Test double grids DiscreteGrowth2" begin
    init = (pop1 = [.5 0. 0.;
                0. 0. 1.;
                0. 0. 0.],
            pop2 = [.5 0. 0.;
                0. 1. 0.;
                0. 0. 0.],)

    exposure= [1. .5 0.;
            .5 0. 0.;
            0. 0. 0.]

    output = ArrayOutput(init; tspan=1:3)
    ruleGrowth =  DiscreteGrowth2{Tuple{:pop1,:pop2},Tuple{:pop1,:pop2}}(intrinsicrate1=1.5, intrinsicrate2=2.0, carrycap = 10.0)

    hood = DispersalKernel{1}(; formulation=ExponentialKernel(15))
    ruleDispersalPop1 = InwardsPopulationDispersal{:pop1,:pop1}(;neighborhood=hood)
    ruleDispersalPop2 = InwardsPopulationDispersal{:pop2,:pop2}(;neighborhood=hood)

    rulesetGrowthDispersal = Ruleset(
     ruleGrowth, ruleDispersalPop1, ruleDispersalPop2;
     overflow=WrapOverflow()
     );

    sim!(output,rulesetGrowthDispersal);
    output[1].pop1
    output[2].pop1

    output[1].pop2
    output[2].pop2
end


text:
function pharmstr2pair(str)
    substr = split(str)
    key = Symbol(substr[length(substr)])
    val = Set{Int}()
    for i = 2:length(substr)-1
        push!(val, parse(Int,substr[i]))
    end
    return Pair{Symbol,Set{Int}}(key, val)
end

function pharmfeatures(mol::UndirectedGraph, nodes::Set{Int}=nodeset(mol))
    strarray = try attributes(mol, :PUBCHEM_PHARMACOPHORE_FEATURES)
    catch e
        ""
    end
    feats = Dict{Symbol,Vector{Set{Int}}}()
    for i=2:length(strarray)
        feat = pharmstr2pair(strarray[i])
        if feat.second ∩ nodes == feat.second
            if haskey(feats, feat.first)
                push!(feats[feat.first], feat.second)
            else
                push!(feats, Pair(feat.first, [feat.second]))
            end
            # rings are not specified as aromatic in PubChem SDF files, so check here for aromaticity
            if feat.first==:rings && sum(isaromatic(mol)[collect(feat.second)]) == length(feat.second)
                feat = Pair(:aromaticrings, feat.second)
                if haskey(feats, feat.first)
                    push!(feats[feat.first], feat.second)
                else
                    push!(feats, Pair(feat.first, [feat.second]))
                end
            end
        end
    end
    return feats
end

pharmfeatures(submol::SubgraphView) = pharmfeatures(submol.graph, nodeset(submol))

pharmfeatures!(mol::UndirectedGraph) = attributes(mol)[:pharmfeatures] = pharmfeatures(mol)

function pharmfeatures!(mol::SubgraphView)
    pmol = mol.graph
    ppf = get!(attributes(pmol), :pharmfeatures, pharmfeatures!(pmol))
    return filter(pr -> pr.first ∈ nodeset(mol), ppf)
end


text:
using DiffEqFlux, OrdinaryDiffEq, Flux, Optim, Plots, AdvancedHMC, MCMCChains
using JLD, StatsPlots
using BSON
using Random

function model_bayes(ode_data::AbstractArray,
                     numwarmup::Int = 500,
                     numsamples::Int = 500,
                     initstepsize::Float64 = 0.45,
                     odesolver = Tsit5;
                     log_dir="")
  u0 = ode_data[:, 1]
  println("u0")
  println(u0)
  datasize = length(ode_data[1,:])
  println("datasize")
  println(datasize)

  tspan = (0.0, 1.0) # (0.0, Float64(datasize) - 1)  # (0.0, Float64(datasize) - 1) # (0.0, 0.5*(Float64(datasize) - 1))
  tsteps = range(tspan[1], tspan[2], length = datasize) # range(tspan[1], tspan[2], length = datasize)

  data_dim = size(ode_data)[1]

  println("tspan")
  println(tspan)

  println("tsteps")
  println(tsteps)

  println("data_dim")
  println(data_dim)

  # ----- define Neural ODE architecture
  dudt2 = FastChain((x, p) -> x.^3,
      FastDense(data_dim, 25, swish), # FastDense(2, 50, tanh),
      FastDense(25, data_dim)) # FastDense(2, 50, tanh),
  prob_neuralode = NeuralODE(dudt2, tspan, odesolver(), saveat = tsteps) #Trapezoid

  # ----- define loss function for Neural ODE
  function predict_neuralode(p)
    Array(prob_neuralode(u0, p))
  end

  function loss_neuralode(p)
    pred = predict_neuralode(p)
    println("pred")
    println(size(pred))
    println("ode_data")
    println(size(ode_data))
    loss = sum(abs2, ode_data .- pred)
    return loss, pred
  end

  # ----- define Hamiltonian log density and gradient
  l(θ) = -sum(abs2, ode_data .- predict_neuralode(θ)) - sum(θ .* θ)

  function dldθ(θ)
    x,lambda = Flux.Zygote.pullback(l, θ)
    grad = first(lambda(1))
    return x, grad
  end

  # ----- define step size adaptor function and sampler
  metric  = DiagEuclideanMetric(length(prob_neuralode.p))

  h = Hamiltonian(metric, l, dldθ)

  integrator = Leapfrog(find_good_stepsize(h, Float64.(prob_neuralode.p)))

  prop = AdvancedHMC.NUTS{MultinomialTS, GeneralisedNoUTurn}(integrator)

  adaptor = StanHMCAdaptor(MassMatrixAdaptor(metric), StepSizeAdaptor(initstepsize, prop.integrator))

  samples, stats = sample(h, prop, Float64.(prob_neuralode.p), numwarmup, adaptor, numsamples; progress=true)

  println("samples")
  println(samples)
  println(size(samples))
  println(length(samples))
  println(size(samples[1]))
  losses = map(x-> x[1],[loss_neuralode(samples[i]) for i in 1:length(samples)])

  if log_dir != ""
    bson(joinpath(log_dir, "bayes_model_no_exo.bson"), nn_model=prob_neuralode, samples=samples, losses=losses)
  end

  (samples, losses, predict_neuralode)
end

function model_bayes_exo(non_exo_data::AbstractArray,
                         exo_data::AbstractArray,
                         numwarmup::Int = 500,
                         numsamples::Int = 500,
                         initstepsize::Float64 = 0.45,
                         odesolver = Trapezoid;
                         log_dir="")

  #@show size(non_exo_data)
  u0 = non_exo_data[:, 1]
  datasize = length(non_exo_data[1,:])
  tspan = (0.0, Float64(datasize) - 1)  # (0.0, 1.0) # (0.0, Float64(datasize) - 1)
  tsteps = range(tspan[1], tspan[2], length = datasize) # range(tspan[1], tspan[2], length = datasize)

  input_data_dim = size(non_exo_data, 1) + size(exo_data, 1)
  output_data_dim = size(non_exo_data, 1)

  # ----- define Neural ODE architecture
  dudt = FastChain((x, p) -> x.^3,
  FastDense(input_data_dim, 50, swish), # FastDense(2, 50, tanh),
  FastDense(50, output_data_dim)) # FastDense(2, 50, tanh),
  p_model = initial_params(dudt)

  function dudt_exo(u::AbstractArray{<:Float64}, p::AbstractArray{<:Float64}, t)
    dudt(vcat(u[:, 1], exo_data[:, t]), p)
  end

  function dudt2(u::AbstractArray{<:Float64}, p::AbstractArray{<:Float64}, t)
    # dudt(vcat(u[:, 1], exo_data[:, t]), p)
    # println(size(u))
    # @show vcat(u, exo_data[:, Int(round(t*(datasize - 1)/(tspan[2] - tspan[1]) + (1 - (datasize - 1)*tspan[1]/(tspan[2] - tspan[1]))))])
    dudt(vcat(u, exo_data[:, Int(round(t*(datasize - 1)/(tspan[2] - tspan[1]) + (1 - (datasize - 1)*tspan[1]/(tspan[2] - tspan[1]))))]), p)
  end

  prob = ODEProblem(dudt2, u0, tspan, nothing)

  # ----- define loss function for Neural ODE
  function predict_neuralode(p)
    _prob = remake(prob, p=p)
    Array(solve(_prob, Tsit5(), saveat=tsteps, abstol = 1e-8, reltol = 1e-6))
  end

  function loss_neuralode(p)
    pred = predict_neuralode(p)
    loss = sum(abs2, non_exo_data .- pred)
    return loss, pred
  end

  # ----- define Hamiltonian log density and gradient
  l(θ) = -sum(abs2, non_exo_data .- predict_neuralode(θ)) - sum(θ .* θ)

  function dldθ(θ)
    x, lambda = Flux.Zygote.pullback(l,θ)
    grad = first(lambda(1))
    return x, grad
  end

  # ----- define step size adaptor function and sampler
  metric = DiagEuclideanMetric(length(p_model))

  h = Hamiltonian(metric, l, dldθ)

  integrator = Leapfrog(find_good_stepsize(h, Float64.(p_model)))

  prop = AdvancedHMC.NUTS{MultinomialTS, GeneralisedNoUTurn}(integrator)

  adaptor = StanHMCAdaptor(MassMatrixAdaptor(metric), StepSizeAdaptor(initstepsize, prop.integrator))

  samples, stats = sample(h, prop, Float64.(p_model), numwarmup, adaptor, numsamples; progress=true)

  losses = map(x-> x[1],[loss_neuralode(samples[i]) for i in 1:length(samples)])

  if log_dir != ""
    bson(joinpath(log_dir, "3_bayes_model_exo$(rand).bson"), nn_model=p_model, samples=samples, losses=losses)
  end

  (samples, losses, predict_neuralode)
end

function model_bayes_exo_with_init(non_exo_data::AbstractArray,
                                   exo_data::AbstractArray,
                                   numwarmup::Int = 500,
                                   numsamples::Int = 500,
                                   initstepsize::Float64 = 0.45,
                                   odesolver = Trapezoid;
				                           initid=81,
                                   log_dir="")

  u0 = vcat(non_exo_data[:, 1], exo_data[:, 1])
  datasize = length(non_exo_data[1,:])
  tspan = (0.0, 1.0) # (0.0, Float64(datasize) - 1) # (0.0, 1.0)  # (0.0, 1.0) # (0.0, Float64(datasize) - 1)
  tsteps = range(tspan[1], tspan[2], length = datasize) # range(tspan[1], tspan[2], length = datasize)

  input_data_dim = size(non_exo_data, 1) + size(exo_data, 1)
  output_data_dim = size(non_exo_data, 1)

  # ----- define Neural ODE architecture
  # dudt = FastChain((x, p) -> x.^3,
  # FastDense(input_data_dim, 50, swish), # FastDense(2, 50, tanh),
  # FastDense(50, input_data_dim)) # FastDense(2, 50, tanh),

  # function dudt_exo(u::AbstractArray{<:Float64}, p::AbstractArray{<:Float64})
  #   dudt(vcat(u[:, 1], exo_data[:, 1]), p)
  # end

  # prob_neuralode = NeuralODE(dudt, tspan, odesolver(), saveat = tsteps) # Trapezoid

  d = BSON.load("/path/to/model")
  prob_neuralode = d[:nn_model]
  init_theta = d[:theta]

  # ----- define loss function for Neural ODE
  function predict_neuralode(p)
    Array(prob_neuralode(u0, p))
  end

  function loss_neuralode(p)
    pred = predict_neuralode(p)
    loss = sum(abs2, non_exo_data .- pred[1:output_data_dim,:])
    return loss, pred
  end

  # ----- define Hamiltonian log density and gradient
  l(θ) = -sum(abs2, non_exo_data .- predict_neuralode(θ)[1:output_data_dim]) - sum(θ .* θ)

  function dldθ(θ)
    x, lambda = Flux.Zygote.pullback(l,θ)
    grad = first(lambda(1))
    return x, grad
  end

  # ----- define step size adaptor function and sampler
  metric = DiagEuclideanMetric(length(init_theta))

  h = Hamiltonian(metric, l, dldθ)

  integrator = Leapfrog(find_good_stepsize(h, Float64.(init_theta)))

  prop = AdvancedHMC.NUTS{MultinomialTS, GeneralisedNoUTurn}(integrator)

  adaptor = StanHMCAdaptor(MassMatrixAdaptor(metric), StepSizeAdaptor(initstepsize, prop.integrator))

  samples, stats = sample(h, prop, Float64.(init_theta), numwarmup, adaptor, numsamples; progress=true)

  losses = map(x-> x[1],[loss_neuralode(samples[i]) for i in 1:length(samples)])

  num = rand(1:100)
  if log_dir != ""
    bson(joinpath(log_dir, "2_bayes_model_exo$(num).bson"), nn_model=prob_neuralode, samples=samples, losses=losses, theta=init_theta)
  end

  (samples, losses, predict_neuralode, num)
end

function find_model(index::Int)
  Core.eval(Main, :(import NNlib))

  data_directory = "/path/to/model"
  folder_names = readdir(data_directory)
  println(folder_names)
  output_dictionaries = []
  for folder in folder_names
    if "model.bson" in readdir(joinpath(data_directory, folder)) # if run has finished
      d = BSON.load(joinpath(data_directory, folder, "model.bson"))
      push!(output_dictionaries, d)
    end
  end
  sorted_dictionaries = sort(output_dictionaries, by=(x -> x[:loss]))
  sorted_dictionaries[index][:nn_model]
end

function plot_model_results(samples, losses, var::String)
  ##################### PLOTS: LOSSES ###############
  scatter(losses, ylabel = "Loss",  yscale= :log, label = "Architecture1: 500 warmup, 500 sample")

  ################### RETRODICTED PLOTS: TIME SERIES #################
  pl = scatter(tsteps, ode_data[1,:], color = :red, label = "Data: CGM", xlabel = "t", title = "CGM & $(var)")
  scatter!(tsteps, ode_data[2,:], color = :blue, label = "Data: $(var)")

  for k in 1:300
    resol = predict_neuralode(samples[100:end][rand(1:400)])
    plot!(tsteps,resol[1,:], alpha=0.1, color = :red, label = "")
    plot!(tsteps,resol[2,:], alpha=0.1, color = :blue, label = "")
  end

  idx = findmin(losses)[2]
  prediction = predict_neuralode(samples[idx])

  plot!(tsteps,prediction[1,:], color = :black, w = 2, label = "")
  plot!(tsteps,prediction[2,:], color = :black, w = 2, label = "")
  savefig(pl, "glucose_$(var)_bin_$(string(bin_size)).png")
end

function plot_model_results(samples, losses)
  ##################### PLOTS: LOSSES ###############
  scatter(losses, ylabel = "Loss",  yscale= :log, label = "Architecture1: 500 warmup, 500 sample")

  ################### RETRODICTED PLOTS: TIME SERIES #################
  pl = plot(tsteps, ode_data[1,:], color = :red, label = "Data: CGM", xlabel = "t", title = "CGM, Steps, and Bolus")
  plot!(tsteps, ode_data[2,:], color = :blue, label = "Data: Steps")
  plot!(tsteps, ode_data[3,:], color = :blue, label = "Data: Bolus")

  for k in 1:300
    resol = predict_neuralode(samples[800:end][rand(1:200)])
    plot!(tsteps,resol[1,:], alpha=0.1, color = :red, label = "")
    plot!(tsteps,resol[2,:], alpha=0.1, color = :blue, label = "")
    plot!(tsteps,resol[3,:], alpha=0.1, color = :blue, label = "")
  end

  idx = findmin(losses)[2]
  prediction = predict_neuralode(samples[idx])

  plot!(tsteps,prediction[1,:], color = :black, w = 2, label = "")
  plot!(tsteps,prediction[2,:], color = :black, w = 2, label = "")
  savefig(pl, "glucose_$(var)_bin_$(string(bin_size)).png")
end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
# @enter landmarkmatching(landmarks0, landmarksT)
# Juno.@enter landmarkmatching(landmarks0, landmarksT)

struct Info{TLT,TΣT,TμT,TL0,TΣ0,y}
     LT::TLT
     ΣT::TΣT
     μT::TμT
     L0::TL0
     Σ0::TΣ0

    function ObsInfo(LT,ΣT,μT,L0,Σ0, y)
         new{typeof(LT),typeof(ΣT),typeof(μT),typeof(L0),typeof(Σ0)}(LT,ΣT,μT,L0,Σ0*y)
    end
end

Info(1,2,3,4,5,10)

if false
    using BridgeLandmarks
    F = [(i==j) * one(UncF) for i in 1:5, j in 1:3]  # pick position indices
    F
    struct Test
        F
    end

    obj = Test(F)
    obj
    show(obj.F)
    show(obj)
    import Base.show
    function show(io::IO, mime::MIME"text/plain",obj::Test)
        print(io,mime, obj.F)
    end

    function show(io::IO, obj::Test)
        show(io,obj.F)
    end
    show(obj)
end

mutable struct MyStruct{tX, tW}
	X::tX #Vector{Bridge.SamplePath}
	W::tW
end

X = rand(2,4)
W = [rand(3)]
A = MyStruct(X,W)

function change!(A)
	A.W = [rand(10)]
end

change!(A)
A

A1 = [[1]]
A2 = copy(A1)
A1 === A2
A1[1] === A2[1]
A2[1]  =[6.0]
print(A1)

A1
ismutable(A1)               # true
A3 = deepcopy(A1)
A3 === A1
A1[1] === A3[1]

nshapes=3
xinit=BL.NState(rand(PointF,5),rand(PointF,5))
t = collect(0:.01:1.0)

X = [BL.initSamplePath(t, xinit) for _ in 1:nshapes]
W = [BL.initSamplePath(t,  zeros(PointF, 2)) for _ in 1:nshapes]
for k in 1:nshapes
	sample!(W[k], BL.Wiener{Vector{PointF}}())
end
x = deepvec(xinit)
∇x = deepcopy(x)

Xk = BL.initSamplePath(t, 0.0*xinit)
X[2] = Xk

# memory allocations, actual state at each iteration is (X,W,Q,x,∇x) (x, ∇x are initial state and its gradient)
Xᵒ = deepcopy(X)
Qᵒ = deepcopy(Q)
Wᵒ = initSamplePath(t,  zeros(StateW, dwiener))

struct Test{S,T}
	a::S
	b::T
	function TestS(v,w)
		new{typeof(v),typeof(w)}(v,w)
	end
end

mutable struct MTest{S,T}
	a::S
	b::T
	function MTestS(v,w)
		new{typeof(v),typeof(w)}(v,w)
	end
end

adjust = function(x::Test,anew)
	TestS(anew,x.b)
end

adjustm = function(x::MTestSanew)
	x.a = anew
	x
end

using BenchmarkTools

x = TestS(rand(4),[3.,4.])
xm = MTestS(rand(4),[3.,4.])

@benchmark adjust(x,rand(10_000))

@benchmark adjustm(xm,rand(10_000))

d = [[:innov, 1.0, 1], [:parameter, 1, 3], [:parameter, 0, 10], [:innov, 1.0, 6],[:rmmala_mom, 0, 4]]

updatescheme = [:innov, :mala_mom]

dfacc_names = push!(updatescheme, :iteration)
dfacc = DataFrame(fill([], length(updatescheme)),updatescheme)
push!(df, [1,0,3])

DataFrames.names!(df, updatescheme)

updatescheme = [:innov, :mala_mom]
accinfo = DataFrame(fill([], length(updatescheme)+1),push!(updatescheme,:iteration))

function adaptpcnstep!(ρ, n, accinfo, nshapes, η; adaptskip = 15, targetaccept=0.5)

    recentmean = mean(accinfo[end-adaptskip+1:end])
        # ind1 =  findall(first.(accinfo).==:innov)[end-adaptskip*nshapes+1:end]
        # recent_mean = mean(last.(accinfo)[ind1])
        if recentmean > targetaccept
            u = sigmoid(invsigmoid(ρ) - η(n))
        else
            u = sigmoid(invsigmoid(ρ) + η(n))
        end
		ρ = u
    nothing
end

sigmoid(z::Real) = 1.0 / (1.0 + exp(-z))
invsigmoid(z::Real) = log(z/(1-z))

adaptskip = 15; targetaccept=0.5
η = n -> min(0.2, 10/n)                          # stepsize cooling function
nshapes = 1
accinfo = rand(Bernoulli(0.2),30)
n = 5

ρ = 0.9
adaptpcnstep!(ρ, n, accinfo, nshapes, η)
ρ

struct MyS{T}
	x::T
	y

	function MyS(ξ)
		y = sum(ξ)
		new{typeof(ξ)}(ξ,y)
	end
end
x = State(rand(PointF,3),zeros(PointF,3))
deepx = deepvec(x)
BL.deepvec2state(deepx)-x # check
x.p
x.q

MyS(1:10)
MyS(rand(6))

x = State(rand(PointF,3),zeros(PointF,3))
deepx = deepvec(x)
BL.deepvec2state(deepx)-x # check
x.p
x.q
BL.deepvec(x.q)
vec(x)

testfun = function(p,q)
	norm(p+q)
end

testfun2(q)  = (p) -> testfun(p,q)
x.x[1,:]
ForwardDiff.gradient(testfun2(x.q), x.p)

using BridgeLandmarks

function split_state(x::BL.NState)
	q =  vec(reinterpret(BL.deepeltype(x), x.x[1,:]))
	p =  vec(reinterpret(BL.deepeltype(x), x.x[2,:]))
	q, p
end

y = State(rand(PointF,3),rand(PointF,3))
q, p = split_state(y)
yy = merge_state(q,p)
yy-y

function merge_state(q,p)
	BL.NState(reinterpret(Point{eltype(q)}, q),reinterpret(Point{eltype(p)}, p))
end


text:
using LinearAlgebra

function decomp_qr(A)

    (m,n) = size(A);
    R = zeros(n,n);
    Q = zeros(m,n);

    for j = 1:n

        soma = zeros(m,1);

        for i = 1 : j-1
            R[i,j] = Q[:,i]' * A[:,j];
            soma = soma + R[i,j]*Q[:,i];
        end

        v = A[:,j] - soma;
        R[j,j] = norm(v);
        Q[:,j] = v/R[j,j];

    end

    return (Q,R);

end


text:
import HarmonicBalance.is_harmonic

@variables a,b,c,t,x(t),f, y(t)

@test is_harmonic(cos(f*t), t)
@test is_harmonic(1, t)
@test !is_harmonic(cos(f*t^2 + a), t)
@test !is_harmonic(a + t, t)

dEOM = DifferentialEquation([a + x, t^2 + cos(t)], [x, y])
@test !is_harmonic(dEOM, t)


text:
import Flamenco
import HTTP
import JSON

include("./request.jl")

function make_and_test_request_delete(headers::Vector{Pair{String,String}}) :: Nothing
    r = HTTP.request("DELETE", "http://127.0.0.1:8000/", headers)
    parse_and_test_request(r, 200, headers, 15, true, false, Dict("Success" => 200))
end

function make_and_test_request_get(headers::Vector{Pair{String,String}}) :: Nothing
    r = HTTP.request("GET", "http://127.0.0.1:8000/", headers)
    parse_and_test_request(r, 200, headers, 230, false, true)
end

function make_and_test_request_options(headers::Vector{Pair{String,String}}) :: Nothing
    r = HTTP.request("OPTIONS", "http://127.0.0.1:8000/", headers)
    parse_and_test_request(r, 200, headers, 2, true, false, Dict())
    compare_http_header(r.headers, "Allow", "POST")
    compare_http_header(r.headers, "Access-Control-Allow-Methods", "POST")
    compare_http_header(r.headers, "Access-Control-Allow-Headers", "X-PINGOTHER, Content-Type")
end

function make_and_test_request_put(headers::Vector{Pair{String,String}}, payload::Dict, content_length::Int64) :: Nothing
    r = HTTP.request("PUT", "http://127.0.0.1:8000/", headers, JSON.json(payload))
    parse_and_test_request(r, 200, headers, content_length, true, false, payload)
end

function make_and_test_request_post(headers::Vector{Pair{String,String}}, payload::Dict, content_length::Int64) :: Nothing
    r = HTTP.request("POST", "http://127.0.0.1:8000/", headers, JSON.json(payload))
    parse_and_test_request(r, 200, headers, content_length, true, false, payload)
end

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

function routes(; request_headers::Dict{String, String}, request_method::String, request_path::String, request_payload::String) :: String
    body::String = "a"
    headers::Dict{String, String} = Dict()
    status_code::Int16 = 404

    if (request_path=="/")
        status_code = 200
        delete!(request_headers, "Content-Length")
        delete!(request_headers, "Date")

        if (request_method=="DELETE")
            headers = merge(
                request_headers,
                Dict("Content-Type" => "application/json")
            )
            body = JSON.json(Dict(
                "Success" => 200
            ))
        elseif (request_method=="GET")
            headers = merge(
                request_headers,
                Dict("Content-Type" => "text/html; charset=UTF-8")
            )
            body = """<!DOCTYPE html>

                    <meta http-equiv="Content-type" content="text/html; charset=utf-8">
                    <title>Flamenco</title>

Content ...

            """
        elseif (request_method=="OPTIONS")
            headers = merge(
                request_headers,
                Dict(
                    "Allow" => "POST",
                    "Access-Control-Allow-Methods" => "POST",
                    "Access-Control-Allow-Headers" => "X-PINGOTHER, Content-Type",
                    "Content-Type" => "application/json"
                )
            )
            body = JSON.json(Dict())
        elseif (request_method=="PUT" || request_method=="POST")
            headers = merge(
                request_headers,
                Dict(
                    "Content-Type" => "application/json"
                )
            )
            body = request_payload
        end
    end

    return Flamenco.write_response(status_code, headers; body=body)
end

function server_close() :: Nothing
    Flamenco.close_server()
end

function server_start()
    @async Flamenco.start_server("127.0.0.1", 8000) do request::Flamenco.Server.Request
        routes(; request_headers=request.headers, request_method=request.method, request_path=request.target, request_payload=request.body)
    end
end


text:
"""
Computes T-Distribution PDF and CDF stably using built-in function.
Returns:
    p(z0 | theta, lambda, z)
    P(z0 | theta, lambda, z)
"""
function likelihood(θ, λ, train, test, theta_params::Union{θ_params{Array{Float64, 2},
    Cholesky{Float64,Array{Float64, 2}}}, Nothing}=nothing, type = "Gaussian")
    s = train.s; s0 = test.s0; X = train.X; X0 = test.X0; z = train.z; n = size(X, 1); p = size(X, 2); k = size(X0, 1)  #unpack
    if theta_params==nothing; theta_params = funcθ(θ, train, test, type); end
    g = boxCox #boxCox by default
    Bθ = theta_params.Bθ
    Cθ = theta_params.Cθ
    Hθ = theta_params.Hθ
    Σθ_inv_X = theta_params.Σθ_inv_X
    choleskyΣθ = theta_params.choleskyΣθ
    βhat = (X'*(choleskyΣθ\X))\(X'*(choleskyΣθ\g(z, λ)))
    qtilde = (expr = g(z, λ)-X*βhat; expr'*(choleskyΣθ\expr))
    m = Bθ*(choleskyΣθ\g(z, λ)) + Hθ*βhat
    t = LocationScale(m[1], sqrt(qtilde[1]*Cθ[1]/(n-p)), TDist(n-p))
    return (z0 -> Distributions.pdf(t, g(z0, λ)), z0 -> Distributions.cdf(t, g(z0, λ)))
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
"""
    kron!(A::GBMatrix, B::GBMatrix, op = BinaryOps.TIMES; kwargs...)::GBMatrix

In-place version of [kron](@ref).
"""
function LinearAlgebra.kron!(
    C::GBMatOrTranspose,
    A::GBMatOrTranspose,
    B::GBMatOrTranspose,
    op::BinaryUnion = BinaryOps.TIMES;
    mask = nothing,
    accum = nothing,
    desc = nothing
)
    mask, accum = _handlenothings(mask, accum)
    desc === nothing && (desc = DEFAULTDESC)
    op = getoperator(op, optype(A, B))
    A, desc, B = _handletranspose(A, desc, B)
    accum = getaccum(accum, eltype(C))
    if op isa TypedBinaryOperator
        libgb.GxB_kron(C, mask, accum, op, A, B, desc)
    elseif op isa TypedMonoid
        libgb.GrB_Matrix_kronecker_Monoid(C, mask, accum, op, A, B, desc)
    elseif op isa TypedSemiring
        libgb.GrB_Matrix_kronecker_Semiring(C, mask, accum, op, A, B, desc)
    else
        throw(ArgumentError("$op is not a valid monoid binary op or semiring."))
    end
    return C
end

function LinearAlgebra.kron!(
    C::GBMatOrTranspose,
    A::GBMatOrTranspose,
    B::GBMatOrTranspose,
    op::Function;
    mask = nothing,
    accum = nothing,
    desc = nothing
)
    return kron!(C, A, B, BinaryOp(op); mask, accum, desc)
end

"""
   kron(A::GBMatrix, B::GBMatrix, op = BinaryOps.TIMES; kwargs...)::GBMatrix

Kronecker product of two matrices using `op` as the multiplication operator.
Does not support `GBVector`s at this time.

# Arguments
- `A::GBMatrix`: optionally transposed.
- `B::GBMatrix`: optionally transposed.
- `op::MonoidBinaryOrRig = BinaryOps.TIMES`: the binary operation which replaces the arithmetic
    multiplication operation from the usual kron function.

# Keywords
- `mask::Union{Nothing, GBMatrix} = nothing`: optional mask.
- `accum::Union{Nothing, AbstractBinaryOp} = nothing`: binary accumulator operation
    where `C[i,j] = accum(C[i,j], T[i,j])` where T is the result of this function before accum is applied.
- `desc = nothing`
"""
function LinearAlgebra.kron(
    A::GBMatOrTranspose,
    B::GBMatOrTranspose,
    op::BinaryUnion = BinaryOps.TIMES;
    mask = nothing,
    accum = nothing,
    desc = nothing
)
    t = inferoutputtype(A, B, op)
    C = GBMatrix{t}(size(A,1) * size(B, 1), size(A, 2) * size(B, 2))
    kron!(C, A, B, op; mask, accum, desc)
    return C
end

function LinearAlgebra.kron(
    A::GBMatOrTranspose,
    B::GBMatOrTranspose,
    op::Function;
    mask = nothing,
    accum = nothing,
    desc = nothing
)
    return kron(A, B, BinaryOp(op); mask, accum, desc)
end


text:
@DeepShiba_Func mutable struct Reshape <: Func
    shape::Tuple{Int}
end

@DeepShiba_Func mutable struct Transpose <: Func end

_Reshape(x, shape) = Reshape([x, ], nothing, nothing, shape)(x)
_Transpose(x) = Transpose([x, ], nothing, nothing)(x)

forward(f::Reshape, x)  = reshape(x, f.shape)
forward(f::Transpose, x) = transpose(x)

function backward(f::Reshape)::AbstractArray{Variable}
    gys = get_gys(f.outputs)
    x_shape = size(f.inputs[1])
    return reshape.(gys, Ref(x_shape))
end

function backward(f::Transpose)::AbstractArray{Variable}
    gys = get_gys(f.outputs)
    return transpose.(gys)
end

Base.sin(x::Variable) = _Sin(x)
Base.cos(x::Variable) = _Cos(x)
Base.tan(x::Variable) = _Tan(x)
Base.tanh(x::Variable) = _Tanh(x)

Base.reshape(x::Variable, shape::Variable) = _Reshape(x, shape.data)
Base.reshape(x::Variable, shape) = _Reshape(x, shape)
Base.reshape(x, shape::Variable) = _Reshape(variable(x), shape.data)
Base.transpose(x::Variable) = _Transpose(x)


text: download("https://cdn.plot.ly/plotly-2.8.3.min.js", joinpath(@__DIR__, "plotly-2.8.3.min.js"))


text:
"""Main DwsimOpt module.
"""
module DwsimOpt

# importing packages
using Test

# export packages
# export OptProblemDef

# includes
# include("SimOpt.jl")

end # module