text:
include("twihesexp.jl")

############################
#
# Parameters
#
############################

#Searches for a suitable f
function findf(la, ea, bound)
    for i in 1:bound
        p = la^ea*i*3-1
        if p%9 == 5 && isprime(p)
            return i*3
        end
    end
end

#Prime and field
const la = big(2)
const ea = 63
const eb = 41
const f = big(165)
const p = la^ea*f-1
const F = Fp2(p)

#cuberoot of One
const w = F(big(0),big(1))
if !isone(w^3)
    error("w is not a cuberoot...")
end

#Starting curve: a = 1 for speed, solve for d to get j(E) = some supersingular j-invariant. See sage-script
#Starting j-invariant = 681983795137470565003 + 525699670422529264746*w
#See sage script for construction of d0 from j-Invariant and a0
const a0 = one(F)
const d0 = F(big(51849466119744637087), big(763718193217644493204))

#Generates torsion basises from random points
function genTorsionBasis(alicePX, alicePY, aliceQX, aliceQY)
    Z = one(F)
    aliceP = EC(alicePX, alicePY, Z)
    aliceQ = EC(aliceQX, aliceQY, Z)

    if !isValid(aliceP, a0, d0)
        error("Alice' P is not a valid point!")
    end

    if !isValid(aliceQ, a0, d0)
        error("Alice' Q is not a valid point!")
    end

    Pa = normalized(mul(f, aliceP, a0))
    Qa = normalized(mul(f, aliceQ, a0))

    Patest = doublePow(Pa, a0, ea-1)
    Qatest = doublePow(Qa, a0, ea-1)

    if !isIdentity(double(Patest, a0))
        error("Pa : Something does not look right with the order of the curve...")
    end

    if !isIdentity(double(Qatest, a0))
        error("Qa : Something does not look right with the order of the curve...")
    end

    if isIdentity(Patest)
        error("Alice' point Pa does not have the right order!")
    end
    if isIdentity(Qatest)
        error("Alice' point Qa does not have the right order!")
    end

    if Patest == Qatest
        error("Alice' points are not independent!")
    end

    return Pa, Qa
end

#See sage script for generating random points
PX1 = F(big(177014855531225039999), big(1202619811203149879755))
PY1 = F(big(312237303984514937680), big(1137280334077319004711))
QX1 = F(big(1328653618852880321918), big(671857271963059083134))
QY1 = F(big(905775967463901689844), big(260557785141733845323))
const Pa, Qa = genTorsionBasis(PX1, PY1, QX1, QY1)

# Compute the optimal strategy of size n, given weights p and q
function computeStrat(n, p, q)
    S = Dict{Int64, Array{Int64, 1}}()
    C = Dict{Int64, Int64}()
    S[1] = []
    C[1] = 0
    for i in 2:n
        b = argmin([C[i-b] + C[b] + b*p + (i-b)*q for b in 1:i-1])
        S[i] = append!(append!([b], S[i-b]), S[b])
        C[i] = C[i-b] + C[b] + b*p + (i-b)*q
    end
    return S[n]
end

#The weights are estimated using @benchmark
const strat_A = Vector(computeStrat(ea, 23358, 32132))


text:
K = 0.155 # [1/s]
T = 80.5 # [s]

start_time_second = 0.0
time_second_interval = 0.01
end_time_second = 200.0
time_list = start_time_second:time_second_interval:end_time_second

@testset "kt.jl" begin
    Ts = 50.0
    δ_list = 10.0 * pi / 180.0 * sin.(2.0 * pi / Ts * time_list) # [rad]
    kt_results = kt_simulate(K, T, time_list, δ_list)
end

@testset "kt_zigzag_test" begin
    target_δ_rad = 10.0 * π / 180.0
    target_ψ_rad_deviation = 10.0 * π / 180.0
    r_list, ψ_list, δ_list =
        kt_zigzag_test(K, T, time_list, target_δ_rad, target_ψ_rad_deviation)
end

@testset "kt least square method func" begin
    Ts = 50.0
    δ_list = 10.0 * pi / 180.0 * sin.(2.0 * pi / Ts * time_list) # [rad]
    r, δ = kt_simulate(K, T, time_list, δ_list)

    data = ShipData(time_list, 0, 0, r, 0, 0, 0, δ, 0)
    K_est, T_est = estimate_kt_lsm(data)
    @test abs(K - K_est) < 1.0
    @test abs(T - T_est) < 10.0
end

@testset "kt bootstrap least square method func" begin
    Ts = 50.0
    δ_list = 10.0 * pi / 180.0 * sin.(2.0 * pi / Ts * time_list) # [rad]
    r, δ = kt_simulate(K, T, time_list, δ_list)
    noize_dist = Normal(0.0, 0.0005)
    r_obs = r + rand(noize_dist, size(r))
    data = ShipData(time_list, 0, 0, r_obs, 0, 0, 0, δ, 0)
    one_sample_size = 5000
    K_est_samples, T_est_samples =
        estimate_kt_lsm_time_window_sampling(data, one_sample_size)
end

@testset "kt mcmc sampling func" begin
    Ts = 50.0
    δ_list = 10.0 * pi / 180.0 * sin.(2.0 * pi / Ts * time_list) # [rad]
    r, δ = kt_simulate(K, T, time_list, δ_list)
    sampling_rate = 10
    time_obs = time_list[1:sampling_rate:end]
    noize_dist = Normal(0.0, 0.005)
    r_w_noize = r + rand(noize_dist, size(r))
    r_obs = r_w_noize[1:sampling_rate:end]
    δ_obs = δ[1:sampling_rate:end]
    data = ShipData(time_obs, 0, 0, r_obs, 0, 0, 0, δ_obs, 0)
    n_samples = 10
    n_chains = 1
    model = create_model_for_mcmc_sample_kt(data)
    chain = nuts_sampling_single_thread(model, n_samples, n_chains)
end

ext: jl

max_stars_repo_name: mitsuyukiLab/JuliaShipMMG

max_issues_repo_name: mitsuyukiLab/JuliaShipMMG


text:
using Pkg
using Pollen

using FluxTraining
const PACKAGE = FluxTraining

LAZY = get(ENV, "POLLEN_LAZY", "false") == "true"

# Create Project
m = PACKAGE
ms = [m, m.Events, m.Phases, m.Loggables]

@info "Creating project..."
project = Project(
    Pollen.Rewriter[
        Pollen.DocumentFolder(pkgdir(m), prefix = "documents"),
        Pollen.ParseCode(),
        Pollen.ExecuteCode(),
        Pollen.PackageDocumentation(ms),
        Pollen.DocumentGraph(),
        Pollen.SearchIndex(),
        Pollen.SaveAttributes((:title,)),
        Pollen.LoadFrontendConfig(pkgdir(m))
    ],
)

DIR = mktempdir()
@info "Serving from directory $DIR"
Pollen.serve(project, DIR, lazy=LAZY, format = Pollen.JSON())


text:
## Indices

function Base.filter(pred, inds::AbstractIndices)
    out = empty(inds)
    for i in inds
        if pred(i)
            insert!(out, i)
        end
    end
    return out
end

# Occasionally it might be faster to copy all and remove just a few indices...
# TODO: export `unsafe_filter`, or hook up to certain datatypes or something
function unsafe_filter(pred, inds::AbstractIndices)
    out = copy(inds)
    Base.unsafe_filter!(pred, out)
    return out
end

struct FilteredIndices{I, Parent <: AbstractIndices{I}, Pred} <: AbstractIndices{I}
    parent::Parent
    pred::Pred
end

Base.IteratorSize(::FilteredIndices) = Base.SizeUnknown # You can still call `length`, but it is slower
Base.length(inds::FilteredIndices) = count(inds.pred, inds.parent)

function Base.in(i::I, inds::FilteredIndices{I}) where {I}
    return i in inds.parent && inds.pred(i)
end

function Base.iterate(inds::FilteredIndices{I}, s...) where {I}
    tmp = iterate(inds.parent, s...)
    tmp === nothing && return nothing
    (i::I, s2) = tmp
    if inds.pred(i)
        return (i, s2)
    else
        return iterate(inds, s2)
    end
end

function filterview(pred, inds::AbstractIndices{I}) where {I}
    return FilteredIndices{I, typeof(inds), typeof(pred)}(inds, pred)
end

## Dictionaries

function Base.filter(pred, dict::AbstractDictionary)
    out = empty(dict)
    for (i, v) in pairs(dict)
        if pred(v)
            insert!(out, i, v)
        end
    end
    return out
end

# Occasionally it might be faster to copy all and remove just a few indices...
# TODO: export `unsafe_filter`, or hook up to certain datatypes or something
function unsafe_filter(pred, dict::AbstractDictionary)
    out = copy(dict)
    Base.unsafe_filter!(pred, dict)
    return out
end

struct FilteredDictionary{I, V, Parent <: AbstractDictionary{I, V}, Pred} <: AbstractDictionary{I, V}
    parent::Parent
    pred::Pred
end

Base.IteratorSize(::FilteredDictionary) = Base.SizeUnknown # You can still call `length`, but it is slower

@inline function Base.keys(dict::FilteredDictionary)
    return filterview(i -> dict.pred(dict.parent[i]), keys(dict.parent))
end

@propagate_inbounds function Base.getindex(dict::FilteredDictionary{I}, i::I) where {I}
    parent = dict.parent
    (hastoken, token) = gettoken(parent, i)
    @boundscheck if !hastoken
        throw(IndexError("Dictionary does not contain index: $i"))
    end
    out = gettokenvalue(parent, token)
    @boundscheck if !dict.pred(out)
        throw(IndexError("Dictionary does not contain index: $i"))
    end
    return out
end

function filterview(pred, inds::AbstractDictionary{I, T}) where {I, T}
    return FilteredDictionary{I, T, typeof(inds), typeof(pred)}(inds, pred)
end

Base.similar(dict::FilteredDictionary, ::Type{T}, indices) where {T} = similar(dict.parent, T, indices)
Base.empty(dict::FilteredDictionary, ::Type{T}) where {T} = similar(dict.parent, T)


text:
function hyperbolic_cross_polynomial(node::AbstractArray{R,1},multi_index::Array{S,2}) where {T<:AbstractFloat,R<:Number,S<:Integer}

  n = 2*maximum(multi_index,dims=1).+1
  d = length(n) # d is the number of dimensions

  base_polynomials        = Array{Array{R,2},1}(undef,d)
  unique_base_polynomials = Array{Array{R,1},2}(undef,size(multi_index))

  # Construct the base polynomials

  @inbounds for i = 1:d
    base_polynomials[i] = chebyshev_polynomial(n[i]-1,node[i])
  end

  # Compute the unique polynomial terms from the base polynomials

  @inbounds for i in axes(multi_index,1)
    for j = 1:d
      if multi_index[i,j] == 0
        unique_base_polynomials[i,j] = [base_polynomials[j][1]]
      else
        unique_base_polynomials[i,j] = [base_polynomials[j][2*multi_index[i,j]], base_polynomials[j][2*multi_index[i,j]+1]]
      end
    end
  end

  # Construct the first row of the interplation matrix

  n = determine_grid_size(multi_index)
  polynomial = Array{R,1}(undef,n[1])
  l = 1
  @inbounds for j in axes(multi_index,1)
    new_polynomials = unique_base_polynomials[j,1]
    for i = 2:d
      new_polynomials = kron(new_polynomials,unique_base_polynomials[j,i])
    end
    m = length(new_polynomials)
    polynomial[l:l+m-1] = new_polynomials
    l += m
  end

  return polynomial

end

function hyperbolic_cross_polynomial(node::AbstractArray{R,1},multi_index::Array{S,2},domain::Union{Array{T,1},Array{T,2}}) where {T<:AbstractFloat,R<:Number,S<:Integer}

  if size(domain,2) != length(node)
    error("domain is inconsistent with the number of dimensions")
  end

  node = copy(node)

  d = length(node)
  for i = 1:d
    node[i] = normalize_node(node[i],domain[:,i])
  end

  polynomial = hyperbolic_cross_polynomial(node,multi_index)

  return polynomial

end

function hyperbolic_cross_evaluate(weights::Array{T,1},node::AbstractArray{R,1},multi_index::Array{S,2}) where {T<:AbstractFloat,R<:Number,S<:Integer}

  poly = hyperbolic_cross_polynomial(node,multi_index)

  estimate = weights'poly

  return estimate

end

function hyperbolic_cross_evaluate(weights::Array{T,1},node::AbstractArray{R,1},multi_index::Array{S,2},domain::Union{Array{T,1},Array{T,2}}) where {T<:AbstractFloat,R<:Number,S<:Integer}

if size(domain,2) != length(node)
  error("domain is inconsistent with the number of dimensions")
end

node = copy(node)

d = length(node)
for i = 1:d
  node[i] = normalize_node(node[i],domain[:,i])
end

estimate = hyperbolic_cross_evaluate(weights,node,multi_index)

return estimate

end

function hyperbolic_cross_evaluate(weights::Array{T,1},polynomial::Array{R,1}) where {T<:AbstractFloat,R<:Number}

  estimate = weights'polynomial

  return estimate

end

function hyperbolic_cross_evaluate(weights::Array{T,1},multi_index::Array{S,2}) where {T<:AbstractFloat,S<:Integer}

function goo(x::Array{R,1}) where {R<:Number}

  return hyperbolic_cross_evaluate(weights,x,multi_index)

end

return goo

end

function hyperbolic_cross_evaluate(weights::Array{T,1},multi_index::Array{S,2},domain::Union{Array{T,1},Array{T,2}}) where {T<:AbstractFloat,S<:Integer}

function goo(x::Array{R,1}) where {R<:Number}

  return hyperbolic_cross_evaluate(weights,x,multi_index,domain)

end

return goo

end


text:
using Thrift
import Thrift.process, Thrift.meta

# include the generated module, which in-turn includes our implementation code in `hello_impl.jl`
include("gen-jl/hello/hello.jl");
import hello.SayHelloProcessor

# create a server instance with our choice of protocol and transport
srvr_processor = SayHelloProcessor()
srvr_transport = TServerSocket(19999)
srvr = TSimpleServer(srvr_transport, srvr_processor, x->x, x->TBinaryProtocol(x), x->x, x->TBinaryProtocol(x))

# start serving client requests
println("starting to serve requests...")
serve(srvr)


text:
function thres_show_cb(widget::Ptr,user_data::Tuple{SpikeSorting.Single_Channel})

    sc, = user_data
    mywidget = convert(CheckButton, widget)
    sc.show_thres=getproperty(mywidget,:active,Bool)
    sc.old_thres=getproperty(sc.adj_thres,:value,Int)
    sc.thres=getproperty(sc.adj_thres,:value,Int)

    nothing
end

function thres_cb(widget::Ptr,user_data::Tuple{SpikeSorting.Single_Channel})

    sc,  = user_data

    mythres=getproperty(sc.adj_thres,:value,Int)
    setproperty!(sc.thres_widgets.sb,:label,string(mythres))
    sc.thres_changed=true

    nothing
end

function gain_check_cb(widget::Ptr,user_data::Tuple{SpikeSorting.Single_Channel})

    sc, = user_data

    mygain=getproperty(sc.gain_widgets.multiply,:active,Bool)

    if mygain
        Gtk.GAccessor.increments(sc.gain_widgets.gainbox,10,10)
    else
        Gtk.GAccessor.increments(sc.gain_widgets.gainbox,1,1)
    end

    nothing
end

function sb2_cb(widget::Ptr,user_data::Tuple{Single_Channel})

    sc, = user_data

    mygain=getproperty(sc.gain_widgets.all,:active,Bool)

    gainval=getproperty(sc.gain_widgets.gainbox,:value,Int)
    mythres=getproperty(sc.adj_thres,:value,Int)

    if mygain==true
    else
    end

    sc.s = -1 .* gainval/1000

    sc.thres_changed=true

    nothing
end


text:
#=

    Model.jl

    A Julia module containing code to solve and simulate the model.

    Author:     Alexander P. Browning
                ======================
                School of Mathematical Sciences
                Queensland University of Technology
                ======================
                ap.browning@icloud.com
                alexbrowning.me

=#
module Model

    using ColorSchemes
    using Distributions
    using Interpolations
    using KernelDensity
    using LinearAlgebra
    using MCMCChains
    using Plots
    using Random
    using SpecialFunctions
    using Statistics
    using StatsBase
    using StatsFuns
    using StatsPlots

    export solve_ode_model, solve_ode_model_all_vars, solve_alt_ode_model
    export GammaAlt, LogNormalAlt, GaussianCopula, TCopula, Copula, MvDependent, quantile_interpolated, vec_to_cor, marginalize
    export plot_fit_univariate, plot_fit_bivariate, plot_distribution, plot_distributions, plot_biv_distribution, plot_biv_distributions
    export simulate_model, simulate_model_noiseless, parameters_to_antibody
    export col_Q̄,col_Ū,col_Q,col_U,col_Q̄Ū,col_QU
    export pdf_ci

    col_Q̄ = "#F59640"  # Quenchable, quenched sample
    col_Ū = "#1BA7AD"  # Unquenchable, quenched sample
    col_Q = "#F54D3A"  # Quenchable, unquenched sample
    col_U = "#2339AC"  # Unquenchable, unquenched sample
    col_Q̄Ū = "#c56f88"
    col_QU = "#a8216b"

    include("deterministic.jl")
    include("distributions.jl")
    include("plotting.jl")
    include("statistical.jl")

end

max_stars_repo_name: ap-browning/internalisation

max_issues_repo_name: ap-browning/internalisation

max_forks_repo_name: ap-browning/internalisation


text:
"""
Functionality for creating and working with Genie plugins.
"""
module Plugins

import Genie
import Pkg, Markdown, Logging

const FILES_FOLDER = "files"
const PLUGINS_FOLDER = Genie.config.path_plugins
const TASKS_FOLDER = Genie.config.path_tasks
const APP_FOLDER = Genie.config.path_app

const path_prefix = joinpath(@__DIR__, "..", FILES_FOLDER, "new_app") |> normpath |> relpath
const FOLDERS = [ joinpath(path_prefix, APP_FOLDER),
                  joinpath(path_prefix, "db"),
                  joinpath(path_prefix, PLUGINS_FOLDER),
                  joinpath(path_prefix, Genie.config.server_document_root) ]

"""
    recursive_copy(path::String, dest::String; only_hidden = true, force = false)

Utility function to copy plugin files from package dir to app.
"""
function recursive_copy(path::String, dest::String; only_hidden = true, force = false)
  for (root, dirs, files) in walkdir(path)
    dest_path = joinpath(dest, replace(root, (path_prefix * '/')=>""))

    try
      mkpath(dest_path)
      @info "Created dir $dest_path"
    catch ex
      @error "Failed to create dir $dest_path"
    end

    for f in files
      (only_hidden && startswith(f, ".")) || continue # only copy hidden files, especially .gitkeep
      try
        cp(joinpath(root, f), joinpath(dest_path, f), force = force)
        @info "Copied $(joinpath(root, f)) to $(joinpath(dest_path, f))"
      catch ex
        @error "Failed to copy $(joinpath(root, f)) to $(joinpath(dest_path, f))"
      end
    end
  end
end

"""
    congrats()

Shows success message and instructions when scaffolding a plugin.
"""
function congrats()
  message = """
  Congratulations, your plugin is ready!
  You can use this default installation function in your plugin's module:
  """
  print(message)

  Markdown.doc"""
  ```julia
  function install(dest::String; force = false) :: Nothing
    src = abspath(normpath(joinpath(@__DIR__, "..", Genie.Plugins.FILES_FOLDER)))

    for f in readdir(src)
      isfile(f) && continue
      isdir(f) || mkpath(joinpath(src, f))

      Genie.Plugins.install(joinpath(src, f), dest, force = force)
    end

    nothing
  end
  ```
  """
end

"""
    scaffold(plugin_name::String, dest::String = "."; force = false)

Scaffolds a new plugin as a Julia project
"""
function scaffold(plugin_name::String, dest::String = "."; force = false)
  plugin_name = replace(plugin_name, " "=>"") |> strip |> string
  dest = normpath(dest) |> abspath
  ispath(dest) || mkpath(dest)

  @info "Generating project file"
  cd(dest)
  Pkg.generate(plugin_name)
  dest = joinpath(dest, plugin_name)

  @info "Scaffolding file structure"
  mkpath(joinpath(dest, FILES_FOLDER))

  for path in FOLDERS
    recursive_copy(path, joinpath(dest, FILES_FOLDER), force = force)
  end

  initializer_path = joinpath(dest, FILES_FOLDER, PLUGINS_FOLDER, lowercase(plugin_name) * ".jl")
  @info "Creating plugin initializer at $initializer_path"
  touch(initializer_path)

  @info "Adding dependencies"

  cd(dest)
  Pkg.activate(".")
  Pkg.add("Genie")

  run(`git init`)
  run(`git add .`)
  run(`git commit -am "initial commit"`)

  congrats()
end

"""
    install(path::String, dest::String; force = false)

Utility to allow users to install a plugin
"""
function install(path::String, dest::String; force = false)
  isdir(dest) || mkpath(dest)
  cd(dest)
  isdir(Genie.config.path_plugins) || mkpath(Genie.config.path_plugins)

  root_length = splitpath(path) |> length
  depth = 0

  for (root, dirs, files) in walkdir(path)
    depth = length(splitpath(root)) - root_length
    dest_path = joinpath(abspath(dest), splitpath(root)[end-depth:end]...)

    try
      mkpath(dest_path)
      @info "Created dir $dest_path"
    catch ex
      @error "Did not create dir $dest_path"
    end

    for f in files
      try
        cp(joinpath(root, f), joinpath(dest_path, f), force = force)
        @info "Copied $(joinpath(root, f)) to $(joinpath(dest_path, f))"
      catch ex
        @error "Did not copy $(joinpath(root, f)) to $(joinpath(dest_path, f))"
      end
    end
  end
end

end


text:
module CellQuadraturesTests

using Test
using FillArrays
using Gridap.Helpers
using Gridap.Fields
using Gridap.ReferenceFEs
using Gridap.Arrays
using Gridap.CellData
using Gridap.TensorValues
using Gridap.Geometry

function test_old_constructor(trian::Triangulation,degree::Int)
  cell_quad = CellQuadrature(trian,degree)
  CellQuadrature(cell_quad.cell_quad,
             cell_quad.cell_point,
             cell_quad.cell_weight,
             cell_quad.trian,
             cell_quad.data_domain_style)
end

domain = (0,1,0,1,0,1)
cells = (4,4,4)
model = simplexify(CartesianDiscreteModel(domain,cells))

trian = Triangulation(model)

degree = 2
quad = test_old_constructor(trian,degree)

x = get_cell_points(quad)
@test isa(x,CellPoint)

v = GenericCellField(get_cell_shapefuns(trian),trian,ReferenceDomain())
u = GenericCellField(lazy_map(transpose,get_data(v)),v.trian,v.domain_style)

m = ∇(v)⋅∇(u)
test_array(m(x),collect(m(x)))

s = integrate(m,quad)
test_array(s,collect(s))

s = ∫( ∇(v)⋅∇(u) )*quad
test_array(s,collect(s))

s = ∫(1)*quad
@test sum(s) ≈ 1
@test ∑(s) ≈ 1

trian_N = BoundaryTriangulation(model)
quad_N = CellQuadrature(trian_N,degree)

s = ∫( ∇(v)⋅∇(u) )*quad_N
test_array(s,collect(s))

s = ∫(1)*quad_N
@test sum(s) ≈ 6
@test ∑(s) ≈ 6

s = ∫( x->1 )*quad_N
@test sum(s) ≈ 6
@test ∑(s) ≈ 6

cell_measure = get_cell_measure(trian)
cell_measure_N = get_cell_measure(trian_N,trian)
@test length(cell_measure) == num_cells(model)
@test length(cell_measure_N) == num_cells(model)
@test sum(cell_measure) ≈ 1
@test sum(cell_measure_N) ≈ 6

quad = CellQuadrature(trian,Quadrature(duffy,2))
s = ∫(1)*quad
@test sum(s) ≈ 1
@test ∑(s) ≈ 1

quad = CellQuadrature(trian,Quadrature(TET,duffy,2))
s = ∫(1)*quad
@test sum(s) ≈ 1
@test ∑(s) ≈ 1

#using Gridap.Visualization
#writevtk(trian,"trian",celldata=["a"=>cell_measure,"b"=>cell_measure_N])

#using Gridap.Fields: MockField, MockBasis, OtherMockBasis
#
#degree = 3
#quad = Quadrature(HEX,degree)
#test_quadrature(quad)
#@test sum(get_weights(quad)) ≈ 1
#
#degree = (1,2,3)
#quad = Quadrature(HEX,degree)
#test_quadrature(quad)
#@test sum(get_weights(quad)) ≈ 1
#
#degree = 3
#quad = Quadrature(TET,degree)
#test_quadrature(quad)
#@test sum(get_weights(quad)) ≈ 0.5*1/3
#
#l = 10
#quad = CellQuadrature(degree,[QUAD,],Fill(1,l))
#q = get_coordinates(quad)
#@test isa(q,Fill)
#w = get_weights(quad)
#@test isa(w,Fill)
#
#quad = CellQuadrature(degree,[QUAD,],ones(Int,l))
#q = get_coordinates(quad)
#@test isa(q,CompressedArray)
#w = get_weights(quad)
#@test isa(w,CompressedArray)
#
#ndofs = 4
#
#z = 2.0
#v = VectorValue(3.0,1.5)
#w = VectorValue(3.4,3.5)
#a = MockBasis{2}(v,ndofs)
#b = MockBasis{2}(w,ndofs)
#c = fill(1.0,ndofs)
#f = OtherMockBasis{2}(ndofs)
#g = MockField{2}(v)
#
#l = 10
#zl = [ z for  i in 1:l]
#cl = fill(c,l)
#fl = Fill(f,l)
#ϕl = lincomb(fl,cl)
#gl = fill(g,l)
#al = Fill(a,l)
#bl = fill(b,l)
#
#gf = GenericCellField(gl,ϕl,Val(true))
#af = GenericCellField(al,ϕl,Val(true),Fill((Base.OneTo(ndofs),),l),Val((:,)))
#bf = GenericCellField(bl,ϕl,Val(true),Fill((Base.OneTo(ndofs),),l),Val((:,)))
#zf = convert_to_cell_field(zl,ϕl)
#df = af*zf
#dft = trialize_cell_basis(df)
#
#@test sum(integrate(1,ϕl,quad)) ≈ 213.33333333333323
#@test sum(integrate(zl,ϕl,quad)) ≈  426.66666666666646
#
#cm = integrate(df⋅dft,ϕl,quad)
#@test cm[1] ≈ fill(560.0,4,4)
#
#cv = integrate(df⋅v,ϕl,quad)
#@test cv[1] ≈ fill(359.9999999999998,4)
#
#quad1 = CellQuadrature(degree,[QUAD,],Fill(1,4))
#quad2 = CellQuadrature(degree,[TRI,],Fill(1,6))
#quad = lazy_append(quad1,quad2)
#@test isa(get_coordinates(quad),AppendedArray)
#@test isa(get_weights(quad),AppendedArray)
#@test isa(get_array(quad),AbstractArray{<:Quadrature})
#@test isa(get_array(quad),AppendedArray)

end # module


text:
module Models

export GANHyperParams
export GAN

using Flux
using Flux.Optimise: AbstractOptimiser
using Parameters

#################
# model structs #
#################

# hyperparameters

@with_kw struct GANHyperParams
    latent_dim::Int = 100
    dscr_loops::Int = 1
    minibatch::Int  = 50
    η_gen::Float32  = 0.0001f0
    η_dscr::Float32 = 0.0001f0
    img_size::Tuple = undef
end

# simple GAN model (GAN)

abstract type AbstractGAN end

mutable struct GAN <: AbstractGAN
    G::Chain{T} where T
    D::Chain{H} where H
    G_opt::AbstractOptimiser
    D_opt::AbstractOptimiser
    hparams::GANHyperParams
end

function GAN(G::Chain{T}, D::Chain{H}; kws...) where {T, H}
    hparams = GANHyperParams(; kws...)
    G_opt = OADAM(hparams.η_gen)
    D_opt = OADAM(hparams.η_dscr)
    return GAN(G, D, G_opt, D_opt, hparams)
end

function GAN(G::Chain{T}, D::Chain{H}, opt::AbstractOptimiser; kws...) where {T,H}
    hparams = GANHyperParams(; kws...)
    G_opt = opt(hparams.η_gen)
    D_opt = opt(hparams.η_dscr)
    return GAN(G, D, G_opt, D_opt, hparams)
end

function GAN(G::Chain{T}, D::Chain{H},
             G_opt::AbstractOptimiser, D_opt::AbstractOptimiser; kws...) where {T,H}
    hparams = GANHyperParams(; kws...)
    G_opt = G_opt(hparams.η_gen)
    D_opt = D_opt(hparams.η_dscr)
    return GAN(G, D, G_opt, D_opt, hparams)
end

Flux.@functor GAN
# TODO: add conditional GAN model (CondGAN)
# TODO: add art GAN model (ArtGAN)

end

max_stars_repo_name: aarontrowbridge/GAN

max_issues_repo_name: aarontrowbridge/GAN

max_forks_repo_name: aarontrowbridge/GAN


text:
module RxNav

using EzXML
using HTTP

export rcui, drugs, interaction, interaction_within_list, prescribable, interact

include("util.jl")
include("API/RxClassAPI.jl")
include("API/RxNormAPI.jl")
include("API/RxTermsAPI.jl")
include("API/DrugInteractionAPI.jl")

"""
    rcui(name)

Take a name of an NDC drug, return its rxcui as String.
"""
function rcui(name)
    try
        doc = getdoc("rcui", HTTP.URIs.escapeuri(name))
        idstring = nodecontent(findfirst("//idGroup/rxnormId", doc))
        return idstring
    catch y
        @warn y
        return ""
    end
end

"""
    drugs(name)

Given a drug name, return a list of all available dosing forms of the drug.
"""
function drugs(name)
    try
        doc = getdoc("drugs", HTTP.URIs.escapeuri(name))
        nameelements = findall("//drugGroup/conceptGroup/conceptProperties/name", doc)
        return nodecontent.(nameelements)
    catch y
        @warn y
        return String[]
    end
end

"""
    interact(list::Vector)
    interact(s1::String, severeonly::Bool=true)
    interact(s1::String, s2::String, args...)

Get a list of interactions for a single drug (or rxcui drug id) or pairwise interactions for more than one drug (or rxcuid).
"""
interact(list::Vector) = if length(list) > 1 interaction_within_list(list); else interact(first(list)) end
interact(s1::String, severeonly::Bool=true) = interaction(s1; ONCHigh=severeonly)
interact(s1::String, s2::String, args...) = interact([[s1, s2]; [x for x in args]])

"""
    interaction(id; ONCHigh = true)

Given a drug name or rxcui id string, return known drug interations for that drug.
If ONCHigh is true only return the ONCHigh database entries, which returns fewer
entries, tending to list only the more significant interactions. Set ONCHigh
to false to get all known interactions, which can be multiple and sometimes redundant.
Returns a `Vector` of `NamedTuple`s as in (drug1, drug2, severity, description)
"""
function interaction(id; ONCHigh = true)
    if !is_in_rxcui_format(id)
        id = rcui(id)
    end
    interactions = NamedTuple[]
    try
        tail = HTTP.URIs.escapeuri(id) * (ONCHigh ? "&sources=ONCHigh" : "")
        doc = getdoc("interaction", tail)
        pairs = findall("//interactionTypeGroup/interactionType/interactionPair", doc)
        for p in pairs
            sev = nodecontent(findfirst("severity", p))
            desc = nodecontent(findfirst("description", p))
            enames = findall("interactionConcept/minConceptItem/name", p)
            if !isempty(enames)
                names = nodecontent.(enames)
                push!(interactions, (drug1=names[1], drug2=names[2], severity=sev, description=desc))
            end
        end
    catch y
        @warn y
    end
    return interactions
end

"""
    interaction_within_list(idlist::Vector{String})

Given a list of drug names or rxcui id strings, return known drug interations for
that combination of drugs. Results are organized pairwise, so if A, B, and C have
mutual interactions this will be reported for example as A with B, A with C, B with C.
Returns a `Vector` of `NamedTuple`s as in (drug1, drug2, severity, description)
"""
function interaction_within_list(idlist::Vector{String})
    for (i, id) in enumerate(idlist)
        if !is_in_rxcui_format(id)
            idlist[i] = rcui(id)
        end
    end
    interactions = NamedTuple[]
    try
        arg = join(map(x -> HTTP.URIs.escapeuri(x), idlist), "+")
        doc = getdoc("interactionlist", arg)
        pairs = findall("//fullInteractionTypeGroup/fullInteractionType/interactionPair", doc)
        for p in pairs
            sev = nodecontent(findfirst("severity", p))
            desc = nodecontent(findfirst("description", p))
            enames = findall("interactionConcept/minConceptItem/name", p)
            if !isempty(enames)
                names = nodecontent.(enames)
                push!(interactions, (drug1=names[1], drug2=names[2], severity=sev, description=desc))
            end
        end
        catch y
        @warn y
    end
    return interactions
end

end  # module


text:
using AddLatest
using Test

import Documenter
import Pkg
import UUIDs

@testset "AddLatest.jl" begin
    cache_1 = RegistryParsingCache()

    @testset "Unit tests" begin
        @testset "assert.jl" begin
            @test_throws AddLatest.PackageNotFoundException AddLatest.assert_package_found(false, "PkgA", nothing)
            @test_throws AddLatest.PackageNotFoundException AddLatest.assert_package_found(false, "PkgA", UUIDs.UUID(0))
            @test_throws AddLatest.PackageAmbiguousException AddLatest.assert_no_ambiguity(true, "PkgA", nothing)
            @test_throws AddLatest.PackageAmbiguousException AddLatest.assert_no_ambiguity(true, "PkgA", UUIDs.UUID(0))
            let
                a = AddLatest.VersionAndTree("1.0.0", "foo")
                b = AddLatest.VersionAndTree("1.0.0", "bar")
                @test_throws AddLatest.RegistriesDisagreeException AddLatest.assert_no_disagreement(a, b)
            end
        end

        @testset "parse-registry.jl" begin
            let
                versions = Dict()
                versions["1.0.0"] = Dict()
                versions["1.0.0"]["yanked"] = true
                AddLatest._delete_yanked_versions!(versions)
            end
            let
                versions = Dict()
                versions["1.0.0"] = Dict()
                versions["1.0.0"]["yanked"] = false
                AddLatest._delete_yanked_versions!(versions)
            end
        end

        @testset "registry-list.jl" begin
            @test AbstractRegistry(LocalFolderRegistry("foo")) == LocalFolderRegistry("foo")
        end
    end

    @testset "Integration tests" begin
        @testset "Integration with the Pkg API" begin
            @testset "without specifying a cache" begin
                Pkg.add(latest("Example"))
            end
            @testset "specify a cache" begin

                @testset "name only" begin
                    Pkg.add(latest("Example"; cache = cache_1))
                    Pkg.add(latest("Example"; cache = cache_1, registries = ["https://github.com/JuliaRegistries/General.git"]))
                    Pkg.add([latest("Example"; cache = cache_1)])
                    Pkg.add(AddLatest._convert_to_package_spec(latest("Example"; cache = cache_1)))
                    Pkg.add(Pkg.Types.PackageSpec(latest("Example"; cache = cache_1)))
                    Pkg.add(Pkg.API.Package(latest("Example"; cache = cache_1)))
                    Pkg.add(AddLatest.convert(Pkg.Types.PackageSpec, latest("Example"; cache = cache_1)))

                    Pkg.add(latest("Example"; cache = cache_1))
                    Pkg.pin(latest("Example"; cache = cache_1))
                    Pkg.rm("Example")
                    Pkg.add([latest("Example"; cache = cache_1)])
                    Pkg.pin([latest("Example"; cache = cache_1)])
                    Pkg.rm("Example")
                end
                @testset "name and uuid" begin
                    Pkg.add(latest("Example", UUIDs.UUID("7876af07-990d-54b4-ab0e-23690620f79a"); cache = cache_1))
                end
            end
        end
        @testset "test a package with a yanked version" begin
            latest("MLJBase"; cache = cache_1)
        end
    end

    Base.empty!(cache_1)

    @testset "Doctests" begin
        Documenter.doctest(AddLatest)
    end
end


text:
using Base.Threads
using LinearAlgebra

#local ttvec rank increase function with noise ϵ_wn
function tt_up_rks_noise(tt_vec,tt_ot_i,rkm,rk,ϵ_wn)
	vec_out = zeros(eltype(tt_vec),size(tt_vec,1),rkm,rk)
	vec_out[:,1:size(tt_vec,2),1:size(tt_vec,3)] = tt_vec
	if !iszero(ϵ_wn)
		if rkm == size(tt_vec,2) && rk>size(tt_vec,3)
			Q = rand_orthogonal(size(tt_vec,1)*rkm,rk-size(tt_vec,3))
			vec_out[:,:,size(tt_vec,3)+1:rk] = ϵ_wn*reshape(Q,size(tt_vec,1),rkm,rk-size(tt_vec,3))
			tt_ot_i =0
		elseif rk == size(tt_vec,3) && rkm>size(tt_vec,2)
			Q = rand_orthogonal(rkm-size(tt_vec,2),size(tt_vec,1)*rk)
			vec_out[:,size(tt_vec,2)+1:rkm,:] = ϵ_wn*reshape(Q,size(tt_vec,1),rkm-size(tt_vec,2),rk)
			tt_ot_i =0
		elseif rk>size(tt_vec,3) && rkm>size(tt_vec,2)
			Q = rand_orthogonal((rkm-size(tt_vec,2))*size(tt_vec,1),(rk-size(tt_vec,3)))
			vec_out[:,size(tt_vec,2)+1:rkm,size(tt_vec,3)+1:rk] = ϵ_wn*reshape(Q,size(tt_vec,1),rkm-size(tt_vec,2),rk-size(tt_vec,3))
		end
	end
	return vec_out
end

"""
returns the ttvector with ranks rks and noise ϵ_wn for the updated ranks
"""
function tt_up_rks(x_tt::ttvector{T},rk_max::Int;rks=vcat(1,rk_max*ones(Int,length(x_tt.ttv_dims)-1),1),ϵ_wn=0.0) where T<:Number
	d = length(x_tt.ttv_dims)
	vec_out = Array{Array{T}}(undef,d)
	out_ot = zeros(Int64,d)
	@assert(rk_max > maximum(x_tt.ttv_rks),"New bond dimension too low")
	n_in = 1
	n_out = prod(x_tt.ttv_dims)
	for i in 1:d
		n_in *= x_tt.ttv_dims[i]
		n_out = Int(n_out/x_tt.ttv_dims[i])
		rks[i+1] = min(rks[i+1],n_in,n_out)
		vec_out[i] = tt_up_rks_noise(x_tt.ttv_vec[i],x_tt.ttv_ot[i],rks[i],rks[i+1],ϵ_wn)
	end
	return ttvector{T}(vec_out,x_tt.ttv_dims,rks,x_tt.ttv_ot)
end

"""
returns the orthogonalized ttvector with root i
"""
function orthogonalize(x_tt::ttvector{T};i=1::Int) where T<:Number
	d = length(x_tt.ttv_dims)
	@assert(1≤i≤d, DimensionMismatch("Impossible orthogonalization"))
	y_tt = copy(x_tt)
	y_tt.ttv_ot[i]=0
	for j in 1:i-1
		y_tt.ttv_ot[j]=1
		y_vectemp = reshape(y_tt.ttv_vec[j],y_tt.ttv_dims[j]*y_tt.ttv_rks[j],y_tt.ttv_rks[j+1])
		F = qr(y_vectemp)
		y_tt.ttv_vec[j] = reshape(Matrix(F.Q)[:,1:y_tt.ttv_rks[j+1]],x_tt.ttv_dims[j],y_tt.ttv_rks[j],y_tt.ttv_rks[j+1])
		@threads for k in 1:x_tt.ttv_dims[j]
			y_tt.ttv_vec[j+1][k,:,:] = F.R[1:y_tt.ttv_rks[j+1],1:y_tt.ttv_rks[j+1]]*y_tt.ttv_vec[j+1][k,:,:]
		end
	end
	for j in d:-1:i+1
		y_tt.ttv_ot[j]=-1
		y_vectemp = reshape(permutedims(y_tt.ttv_vec[j],[2,1,3]),y_tt.ttv_rks[j],y_tt.ttv_dims[j]*y_tt.ttv_rks[j+1])
		F = lq(y_vectemp)
		y_tt.ttv_vec[j] = permutedims(reshape(Matrix(F.Q)[1:y_tt.ttv_rks[j],:],y_tt.ttv_rks[j],x_tt.ttv_dims[j],y_tt.ttv_rks[j+1]),[2 1 3])
		@threads for k in 1:x_tt.ttv_dims[j]
			y_tt.ttv_vec[j-1][k,:,:] = y_tt.ttv_vec[j-1][k,:,:]*F.L[1:y_tt.ttv_rks[j],1:y_tt.ttv_rks[j]]
		end
	end
	return y_tt
end

"""
returns a TT representation where the singular values lower than tol are discarded
"""
function tt_rounding(x_tt::ttvector{T};tol=1e-12) where T<:Number
	d = length(x_tt.ttv_dims)
	y_rks = copy(x_tt.ttv_rks)
	y_vec = copy(x_tt.ttv_vec)
	for j in 1:d-1
		A = zeros(T,x_tt.ttv_dims[j],y_rks[j],x_tt.ttv_dims[j+1],y_rks[j+2])
		@tensor A[a,b,c,d] = y_vec[j][a,b,z]*y_vec[j+1][c,z,d]
		u,s,v = svd(reshape(A,size(A,1)*size(A,2),:);alg=LinearAlgebra.QRIteration())
		rtol = norm(s)*tol
		Σ = s[s.>rtol]
		y_rks[j+1] = length(Σ)
		y_vec[j] = reshape(u[:,s.>rtol],x_tt.ttv_dims[j],y_rks[j],:)
		y_vec[j+1] = permutedims(reshape(Diagonal(Σ)*v'[s.>rtol,:],:,x_tt.ttv_dims[j+1],y_rks[j+2]),[2 1 3])
	end
	for j in d:-1:2
		A = zeros(T,x_tt.ttv_dims[j-1],y_rks[j-1],x_tt.ttv_dims[j],y_rks[j+1])
		@tensor A[a,b,c,d] = y_vec[j-1][a,b,z]*y_vec[j][c,z,d]
		F = svd(reshape(A,size(A,1)*size(A,2),:),alg=LinearAlgebra.QRIteration())
		rtol = norm(F.S)*tol
		Σ = F.S[F.S .>rtol]
		y_rks[j] = length(Σ)
		y_vec[j] = permutedims(reshape(F.Vt[F.S.>rtol,:],:,x_tt.ttv_dims[j],y_rks[j+1]),[2 1 3])
		y_vec[j-1] = reshape(F.U[:,F.S.>rtol]*Diagonal(Σ),x_tt.ttv_dims[j-1],y_rks[j-1],:)
	end
	return ttvector{T}(y_vec,x_tt.ttv_dims,y_rks,vcat(0,ones(Int64,d-1)))
end

"""
returns the rounding of the TT operator
"""
function tt_rounding(A_tto::ttoperator;tol=1e-12)
	return ttv_to_tto(tt_rounding(tto_to_ttv(A_tto);tol=tol))
end

"""
returns the singular values of the reshaped tensor x[μ_1⋯μ_k;μ_{k+1}⋯μ_d] for all 1≤ k ≤ d
"""
function tt_svdvals(x_tt::ttvector;tol=1e-14)
	d = length(x_tt.ttv_dims)
	Σ = Array{Array{Float64,1},1}(undef,d-1)
	y_tt = orthogonalize(x_tt)
	y_rks = y_tt.ttv_rks
	for j in 1:d-1
		A = zeros(y_tt.ttv_dims[j],y_rks[j],y_tt.ttv_dims[j+1],y_rks[j+2])
		@tensor A[a,b,c,d] = y_tt.ttv_vec[j][a,b,z]*y_tt.ttv_vec[j+1][c,z,d]
		u,s,v = svd(reshape(A,size(A,1)*size(A,2),:),alg=LinearAlgebra.QRIteration())
		Σ[j] = s[s.>tol]
		y_rks[j+1] = length(Σ[j])
		y_tt.ttv_vec[j+1] = permutedims(reshape(Diagonal(Σ[j])*v'[s.>tol,:],:,y_tt.ttv_dims[j+1],y_rks[j+2]),[2 1 3])
	end
	return Σ
end

function sv_trunc(s::Array{Float64},tol)
	if tol==0.0
		return s
	else
		d = length(s)
		i=0
		weight = 0.0
		norm2 = dot(s,s)
		while (i<d) && weight<tol*norm2
			weight+=s[d-i]^2
			i+=1
		end
		return s[1:(d-i+1)]
	end
end

function left_compression(A,B;tol=1e-12)
    dim_A = [i for i in size(A)]
    dim_B = [i for i in size(B)]

    B = permutedims(B, [2,1,3]) #B r_1 x n_2 x r_2
    U = reshape(A,:,dim_A[3])
    u,s,v = svd(U*reshape(B, dim_B[2],:),full=false) #u is the new A, dim(u) = n_1r_0 x tilde(r)_1
    s_trunc = sv_trunc(s,tol)
    dim_B[2] = length(s_trunc)
    U = reshape(u[:,1:dim_B[2]],dim_A[1],dim_A[2],dim_B[2])
    B = reshape(Diagonal(s_trunc)*v[:,1:dim_B[2]]',dim_B[2],dim_B[1],dim_B[3])
    return U, permutedims(B,[2,1,3])
end

"""
parallel compression of the ttvector
TODO refactoring
"""
function tt_compression_par(X::ttvector;tol=1e-14,Imax=2)
    Y = deepcopy(X.ttv_vec) :: Array{Array{Float64,3},1}
    rks = deepcopy(X.ttv_rks) :: Array{Int64}
    d = length(X.ttv_dims)
    rks_prev = zeros(Integer,d)
    i=0
    while norm(rks-rks_prev)>0.1 && i<Imax
        i+=1
        rks_prev = deepcopy(rks) :: Array{Int64}
        if mod(i,2) == 1
            @threads for k in 1:floor(Integer,d/2)
                Y[2k-1], Y[2k] = left_compression(Y[2k-1], Y[2k], tol=tol)
                rks[2k-1] = size(Y[2k-1],3)
            end
        else
            @threads for k in 1:floor(Integer,(d-1)/2)
                Y[2k], Y[2k+1] = left_compression(Y[2k], Y[2k+1], tol=tol)
                rks[2k] = size(Y[2k],3)
            end
        end
    end
    return ttvector(Y,X.ttv_dims,rks,zeros(Integer,d))
end

function tt_compression_par(A::ttoperator;tol=1e-14,Imax=2)
	return ttv_to_tto(tt_compression_par(tto_to_ttv(A);tol=tol,Imax=Imax))
end

max_stars_repo_name: msdupuy/Tensor-Train-Julia

max_issues_repo_name: msdupuy/Tensor-Train-Julia

max_forks_repo_name: msdupuy/Tensor-Train-Julia


text:
using IRTools: blocks, returnvalue, branches, isreturn, predecessors

isvar(v) = typeof(v) == Variable

# zt -- this is far too long for what it does
function returnvars(b)
  if isreturn(b)
    ret = returnvalue(b)
    if ret isa Variable
      return [ret]
    else
      return Variable[]
    end
  else
    return Variable[]
  end
end

"Vars that are used in branching in block `b`"
branchvars(b) =
  reduce(vcat, [filter(isvar, [br.args; br.condition]) for br in  branches(b)], init = Variable[])

stmtvars(s::Statement) =
  filter(isvar, s.expr.args)

"vars used in statements of `b`"
stmts(b) =
  reduce(vcat, [filter(isvar, b[v].expr.args) for v in keys(b)], init = Variable[])

"vars used in some form by block `b`"
usedvars(b) = stmts(b) ∪ branchvars(b) ∪ returnvars(b)

"All"
branchpreds(block) =
  filter(br -> br.block == block, branches(predecessors(block)))

"Variables that are produced by `b`: inputs or lhs of statements"
producedvars(b) = arguments(b) ∪ keys(b)

"Variables that are used in block `b` but neither defined in `b` nor inputs to `b`"
locallyundefinedvars(b) = setdiff(usedvars(b), producedvars(b))

"Repalce all occurances of `x` with `y` variable `v1` with variable `v2` in block `b`"
function varreplace!(b, v1, v2)
  IRTools.prewalk!(b) do x
    if x isa Variable && x == v1
      v2
    else
      x
    end
  end
end

"""
`passvars!(ir)`

Removes implicit use of variable.
Updates `ir` such that if a block uses some variable `v` then `v` is an input to that block.
"""
function passvars!(ir)
  blockidss_ = Set{Int}([b.id for b in IRTools.blocks(ir)]) ## zt: do we want to be putting ir into
  m = Dict{Variable, Set{Variable}}()  # Mapping from vars to equivalence class
  while !isempty(blockidss_)
    bid = pop!(blockidss_)
    b = IRTools.block(ir, bid)
    undefvars = locallyundefinedvars(b)

    # Nothing to do if no undefined vars, skip
    if !isempty(undefvars)
      newargs = Variable[]
      for v in undefvars
        # Is `v` truly undefined (locally), or is it equal to one of the inputs
        res = findfirst(arg -> arg in keys(m) && v in m[arg], arguments(b))
        # If so : add a new argument
        if isnothing(res)
          arg = argument!(b)      # add the enecessary inputs
          push!(newargs, v)
          # zt -- doing multiple passes, could do more efficiently with one
          varreplace!(b, v, arg)  #
          push!(get!(m, arg, Set{Variable}()), v)
        else
          varreplace!(b, v, arguments(b)[res])
        end
      end

      isempty(newargs) && continue
      # If we have added new arguments we need to update all the branch points that lead here
      for pa in predecessors(b)
        # For each block which branches to `b`, need to update the branch poins
        for br in branches(pa)
          if br.block == b.id
            offset = length(br.args) - length(newargs)
            for i = 1:length(newargs)
              br.args[i + offset] = newargs[i]
            end
            @assert !any(isnothing, br.args)
          end
        end
        push!(blockidss_, pa.id)    # Add to queue (since we've updated it and it may now have undefvars vars)
      end
    end
  end
  ir
end


text:
content = open(metadatapath, "r+") do file
   read(file, String)
end

content = replace(content, r"[àâ]" => 'a')
content = replace(content, r"[èéê]" => 'e')
content = replace(content, 'ü' => 'u')
content = replace(content, '’' => '\'')
content = replace(content, r"[“”]" => '\"')

open(metadatapath, "w") do file
   write(file, content)
end


text:
#### 03b - Random Forests predictions ####
using RCall
R"source(file.path('src', 'required.R'))" # bug with `velox` if not called here
include("required.jl")

## Conditional arguments
# save_data = true

## Perform RandomForests
begin
    R"""
    ## 0. Load packages ####
    source(file.path("src", "required.R"))

    # Conditional evaluations
    # subset_qc <- TRUE # subset to QC data (optional)
    # create_models <- TRUE # train models
    # save_models <- TRUE # save & overwrite models

    ## 1. Load data ####

    # Prepare data
    source(here("src", "02a_training_data-preparation.R"))

    # Remove sites with NA values
    inds_na <- map(env_full, ~ which(is.na(.x)))
    (inds_na <- sort(unique(unlist(inds_na))))
    vars_nona <- vars_full[-inds_na,]

    # Load model
    load(here("data", "rdata", "rf_models.RData"))

    # Make prediction
    system.time(
        rf_pred <- pblapply(ranger_models, function(x) predict(x, vars_nona))
    )

    # Extract predictions
    predictions <- map(rf_pred, "predictions") %>%
        map_df(~ as.numeric(levels(.x))[.x])

    # Add sites with NAs
    predictions_full <- matrix(
        NA,
        nrow = nrow(vars_full),
        ncol = ncol(predictions)
    )
    colnames(predictions_full) <- colnames(predictions)
    predictions_full[-inds_na,] <- as.matrix(predictions)
    """
end
@rget predictions predictions_full inds_na

## Create Y matrices

# Get matrix Y
Y = replace(predictions_full, missing => nothing) |> Array{Union{Nothing,Float32}}
# Set values to nothing if no species present
inds_zeros = _indsnotobs(Y)
Y[inds_zeros, :] .= nothing

## Create distributions

# Load raw distributions (for grid size)
glossary = CSV.read(joinpath("data", "proc", "glossary.csv"), DataFrame)
spenames = filter(:type => ==("species"), glossary).full_name
distributions = [
    geotiff(SimpleSDMPredictor, joinpath("data", "proc", "distributions_raw.tif"), i) for
    i in eachindex(spenames)
]
raw_distributions = distributions
# Get layer dimensions & limits
dims = size(raw_distributions[1].grid)
lims = (
    left=raw_distributions[1].left,
    right=raw_distributions[1].right,
    bottom=raw_distributions[1].bottom,
    top=raw_distributions[1].top,
)

# Create RF distribution layers
Ydistrib = replace(Y, 0.0 => nothing)
Ygrids = [Ydistrib[:, col] for col in 1:size(Ydistrib, 2)]
Ygrids = reshape.(Ygrids, dims...) .|> Array
distributions = SimpleSDMResponse.(Ygrids, lims...)

## Export results
# save_data = true
if (@isdefined save_data) && save_data == true
    @save joinpath("data", "jld2", "rf-distributions.jld2") distributions
    _zip_jld2(
        joinpath("data", "jld2", "rf-distributions.zip"),
        joinpath("data", "jld2", "rf-distributions.jld2"),
    )
    touch(joinpath("data", "jld2", "rf-distributions.jld2"))
end

## Get richness & LCBD

richness = calculate_richness(Y, distributions[1])
lcbd = calculate_lcbd(Y, distributions[1])

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
function D3Trees.D3Tree(D::AdaOPSTree; title="AdaOPS Tree", kwargs...)
    lenb = D.b
    lenba = D.ba
    len = lenb + lenba
    children = Vector{Vector{Int}}(undef, len)
    text = Vector{String}(undef, len)
    tt = fill("", len)
    link_style = fill("", len)
    for b in 1:lenb
        children[b] = D.children[b] .+ lenb
        text[b] = @sprintf("""
                           o:%s prob:%6.2f
                           l:%6.2f, u:%6.2f""",
                           b==1 ? "<root>" : string(D.obs[b]),
                           D.obs_prob[b],
                           D.l[b],
                           D.u[b],
                          )
        tt[b] = """
                o: $(b==1 ? "<root>" : string(D.obs[b]))
                prob: $(D.obs_prob[b])
                l: $(D.l[b])
                u: $(D.u[b])
                $(length(D.children[b])) children
                """
        link_width = 2.0
        link_style[b] = "stroke-width:$link_width"
        for ba in D.children[b]
            link_style[ba+lenb] = "stroke-width:$link_width"
        end

        for ba in D.children[b]
            children[ba+lenb] = D.ba_children[ba]
            text[ba+lenb] = @sprintf("""
                                     a:%s r:%6.2f |ϕ|:%2d
                                     l:%6.2f, u:%6.2f""",
                                     D.ba_action[ba], D.ba_r[ba], length(D.ba_particles[ba]),
                                     D.ba_l[ba], D.ba_u[ba],
                                     )
            tt[ba+lenb] = """
                          a: $(D.ba_action[ba])
                          r: $(D.ba_r[ba])
                          |ϕ|:$(length(D.ba_particles[ba]))
                          l: $(D.ba_l[ba])
                          u: $(D.ba_u[ba])
                          $(length(D.ba_children[ba])) children
                          """
        end

    end
    return D3Tree(children;
                  text=text,
                  tooltip=tt,
                  link_style=link_style,
                  title=title,
                  kwargs...
                 )
end

Base.show(io::IO, mime::MIME"text/html", D::AdaOPSTree) = show(io, mime, D3Tree(D))
Base.show(io::IO, mime::MIME"text/plain", D::AdaOPSTree) = show(io, mime, D3Tree(D))


text:
# Curve methods implementation

_zero_rate{METHOD<:RateInterpolation}(method::METHOD, curve::AbstractIRCurve, maturity::Date) = _zero_rate(method, curve_get_dtm(curve), curve_get_zero_rates(curve), days_to_maturity(curve, maturity))

function _zero_rate(comp::CompositeInterpolation, curve::AbstractIRCurve, maturity::Date)
	dtm = days_to_maturity(curve, maturity)

	if dtm < curve_get_dtm(curve)[1]
		return _zero_rate(comp.before_first, curve, maturity)
	elseif dtm > curve_get_dtm(curve)[end]
		return _zero_rate(comp.after_last, curve, maturity)
	else
		return _zero_rate(comp.inner, curve, maturity)
	end
end

# Returns tuple (index_a, index_b) for input vector x
# for interpolands on linear interpolation on point x_out
function _interpolationpoints{T}(x::Vector{T}, x_out::T)
	local index_a::Int
	local index_b::Int

	if x_out <= x[1]
		# Interpolation point is before first vertice
		# Slope will be determined by the 1st and 2nd vertices
		index_a = 1
		index_b = 2
	elseif x_out >= x[end]
		# Interpolation point is after last vertice
		# Slope will be determined by the last and last-1 vertices
		index_b = length(x)
		index_a = index_b - 1
	else
		# Inner point
		index_a =  findlast(a -> a < x_out, x) # last element before x_out on x
		index_b = findfirst(b -> b >=  x_out, x) # first element after x_out on x
	end

	return index_a, index_b
end

# Perform Linear interpolation. Slope is determined by points (Xa, Ya) and (Xb, Yb).
# Interpolation occurs on point (x_out, returnvalue)
_linearinterp{TX, TY}(Xa::TX, Ya::TY, Xb::TX, Yb::TY, x_out::TX) = (x_out - Xa) * (Yb - Ya) / (Xb - Xa) + Ya

# Linear interpolation of zero_rates
function _zero_rate(::Linear, x::Vector{Int}, y::Vector{Float64}, x_out::Int)
	# If this curve has only 1 vertice, this will be a flat curve
	if length(x) == 1
		return y[1]
	end

	index_a, index_b = _interpolationpoints(x, x_out)
	return _linearinterp(x[index_a], y[index_a], x[index_b], y[index_b], x_out)
end

# Step Function
function _zero_rate(::StepFunction, x::Vector{Int}, y::Vector{Float64}, x_out::Int)
	# If this curve has only 1 vertice, this will be a flat curve
	if length(x) == 1
		return y[1]
	end

	if x_out <= x[1]
		# Interpolation point is before first vertice
		# The result will be extrapolated using the first vertice zero_rate
		return y[1]
	elseif x_out >= x[end]
		# Interpolation point is after last vertice
		# The result will be extrapolated using the last vertice zero_rate
		return y[end]
	else
		# Inner point
		return y[findlast(a -> a <= x_out, x)] # last element before x_out on x
	end
end

# Flat Forward is linear interpolation on the log of discountfactors
# Maybe not useful for SimpleCompounding curves.
function _zero_rate(::FlatForward, curve::AbstractIRCurve, maturity::Date)
	# If this curve has only 1 vertice, this will be a flat curve
	if length(curve_get_zero_rates(curve)) == 1
		return curve_get_zero_rates(curve)[1]
	end

	x_out = days_to_maturity(curve, maturity)
	curve_dtm = curve_get_dtm(curve)
	curve_zero_rates = curve_get_zero_rates(curve)
	index_a, index_b = _interpolationpoints(curve_dtm, x_out)
	Xa = curve_dtm[index_a]
	Ya = curve_zero_rates[index_a]
	Xb = curve_dtm[index_b]
	Yb = curve_zero_rates[index_b]

	_daysperyear_ = daysperyear(curve_get_daycount(curve))
	year_fraction_a = Xa / _daysperyear_
	logPa = log(_discountfactor(curve_get_compounding(curve), Ya, year_fraction_a))

	year_fraction_b = Xb / _daysperyear_
	logPb = log(_discountfactor(curve_get_compounding(curve), Yb, year_fraction_b))

	year_fraction_x = x_out / _daysperyear_
	logPx = _linearinterp(year_fraction_a, logPa, year_fraction_b, logPb, year_fraction_x)

	return discountfactor_to_rate(curve_get_compounding(curve), exp(logPx), year_fraction_x)
end

function _zero_rate(::NelsonSiegel, curve::AbstractIRCurve, maturity::Date)

	# beta1 = param[1]
	# beta2 = param[2]
	# beta3 = param[3]
	# lambda = param[4]

	param = curve_get_model_parameters(curve)
	t = yearfraction(curve, maturity)
	_exp_lambda_t_ = exp(-param[4]*t)
	F_beta2 = (1.0 - _exp_lambda_t_) / (param[4]*t)

	return param[1] + param[2]*F_beta2 + param[3]*(F_beta2 - _exp_lambda_t_)
end

function _zero_rate(::Svensson, curve::AbstractIRCurve, maturity::Date)

	# beta1 = param[1]
	# beta2 = param[2]
	# beta3 = param[3]
	# beta4 = param[4]
	# lambda1 = param[5]
	# lambda2 = param[6]

	param = curve_get_model_parameters(curve)
	t = yearfraction(curve, maturity)
	_exp_lambda1_t_ = exp(-param[5]*t)
	_exp_lambda2_t_ = exp(-param[6]*t)
	F_beta2 = (1.0 - _exp_lambda1_t_) / (param[5]*t)

	return param[1] + param[2]*F_beta2 + param[3]*(F_beta2 - _exp_lambda1_t_) +
			param[4]*( (1.0 - _exp_lambda2_t_)/(param[6]*t) - _exp_lambda2_t_)
end

include("splines.jl")

function _zero_rate(::CubicSplineOnRates, x::Vector{Int}, y::Vector{Float64}, x_out::Int)
	sp = splinefit(x, y)
	return splineint(sp, x_out)
end

function _zero_rate(::CubicSplineOnRates, curve::IRCurve, maturity_vec::Vector{Date})
	sp = splinefit(curve_get_dtm(curve), curve_get_zero_rates(curve))

	l = length(maturity_vec)
	rates = Array(Float64, l)

	for i in 1:l
		dtm = days_to_maturity(curve, maturity_vec[i])
		rates[i] = splineint(sp, dtm)
	end
	return rates
end

# Aux function for _zero_rate(::CubicSplinesOnDiscountFactors, ...) methods
function _splinefit_discountfactors(curve::AbstractIRCurve)
	dtm_vec = curve_get_dtm(curve)
	curve_rates_vec = curve_get_zero_rates(curve)
	l = length(dtm_vec)
	yf_vec = Array(Float64, l)
	discount_vec = Array(Float64, l)

	for i = 1:l
		yf_vec[i] = dtm_vec[i] / daysperyear(curve_get_daycount(curve))
		discount_vec[i] = _discountfactor(curve_get_compounding(curve), curve_rates_vec[i], yf_vec[i])
	end

	return splinefit(yf_vec, discount_vec)
end

function _zero_rate(::CubicSplineOnDiscountFactors, curve::AbstractIRCurve, maturity::Date)
	sp = _splinefit_discountfactors(curve)
	yf_maturity = yearfraction(curve_get_daycount(curve), curve_get_date(curve), maturity)
	result_discount_factor = splineint(sp, yf_maturity)
	return discountfactor_to_rate(curve_get_compounding(curve), result_discount_factor, yf_maturity)
end

function _zero_rate(::CubicSplineOnDiscountFactors, curve::AbstractIRCurve, maturity_vec::Vector{Date})
	sp = _splinefit_discountfactors(curve)
	mat_vec_len = length(maturity_vec)

	yf_maturity_vec = Array(Float64, mat_vec_len)
	for i in 1:mat_vec_len
		yf_maturity_vec[i] = yearfraction(curve_get_daycount(curve), curve_get_date(curve), maturity_vec[i])
	end

	return discountfactor_to_rate(curve.compounding, splineint(sp, yf_maturity_vec), yf_maturity_vec)
end

# Generate vector functions
for elty in (:FlatForward, :CompositeInterpolation, :StepFunction, :Linear, :NelsonSiegel, :Svensson)
	@eval begin
		function _zero_rate(m::$elty, curve::AbstractIRCurve, maturity_vec::Vector{Date})
			l = length(maturity_vec)
			rates = Array(Float64, l)
			for i = 1:l
				rates[i] = _zero_rate(m, curve, maturity_vec[i])
			end
		return rates
		end
	end
end


text:
using Gadfly, RDatasets, Distributions

set_default_plot_size(6inch, 3inch)

xs = dataset("datasets", "faithful").Waiting

quantiles = [0.1, 0.9]
p1 = plot(x=xs, Geom.density, Guide.title("quantiles=$quantiles"),
        layer(x=xs, Stat.quantile_bars(quantiles=quantiles), Geom.segment))

bandwidth = 0.6
p2 = plot(x=xs, Stat.density(bandwidth=bandwidth), Geom.line, Guide.title("bandwidth=$bandwidth"),
    layer(x=xs, Stat.quantile_bars(bandwidth=bandwidth), Geom.segment))

p = hstack(p1, p2)


text:
###################################
##General ABC SMC algorithm
###################################

"
Details of an application, mainly the data and model of interest.
The following related methods should be defined.

- `getDistances(a::ApplicationDetails, particles::Particles)`
- `report(a::ApplicationDetails, p::Particles, targetε::Float64)`
- `updateParticles!(a::ApplicationDetails, p::Particles, θ::Array{Float64, 1})`
"
abstract ApplicationDetails

"
Details of the current weighted particles during an SMC run.
This should include fields:

- `nparticles` Number of particles
- `llest` Current log likelihood estimate
- `ε` Current value of epsilon
- `iteration` How many SMC iterations have been performed
- `reachedBound` True if ABChd has terminated early because a lower bound specified on the log likelihood has been reached.
- `schedule` Schedule of epsilon values used (or, for an adaptive algorithm, schedule so far)
- `maxz` Maximum slice sampling value of |z| in previous SMC iteration

The following related methods should be defined:

- `Particles(a::AlgorithmDetails, nparticles::Int, parameters::Array{Float64, 1})`, a constructor. This should initialise each particle from the model specified by `a` under specified `parameters` value.
- `shuffle!(particles::Particles, indices::Array{Int, 1})` Updates particles so ith new particle is a copy of indices[i]th old particle.
"
abstract Particles

"
Performs or continues an SMC run.

## Arguments

- `a` Details of the application.
- `parameters` Vector of parameters whose likelihood is to be estimated.
- `schedule` Schedule of ε values. (Only required in non-adaptive version.)
- `targetε` Target ε value. (Only required in adaptive version.)
- `nparticles` Number of particles to use
- `k` Number of particles to accept at each iteration. (Only required in adaptive version.)
- `llBound` Lower bound on the log-likelihood. The algorithm terminates if the estimate is guaranteed to be lower than this.
- `d` (Effective) data dimension. Used in the normalising constant of the likelihood estimate, so optional.
- `quiet` When `false` progress reports are shown.

## Returns a `Particles` object with all the final information
"
##Non-adaptive version:
function doSMC(a::ApplicationDetails, parameters::Array{Float64, 1}, schedule::Array{Float64, 1}; nparticles=200, llBound=-Inf, d=0, quiet=false)
    ##Initialise SMC, including sampling from prior
    particles::Particles = Particles(a, nparticles, parameters)
    particles.schedule = schedule
    particles.llest = lgamma(1+d/2)-(d/2)*log(pi)-d*log(schedule[end]) ##log(1/V) where V is Lebesgue measure of a radius schedule[end] and dimension d
    particles.iteration = 0 ##Number of *completed* iterations

    ##Main loop
    nits = length(schedule)
    for ε in schedule
        particles.ε = ε
        dists::Array{Float64, 1} = getDistances(a, particles)
        acc = (dists .<= ε)
        particles.llest += log(mean(acc))
        particles.iteration += 1
        quiet || report(a, particles, schedule[end], d, false)
        if particles.llest <= llBound
            particles.reachedBound = true
            break
        end
        if particles.iteration < nits
            ##These steps can be skipped on the last iteration
            resample!(particles, acc)
            updateParticles!(a, particles, parameters)
        end
    end
    particles
end

##Adaptive version:
function doSMC(a::ApplicationDetails, parameters::Array{Float64, 1}, targetε::Float64; nparticles=200, k=ceil(Int, nparticles/2), llBound=-Inf, d=0, quiet=false)
    ##Initialise SMC by doing first iteration
    particles::Particles = Particles(a, nparticles, parameters)
    particles.ε = Inf
    particles.llest = lgamma(1+d/2)-(d/2)*log(pi)-d*log(targetε) ##log(1/V) where V is Lebesgue measure of a radius targetε and dimension d
    particles.iteration = 1 ##Number of *completed* iterations
    schedule = Float64[Inf]
    quiet || report(a, particles, targetε, d, true)

    ##Main loop
    while particles.ε > targetε
        dists::Array{Float64, 1} = getDistances(a, particles)
        particles.ε = max(targetε, select(dists, k))
        push!(schedule, particles.ε)
        acc = (dists .<= particles.ε)
        particles.llest += log(mean(acc))
        particles.iteration += 1
        quiet || report(a, particles, targetε, d, true)
        if (particles.llest <= llBound)
            particles.reachedBound = true
            break
        end
        if particles.ε == targetε
            break
        end
        resample!(particles, acc)
        updateParticles!(a, particles, parameters)
    end
    particles.schedule = schedule
    particles
end

##DEFAULT SMC FUNCTIONS

"
Perform a resampling step.
"
function resample!(particles::Particles, acc::BitArray)
    if sum(acc) == 0
        error("Resampling failed: no particles accepted")
    end
    N = particles.nparticles
    ind_alive = (1:N)[acc] ##Indices of alive particles
    ##Multinomial resampling
    indices = rand(ind_alive, N)

    ##Stratified resampling
    ##nAlive = sum(acc) ##Number of alive particles
    ##u = (1:N) - rand(N) ##Stratified uniform draws
    ##i = ceil(Int, u * nAlive/N) ##Random indices in 1:nAlive
    ##indices = ind_alive[i]

    ##Update the particles
    shuffle!(particles, indices)
    nothing
end

"
Report details of current iteration
"
function report(a::ApplicationDetails, p::Particles, targetε::Float64, d::Int, adaptive::Bool)
    print("Iteration $(p.iteration), eps=$(p.ε)")
    if !adaptive
        llest = p.llest
        if d>0 ##Avoids evaluating 0*Inf
            llest += d*(log(targetε) - log(p.ε))
        end
        print(", log of likelihood estimate=$llest, upper bound on final log likelihood estimate=$(p.llest)")
    end
    print("\n")
    nothing
end


text:
using MathOptFormat
using Test

const MOI = MathOptFormat.MOI
const MOIU = MOI.Utilities

@testset "MathOptFormat tests" begin
    @testset "$(file)" for file in ["CBF", "LP", "MOF", "MPS", "SDPA"]
        include("$(file)/$(file).jl")
    end

    @testset "Copying options" begin
        models = [
            MathOptFormat.CBF.Model,
            MathOptFormat.LP.Model,
            MathOptFormat.MOF.Model,
            MathOptFormat.MPS.Model,
            MathOptFormat.SDPA.Model,
        ]
        for src in models
            model_src = src()
            for dest in models
                model_dest = dest()
                MOI.copy_to(model_dest, model_src)
                @test !isempty(sprint(write, model_dest))
            end
            model_dest = MOIU.MockOptimizer(MOIU.Model{Float64}())
            MOI.copy_to(model_dest, model_src)
        end
    end

    @testset "Calling MOF.compressed_open" begin
        for cs in [MathOptFormat.Bzip2(), MathOptFormat.Gzip()]
            for open_type in ["a", "r+", "w+", "a+"]
                @test_throws ArgumentError MathOptFormat.compressed_open(
                    (x) -> nothing, "dummy.gz", open_type, cs
                )
            end
        end
    end

    @testset "Provided compression schemes" begin
        filename = joinpath(@__DIR__, "MPS", "free_integer.mps")
        model = MathOptFormat.Model(filename = filename)
        MOI.read_from_file(model, filename)
        filename = joinpath(@__DIR__, "free_integer.mps")
        MOI.write_to_file(model, filename * ".garbage")
        for ext in ["", ".bz2", ".gz"]
            MOI.write_to_file(model, filename * ext)
            model2 = MathOptFormat.Model(filename = filename * ext)
            MOI.read_from_file(model2, filename)
        end

        sleep(1.0)  # Allow time for unlink to happen.
        for ext in ["", ".garbage", ".bz2", ".gz"]
            rm(filename * ext, force = true)
        end
    end

    @testset "Model" begin
        for (format, model) in [
            (MathOptFormat.FORMAT_CBF, MathOptFormat.CBF.Model()),
            (MathOptFormat.FORMAT_LP, MathOptFormat.LP.Model()),
            (MathOptFormat.FORMAT_MOF, MathOptFormat.MOF.Model()),
            (MathOptFormat.FORMAT_MPS, MathOptFormat.MPS.Model()),
            (MathOptFormat.FORMAT_SDPA, MathOptFormat.SDPA.Model()),
        ]
            @test typeof(
                MathOptFormat.Model(format = format, filename = "foo.bar")
            ) == typeof(model)
        end
        @test_throws(
            ErrorException(
                "When `format==FORMAT_AUTOMATIC` you must pass a `filename`."
            ),
            MathOptFormat.Model(format = MathOptFormat.FORMAT_AUTOMATIC)
        )
        for (ext, model) in [
            (".cbf", MathOptFormat.CBF.Model()),
            (".lp", MathOptFormat.LP.Model()),
            (".mof.json", MathOptFormat.MOF.Model()),
            (".mps", MathOptFormat.MPS.Model()),
            (".sdpa", MathOptFormat.SDPA.Model()),
        ]
            @test typeof(MathOptFormat.Model(filename = "a$(ext)")) ==
                typeof(model)
            @test typeof(MathOptFormat.Model(filename = "a$(ext).gz")) ==
                typeof(model)
        end
        @test_throws(
            ErrorException("Unable to automatically detect format of a.b."),
            MathOptFormat.Model(filename = "a.b")
        )
    end
end


text:
function step(model::AbstractSDE, scheme::EulerMaruyama, t0, x0, Δw)
  μ = drift(model, t0, x0)
  σ = diffusion(model, t0, x0)
  x = x0 + μ * scheme.Δt + σ * Δw
  x
end

function _normal_transition_params(model, scheme::EulerMaruyama, t0, x0, x1)
  μ = x0 + drift(model, t0, x0) * scheme.Δt
  σ = diffusion(model, t0, x0)
  Σ = scheme.Δt * σ * σ'
  μ, Σ
end


text:
using PowerSystems
using NLsolve
const PSY = PowerSystems

raw_file = joinpath(TEST_FILES_DIR, "benchmarks/psse/EXAC1/TVC_System_32.raw")
dyr_file = joinpath(TEST_FILES_DIR, "benchmarks/psse/EXAC1/TVC_System.dyr")

sys = System(raw_file, dyr_file)
source = first(get_components(Source, sys))
PSY.set_X_th!(source, 0.01)

for l in get_components(PSY.PowerLoad, sys)
    PSY.set_model!(l, PSY.LoadModels.ConstantImpedance)
end


text:
using Distances

function weights(S::AbstractArray, R::AbstractArray, metric::Function = Distances.euclidean)
    ws = 0.0
    for s in S
        for r in R
            ws += metric(s,r)
        end
    end
    ws
end

function weights_pairwise(S::AbstractMatrix, R::AbstractMatrix, metric::Distances.PreMetric = Distances.Euclidean())
    sum(Distances.pairwise(metric,S,R))
end

function weights_half(S::AbstractArray, R::AbstractArray, metric::Function = Distances.euclidean)
    ws = 0.0
    for (i,s) in enumerate(S)
        for (j,r) in enumerate(R)
            if i > j
                ws += metric(s,r)
            end
        end
    end
    ws
end

# https://www.coursera.org/learn/cluster-analysis/lecture/jDuBD/6-7-internal-measures-for-clustering-validation
function intra_cluster_weights(C::AbstractArray)
    W_in = 0.0
    N_in = 0
    for (i,S) in enumerate(C)
        W_in += weights_half(S,S)
        N_in += binomial(length(S),2)
    end
    W_in, N_in
end

function intra_cluster_weights_pairwise(C::AbstractArray)
    W_in = 0.0
    N_in = 0
    for (i,S) in enumerate(C)
        # NOTE: problem 1: setindex! is not differentiable
        # NOTE: problem 2: Tracked Float64 can be convert to Float64
        W_in = W_in + weights_pairwise(S,S)
        N_in += binomial(size(S)[1],2)
    end
    0.5*W_in, N_in
end

function inter_cluster_weights(C::AbstractArray)
    W_out = 0.0
    N_out = 0
    for (i,S) in enumerate(C)
        for (j,R) in enumerate(C)
            if j > i
                W_out += weights(S,R)
                N_out += size(S)[1] * size(R)[1]
            end
        end
    end
    W_out, N_out
end

function inter_cluster_weights_pairwise(C::AbstractArray)
    W_out = 0.0
    N_out = 0
    for (i,S) in enumerate(C)
        for (j,R) in enumerate(C)
            if j > i
                W_out = W_out + weights_pairwise(S,R)
                N_out += size(S)[1] * size(R)[1]
            end
        end
    end
    W_out, N_out
end

function betacv(C::AbstractArray)
    W_in, N_in = intra_cluster_weights(C)
    W_out, N_out = inter_cluster_weights(C)

    (W_in / N_in) / (W_out / N_out)
end

function betacv_fused(C::AbstractArray; metric::Function = Distances.euclidean)
    K = length(C)
    W_in, W_out = 0.0, 0.0
    N_in, N_out = 0,   0

    for i in 1:K
        # intra cluster weights
        for k in 1:length(C[i])
            for l in 1:length(C[i])
                if k > l
                    W_in += metric(C[i][k],C[i][l])
                end
            end
        end
        N_in += binomial(length(C[i]),2)
        # inter cluster weights
        for j in 1:K
            if j > i
                for s in C[i]
                    for r in C[j]
                        W_out += metric(s,r)
                    end
                end
                N_out += size(C[i])[1] * size(C[j])[1]
            end
        end
    end

    (W_in / N_in) / (W_out / N_out)
end

function betacv_pairwise(C::AbstractArray)
    W_in, N_in = intra_cluster_weights_pairwise(C)
    W_out, N_out = inter_cluster_weights_pairwise(C)

    (W_in / N_in) / (W_out / N_out)
end

max_stars_repo_name: swiesend/FluxMO

max_issues_repo_name: swiesend/FluxMO

max_forks_repo_name: swiesend/FluxMO


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

using Markdown
using InteractiveUtils

# ╔═╡ 5248e43c-3761-11eb-1d42-bd8f8b9f19a7
begin
	import Pkg
	Pkg.add("FFTW")
	Pkg.add("Plots")
	Pkg.add("DSP")
	Pkg.add("PlutoUI")
	Pkg.add("WAV")
	Pkg.add("FileIO")
	Pkg.add("WebIO")
	Pkg.add("OffsetArrays")
    using FFTW
    using Plots
    using DSP
    using PlutoUI
    using WAV
    using FileIO
    using WebIO
	using OffsetArrays
end

# ╔═╡ 49d931b2-3761-11eb-1bc2-25a1d3b38eba
md"# Fast Fourier Transforms"

# ╔═╡ 556bc8c6-3761-11eb-0609-059485893468
md"### Playing with sound"

# ╔═╡ de377bcc-3971-11eb-1920-01ae73829171
md"This downloads five simple two second clips.  The first four are individual notes being hummed, and the last is the result of overlaying these four tracks."

# ╔═╡ 6374807c-3761-11eb-03e6-0f45dd0fd4e4
begin
	files = download.([
		"https://www.dropbox.com/s/3zygy6n72prdv9x/Sound1.wav?dl=0",
		"https://www.dropbox.com/s/qxb79tw6jhillym/Sound2.wav?dl=0",
		"https://www.dropbox.com/s/evo14tvdujkcssc/Sound3.wav?dl=0",
		"https://www.dropbox.com/s/sfxpt3erqelc7jt/Sound4.wav?dl=0",
		"https://www.dropbox.com/s/xlavzvtoaxe0oe7/Sound1234.wav?dl=0",
	])
end

# ╔═╡ af9a9c96-3803-11eb-1ec1-6315d95b164a
sounds = wavread.(files)

# ╔═╡ ffd94950-3803-11eb-36f9-93f943d7af8f
signals = [sound[1][:, 1] for sound in sounds]

# ╔═╡ d26c7384-396b-11eb-03f1-59108a815886
md"When you zoom in on one of the first four sounds, it looks like a simple periodic pattern with some frequency."

# ╔═╡ c6bf9c00-3803-11eb-12f6-87afda3aea9f
plot(signals[1], x_lims=(20000, 22000))

# ╔═╡ 0d9bca86-396c-11eb-1b96-9987e8e188bc
md"The final sound is the sum of these first four, and gives a more complicated wave"

# ╔═╡ 1140ee78-3804-11eb-34b5-6354d4db9642
plot(signals[5], x_lims=(20000, 22000))

# ╔═╡ 2291650e-396c-11eb-26e4-b791a57b87e4
md"The FFT pulls out the frequency components of each sound, which for musical notes like these looks like a spike above the base frequency itself, with smaller spikes at whole number multiples of that frequency, which are the overtones.

In this case, when we say \"frequency\", unless we rescale the x-axis it won't be cycles per unit time, but instead cycles per T units of time where T is the duration of the clip.  The FFT just sees a list of sample values, it has no way of knowing how many samples per second there are, or that seconds are the unit of time we might care about, so the only cononical unit of time would be the clip's duration.
"

# ╔═╡ 605fc966-3804-11eb-2aef-f1f934d49e36
signal_ffts = fft.(signals)

# ╔═╡ a3aab82a-3804-11eb-3d2c-6f8e9eb9a12f
plot(
	[abs.(signal_ffts[i]) for i in 1:4],
	x_lims=(1, 1000)
)

# ╔═╡ a57789ae-3805-11eb-21fc-c7b776fa2627
plot(
	abs.(signal_ffts[5]),
	x_lims=(1, 1000)
)

# ╔═╡ 66a98e40-3761-11eb-1b60-c1a3af8a062d
md"### Writing our own DFT function"

# ╔═╡ f15076fe-396d-11eb-0de9-1380a93497b9
md"By reading off the mathematical definition directly, we can implement the discrete fourier transform ourselves fairly simply."

# ╔═╡ 2839a190-3808-11eb-2554-99897f4542ca
function dft(signal)
	N = length(signal)
	# Precompute all the roots of unity we'll need.
	# Zeta here refers to exp(-2π * im / N).
	# OffsetArray lets us treat this like a zero-indexed array,
	# so that zeta_powers[n] = zeta^n
	zeta_powers = OffsetArray(
		[exp(-2π * im * n / N) for n in 0:(N - 1)],
		0:(N-1)
	)
	return [
		sum(
			signal[n + 1] * zeta_powers[(n * f) % N]
			for n in 0:(N - 1)
		)
		for f in 0:(N - 1)
	]
end

# ╔═╡ 429fe2ba-396e-11eb-0e01-13a0db7f9f03
md"We can verify that this gives the same output as the fft function, however by noting the runtime displayed in the lower right of each julia cell, you can see that our implmentation is much slower.  Fundeamentally this is because our naieve implementation has O(N^2) complexity, whereas the FFT using a more clever algorithm which is O(N*log(N))"

# ╔═╡ 93caa8dc-396e-11eb-36ea-2fe87a5b7a16
N = 10000 # Change this and note the runtimes below

# ╔═╡ 114c966a-3809-11eb-2b84-f1323e9b720d
signal = signals[1][1:N]

# ╔═╡ ac0cb6e2-396e-11eb-1f46-1d8b2100bb44
signal_fft = fft(signal)

# ╔═╡ 156b7572-3809-11eb-0d63-c9be4201e6ba
signal_dft = dft(signal)

# ╔═╡ 41f22476-3809-11eb-305a-51caf3dde226
isapprox(signal_dft, signal_fft)

# ╔═╡ Cell order:
# ╟─49d931b2-3761-11eb-1bc2-25a1d3b38eba
# ╟─5248e43c-3761-11eb-1d42-bd8f8b9f19a7
# ╟─556bc8c6-3761-11eb-0609-059485893468
# ╟─de377bcc-3971-11eb-1920-01ae73829171
# ╠═6374807c-3761-11eb-03e6-0f45dd0fd4e4
# ╠═af9a9c96-3803-11eb-1ec1-6315d95b164a
# ╠═ffd94950-3803-11eb-36f9-93f943d7af8f
# ╟─d26c7384-396b-11eb-03f1-59108a815886
# ╠═c6bf9c00-3803-11eb-12f6-87afda3aea9f
# ╟─0d9bca86-396c-11eb-1b96-9987e8e188bc
# ╠═1140ee78-3804-11eb-34b5-6354d4db9642
# ╟─2291650e-396c-11eb-26e4-b791a57b87e4
# ╠═605fc966-3804-11eb-2aef-f1f934d49e36
# ╠═a3aab82a-3804-11eb-3d2c-6f8e9eb9a12f
# ╠═a57789ae-3805-11eb-21fc-c7b776fa2627
# ╟─66a98e40-3761-11eb-1b60-c1a3af8a062d
# ╟─f15076fe-396d-11eb-0de9-1380a93497b9
# ╠═2839a190-3808-11eb-2554-99897f4542ca
# ╟─429fe2ba-396e-11eb-0e01-13a0db7f9f03
# ╠═93caa8dc-396e-11eb-36ea-2fe87a5b7a16
# ╠═114c966a-3809-11eb-2b84-f1323e9b720d
# ╠═ac0cb6e2-396e-11eb-1f46-1d8b2100bb44
# ╠═156b7572-3809-11eb-0d63-c9be4201e6ba
# ╠═41f22476-3809-11eb-305a-51caf3dde226

max_stars_repo_name: mathijsvdv/ComputationalThinking

max_issues_repo_name: mathijsvdv/ComputationalThinking

max_forks_repo_name: mathijsvdv/ComputationalThinking


text:
module Fermionic
using SparseArrays
using LinearAlgebra

include("base_functions.jl")
include("operators.jl")
include("operators_fixed.jl")
include("states.jl")
include("states_fixed.jl")
include("correlations.jl")
include("logic_gates.jl")
include("mixed.jl")

export Op, dim, basis, a, ad, ada, aad, aa, adad, adf, af, vacuum
export fixed, basis_m,  Op_fixed
export State, st, ope, typ, rhosp, rhoqsp, non_zero
export State_fixed, nume, fixed_state, unfixed_state
export eigensp, ssp, eigenqsp, sqsp, majorization_sp, majorization_qsp, n_avg, coef, rhom, rhomd, rhom2, trp
export sigma_x, sigma_y, sigma_z, phase, hadamard, ucnot, swap
export rhosp_mixed, eigensp_mixed
end # module


text:
@async function nested()
    return 42
end
@async function async_main()
    @assert(wait(nested()) == 42)
    println("OK")
end
@async function echo_server()
    server = listen(2001)
    while true
        sock = wait(accept(server))
        @async function writer()
            while isopen(sock)
                data = wait(readline(sock))
                wait(write(sock, upper(data)))
            end
        end
        wait(writer())
    end
end
if abspath(PROGRAM_FILE) == @__FILE__
    run(asyncio, async_main())
    run(asyncio, echo_server())
end


text:
include("../src/SFFM.jl")
using LinearAlgebra, Plots

## define the model(s)
@testset "$modelfile" for modelfile in ["testModel1.jl"; "testModel2.jl"; "testModel3.jl"]
    include(modelfile)

    ## section 4.3: the marginal stationary distribution of X
    ## mesh
    Δ = 0.4
    Nodes = collect(approxBounds[1, 1]:Δ:approxBounds[1, 2])
    NBases = 2
    Basis = "lagrange"
    mesh = SFFM.DGMesh(
        model,
        Nodes,
        NBases,
        Basis = Basis,
    )
    Ψₓ = SFFM.PsiFunX(model)
    ξₓ = SFFM.MakeXiX(model, Ψₓ)

    pₓ, πₓ, Πₓ, Kₓ = SFFM.StationaryDistributionX(model, Ψₓ, ξₓ)

    analyticX = (
            pm = [pₓ[:];0;0],
            distribution = πₓ(SFFM.CellNodes(mesh)),
            x = SFFM.CellNodes(mesh),
            type = "density"
        )

    @test sum(Πₓ(40))-1 < 1e-6

    @testset "DG witn NBases $NBases" for NBases in 1:2
            mesh = SFFM.DGMesh(
                model,
                Nodes,
                NBases,
                Basis = Basis,
            )

            # compute the marginal via DG
            All = SFFM.MakeAll( model, mesh, approxType = "projection")
            @test isapprox(sum(All.B.B,dims=2), zeros(size(All.B.B,1)), atol = sqrt(eps()))
            @test -sum(All.B.BDict["++"],dims=2) ≈ sum(All.B.BDict["+-"],dims=2) + sum(All.B.BDict["+0"],dims=2)
            @test -sum(All.B.BDict["--"],dims=2) ≈ sum(All.B.BDict["-+"],dims=2) + sum(All.B.BDict["-0"],dims=2)
            @test -sum(All.B.BDict["00"],dims=2) ≈ sum(All.B.BDict["0-"],dims=2) + sum(All.B.BDict["0+"],dims=2)
            @test -sum(All.D["++"](s=0),dims=2) ≈ sum(All.D["+-"](s=0),dims=2)
            @test sum(All.D["-+"](s=0),dims=2) ≈ -sum(All.D["--"](s=0),dims=2)

            Ψ = SFFM.PsiFun(All.D)
            @test sum(Ψ,dims=2) ≈ ones(size(Ψ,1))

            # the distribution of X when Y first returns to 0
            ξ = SFFM.MakeXi(All.B.BDict, Ψ)
            @test sum(ξ) ≈ 1

            marginalX, p, K = SFFM.MakeLimitDistMatrices(
                All.B.BDict,
                All.D,
                All.R.RDict,
                Ψ,
                ξ,
                mesh,
                model,
            )

            # convert marginalX to a distribution for plotting
            Dist = SFFM.Coeffs2Dist(
                    model,
                    mesh,
                    marginalX,
                    type="probability",
                )
            @test sum(Dist.pm) + sum(Dist.distribution) ≈ 1
    end
end

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
"""
```
default_settings!(m::AbstractModel)
```

Default Settings are constructed, initialized and added to `m.settings`.
"""
function default_settings!(m::AbstractModel)

    settings = m.settings

    # I/O File locations
    saveroot = normpath(joinpath(dirname(@__FILE__), "..","save"))
    datapath = normpath(joinpath(dirname(@__FILE__), "..","save","input_data"))

    settings[:saveroot] = Setting(:saveroot, saveroot, "Root of data directory structure")
    settings[:dataroot] = Setting(:dataroot, datapath, "Input data directory path")

    # Data settings for released and conditional data. Default behavior is to set vintage
    # of data to today's date.
    vint = Dates.format(now(), "yyyymmdd")
    settings[:data_vintage] = Setting(:data_vintage, vint, true, "vint",
        "Data vintage")
    settings[:data_id] = Setting(:data_id, 3,
        "Dataset identifier")

    return settings
end

"""
```
default_test_settings!(m::AbstractModel)
```

The following Settings are constructed, initialized and added to
`m.test_settings`. Their purposes are identical to those in
`m.settings`, but these values are used to test package.

### I/O Locations and identifiers
- `saveroot::Setting{String}`: A temporary directory in /tmp/
- `dataroot::Setting{String}`: dsgeroot/test/reference/
- `data_vintage::Setting{String}`: \"_REF\"
"""
function default_test_settings!(m::AbstractModel)

    test = m.test_settings

    # I/O
    dataroot = normpath(joinpath(dirname(@__FILE__), "..", "test", "reference", "input_data"))
    saveroot = mktempdir()

    # General
    test[:saveroot] = Setting(:saveroot, saveroot,
        "Where to write files when in test mode")
    test[:dataroot] = Setting(:dataroot, dataroot,
        "Location of input files when in test mode" )
    test[:data_vintage] = Setting(:data_vintage, "REF", true, "vint",
        "Reference data identifier")

    return test
end


text:
"""
    Atom{N, T} <: Predicate{N}

An atomic predicate of arity `N`.

### Fields

- `p` -- Julia predicate
"""
struct Atom{N, T} <: Predicate{N}
    p::T

    function Atom(p::T; N::Int=1) where {T}
        return new{Val{N}, T}(p)
    end
end

function Base.:(==)(a1::Atom, a2::Atom)
    return a1.p == a2.p
end

# function-like evaluation
@inline function (a::Atom)(args...)
    evaluate(a, args...)
end

function evaluate(a::Atom, args...)
    assert_same_length(a, args...)
    return a.p(args...)
end


text:
import CGinJulia.widemult, CGinJulia.det2, CGinJulia.det3, CGinJulia.det4, CGinJulia.det4ck

## I don't want Float64 to widen to BigFloat, but to Double64
@testset "wide multiplication and DoubleFloats" begin
	f = 1f0; g = 2f0; h = 2f0^22-1
	df = 1.0; dg = 2.0; dh = 2.0^52-1
	@test(typeof(h) == Float32)
	@test(typeof(widemult(f,g))==Float64)
	@test(typeof(widemult(df,dg))==Double64)
	@test(widemult(g,h)+2 == 2.0^23)
	@test(widemult(dg,dh)+2 == 2.0^53)
end
##

@testset "determinants" begin
	e1 = (1,0,0,0); e2 = (0,1,0,0); e3 = (0,0,1,0); e4 = (0,0,0,1);
	r = (1,2,3,5); s = (7,11,13,17); t = (1, 2, 4, 8); u = (1, 10, 100, 1000)
	@test(det2(e1,e2, 1,2) == 1)
	@test(det3(e1,e2,e3, 1,2,3) == 1)
	@test(det4(e1,e2,e3,e4, 1,2,3,4) == 1)
	@test(det2(r,s, 1,2) == -3)
	@test(det3(r,s,t, 1,2,3) == -3)
	@test(det4(r,s,t,u, 1,2,3,4) == -2160)
	@test(det4(r,s,t,u, 1,2,3,4) == det4ck(r,s,t,u, 1,2,3,4))
	@test(det4(r,s,t,u, 1,2,3,4) == det4dbl(r,s,t,u, 1,2,3,4))
end

max_stars_repo_name: snoeyink/CGinJulia

max_issues_repo_name: snoeyink/CGinJulia

max_forks_repo_name: snoeyink/CGinJulia


text:
const opts = Base.JLOptions()
const inline_flag = opts.can_inline == 1 ? `` : `--inline=no`
const cov_flag = (opts.code_coverage == 1) ? `--code-coverage=user` :
                 (opts.code_coverage == 2) ? `--code-coverage=all` :
                 ``

function run_test(script)
    srvrscript = joinpath(dirname(@__FILE__), script)
    srvrcmd = `$(joinpath(Sys.BINDIR, "julia")) $cov_flag $inline_flag $script`
    println("Running tests from ", script, "\n", "="^60)
    ret = run(srvrcmd)
    println("Finished ", script, "\n", "="^60)
    nothing
end

run_test("runtests_capi.jl")
run_test("runtests_base.jl")


text:
# Delimited files with a single column

Base.length(d :: DataFrameX) = size(d.df, 1);
row_indices(d :: DataFrameX) = 1 : (length(d) - 1);
# This is where the actual data are expected
col_index(d :: DataFrameX) = 2;
data_vector(d :: DataFrameX) =
    Vector{Double}(d.df[row_indices(d), col_index(d)]);

## -------------  Read delimited file

function read_by_x(fPath :: AbstractString)
    df = read_delim_file_to_df(fPath);
    # No of rows does not include header
    nr, nc = size(df);
    @check nc == 2
    @check nr > 3
    return DataFrameX(df);
end

read_vector_by_x(fPath) = data_vector(read_by_x(fPath));

# ----------------------


text:
using CuArrays
using NNlib: logsoftmax, logσ

# Cost functions

mse(ŷ, y) = sum((ŷ .- y).^2) * 1 // length(y)

function _crossentropy(ŷ::AbstractVecOrMat, y::AbstractVecOrMat, weight::Nothing)
  return -sum(y .* log.(ŷ)) * 1 // size(y, 2)
end

function _crossentropy(ŷ::AbstractVecOrMat, y::AbstractVecOrMat, weight::Number)
  return -sum(y .* log.(ŷ)) .* weight * 1 // size(y, 2)
end

function _crossentropy(ŷ::AbstractVecOrMat, y::AbstractVecOrMat, weight::AbstractVector)
  return -sum(y .* log.(ŷ) .* weight) * 1 // size(y, 2)
end

crossentropy(ŷ::AbstractVecOrMat, y::AbstractVecOrMat; weight=nothing) = _crossentropy(ŷ, y, weight)

function logitcrossentropy(logŷ::AbstractVecOrMat, y::AbstractVecOrMat; weight = 1)
  return -sum(y .* logsoftmax(logŷ) .* weight) * 1 // size(y, 2)
end

"""
    binarycrossentropy(ŷ, y; ϵ=eps(ŷ))

Return `-y*log(ŷ + ϵ) - (1-y)*log(1-ŷ + ϵ)`. The ϵ term provides numerical stability.

    julia> binarycrossentropy.(σ.([-1.1491, 0.8619, 0.3127]), [1, 1, 0.])
    3-element Array{Float64,1}:
    1.4244
    0.352317
    0.86167
"""
binarycrossentropy(ŷ, y; ϵ=eps(ŷ)) = -y*log(ŷ + ϵ) - (1 - y)*log(1 - ŷ + ϵ)

# Re-definition to fix interaction with CuArrays.
CuArrays.@cufunc binarycrossentropy(ŷ, y; ϵ=eps(ŷ)) = -y*log(ŷ + ϵ) - (1 - y)*log(1 - ŷ + ϵ)

"""
    logitbinarycrossentropy(logŷ, y)

`logitbinarycrossentropy(logŷ, y)` is mathematically equivalent to `binarycrossentropy(σ(logŷ), y)`
but it is more numerically stable.

    julia> logitbinarycrossentropy.([-1.1491, 0.8619, 0.3127], [1, 1, 0.])
    3-element Array{Float64,1}:
     1.4244
     0.352317
     0.86167
"""
logitbinarycrossentropy(logŷ, y) = (1 - y)*logŷ - logσ(logŷ)

# Re-definition to fix interaction with CuArrays.
CuArrays.@cufunc logitbinarycrossentropy(logŷ, y) = (1 - y)*logŷ - logσ(logŷ)

"""
    normalise(x::AbstractArray; dims=1)

Normalises `x` to mean 0 and standard deviation 1, across the dimensions given by `dims`. Defaults to normalising over columns.

    julia> a = reshape(collect(1:9), 3, 3)
    3×3 Array{Int64,2}:
     1  4  7
     2  5  8
     3  6  9

    julia> normalise(a)
    3×3 Array{Float64,2}:
     -1.22474  -1.22474  -1.22474
      0.0       0.0       0.0
      1.22474   1.22474   1.22474

    julia> normalise(a, dims=2)
    3×3 Array{Float64,2}:
     -1.22474  0.0  1.22474
     -1.22474  0.0  1.22474
     -1.22474  0.0  1.22474
"""
function normalise(x::AbstractArray; dims=1)
  μ′ = mean(x, dims = dims)
  σ′ = std(x, dims = dims, mean = μ′, corrected=false)
  return (x .- μ′) ./ σ′
end


text:
#= samin.jl. Translation of samin.cc code (in the optim package
   for GNU Octave) to Julia

  Copyright (C) 2016   Michael Creel   <michael.creel@uab.es>

  This program is free software you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation either version 2 of the License, or
  (at your option) any later version.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA

  References:

  The code follows the article
  Goffe, William L. (1996) "SIMANN: A Global Optimization Algorithm
	using Simulated Annealing " Studies in Nonlinear Dynamics & Econometrics
  	Oct96, Vol. 1 Issue 3.

  The code uses the same names for control variables,
  for the most part. A notable difference is that the initial
  temperature is found automatically to ensure that the active
  bounds when the temperature begins to reduce cover the entire
  parameter space (defined as a n-dimensional
  rectangle that is the Cartesian product of the
  (lb_i, ub_i), i = 1,2,..n

  Also of note:
  Corana et. al., (1987) "Minimizing Multimodal Functions of Continuous
	Variables with the "Simulated Annealing" Algorithm",
 	ACM Transactions on Mathematical Software, V. 13, N. 3.

  Goffe, et. al. (1994) "Global Optimization of Statistical Functions
 	with Simulated Annealing", Journal of Econometrics,
 	V. 60, N. 1/2.

usage: x, obj, convergence, details = samin(f,
                                            x_init,
                                            lb,
                                            ub,
                                            nt,
                                            ns,
                                            rt,
                                            maxevals,
                                            neps,
                                            functol,
                                            paramtol,
                                            verbosity)
Arguments:
REQUIRED
* f: objective function
* x_init: starting value
* lb:  vector of lower bounds
* ub: vector of upper bounds
KEYWORDS
* nt:  integer: # reduce temperature every nt*ns*dim(x_init) evaluations
* ns:  integer: # adjust bounds every ns*dim(x_init) evaluations
* rt:  (0 < rt <1): geometric temperature reduction factor: when temp changes, new temp is t=rt*t
* maxevals:  integer: limit on function evaluations
* neps:  integer:  number of previous best values the final result is compared to
* functol:    (> 0): the required tolerance level for function value
                   comparisons
* paramtol:   (> 0): the required tolerance level for parameters
* verbosity - scalar: 0, 1, or 2.
	* 0 = no screen output
	* 1 = only final results to screen
	* 2 = summary every temperature change

Returns:
* x: the minimizer
* obj: the value of f() at x
* convergence:
	0 if no convergence within maxevals function evaluations
	1 if normal convergence to a point interior to the parameter space
	2 if convergence to point very near bounds of parameter space
	  (suggest re-running with looser bounds)
* details: a px3 matrix. p is the number of times improvements were found.
           The columns record information at the time an improvement was found
           * first: cumulative number of function evaluations
           * second: temperature
           * third: function value

Example: see samin_test
=#

function samin(obj_fn, x, lb, ub; nt=5, ns=5, rt=0.5, maxevals=1000000, neps=5, functol=1e-8, paramtol=1e-5, verbosity=2)
    n = size(x,1) # dimension of parameter
    #  Set initial values
	nacc = 0 # total accepted trials
	t = 2.0 # temperature - will initially rise or fall to cover parameter space. Then it will fall
	converge = 0 # convergence indicator 0 (failure), 1 (normal success), or 2 (convergence but near bounds)
	coverage_ok = 0 # has parameter space been covered? When turns to 1, temperature starts to fall
	# most recent values, to compare to when checking convergend
	fstar = typemax(Float64)*ones(neps)
	# Initial obj_value
	xopt = copy(x)
    f = obj_fn(x)
	fopt = copy(f) # give it something to compare to
	func_evals = 0 # total function evaluations (limited by maxeval)
	details = [func_evals t fopt]
    bounds = ub - lb
	# check for out-of-bounds starting values
	for i = 1:n
		if(( x[i] > ub[i]) || (x[i] < lb[i]))
            @printf("samin: initial parameter %d out of bounds\n", i)
            converge = 0
			return xopt, fopt, converge, details
		end
    end
	# main loop, first increase temperature until parameter space covered, then reduce until convergence
	while converge==0
		# statistics to report at each temp change, set back to zero
		nup = 0
		nrej = 0
		nnew = 0
		ndown = 0
		lnobds = 0

		# repeat nt times then adjust temperature
		for m = 1:nt
			# repeat ns times, then adjust bounds
            nacp = zeros(n)
			for j = 1:ns
				# generate new point by taking last and adding a random value
				# to each of elements, in turn
				for h = 1:n
					# new Sept 2011, if bounds are same, skip the search for that vbl.
                    # Allows restrictions without complicated programming
					if (lb[h] != ub[h])
						xp = copy(x)
						xp[h] += (2.0 * rand() - 1.0) * bounds[h]
						if((xp[h] < lb[h]) || (xp[h] > ub[h]))
							xp[h] = lb[h] + (ub[h] - lb[h]) * rand()
							lnobds += 1
                        end
						# Evaluate function at new point
						fp = obj_fn(xp)
						func_evals += 1
						#  Accept the new point if the function value decreases
						if (fp <= f)
							x = copy(xp)
							f = copy(fp)
							nacc += 1 # total number of acceptances
							nacp[h] += 1 # acceptances for this parameter
							nup += 1
							#  If lower than any other point, record as new optimum
							if(fp < fopt)
								xopt = copy(xp)
								fopt = copy(fp)
								nnew +=1
								details = [details; [func_evals t fp]]
                            end
						# If the point is higher, use the Metropolis criteria to decide on
						# acceptance or rejection.
						else
							p = exp(-(fp - f) / t)
							if(rand() < p)
								x = copy(xp)
								f = copy(fp)
								nacc += 1
								nacp[h] += 1
								ndown += 1
							else
                                nrej += 1
                            end
                        end
					    # If maxevals exceeded, terminate the algorithm
					    if(func_evals >= maxevals)
						    if(verbosity >= 1)
							    println("================================================")
					            println("SAMIN results")
							    println("NO CONVERGENCE: MAXEVALS exceeded")
					            @printf("\n     Obj. value:  %16.5f\n\n", fopt)
        			            println("       parameter      search width")
                                for i=1:n
                                    @printf("%16.5f  %16.5f \n", xopt[i], bounds[i])
                                end
					            println("================================================")
                            end
                            converge = 0
						    return xopt, fopt, converge, details
                        end
                    end
                end
            end
            #  Adjust bounds so that approximately half of all evaluations are accepted
			test = 0
			for i = 1:n
				if (lb[i] != ub[i])
					ratio = nacp[i] / ns
                    if(ratio > 0.6) bounds[i] = bounds[i] * (1.0 + 2.0 * (ratio - 0.6) / 0.4) end
                    if(ratio < .4) bounds[i] = bounds[i] / (1.0 + 2.0 * ((0.4 - ratio) / 0.4)) end
                    # keep within initial bounds
					if(bounds[i] >= (ub[i] - lb[i]))
						bounds[i] = ub[i] - lb[i]
						test += 1
                    end
				else
                    test = test + 1 # make sure coverage check passes for the fixed parameters
                end
            end
            nacp = nacp - nacp # set back to zero
			# check if we cover parameter space, if we have yet to do so
            if (coverage_ok != 1) coverage_ok = (test == n) end
        end

        # intermediate output, if desired
		if(verbosity == 2)
			println("samin: intermediate results before next temperature change")
			println("temperature: ", t)
			println("current best function value: ", fopt)
			println("total evaluations so far: ", func_evals)
			println("total moves since last temperature reduction: ", nup + ndown + nrej)
			println("downhill: ", nup)
			println("accepted uphill: ", ndown)
			println("rejected uphill: ", nrej)
			println("out of bounds trials: ", lnobds)
			println("new minima this temperature: ", nnew)
            println()
			println("       parameter      search width")
            for i=1:n
                @printf("%16.5f  %16.5f \n", xopt[i], bounds[i])
            end
			println()
        end
        # Check for convergence, if we have covered the parameter space
		if (coverage_ok)

            # last value close enough to last neps values?
			fstar[1] = f
			test = 0
            for i=1:neps
                test += (abs(f - fstar[i]) > functol)
            end
			test = (test > 0) # if different from zero, function conv. has failed

			# last value close enough to overall best?
			if (((fopt - f) <= functol) && (!test))
				# check for bound narrow enough for parameter convergence
				for i = 1:n
					if (bounds[i] > paramtol)
						converge = 0 # no conv. if bounds too wide
						break
					else
                        converge = 1
                    end
                end
            end
            # check if optimal point is near boundary of parameter space, and change message if so
            if (converge == 1) && (lnobds > 0)
                converge = 2
            end
            # Like to see the final results?
			if (converge > 0)
				if (verbosity >= 1)
					println("================================================")
					println("SAMIN results")
                    if (converge == 1)
                        println("==> Normal convergence <==")
                    end
					if (converge == 2)
						println("==> WARNING <==: Last point satisfies convergence criteria,")
						println("but is near boundary of parameter space.")
						println(lnobds, " out of  ", (nup+ndown+nrej), " evaluations were out of bounds in the last round.")
						println("Expand bounds and re-run, unless this is a constrained minimization.")
                    end
					@printf("\n     Obj. value:  %16.5f\n\n", fopt)
        			println("       parameter      search width")
                    for i=1:n
                        @printf("%16.5f  %16.5f \n", xopt[i], bounds[i])
                    end
					println("================================================")
                end
            end
            # Reduce temperature, record current function value in the
			# list of last "neps" values, and loop again
			t = rt * t
            for i = neps:-1:2
                fstar[i] = fstar[i-1]
            end
			f = copy(fopt)
			x = copy(xopt)
		else  # coverage not ok - increase temperature quickly to expand search area
            t = t*t
            for i = neps:-1:2
                fstar[i] = fstar[i-1]
            end
			f = fopt
			x = xopt
        end
    end
    return xopt, fopt, converge, details
end