text:
# Convenience wrapper
save_subdomain_temperature!(m, state, aux) =
    save_subdomain_temperature!(m, m.moisture, state, aux)

using KernelAbstractions: @print

"""
    save_subdomain_temperature!(
        m::AtmosModel,
        moist::EquilMoist,
        state::Vars,
        aux::Vars,
    )

Updates the subdomain sensible temperature, given:
 - `m`, an `AtmosModel`
 - `moist`, an `EquilMoist` model
 - `state`, state variables
 - `aux`, auxiliary variables
"""
function save_subdomain_temperature!(
    m::AtmosModel,
    moist::EquilMoist,
    state::Vars,
    aux::Vars,
)
    N_up = n_updrafts(m.turbconv)
    ts = recover_thermo_state(m, state, aux)
    ts_up = new_thermo_state_up(m, state, aux, ts)
    ts_en = new_thermo_state_en(m, state, aux, ts)

    @unroll_map(N_up) do i
        aux.turbconv.updraft[i].T = air_temperature(ts_up[i])
    end
    aux.turbconv.environment.T = air_temperature(ts_en)
    return nothing
end

# No need to save temperature for DryModel.
function save_subdomain_temperature!(
    m::AtmosModel,
    moist::DryModel,
    state::Vars,
    aux::Vars,
) end

max_forks_repo_name: jm-c/CLIMA


text:
@testset "transformed_parameters.jl" begin

    @testset "Two Gaussians" begin

        # Two Gaussian peaks with σ(μ) an affine function
        #
        model = Gaussian_Peak() + Gaussian_Peak()

        # model θ is:
        # idx:  1   2  3   4   5   6
        #  θ : h1, μ1, σ1, h2, μ2, σ2
        #
        # hence:
        #
        src  = [2,5]
        dest = [3,6]

        # The σ(μ) relation: σ(μ=1)=1 and σ(μ=100)=10 for [θ̂A, θ̂B]=[1, 1]
        #
        f_σ_μ = Map_Affine(1.0=>1.0,100.0=>10.0)

        # Now the model_σ_μ parameters are:
        #
        # [ h1, μ1, h2, μ2, θ̂A, θ̂B ]
        #
        # the initial model is called with:
        #
        # [ h1, μ1, σ1 = affine(μ1, θ̂A, θ̂B), h2, μ2, , σ2 = affine(μ2, θ̂A, θ̂B) ]
        #
        model_σ_μ = model_σ_μ = create_model_transform_src_insert_dest(model, f_σ_μ, src=>dest)

        # Test some computations
        #
        model_σ_μ_θ = Float64[1,5, 2, 95, 1, 1]
        model_θ = get_model_θ(model_σ_μ, model_σ_μ_θ)

        @test model_θ ≈ [1.0, 5.0, 1.3636363636363638, 2.0, 95.0,
        9.545454545454547]

        X = [1:10.0:101;]
        Y = eval_y(model_σ_μ, X, model_σ_μ_θ)

        @test Y ≈ [0.013538440136652369, 6.252150377485093e-5,
        1.7807741296584817e-13, 3.4628155839082785e-10,
        2.2470318744914768e-7, 4.865734936177968e-5,
        0.003515984927922131, 0.08478224390932892, 0.6822163961931842,
        1.8318872636421548, 1.6414765303163972]

    end

end

ext: jl


text:
function main!(r)
    while (!isborder(r, Nord) && !isborder(r, West))==true
        move!(r, Nord)
        move!(r, West)
        putmarker!(r)
    end
    while ismarker(r)
        move!(r,Ost)
        move!(r,Sud)
    end
    while (!isborder(r, Sud) && !isborder(r, West))==true
        move!(r, Sud)
        move!(r, West)
        putmarker!(r)
    end
    while ismarker(r)
        move!(r,Ost)
        move!(r,Nord)
    end
    while (!isborder(r, Nord) && !isborder(r, Ost))==true
        move!(r, Nord)
        move!(r, Ost)
        putmarker!(r)
    end
    while ismarker(r)
        move!(r,West)
        move!(r,Sud)
    end
    while (!isborder(r, Sud) && !isborder(r, Ost))==true
        move!(r, Sud)
        move!(r, Ost)
        putmarker!(r)
    end
    while ismarker(r)
        move!(r,Nord)
        move!(r,West)
    end
    putmarker!(r)
end


text:
using BEAST
using CompScienceMeshes

o, x, y, z = euclidianbasis(3)

Γ = readmesh(Pkg.dir("BEAST","examples","sphere2.in"))
RT = raviartthomas(Γ)

κ = 1.0
t = Maxwell3D.singlelayer(gamma=im*κ)
E = Maxwell3D.planewave(direction=z, polarization=x, wavenumber=κ)
e = (n × E) × n

@hilbertspace j
@hilbertspace k

efie = @discretise t[k,j]==e[k]  j∈RT k∈RT

u = gmres(efie)
println("Solution computed")

isdefined(:plotresults) || (plotresults = false)
plotresults && (postproc = true)
isdefined(:postproc) || (postproc = false)

if postproc

    Φ, Θ = [0.0], linspace(0,π,100)
    pts = [point(cos(ϕ)*sin(θ), sin(ϕ)*sin(θ), cos(θ)) for ϕ in [0.0] for θ in Θ]
    ffd = potential(MWFarField3D(im*κ), pts, u, RT)
    println("far field computed")

    fcr, geo = facecurrents(u, RT)
    println("Face currents computed")

    nx, nz = 50, 100
    X, Z = linspace(-2,2,nx), linspace(-4,4,nz)
    grid = [point(x,0,z) for x in X, z in Z]
    nfd = potential(MWSingleLayerField3D(κ), grid, u, RT)
    nfd = reshape(nfd, (nx,nz))
    nfd .-= E.(grid)
    println("near field computed.")

end

if plotresults
    @eval begin
        using PlotlyJS
        include(Pkg.dir("CompScienceMeshes","examples","plotlyjs_patches.jl"))
        t1 = scatter(x=Θ, y=real.(norm.(ffd)))
        t2 = patch(geo, real.(norm.(fcr)))
        t3 = heatmap(x = X, y = Z, z = clamp.(real.(norm.(nfd)), 0.0, 2.0))
    end
end

#plot(t1) # uncomment to plot the far field (line2d)
#plot(t2) # uncomment to plot the induced current (mesh3d)
#plot(t3) # uncomment to plot the near field (heatmap)

nothing;


text:
using Parameters
using DifferentialEquations
using Blink
using Interact
using Plots
gr()

# Bring in the model
include("lv_model.jl")

# Build a Graphical User Interface (GUI) to explore the parameters
let
    u0 = [0.5, 0.5]
    t_span = (0.0, 100.0)

    ui = @manipulate for α12 in slider(0.1:0.01:2.0; label = "α12", value = 0.8),
                         α21 in slider(0.1:0.01:2.0; label = "α21", value = 0.3),
                         K1 in slider(0.1:0.01:3.0; label = "K1", value = 1.5),
                         K2 in slider(0.1:0.01:3.0; label = "K2", value = 1.0),
                         r1 in slider(0.1:0.01:10.0; label = "r1", value = 1.0),
                         r2 in slider(0.1:0.01:10.0; label = "r2", value = 2.0)
        p = LVPar(α12 = α12, α21 = α21, K1 = K1, K2 = K2, r1 = r1, r2 = r2)
        prob = ODEProblem(lv_comp, u0, t_span, p)
        sol = solve(prob, reltol = 1e-8)
        plot(sol)
        xlabel!("time")
        ylabel!("Density")
    end

    w = Blink.AtomShell.Window()
    body!(w, ui)
end

max_stars_repo_name: McCannLab/TheoryClub

max_issues_repo_name: McCannLab/TheoryClub

max_forks_repo_name: McCannLab/TheoryClub


text:
# addprocs(7)
# ask first for user convenience
using Caesar

include("/home/dehann/Documents/blandauthlocal.jl")
cloudGraph, addrdict = standardcloudgraphsetup(nparticles=true, addrdict=addrdict)

using RoME
using CloudGraphs, Neo4j

# using IncrementalInference

# session = addrdict["session"]
session = "TestHexagonalDrive11"
addrdict["num particles"] = "100"
Nparticles = parse(Int, addrdict["num particles"])
println("Attempting to solve session $(session) with $(Nparticles) particles per marginal...")

fg = initfg(sessionname=session, cloudgraph=cloudGraph)

setBackendWorkingSet!(fg.cg.neo4j.connection, session)

# TODO -- incremental graph and subgraphs are works in progress
while true
  println("=================================================")
  fg = initfg(sessionname=session, cloudgraph=cloudGraph)

  setBackendWorkingSet!(fg.cg.neo4j.connection, session)

  println("get local copy of graph")

  # removeGenericMarginals!(conn) # function should not be necessary, but fixes a minor bug following elimination algorithm
  if fullLocalGraphCopy!(fg)
    tree = wipeBuildNewTree!(fg,drawpdf=true)
    # removeGenericMarginals!(conn)

    # while true # repeat while graph unchanged
      # okay now do the solve
      inferOverTree!(fg, tree, N=Nparticles)
      # if true # current hack till test is inserted
      #   break;
      # end
    # end
  else
    sleep(0.2)
  end
end

# fg = initfg(sessionname=session, cloudgraph=cloudGraph)
# fullLocalGraphCopy!(fg, conn)
#
# ls(fg)
#
# tree = wipeBuildNewTree!(fg,drawpdf=true)
#
# inferOverTreeR!(fg, tree, N=100)

  #


text:
import ..update!, ..select_action, ..reset!

using Random
using StaticArrays

using Ramnet
using Ramnet.Models.AltDiscriminators: FunctionalDiscriminator
using Ramnet.Optimizers
import Ramnet.Loss

struct BinaryActionPolicyLoss <: Loss.AbstractLoss end

# Policy approximator with only two admissible actions: {-1, 1}.
struct BinaryActionPolicy{OS,T<:Real,O<:StaticArray{Tuple{OS},T,1},E<:AbstractEncoder{T}}
    opt::FunctionalOptimizer{BinaryActionPolicyLoss}
    f::FunctionalDiscriminator{1,T,E}
    rng::MersenneTwister
end

function Loss.grad(::BinaryActionPolicyLoss, p::BinaryActionPolicy{OS,T,O,E}, observation::O, action::Int, G) where {OS,T,O,E}
    sign(action) * (1 - probability(p, action, observation)) * G
end

function BinaryActionPolicy(::Type{O}, n::Int, encoder::E; start_learning_rate::Float64 = 0.1, end_learning_rate::Float64 = 1e5, learning_rate_decay::Int = 1000, epochs::Int = 1, partitioner::Symbol = :uniform_random, seed::Union{Nothing,UInt} = nothing) where {OS,T<:Real,O<:StaticArray{Tuple{OS},T,1},E<:AbstractEncoder{T}}
    BinaryActionPolicy{OS,T,O,E}(
        FunctionalOptimizer(BinaryActionPolicyLoss(); start_learning_rate, end_learning_rate, learning_rate_decay, epochs),
        FunctionalDiscriminator{1}(OS, n, encoder, partitioner; seed),
        isnothing(seed) ? MersenneTwister() : MersenneTwister(seed)
    )
end

# Compute π(action | observation)
function probability(p::BinaryActionPolicy{OS,T,O,E}, action::Int, observation::O) where {OS,T<:Real,O<:StaticArray{Tuple{OS},T,1},E<:AbstractEncoder{T}}
    return 1 / (1 + exp(sign(action) * predict(p.f, observation)[1]))
end

# Compute the mode of π(action | observation)
function mode(p::BinaryActionPolicy{OS,T,O,E}, observation::O) where {OS,T<:Real,O<:StaticArray{Tuple{OS},T,1},E<:AbstractEncoder{T}}
    return signbit(predict(p.f, observation)[1]) ? 1 : -1
end

# TODO: Validate action (must be ether -1 or +1)
function update!(p::BinaryActionPolicy{OS,T,O,E}, observation::O, action::Int, G) where {OS,T<:Real,O<:StaticArray{Tuple{OS},T,1},E<:AbstractEncoder{T}}
    gradient = Optimizers.grad(p.opt, p, observation, action, G)

    train!(p.f, observation, learning_rate(p.opt) * gradient)

    return nothing
end

function select_action(p::BinaryActionPolicy{OS,T,O,E}, observation::O; deterministic::Bool = false) where {OS,T<:Real,O<:StaticArray{Tuple{OS},T,1},E<:AbstractEncoder{T}}
    if deterministic
        return mode(p, observation)
    end

    return rand(p.rng) < probability(p, 1, observation) ? 1 : -1
end

function reset!(p::BinaryActionPolicy; seed::Union{Nothing,UInt} = nothing)
    Ramnet.reset!(p.opt)
    Ramnet.reset!(p.f; seed)
    Random.seed!(p.rng, seed)

    return nothing
end


text:
module JuliaFFT

export fft

"""
    fft(xs)

Computes the discrete Fourier transform of the array `xs`.

``y_k = \\sum_{j=1}^N \\exp\\left[ - 2 \\pi i \\frac{(j-1)(k-1)}{N} \\right] x_j``

The implementation is fully generic; as long as the input types have
exponentiation, addition, and multiplication (by complex numbers) defined, the
algorithm will work.  Not particularly fast (about 2ms for FFT of 1024 element
array on my modern-ish MacBook).

Currently only works for 1D arrays.

"""
function fft(xs)
    n = size(xs,1)

    if n == 1
        return xs
    else
        xso = xs[1:2:end]
        xse = xs[2:2:end]

        ffo = fft(xso)
        ffe = fft(xse)

        runity = exp.(-2im*pi*range(0, stop=n/2-1)/n)

        return vcat(ffo .+ runity .* ffe, ffo .- runity.*ffe)
    end
end

end # module


text:
include("../../../iron_line_profile.jl")

verticle_line_x = fill(6.4, 100)
verticle_line_y = LinRange(0, 1.2, 100)

mass_vars = [10, 100, 10e6]
spin_vars = [0.0, 0.5, 0.998]
angle_vars = [15.0, 40.0, 85.0]

mass_fixed = 10
spin_fixed = 0.998
angle_fixed = 40.0

x_vals_vals = []
bins_vals = []
labels = []
for spin in spin_vars
    x_vals, bins = iron_line_profile(;
                                    spin=spin,
                                    mass=mass_fixed,
                                    obs_angle=angle_fixed,
                                    tolerance=1e-12,
                                    dtmax=0.5,
                                    size_multiplier=10,
                                    resolution=3000,
                                    output="data",
                                    nbins=200,
                                    fov=10,
                                    )
    push!(x_vals_vals, x_vals)
    push!(bins_vals, bins)
    push!(labels, "a = $spin")
end
fixed_labels=permutedims(labels)
plt1 = plot(x_vals_vals, bins_vals, label=fixed_labels, grid=false, framestyle=:box)
plot!(
    verticle_line_x,
    verticle_line_y,
    linestyle=:dash,
    linecolor=:black,
    label=false,
    ylims=(0,1.2),
    lw=2
    )
png(plt1, "iron_line_comp_spin.png")

x_vals_vals = []
bins_vals = []
labels = []
for angle in angle_vars
    x_vals, bins = iron_line_profile(;
                                    obs_angle=angle,
                                    spin=spin_fixed,
                                    mass=mass_fixed,
                                    tolerance=1e-12,
                                    dtmax=0.5,
                                    size_multiplier=10,
                                    resolution=3000,
                                    output="data",
                                    nbins=200,
                                    fov=10,
                                    )
    push!(x_vals_vals, x_vals)
    push!(bins_vals, bins)
    push!(labels, "i = $angle")
end
fixed_labels=permutedims(labels)
plt1 = plot(x_vals_vals, bins_vals, label=fixed_labels, grid=false, framestyle=:box)
plot!(
    verticle_line_x,
    verticle_line_y,
    linestyle=:dash,
    linecolor=:black,
    label=false,
    ylims=(0,1.2),
    lw=2
    )
png(plt1, "iron_line_comp_angle.png")


text:
using Stheno: UnaryMean, BinaryMean, BinaryKernel, BinaryCrossKernel, LhsCross, RhsCross,
    OuterCross, OuterKernel

@benchset "compose" begin
    let
        m1, m2 = CustomMean(sin), CustomMean(cos)
        @benchset "UnaryMean CPU" create_benchmarks(UnaryMean(exp, m1))
        @benchset "BinaryMean CPU" create_benchmarks(BinaryMean(+, m1, m2))
    end
    let
        k1, k2 = EQ(), OneKernel()
        @benchset "BinaryKernel EQ CPU" create_benchmarks(BinaryKernel(*, k1, k1))
        @benchset "BinaryKernel Const. CPU" create_benchmarks(BinaryKernel(*, k2, k2))
    end
    let
        k, f = EQ(), abs2
        @benchset "LhsCross abs2 * EQ CPU" create_benchmarks(LhsCross(f, k))
        @benchset "RhsCross EQ * abs2 CPU" create_benchmarks(RhsCross(k, f))
        @benchset "OuterCross f * EQ * f CPU" create_benchmarks(OuterCross(f, k))
        @benchset "OuterKernel f * EQ * f CPU" create_benchmarks(OuterKernel(f, k))
    end
end


text:
"""
    range(length)

A range of the specified `length` starting at one. Alias for `Base.OneTo(length)``.

Equivalent to range(1; length) in Base.
"""
Base.range(length::Integer) = Base.OneTo(length)

"""
    range(start, length)

A [`UnitRange`](@ref) of the specified `length` starting at `start`.

See also `start:stop` for a UnitRange with a specific `stop`.

Equivalent to range(start; length) in Base
"""
Base.range(start::Number, length::Integer) = Base.range(start; length)

"""
    range(start, stop, length)

A range from `start` to `stop` (inclusive) of the specified `length`.

Equivalent to range(start, stop; length) in Base
"""
Base.range(start::Number, stop, length::Integer) = Base.range(start, stop; length)

"""
    range(start, step, stop, length)

This form allows all four arguments to be provided positionally. One of them must be `nothing`.

Equivalent to Base._range
"""
Base.range(start, step, stop, length) = Base._range(start, step, stop, length)


text:
module Sim
include("./name2idx/parameters.jl")
include("./name2idx/species.jl")
include("./set_model.jl")

using .C
using .V
using Printf
using DelayDiffEq

const dt = 0.01

function solvedde(
        diffeq::Function,
        u0::Vector{Float64},
        history::Vector{Float64},
        tspan::Tuple{Float64,Float64},
        p::Vector{Float64},
        tau::Float64)
    h(p,t) = history
    lags = [tau]
    prob = DDEProblem(diffeq,u0,h,tspan,p;constant_lags=lags)
    alg = MethodOfSteps(BS3())
    sol = solve(prob,alg,saveat=dt)
    return sol
end

const tspan = (0.0,48.0)
const t = collect(tspan[1]:dt:tspan[end])

const n_condition = 4

p53_conc = Matrix{Float64}(undef, length(t), n_condition)

for i in 1:n_condition
    local p::Vector{Float64}
    local u0::Vector{Float64}

    if i == 1
        p = param_values_U2OS()
        u0 = initialValues_U2OS()
        p[C.Eto] = 1
    elseif i == 2
        p = param_values_U2OS()
        u0 = initialValues_U2OS()
        p[C.Eto] = 100
    elseif i == 3
        p = param_values_MCF7()
        u0 = initialValues_MCF7()
        p[C.Eto] = 5
    elseif i == 4
        p = param_values_MCF7()
        u0 = initialValues_MCF7()
        p[C.Eto] = 200
    end

    sol = solvedde(diffeq,u0,u0,tspan,p,2.1)
    for j in eachindex(t)
        p53_conc[j,i] = sol.u[j][V.p53]
    end
end

end  # module


text:
#=
parameters:
- Julia version:
- Author: Dan
- Date: 2021-07-29
=#
include("global.jl")


text:
module OceanBoundaryLayerParameterEstimation

pushfirst!(LOAD_PATH, joinpath(@__DIR__, "..", ".."))

using ParameterEstimocean

export
    # lesbrary.jl
    TwoDaySuite, FourDaySuite, SixDaySuite,
    lesbrary_ensemble_simulation,
    estimate_η_covariance,
    SyntheticObservationsBatch,
    two_day_suite_path_4m, four_day_suite_path_4m, six_day_suite_path_4m,
    two_day_suite_path_2m, four_day_suite_path_2m, six_day_suite_path_2m,
    two_day_suite_path_1m, four_day_suite_path_1m, six_day_suite_path_1m,

    # closure_parameters.jl
    bounds,
    default,
    ParameterSet,
    named_tuple_map,
    CATKEParametersRiDependent,
    CATKEParametersRiIndependent,
    CATKEParametersRiDependentConvectiveAdjustment,
    CATKEParametersRiIndependentConvectiveAdjustment,
    RiBasedParameterSet,

    # mcmc.jl
    markov_chain,

    # eki_visuals.jl
    plot_parameter_convergence!,
    plot_pairwise_ensembles!,
    plot_error_convergence!,

    # uq_visuals.jl
    plot_mcmc_densities,
    plot_mcmc_densities!,

    # visualize_profile_predictions.jl
    visualize!,
    open_output_file

include("lesbrary.jl")
include("closure_parameters.jl")
include("mcmc.jl")
include("eki_visuals.jl")
include("visualize_profile_predictions.jl")
include("uq_visuals.jl")

end # module


text:
module Nash

#packages
using Conda
if ! ("sympy" in get.(Conda.parseconda(`list`),"name",nothing))
    Conda.runconda(`install sympy -c anaconda -y`)
end

using Base.Threads
using Distributions
using LightGraphs
using LinearAlgebra
using CDDLib
using Polyhedra
using MeshCat
using Combinatorics
using IterTools
using QuantEcon
using Plots
using Roots
using SymPy
using Manifolds

#utilities.jl
export outer, plot_br
#games.jl
export generate_game, random_2players_game, random_nplayers_game
#equilibria.jl
export get_payoff, best_reply, is_nash_q, iterate_best_reply
#markov.jl
export game2markov, plot_markov
#symmetric.jl
export random_symmetric_2players_game, find_symmetric_nash_equilibrium_2players_game, create_symmetries_graph, check_equality_condition, find_all_equalities
#NEOptim.jl
export vFunction, cart_prod_simplices
#replicator.jl
export create_replicator_eqs

#files
include("utilities.jl")
include("games.jl")
include("equilibria.jl")
include("markov.jl")
include("symmetric.jl")
include("NEOptim.jl")
include("replicator.jl")
end


text:
using CUDA
using MadNLPTests

function _compare_gpu_with_cpu(kkt_system, n, m, ind_fixed)
    madnlp_options = Dict{Symbol, Any}(
        :kkt_system=>MadNLP.DENSE_CONDENSED_KKT_SYSTEM,
        :linear_solver=>MadNLPLapackGPU,
        :print_level=>MadNLP.ERROR,
    )

    nlp = MadNLPTests.DenseDummyQP(; n=n, m=m, fixed_variables=ind_fixed)

    h_ips = MadNLP.InteriorPointSolver(nlp; option_dict=copy(madnlp_options))
    MadNLP.optimize!(h_ips)

    # Init KKT on the GPU
    TKKTGPU = kkt_system{Float64, CuVector{Float64}, CuMatrix{Float64}}
    opt = MadNLP.Options(; madnlp_options...)
    # Instantiate Solver with KKT on the GPU
    d_ips = MadNLP.InteriorPointSolver{TKKTGPU}(nlp, opt; option_linear_solver=copy(madnlp_options))
    MadNLP.optimize!(d_ips)

    # # Check that both results match exactly
    @test h_ips.cnt.k == d_ips.cnt.k
    @test h_ips.obj_val ≈ d_ips.obj_val atol=1e-10
    @test h_ips.x ≈ d_ips.x atol=1e-10
    @test h_ips.l ≈ d_ips.l atol=1e-10
end

@testset "MadNLPGPU ($(kkt_system))" for kkt_system in [
        MadNLP.DenseKKTSystem,
        MadNLP.DenseCondensedKKTSystem,
    ]
    @testset "Size: ($n, $m)" for (n, m) in [(10, 0), (10, 5), (50, 10)]
        _compare_gpu_with_cpu(kkt_system, n, m, Int[])
    end
    @testset "Fixed variables" begin
        n, m = 20, 0 # warning: setting m >= 1 does not work in inertia free mode
        _compare_gpu_with_cpu(kkt_system, n, m, Int[1, 2])
    end
end


text:
module RegressionFormulae

using StatsModels
using Combinatorics
using Base.Iterators

using StatsModels: apply_schema

const TermTuple = NTuple{N, AbstractTerm} where N
const Schemas = Union{StatsModels.Schema, StatsModels.FullRank}

include("fulldummy.jl")
include("power.jl")
include("nesting.jl")

end # module


text:
int main() {
  for(boolean i : new int[20]) {
    boolean x = i;
  }
  return 0;
}

max_stars_repo_name: fredefox/javalette

max_issues_repo_name: fredefox/javalette

max_forks_repo_name: fredefox/javalette


text:
using GLVisualize, GLWindow

if !isdefined(:runtests)
	window = glscreen()
	timesignal = bounce(linspace(0f0,1f0,360))
end

function volume_data(N)
	vol 	= Float32[sin(x/15.0)+sin(y/15.0)+sin(z/15.0) for x=1:N, y=1:N, z=1:N]
	max 	= maximum(vol)
	min 	= minimum(vol)
	vol 	= (vol .- min) ./ (max .- min)
end

vol = visualize(volume_data(128), :iso, isovalue=timesignal)
_view(vol, window)

if !isdefined(:runtests)
	renderloop(window)
end


text:
"""

    random_graph(n,p=0.5, has_self_loops=false)

generate a random directed graph `g` with `n` nodes and the probability to include edge is `p`. If `has_self_loops`, `g` will include self loops.
"""
function random_graph(n::Integer, p::Real = 0.5; has_self_loops = false)
    g = DiGraph()
    for i = 1:n
        ind = 1
        for j = ind:n
            add_node!(g, j)
            if rand() <= p && (i != j || has_self_loops)
                add_edge!(g, i, j)
            end
        end
    end
    return g
end

"""
    random_evolving_graph(nn, nt, p=0.5, is_directed=true, has_self_loops=false)
    random_evolving_graph(g, nn, nt, p=0.5, is_directed=true, has_self_loops=false)

Generate a random evolving graph `g` with `nn` nodes, `nt` timestamps. The random evolving graph has integer nodes and timestamps.
 The probability to include each edge is equal to `p`. If `has_self_loops`, we allow `g` to have self loops. If `g` is not given, we generate a random IntAdjacencyList evolving graph.

# Example

```
julia> using EvolvingGraphs

julia> random_evolving_graph(10, 4, 0.1)
Directed IntAdjacencyList (10 nodes, 37 static edges, 4 timestamps)

julia> g = EvolvingGraph()
Directed EvolvingGraph 0 nodes, 0 static edges, 0 timestamps

julia> random_evolving_graph(g, 20, 3, 0.5)
Directed EvolvingGraph 20 nodes, 558 static edges, 3 timestamps
```
"""
function random_evolving_graph(g::AbstractEvolvingGraph, nn::Int, nt::Int, p::Real = 0.5;
                               has_self_loops = false)
    for t in 1:nt
        for i = 1:nn
            is_directed(g) ? ind =1 : ind = i
            for j = ind:nn
                    if (rand() <= p && (i != j || has_self_loops))
                        add_edge!(g, i, j, t)
                    end
            end
        end
    end
    g
end
function random_evolving_graph(nn::Int, nt::Int,p::Real = 0.5;
                               is_directed = true,has_self_loops  = false)
    g = IntAdjacencyList(nn, nt, is_directed = is_directed)
    random_evolving_graph(g, nn, nt, p, has_self_loops = has_self_loops)
end


text:
# Autogenerated wrapper script for LMDB_jll for x86_64-linux-gnu
export liblmdb, mdb_copy, mdb_dump, mdb_load, mdb_stat

JLLWrappers.@generate_wrapper_header("LMDB")
JLLWrappers.@declare_library_product(liblmdb, "liblmdb.so")
JLLWrappers.@declare_executable_product(mdb_copy)
JLLWrappers.@declare_executable_product(mdb_dump)
JLLWrappers.@declare_executable_product(mdb_load)
JLLWrappers.@declare_executable_product(mdb_stat)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        liblmdb,
        "lib/liblmdb.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        mdb_copy,
        "bin/mdb_copy",
    )

    JLLWrappers.@init_executable_product(
        mdb_dump,
        "bin/mdb_dump",
    )

    JLLWrappers.@init_executable_product(
        mdb_load,
        "bin/mdb_load",
    )

    JLLWrappers.@init_executable_product(
        mdb_stat,
        "bin/mdb_stat",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module PrettyEinsum
using TensorOperations

export @Σein!_str, @Σein_str, @Σein!, @Σein

include("parse.jl")
include("macrostr.jl")
include("macro.jl")

end


text:
using MCObservables

function test_simpleObservable(numtest, MCS)
  inform("Start testing SimpleObservable")

  nsigma = 3

  counts = zeros(Int, nsigma)

  for i in 1:numtest
    obs = SimpleObservable()
    for j in 1:MCS
      obs << randn()
    end
    m = abs(mean(obs))
    er = stderror(obs)
    for isigma in 1:nsigma
      if m < isigma*er
        counts[isigma] += 1
      end
    end
  end
  sigma = counts ./ numtest
  @show sigma

  @test 0.68 < sigma[1] < 0.69
  @test 0.95 < sigma[2] < 0.96
  @test 0.99 < sigma[3]
  inform("Success testing SimpleObservable")
end

result_simpleObservable = test_simpleObservable(100_000, 8192)


text:
using CSV
using DataFrames
using Statistics

"""
    CalculateIndex(fold::Integer, len::Integer, index::Integer)
Calcula el índice final de la partición `index`de  de un conjunto de `len` elementos, dividido en `fold` partes.
Si la división no es exacta se cogerán más elementos en las particiones primeras, todo elemento de la lista tiene asociado una partición.

Ejemplos:
# tamaño 7, tomado en 3 folds [3,5, 7], 123 45 67
l1 = [3,5,7]
l1_test = [CalculateIndex(3,7,i) for i in 1:3 ]
# tamaño 8, tomado en 3 folds [3,6,8] 123 456 78
l2 = [3,6,8]
l2_test = [CalculateIndex(3,8,i) for i in 1:3 ]
# tamaño 9, tomado en 4 folds [3,5, 7, 9] 123 45 67 89
l3 = [3,5,7, 9]
l3_test = [CalculateIndex(4,9,i) for i in 1:4 ]
"""
CalculateIndex(fold::Integer, len::Integer, index::Integer)=index*floor(Int, len/fold)+((mod(len,fold) >= (index)) ? index : mod(len,fold))

"""
    VerboseCrossValidation(data, labels, folds_number, learner_algorithm, file_name::String,
    umbral_tasa_reduccion=0.1, a = 0.5)
Divide el conjunto de datos en `folds_number` particiones disjuntas.
Repite `folds_number` veces:
entrena con `folds_number-1` conjuntos y valida con el restante.

La validación será la tasa media de acierto, también se evaluará el tiempo.
"""
function VerboseCrossValidation(data, labels, folds_number, learner_algorithm, file_name::String,
    umbral_tasa_reduccion=0.1, a = 0.5)

    # Datos que vamos a escribir
    tamaño = folds_number+2
    dfParticion = map( x -> "Partición $x", 1:(tamaño))
    dfTime = Array{Float64}(undef, tamaño)
    dfClasificacion = Array{Float64}(undef, tamaño)
    dfReduccion =  Array{Float64}(undef, tamaño)
    dfAgregacion = Array{Float64}(undef, tamaño)
    dfW = Array{Vector{Float64}}(undef, tamaño)
    # calculate folds
    len = length(labels)
    index = [
        CalculateIndex(folds_number,len,i) for i in 0:folds_number
        ]
    lk = ReentrantLock()

    Threads.@threads for i in 1:folds_number
        lock(lk)
            println("Comienza la iteración $i/$folds_number de CV para $file_name")
        unlock(lk)
        # select train data and test
        train_index = filter(
            x-> x<= index[i] || x > index[i+1] , 1:len
        )

        test_index = (index[i]+1):index[i+1]

        train_data = data[train_index , :]
        train_labels = labels[train_index]
        test_data = data[test_index, :]
        test_labels = labels[test_index]

        time = @elapsed begin
            # train
            clasificator, F_w, w = learner_algorithm(train_data, train_labels)
            # test
            estimations =  map(clasificator, eachrow(test_data))
        end
        # Calculamos la precisión
        aciertos = sum(estimations .== test_labels)
        accuracy = (aciertos/ (index[i+1] - index[i]) )* 100
        # Calculamos tasa de reducción valor entre [0,100]
        tasa_reducion = 100*sum(
            map(x-> (x < umbral_tasa_reduccion) ? 1 : 0, w)
            )/length(w)
        # Evaluación combinación de ambas
        evaluacion = a*accuracy + (1-a)*tasa_reducion

        # Guardamos datos que se escribirán en el fichero
       #            s -> x1000 ms
        dfTime[i] = time#*1000
        dfAgregacion[i] = evaluacion
        dfClasificacion[i] = accuracy
        dfReduccion[i] = tasa_reducion
        dfW[i] = w
        lock(lk)
        # habría que hacer esto con un cerrojo
            println("------------------------------------")
            println("Termina iteración $i/$folds_number de CV para $file_name con:
            Tiempo (s): $(dfTime[i])  tasa clasificación: $accuracy  tasa reducción: $tasa_reducion agregación: $evaluacion
            w = $w\n")
            println("------------------------------------")
        unlock(lk)

    end
    # Añadimos medias
    # Calculamos las medias
    tamaño -=1
    dfParticion[tamaño] = "Medias "
    dfTime[tamaño]= mean(dfTime[1:folds_number])
    dfClasificacion[tamaño]= mean(dfClasificacion[1:folds_number])
    dfReduccion[tamaño]= mean(dfReduccion[1:folds_number])
    dfAgregacion[tamaño]= mean(dfAgregacion[1:folds_number])
    dfW[tamaño]= mean(dfW[1:folds_number])
    # Calculamos las desviaciones típicas
    tamaño += 1
    dfParticion[tamaño] = "Desviación típica"
    dfTime[tamaño]= std(dfTime[1:folds_number])
    dfClasificacion[tamaño]= std(dfClasificacion[1:folds_number])
    dfReduccion[tamaño]= std(dfReduccion[1:folds_number])
    dfAgregacion[tamaño]= std(dfAgregacion[1:folds_number])
    dfW[tamaño]= std(dfW[1:folds_number])

    # Escribimos en el fichero
    DF = DataFrame(
        Nombre_Fila = dfParticion,
        Clasificación = dfClasificacion,
        Reducción = dfReduccion,
        Agregación = dfAgregacion,
        Tiempo = dfTime,
        Pesos = dfW
    )
    println("Se proceden a escribir en $file_name los datos $DF")
    CSV.write(file_name, DF)

    return DF
end

"""
    LeaveOneOut(data, labels, learner_algorithm)
"""
function LeaveOneOut(data, labels, learner_algorithm)
    # calculate folds
    len = length(labels)
    accuracy = 0
    for i in 1:len
        index = union(1:(i-1), (i+1):len)
        clasificator = learner_algorithm(data[index , :], labels[index])

        accuracy += ( clasificator(data[i]) == labels[i]) ? 0 : 1
    end
    accuracy = accuracy*100/len
    return accuracy
end

function VerboseCrossValidationS(data, labels, folds_number, learner_algorithm, file_name::String,
    umbral_tasa_reduccion=0.1, a = 0.5)

    # Datos que vamos a escribir
    tamaño = 1+2
    dfParticion = map( x -> "Partición $x", 1:(tamaño))
    dfTime = Array{Float64}(undef, tamaño)
    dfClasificacion = Array{Float64}(undef, tamaño)
    dfReduccion =  Array{Float64}(undef, tamaño)
    dfAgregacion = Array{Float64}(undef, tamaño)
    dfW = Array{Vector{Float64}}(undef, tamaño)
    # calculate folds
    len = length(labels)
    index = [
        CalculateIndex(folds_number,len,i) for i in 0:folds_number
        ]
    lk = ReentrantLock()

    Threads.@threads for i in 1:1
        lock(lk)
            println("Comienza la iteración $i/$folds_number de CV para $file_name")
        unlock(lk)
        # select train data and test
        train_index = filter(
            x-> x<= index[i] || x > index[i+1] , 1:len
        )

        test_index = (index[i]+1):index[i+1]

        train_data = data[train_index , :]
        train_labels = labels[train_index]
        test_data = data[test_index, :]
        test_labels = labels[test_index]

        time = @elapsed begin
            # train
            clasificator, F_w, w = learner_algorithm(train_data, train_labels)
            # test
            estimations =  map(clasificator, eachrow(test_data))
        end
        # Calculamos la precisión
        aciertos = sum(estimations .== test_labels)
        accuracy = (aciertos/ (index[i+1] - index[i]) )* 100
        # Calculamos tasa de reducción valor entre [0,100]
        tasa_reducion = 100*sum(
            map(x-> (x < umbral_tasa_reduccion) ? 1 : 0, w)
            )/length(w)
        # Evaluación combinación de ambas
        evaluacion = a*accuracy + (1-a)*tasa_reducion

        # Guardamos datos que se escribirán en el fichero
       #            s -> x1000 ms
        dfTime[i] = time#*1000
        dfAgregacion[i] = evaluacion
        dfClasificacion[i] = accuracy
        dfReduccion[i] = tasa_reducion
        dfW[i] = w
        lock(lk)
        # habría que hacer esto con un cerrojo
            println("------------------------------------")
            println("Termina iteración $i/$folds_number de CV para $file_name con:
            Tiempo (s): $(dfTime[i])  tasa clasificación: $accuracy  tasa reducción: $tasa_reducion agregación: $evaluacion
            w = $w\n")
            println("------------------------------------")
        unlock(lk)

    end
    # Añadimos medias
    # Calculamos las medias
    folds_number = 1
    tamaño -=1
    dfParticion[tamaño] = "Medias "
    dfTime[tamaño]= mean(dfTime[1:folds_number])
    dfClasificacion[tamaño]= mean(dfClasificacion[1:folds_number])
    dfReduccion[tamaño]= mean(dfReduccion[1:folds_number])
    dfAgregacion[tamaño]= mean(dfAgregacion[1:folds_number])
    dfW[tamaño]= mean(dfW[1:folds_number])
    # Calculamos las desviaciones típicas
    tamaño += 1
    dfParticion[tamaño] = "Desviación típica"
    dfTime[tamaño]= std(dfTime[1:folds_number])
    dfClasificacion[tamaño]= std(dfClasificacion[1:folds_number])
    dfReduccion[tamaño]= std(dfReduccion[1:folds_number])
    dfAgregacion[tamaño]= std(dfAgregacion[1:folds_number])
    dfW[tamaño]= std(dfW[1:folds_number])

    # Escribimos en el fichero
    DF = DataFrame(
        Nombre_Fila = dfParticion,
        Clasificación = dfClasificacion,
        Reducción = dfReduccion,
        Agregación = dfAgregacion,
        Tiempo = dfTime,
        Pesos = dfW
    )
    println("Se proceden a escribir en $file_name los datos $DF")
    CSV.write(file_name, DF)

    return DF
end


text:
function unpackdemo(D::Int = 2, M::Int = 3)

    W_array = collect( Utilities.generaterandomposdefmat(D) for m = 1:M )
    L_array = collect( cholesky(W_array[m]) for m = 1:M )

    b_array = collect( randn(D) for m = 1:M )
    c_array = rand(M)

    v = packbcW(b_array, c_array, W_array)

    b, c, W = unpackbcW(v,M)

    return v, b, c, W, b_array, c_array, W_array
end

max_stars_repo_name: RoyCCWang/RWJuliaUtilities

max_issues_repo_name: RoyCCWang/RWJuliaUtilities

max_forks_repo_name: RoyCCWang/RWJuliaUtilities


text:
@testset "Periodic Kernel" begin
    x = rand()*2
    v1 = rand(3)
    v2 = rand(3)
    r = rand(3)
    k = PeriodicKernel(r = r)
    @test kappa(k, x) ≈ exp(-0.5x)
    @test k(v1, v2) ≈ exp(-0.5 * sum(abs2, sinpi.(v1 - v2) ./ r))
    @test k(v1, v2) == k(v2, v1)
    @test PeriodicKernel(3)(v1, v2) == PeriodicKernel(r = ones(3))(v1, v2)
    @test repr(k) == "Periodic Kernel, length(r) = $(length(r))"

    # Standardised tests.
    TestUtils.test_interface(PeriodicKernel(r=[0.9]), Vector{Float64})
    TestUtils.test_interface(PeriodicKernel(r=[0.9, 0.9]), ColVecs{Float64})
    TestUtils.test_interface(PeriodicKernel(r=[0.8, 0.7]), RowVecs{Float64})

    # test_ADs(r->PeriodicKernel(r =exp.(r)), log.(r), ADs = [:ForwardDiff, :ReverseDiff])
    @test_broken "Undefined adjoint for Sinus metric, and failing randomly for ForwardDiff and ReverseDiff"
    test_params(k, (r,))
end


text:
module StructDatabaseMapping
export select_one, update!, delete!, drop_table!, clean_table!,
       update!, exists, select_all

export DBMapper, register!, create_table, DBId,
       Nullable, analyze_relations, ForeignKey,  Model, getid, Foreign,
       Cascade, SetNull, configure_relation,
       Restrict, database_kind

using Dates
using Requires
using JSON
using DBInterface

import Base.insert!

abstract type Connection end
abstract type DatabaseKind end
abstract type Model end

include("Connection/Pool.jl")

"""
    mutable struct Field

Internal field representation
"""
mutable struct Field
    name::Symbol
    struct_field::Symbol
    type::Type
    nullable::Bool
    default::Union{Any, Nothing}
end
function Field(name::Symbol, type::Type)
    primary_key = false
    nullable = false
    db_field_name = name
    if type <: ForeignKey
        db_field_name = Symbol(string(db_field_name) * "_id")
    elseif type <: AbstractNullable
        nullable = true
    end

    return Field(db_field_name, name, type, nullable, primary_key)
end

mutable struct Key
    field::Array{Field}
    has_auto_value::Bool
    is_primary::Bool
end
Key(f::Field; primary::Bool=false, auto_value::Bool=false) = Key([f], auto_value, primary)

abstract type OnEventAction end
struct Cascade <: OnEventAction end
struct Restrict <: OnEventAction end
struct SetNull <: OnEventAction end
struct SetDefault <: OnEventAction
    value
end
struct Set <: OnEventAction
end
struct DoNothing <: OnEventAction end

mutable struct Relation
    referenced_table::String
    referenced_field::Symbol
    local_field::Symbol
    on_delete::OnEventAction
    on_update::OnEventAction
end
function Relation(referenced_table::String, referenced_field::Symbol, local_field::Symbol)
    return Relation(referenced_table, referenced_field, local_field, DoNothing(), DoNothing())
end

mutable struct Table
    name::String
    data_type
    fields::Dict{Symbol, Field}
    relations::Dict{Field, Relation}
    primary_key::Key
end
function idfield(t::Table)  :: Symbol
    return t.primary_key.field[1].name
end
isprimarykey(f::Field, t::Table) = t.primary_key.field[1] === f
function fieldlist(t::Table) :: Array{Field}
    return sort(collect(values(t.fields)), by=x->x.name)
end

function foreignfields(t::Table) :: Array{Field}
    return sort(collect(keys(t.relations)))
end

mutable struct DBMapper
    tables::Dict{DataType, Table}
    pool::ConnectionPool
    dirty::Bool
    function DBMapper(database_builder::Function)
        return new(Dict{DataType,Table}(), SimplePool(database_builder), false)
    end
end

@enum RELATION_TYPE  ONE_TO_MANY=1 ONE_TO_ONE=2

abstract type AbstractNullable{T} end
"""
    has_value(data::T) where T<:AbstractNullable
"""
has_value(data::T) where T<:AbstractNullable = !isnothing(data.x)

"""
    set!(data::T, elem::J) where T<:AbstractNullable{J}
"""
set!(data::T, elem::J) where T<:AbstractNullable{J} where J = data.x = elem

"""
    mutable struct Nullable{T} <: AbstractNullable{T}

Optional field
"""
mutable struct Nullable{T} <: AbstractNullable{T}
    x::Union{T,Nothing}
end

function Base.convert(::Type{F}, x::T) where F<:AbstractNullable{T}  where T
    return F(x)
end

function Base.convert(::Type{F}, x::Nothing) where F<:AbstractNullable{T}  where T
    return F(nothing)
end

"""
    mutable struct Nullable{T} <: AbstractNullable{T}

Model identifier
"""
mutable struct DBId{T} <: AbstractNullable{T}
    x::Union{T,Nothing}
end
function DBId{T}() where T
    return DBId{T}(nothing)
end

"""
    idfield(mapper::DBMapper, T::DataType)  :: Symbol

Return the identifier field name for the given type
"""
function idfield(mapper::DBMapper, T::DataType)  :: Symbol
    return idfield(mapper.tables[T])
end
"""
    function getid(elem::T, mapper::DBMapper) where T<:Model

Return the identifier value for the given type
"""
function getid(elem::T, mapper::DBMapper) where T<:Model
    getfield(elem, idfield(mapper, T)).x
end

"""
    function setid!(elem::T, mapper::DBMapper, id)  where T<:Model

Set the identifier value for the given type
"""
function setid!(elem::T, mapper::DBMapper, id)  where T<:Model
    id_field = getfield(elem, idfield(mapper, T))
    set!(id_field, id)
end
"""
    idfieldtype(::Type{T}, mapper::DBMapper) where T<:Model

Obtain the general id field type for a given type
for a field with type DBId{Integer} should return DBId{Integer}
"""
idfieldtype(::Type{T}, mapper::DBMapper) where T<:Model =  fieldtype(T, idfield(mapper, T))

"""
    function idtype(::Type{T}, mapper::DBMapper) where T<:Model

Obtain the element type of the identifier field for a given type
For a field with type DBId{Integer} should return Integer
"""
function idtype(::Type{T}, mapper::DBMapper) where T<:Model
    id_field_type = fieldtype(T, idfield(mapper, T))
    return element_type(id_field_type)
end

"""
    mutable struct Nullable{T} <: AbstractNullable{T}

Foreign Key model field

The foreign key field can be constructed with the T element
or with nothing. The `loaded` field is used in the lazy loading of the foreign
element. When the foreign element is loaded the `loaded` attributed is set
to true

```
struct ModelA <: Model ... end
struct ModelB <: Model
    ...
    foreign_field::ForeignKey{ModelA}
end
```
"""
mutable struct ForeignKey{T} <: AbstractNullable{T}
    data::Union{T, Nothing}
    loaded::Bool
end
function ForeignKey{T}(;data::Union{T,Nothing}=nothing, loaded::Bool=false) where T<:Model
    return ForeignKey{T}(data, loaded)
end
function ForeignKey{T}(data::Union{T,Nothing}) where T<:Model
    return ForeignKey{T}(data, true)
end
const Foreign{T} = Union{T, ForeignKey{T}}

function Base.get(v::ForeignKey{T}, mapper::DBMapper) where T
    if v.loaded
        return v.data
    end
    params = Dict{Symbol, Any}(idfield(mapper, T)=>getid(v.data, mapper))
    v.data = select_one(mapper, T;params...)
    v.loaded = true
    return v.data
end

"""
    function register!(mapper::DBMapper, d::Type{T}; table_name::String="") where T <: Model

Register an element of type T into the mapper.
A Table type is created for the given Model and is stored in the dabase. This function must be called
for each type you want to use with this library.
"""
function register!(mapper::DBMapper, d::Type{T}; table_name::String="") where T <: Model
    mapper.dirty = true
    if table_name == ""
        table_name = String(split(lowercase(string(T)), ".")[end])
    end
    fields = Dict{Symbol, Field}()
    primary_key = nothing
    for (field_name, field_type) in zip(fieldnames(d), fieldtypes(d))
        db_field = Field(field_name, field_type)
        if field_type <: DBId
            primary_key = Key(db_field; primary=true, auto_value=element_type(db_field.type) <: Integer)
        end
        fields[field_name] = db_field
    end
    if primary_key === nothing
        throw("The struct $(T) does not have a primary key. Use the type DBId")
    end
    table = Table(table_name,
                  T,
                  fields,
                  Dict{Field, Relation}(),
                  primary_key)
    mapper.tables[T] =  table
end

"""
    function analyze_relations(mapper::DBMapper)

Updates the relations dict of each table.

After calling this function the mapper state is not dirty.
"""
function analyze_relations(mapper::DBMapper)
    for table in values(mapper.tables)
        for field in fieldlist(table)
            if field.type <: ForeignKey
                referenced_table = mapper.tables[element_type(field.type)]
                referenced_field = referenced_table.primary_key.field[1]
                table.relations[field] = Relation(referenced_table.name,
                                                  referenced_field.name,
                                                  field.name)
            end
        end
    end
    mapper.dirty = false
end

"""
    function column_names(mapper::DBMapper, T::DataType) :: Array

Return the table field names for a given struct.
"""
function column_names(mapper::DBMapper, T::DataType) :: Array
    table = mapper.tables[T]
    return map(field->field.name, fieldlist(table))
end

"""
    function struct_fields(mapper::DBMapper, T::DataType) :: Array{Symbol}

Return the field names for a given struct.
Perhaps should be replaced directly with fieldnames
"""
function struct_fields(mapper::DBMapper, T::DataType) :: Array{Symbol}
    table = mapper.tables[T]
    return map(field->field.struct_field, fieldlist(table))
end

normalize(dbtype, x) = x
normalize(dbtype, x::DBId{T}) where T = x.x
normalize(dbtype, x::ForeignKey{T}) where T = normalize(dbtype, x.data.id)
normalize(dbtype, x::AbstractDict) where T = JSON.json(x)
function struct_field_values(mapper, elem::T; ignore_primary_key::Bool=true, fields::Array{Symbol}=Symbol[]) where T
    table = mapper.tables[T]
    column_names = []
    valuelist = []
    for field in fieldlist(table)
        if length(fields) > 0 && !(field.struct_field in fields)
            continue
        end
        if isprimarykey(field, table) && table.primary_key.has_auto_value && ignore_primary_key
            continue
        end
        field_value = getfield(elem, field.struct_field)

        if (field.nullable == true) && !has_value(field_value)
            continue
        end
        push!(column_names, field.name)
        push!(valuelist, normalize(mapper.pool.dbtype, getfield(elem, field.struct_field)))
    end
    return (column_names, valuelist)
end

function check_valid_type(mapper::DBMapper, T::DataType)
    if !haskey(mapper.tables, T)
        throw("Invalid type")
    end
    if mapper.dirty == true
        analyze_relations(mapper)
    end
end

"""
    function create_table(mapper::DBMapper, T::Type{<:Model}; if_not_exists::Bool=true)

Create the table for the given model
"""
function create_table(mapper::DBMapper, T::Type{<:Model}; if_not_exists::Bool=true)
    check_valid_type(mapper, T)
    create_table(mapper, mapper.pool.dbtype, T; if_not_exists=if_not_exists)
end
function create_table(mapper::DBMapper, dbtype::DataType, T::Type{<:Model}; if_not_exists::Bool=true)
    create_table(mapper, database_kind(dbtype), T; if_not_exists=if_not_exists)
end

"""
    function insert!(mapper::DBMapper, elem::T) where T <: Model

Insert the element in the database

# Arguments
- `mapper::DBMapper`: The database mapper
- `elem::T where T<:Model`: Instantied model to insert

```
struct Author <: Model ... end
insert!(mapper, Author(name="some name", age=30))
```
"""
function insert!(mapper::DBMapper, elem::T) where T <: Model
    check_valid_type(mapper, T)
    insert!(mapper, mapper.pool.dbtype, elem)
end
function insert!(mapper::DBMapper, dbtype::DataType, elem::T) where T <: Model
    insert!(mapper, database_kind(dbtype), elem)
end

"""
    function update!(mapper::DBMapper, elem::T; fields::Array{Symbol}=Symbol[]) where T<:Model

Insert the element in the database

# Arguments
- `mapper::DBMapper`: The database mapper
- `elem::T where T<:Model`: Instantied model to insert
- `fields::Array{Symbol}`: Optional. Array of fields to update.

```
struct Author <: Model ... end
author = Author(name="some name", age=30)
update!(mapper, author)
update!(mapper, author, fields=[:age])
```
"""
function update!(mapper::DBMapper, elem::T; fields::Array{Symbol}=Symbol[]) where T<:Model
    check_valid_type(mapper, T)
    update!(mapper, mapper.pool.dbtype, elem; fields=fields)
end
function update!(mapper::DBMapper, dbtype::DataType, elem::T; fields::Array{Symbol}=Symbol[]) where T <:Model
    update!(mapper, database_kind(dbtype), elem; fields=fields)
end

"""
    function select_one(mapper::DBMapper, T::Type{<:Model}; kwargs...)

Select one element from the database

# Arguments
- `mapper::DBMapper`: The database mapper
- `T::DataType`: Datatype of a registered model we want to select
- `kwargs`: fields we want to search for. The param pk cand be used as generic way
to identify the primary key of the struct

```
struct Author <: Model ... end
select_one(mapper, Author, name="Borges")
```
"""
function select_one(mapper::DBMapper, T::Type{<:Model}; kwargs...)
    check_valid_type(mapper, T)
    params = Dict(kwargs...)
    if haskey(params, :pk)
        params[idfield(mapper, T)] = params[:pk]
        pop!(params, :pk)
    end
    return select_one(mapper, mapper.pool.dbtype, T; params...)
end
function select_one(mapper::DBMapper, dbtype, T::Type{<:Model}; kwargs...)
    return select_one(mapper, database_kind(dbtype), T; kwargs...)
end

"""
    clean_table!(mapper::DBMapper, T::Type{<:Model})

Remove all elements of the type T.

# Arguments
- `mapper::DBMapper`: The database mapper
- `T::Type{<:Model}`: Datatype of a registered model

In cases when possible the structure where those elements are stored (tables in relational case)
"""
function clean_table!(mapper::DBMapper, T::Type{<:Model})
    check_valid_type(mapper, T)
    clean_table!(mapper, mapper.pool.dbtype, T)
end
function clean_table!(mapper::DBMapper, dbtype::DataType, T::Type{<:Model})
    clean_table!(mapper, database_kind(dbtype), T)
end

"""
    function drop_table!(mapper::DBMapper, T::DataType)

Eliminates (when possible) the struct data from the DB

# Arguments
- `mapper::DBMapper`: The database mapper
- `T::Type{<:Model}`: Datatype of a registered model

"""
function drop_table!(mapper::DBMapper, T::DataType)
    check_valid_type(mapper, T)
    drop_table!(mapper, mapper.pool.dbtype, T)
end
function drop_table!(mapper::DBMapper, dbtype::DataType, T::Type{<:Model})
    drop_table!(mapper, database_kind(dbtype), T)
end

"""
function exists(mapper::DBMapper, T::Type{T}; kwargs...) where T <: Model

Return wether the element exists in the database

# Arguments
- `mapper::DBMapper`: The database mapper
- `T::DataType`: Datatype of a registered model we want to know their existence
- `kwargs`: fields we want to search for existence. The param pk cand be used as generic way
to identify the primary key of the struct

```
struct Author <: Model ... end
exists(mapper, Author, name="some name", age=30)
```
"""
function exists(mapper::DBMapper, ::Type{T}; kwargs...) where T <: Model
    check_valid_type(mapper, T)
    params = Dict(kwargs...)
    if haskey(params, :pk)
        params[idfield(mapper, T)] = params[:pk]
        pop!(params, :pk)
    end
    exists(mapper, mapper.pool.dbtype, T; params...)
end
"""
    function exists(mapper::DBMapper, elem::T) where T<:Model

Return wether the element exists in the database

# Arguments
- `mapper::DBMapper`: The database mapper
- `elem<:Model`: Element to determine its existence

```
struct Author <: Model ... end
exists(mapper, Author(name="some name", age=30))
```
"""
function exists(mapper::DBMapper, elem::T) where T<:Model
    check_valid_type(mapper, T)
    return exists(mapper, T; pk=getid(elem, mapper))
end
function exists(mapper::DBMapper, dbtype, T::Type{<:Model}; kwargs...)
    exists(mapper, database_kind(dbtype), T; kwargs...)
end

"""
    function delete!(mapper::DBMapper, T::Type{<:Model}; kwargs...)

Remove elements from the database

# Arguments
- `mapper::DBMapper`: The database mapper
- `T::DataType`: Datatype of a registered model we want to delete
- `kwargs`: fields we want to search for existence.

```
struct Author <: Model ... end
delete!(mapper, Author, name="some name", age=30)
```
"""
function Base.delete!(mapper::DBMapper, T::Type{<:Model}; kwargs...)
    check_valid_type(mapper, T)
    check_valid_type(mapper, T)
    params = Dict(kwargs...)
    if haskey(params, :pk)
        params[idfield(mapper, T)] = params[:pk]
        pop!(params, :pk)
    end
    return delete!(mapper, mapper.pool.dbtype, T; params...)
end
function Base.delete!(mapper::DBMapper, dbtype::DataType, T::Type{<:Model}; kwargs...)
    return delete!(mapper, database_kind(dbtype), T; kwargs...)
end
"""
    function delete!(mapper::DBMapper, T::Type{<:Model}; kwargs...)

Remove elements from the database

# Arguments
- `mapper::DBMapper`: The database mapper
- `elem<:Model`: Element to delete

The element needs a valid identifier.

```
struct Author <: Model ... end
delete!(mapper, Author(id=valid_id, name="some name", age=30))
```
"""
function Base.delete!(mapper::DBMapper, elem::T) where T<:Model
    check_valid_type(mapper, T)
    return delete!(mapper, T; pk=getid(elem, mapper))
end

"""
    select_all(mapper::DBMapper, T::Type{<:Model}; ; fields::Array{Symbol}=[], kwargs...)

Select all the elements that meet a criteria

# Arguments
- `mapper::DBMapper`: The database mapper
- `T::DataType`: Datatype of a registered model we want to search for
- `kwargs`: criteria we want to search

```
struct Author <: Model ... end
select_all(mapper, Author, age=30)
```
"""
function select_all(mapper::DBMapper, T::Type{<:Model};  fields::Array{Symbol}=Symbol[], kwargs...)
    check_valid_type(mapper, T)
    return select_all(mapper, mapper.pool.dbtype, T; fields=fields, kwargs...)
end
function select_all(mapper::DBMapper, dbtype::DataType, T::Type{<:Model};  fields::Array{Symbol}=Symbol[], kwargs...)
    return select_all(mapper, database_kind(dbtype), T; fields=fields, kwargs...)
end

"""
    database_kind(c::Type{T})

Return the kind of the database: Relational or NonRelational

The parameter should be the type of the connection of the concrete
database being used.
Every concrete database implementation must override this function.
"""
databasea_kind(c::Type{T}) where T = throw("Unknow database kind")

function configure_relation(mapper::DBMapper, T::Type, field_name::Symbol; on_delete::OnEventAction=DoNothing(), on_update::OnEventAction=DoNothing())
    if mapper.dirty == true
        analyze_relations(mapper)
    end
    table = mapper.tables[T]
    table.relations[table.fields[field_name]].on_delete = on_delete
    table.relations[table.fields[field_name]].on_update = on_update
end

"""
    element_type(x::Type{T}) where T

Return the element type of generic elements or the type itself for a regular type
Also it convert complex types to common one, ex: <:AbstractDict to Dict
"""
element_type(x::Type{T}) where T = x
element_type(x::DataType) = x
element_type(data::Type{<:AbstractNullable{T}}) where T = T
element_type(x::Type{<:AbstractDict}) = Dict

"""
    unmarshal(mapper::DBMapper, dbtype::DataType, dest::Type, orig)

Convert data obtained form the database to the julia representation.
It has three level of dispatch, the database specific, the database kind specific and
the type specific.
Return the orig element converted to the dest type. The dest type is usually obtained
from the struct fields.
"""
unmarshal(mapper::DBMapper, dbtype::DataType, dest::Type, orig) = unmarshal(mapper, database_kind(dbtype), dest, orig)
unmarshal(mapper::DBMapper, dbkind::Type{T}, dest::Type, orig) where T<:DatabaseKind =  unmarshal(mapper, dest, orig)
unmarshal(mapper::DBMapper, dest::Type, orig) = unmarshal(dest, orig)
unmarshal(dest::DataType, orig) = orig
unmarshal(dest::Type{DateTime}, orig::String)  = DateTime(orig)
unmarshal(D::Type{<:Dict}, d::String) = JSON.parse(d, dicttype=D)
unmarshal(mapper::DBMapper, x) = unmarshal(x)
unmarshal(mapper::DBMapper, ttype::Type{T}, x::String) where T = unmarshal(T, x)
unmarshal(x) = x
unmarshal(d::Type{T}, b::String) where T<:Number = parse(T, b)
unmarshal(d::Type{Integer}, b::String) = parse(Int64, b)
unmarshal(d::Type{String}, b::String) = b
unmarshal(::Type{DBId{T}}, x::String) where T<:Integer = parse(UInt64, x)
unmarshal(::Type{DBId{T}}, x::String) where T<:AbstractString = x

include(joinpath(@__DIR__, "Relational", "Relational.jl"))
include(joinpath(@__DIR__, "NonRelational", "NonRelational.jl"))

included_sources = []
function include_once(path::AbstractString)
    if path ∉ included_sources
        push!(included_sources, path)
        include(path)
    end
end

function __init__()
    @require SQLite="0aa819cd-b072-5ff4-a722-6bc24af294d9" begin
        include_once(joinpath(@__DIR__, "Relational", "SQLite.jl"))
    end
    @require LibPQ="194296ae-ab2e-5f79-8cd4-7183a0a5a0d1" begin
        include_once(joinpath(@__DIR__,  "Relational", "PostgreSQL.jl"))
    end
    @require Redis="0cf705f9-a9e2-50d1-a699-2b372a39b750" begin
        include_once(joinpath(@__DIR__, "NonRelational", "Redis.jl"))
    end

end

end # module

max_stars_repo_name: lucianolorenti/StructDatabaseMapping

max_issues_repo_name: lucianolorenti/StructDatabaseMapping

max_forks_repo_name: lucianolorenti/StructDatabaseMapping


text:
# issite methods for naieve site counting
# =======================================
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

"Test whether a nucleic acid in a sequence is ambiguous."
issite{T<:NucleicAcid}(::Type{Ambiguous}, a::T) = isambiguous(a)

"Test whether a nucleotide site of two aligned sequences has ambiguities."
@inline function issite{T<:NucleicAcid}(::Type{Ambiguous}, a::T, b::T)
    return issite(Ambiguous, a) | issite(Ambiguous, b)
end

"Test whether a nucleic acid in a sequence is a gap character."
issite{T<:NucleicAcid}(::Type{Gap}, a::T) = isgap(a)

"Test whether a nucleotide site of two aligned sequences has gap characters."
@inline function issite{T<:NucleicAcid}(::Type{Gap}, a::T, b::T)
    return issite(Gap, a) | issite(Gap, b)
end

"Test whether a nucleic acid in a sequence is certain."
issite{T<:NucleicAcid}(::Type{Certain}, a::T) = iscertain(a)

"Test whether a nucleotide site of two aligned sequences has two certain NucleicAcids."
@inline function issite{T<:NucleicAcid}(::Type{Certain}, a::T, b::T)
    return issite(Certain, a) & issite(Certain, b)
end

"Test whether a nucleotide site of two aligned sequences, constitutes a match."
issite{T<:NucleicAcid}(::Type{Match}, a::T, b::T) = a == b

"Test whether a nucleotide site of two aligned sequences, constitutes a mismatch."
issite{T<:NucleicAcid}(::Type{Mismatch}, a::T, b::T) = a != b

"Test whether a nucleotide site of two aligned seqences certainly constitutes a conserved site."
@inline function issite{T<:NucleicAcid}(::Type{Conserved}, a::T, b::T)
    return issite(Certain, a, b) & issite(Match, a, b)
end

"Test whether a nucleotide site of two aligned seqences certainly constitutes a mutation."
@inline function issite{T<:NucleicAcid}(::Type{Mutated}, a::T, b::T)
    return issite(Certain, a, b) & issite(Mismatch, a, b)
end

"Test whether a nucleotide site of two aligned sequences certainly constitutes a transition mutation."
@inline function issite{T<:NucleicAcid}(::Type{Transition}, a::T, b::T)
    return issite(Mutated, a, b) & ((ispurine(a) & ispurine(b)) | (ispyrimidine(a) & ispyrimidine(b)))
end

"Test whether a nucleotide site of two aligned sequences certainly constitutes a transversion mutation."
@inline function issite{T<:NucleicAcid}(::Type{Transversion}, a::T, b::T)
    return issite(Mutated, a, b) & ((ispurine(a) & ispyrimidine(b)) | (ispyrimidine(a) & ispurine(b)))
end


text:
import Base: ==, +, -, *, /, merge, merge!

#for StatsBase histogram
for op in (:+, :-, :*, :/)
    @eval ($op)(h1::Histogram, h2::Histogram) =
    (==)(h1.edges,h2.edges) ? Histogram(h1.edges, broadcast($op, h1.weights, h2.weights)) : throw(DimensionMismatch("The bins of the two histograms don't match"))
    @eval ($op)(h1::Histogram, n::Real) = Histogram(h1.edges, broadcast($op, h1.weights, n))
end

for T in (:Hist1D,:Hist2D)
    for op in (:+, :-)
        @eval function ($op)(h1::($T), h2::($T))
            h1.hist.edges != h2.hist.edges && throw(DimensionMismatch("Edges don't match"))
            _f(counts) = any(x -> x<0, counts)
            _hist = ($op)(h1.hist,  h2.hist)
            ($T)(_hist, h1.sumw2 + h2.sumw2)
        end
    end

    @eval function *(h1::($T), num::Real)
        _f(counts) = any(x -> x<0, counts)
        _f(h1.hist.weights) && error("Can't do * when bin count is negative")
        _hist = *(h1.hist, num)
        ($T)(_hist, h1.sumw2 * num^2)
    end

    # https://github.com/aminnj/yahist/blob/4a5767f181ec7fdcc4af18cf15ceedd1c2f89019/yahist/hist1d.py#L427-L430
    @eval function /(h1::($T), h2::($T))
        _f(counts) = any(x -> x<0, counts)
        (_f(h1.hist.weights) || _f(h2.hist.weights)) && error("Can't do / when bin count is negative")
        _hist = /(h1.hist, h2.hist)
        _sumw2 =  @. h1.sumw2 / (h2.hist.weights ^ 2) +
                (sqrt(h2.sumw2) * h1.hist.weights / (h2.hist.weights ^ 2)) ^ 2

        ($T)(_hist, _sumw2)
    end

    @eval function ==(h1::$T, h2::$T)
        h1.hist == h2.hist &&
        h1.sumw2 == h2.sumw2
    end

    @eval function merge!(h1::$T, h2::$T)
        h1.hist.edges != h2.hist.edges && throw(DimensionMismatch("The dimension doesn't match"))
        lock(h1)
        h1.hist.weights .+= h2.hist.weights
        h1.sumw2 .+= h2.sumw2
        unlock(h1)
        ($T)(h1.hist)
    end

    @eval merge(h1::$T, h2::$T) = merge!(deepcopy(h1), h2)
end

function merge(hists...)
    h = deepcopy(first(hists))
    length(hists) == 1 && return h
    foreach(x-> merge!(h, x), hists[2:end])
    h
end


text:
########################################################################################################
## Details
# This file implements the sampled-based estimation proceedure for UQ inference problems.
#       We make heavy use of AdvancedHMC.jl and Distributions.jl, many of the functions below are wrappers
#########################################################################################################
function Sample(pdist :: ArbitraryDistribution; num_chains = 1, estimate_map = false, return_chain = false, num_samples = 100, num_adapts = 300, verbose = false)

    if estimate_map == true
        MAP(pdist; verbose = verbose)
    end

    # Define a Hamiltonian system
    x0 = inverse(pdist.t, pdist.x)
    n = length(x0)
    @model function demo()
        x ~ pdist.prior
        Turing.@addlogprob! pdist.distribution(x)
    end

    model = demo()

    # # Start sampling.
    if num_chains == 1
        chain = Turing.sample(model, NUTS(num_adapts, 0.5; max_depth = 15), num_samples, init_theta = x0)
        samples = copy(group(chain, "x").value.data[:,:,1])
    else
        chain = Turing.sample(model, NUTS(num_adapts, 0.5, 15), MCMCThreads(),  num_samples, num_chains, init_theta = x0.*(1+0.1*randn()) )
        samples = copy( reshape(group(chain, "x").value.data, num_samples*num_chains, n) )
    end
    samples = [samples[i, :] for i = 1:size(samples, 1)]
    if verbose
        show(stdout, "text/plain", chain)
        println(" ")
    end

    # Extract Samples
    pdist.samples = transform.(pdist.t, samples)
    values = pdist.distribution.(samples)
    map = samples[argmax( values )]
    pdist.x = transform(pdist.t, map)
    if return_chain
        return chain, samples
    else
        return samples
    end

end


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

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing4.31"

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
end

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
function directorySearch(directory, searchString)
    outfile  = open(datadir("fileList.txt"),"w")
    fileList = filter(x->occursin(searchString, x), readdir(directory))

    for file in fileList
        println(outfile, file)
    end
    close(outfile)
end

# ╔═╡ 0b9bebfe-2525-11eb-3a90-fd98a5bcdace
directorySearch(pwd(),".jl")

# ╔═╡ 8f0d911c-2526-11eb-3c68-d38a7ede16cb
readlines(datadir("fileList.txt"), keep=true)

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing4.31"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═0b9bebfe-2525-11eb-3a90-fd98a5bcdace
# ╠═8f0d911c-2526-11eb-3c68-d38a7ede16cb
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
using Clapeyron, Test

const SOL = Clapeyron.Solvers

quadratic(x) = x*x - 4
rosenbrock(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
golden_number_fixpoint(x) = one(x) + one(x)/x
@testset "Solvers Module" begin
    @testset "ad_newton" begin
        x0 = 2+rand()
        fdf = SOL.f∂f(quadratic,x0)
        @test fdf[1] ≈ quadratic(x0)
        @test fdf[2] ≈ 2*x0
        @test SOL.ad_newton(quadratic,x0) ≈ 2.0
    end

    @testset "box optimize" begin
        #example of FMinBox in Optim.jl
        lower = [1.25, -2.1]
        upper = [Inf, Inf]
        solution = [1.25, 1.5625]
        initial_x = [2.0, 2.0]
        res = SOL.box_optimize(rosenbrock,initial_x,lower,upper)
        @test all(res.info.minimizer .≈ solution)
        @test res.info.minimum ≈ 0.0625
    end

    @testset "fixpoint" begin
        #example of FMinBox in Optim.jl
        x0 = 2+ rand()
        ϕ = (1+sqrt(5))/2
        @test SOL.fixpoint(golden_number_fixpoint,x0) ≈ ϕ
    end

    function f_diffmcp!(fvec, x)
        fvec[1] = (1-x[1])^2-1.01
    end

    # A difficult MCP.
    #
    # Presented in Miranda and Fackler (2002): "Applied Computational Economics and
    # Finance", p. 51
    #does not converge in NLSolve.jl converges here with NLSolvers.jl
    @testset "nlsolve" begin
        res = SOL.nlsolve(f_diffmcp!,[0.0])
        solution = SOL.x_sol(res)
        v = [1.0]
        f_diffmcp!(v,solution)
        @test v[1] <= sqrt(eps(Float64))
    end
end


text:
using Logging
using JuMP, Gurobi
using DecisionProgramming
using CSV, DataFrames, PrettyTables

# Setting subproblem specific parameters
const chosen_risk_level = 13

# Reading tests' technical performance data (dummy data in this case)
data = CSV.read("CHD_preventative_care_data.csv", DataFrame)

# Bayes posterior risk probabilities calculation function
# prior = prior risk level for which the posterior risk distribution is calculated for,
# t = test done
# returns a 100x1 vector with the probabilities of getting CHD given the prior risk level and test result
# for no test done (i.e. t = 3) returns a zero vector
function update_risk_distribution(prior::Int64, t::Int64)
    if t == 1 # the test is TRS
        # P(TRS = result | sick) = P(TRS_if_sick = result) * P(sick) = P(TRS_if_sick = result) * P(prior_risk)
        numerators = data.TRS_if_sick .* data.risk_levels[prior]

        # P(TRS = result) = P(TRS_if_sick = result) * P(sick) + P(TRS_if_healthy = result) * P(healthy)
        denominators = data.TRS_if_sick .* data.risk_levels[prior]  + data.TRS_if_healthy .* (1 - data.risk_levels[prior])

        posterior_risks = numerators./denominators

        # if the denominator is zero, post_risk is NaN, changing those to 0
        for i = 1:101
            if isnan(posterior_risks[i])
                posterior_risks[i] = 0
            end
        end

        return posterior_risks

    elseif t == 2 #the test is GRS
        numerators = (data.GRS_if_sick .* data.risk_levels[prior])
        denominators = data.GRS_if_sick .* data.risk_levels[prior]  + data.GRS_if_healthy .* (1 .- data.risk_levels[prior])

        posterior_risks =  numerators./denominators

        # if the denominator is zero, post_risk is NaN, changing those to 0
        for i = 1:101
            if isnan(posterior_risks[i])
                posterior_risks[i] = 0
            end
        end

        return posterior_risks

    else # no test performed
        risks_unchanged = zeros(100,1)

        return risks_unchanged

    end
end

# State probabilites calculation function
# risk_p = the resulting array from update_risk_distribution
# t = test done
# h = CHD or no CHD
# returns the probability distribution in 101x1 vector for the states of the R node given the prior risk level (must be same as to function update_risk_distribution), test t and health h
function state_probabilities(risk_p::Array{Float64}, t::Int64, h::Int64, prior::Int64)

    #if no test is performed, then the probabilities of moving to states (other than the prior risk level) are 0 and to the prior risk element is 1
    if t == 3
        state_probabilites = zeros(101, 1)
        state_probabilites[prior] = 1.0
        return state_probabilites
    end

    # return vector
    state_probabilites = zeros(101,1)

    # copying the probabilities of the scores for ease of readability
    if h == 1 && t == 1    # CHD and TRS
        p_scores = data.TRS_if_sick
    elseif t ==1    # no CHD and TRS
        p_scores = data.TRS_if_healthy
    elseif h == 1 && t == 2 # CHD and GRS
        p_scores = data.GRS_if_sick
    else # no CHD and GRS
        p_scores = data.GRS_if_healthy
    end

    for i = 1:101 #iterating through all risk levels 0%, 1%, ..., 99%, 100% in data.risk_levels
        for j = 1:101 #iterates through all risk estimates in risk_p
            #finding all risk estimates risk_p[j] within risk level i
            # risk_level[i] <= risk_p < risk_level[i]
            if i < 101 && data.risk_levels[i] <= risk_p[j] && risk_p[j] < data.risk_levels[i+1]
                state_probabilites[i] += p_scores[j]
            elseif i == 101 && data.risk_levels[i] <= risk_p[j] #special case: the highest risk level[101] = 100%
                state_probabilites[i] += p_scores[j]
            end
        end
    end

    return state_probabilites
end

function analysing_results(Z::DecisionStrategy, sprobs::StateProbabilities)

    d = Z.D[1] #taking one of the decision nodes to retrieve the information_set_R
    information_set_R = vec(collect(paths(S[d.I_j])))
    results = DataFrame(Information_set = map( x -> string(x) * "%", [0:1:100;]))
    # T1
    Z_j = Z.Z_j[1]
    probs =  map(x -> x > 0 ? 1 : 0, get(sprobs.probs, 1,0)) #these are zeros and ones
    dec = [Z_j(s_I) for s_I in information_set_R]
    results[!, "T1"] = map(x -> x == 0 ? "" : "$x", probs.*dec)

    # T2
    Z_j = Z.Z_j[2]
    probs = map(x -> x > 0 ? 1 : 0, (get(sprobs.probs, 4,0))) #these are zeros and ones
    dec = [Z_j(s_I) for s_I in information_set_R]
    results[!, "T2"] = map(x -> x == 0 ? "" : "$x", probs.*dec)

    # TD
    Z_j = Z.Z_j[3]
    probs = map(x -> x > 0 ? 1 : 0, (get(sprobs.probs, 6,0))) #these are zeros and ones
    dec = [Z_j(s_I) for s_I in information_set_R]
    results[!, "TD"] = map(x -> x == 0 ? "" : "$x", probs.*dec)

    pretty_table(results)
end

const R0 = 1
const H = 2
const T1 = 3
const R1 = 4
const T2 = 5
const R2 = 6
const TD = 7
const TC = 8
const HB = 9

const H_states = ["CHD", "no CHD"]
const T_states = ["TRS", "GRS", "no test"]
const TD_states = ["treatment", "no treatment"]
const R_states = map( x -> string(x) * "%", [0:1:100;])
const TC_states = ["TRS", "GRS", "TRS & GRS", "no tests"]
const HB_states = ["CHD & treatment", "CHD & no treatment", "no CHD & treatment", "no CHD & no treatment"]

@info("Creating the influence diagram.")
S = States([
    (length(R_states), [R0, R1, R2]),
    (length(H_states), [H]),
    (length(T_states), [T1, T2]),
    (length(TD_states), [TD])
])

C = Vector{ChanceNode}()
D = Vector{DecisionNode}()
V = Vector{ValueNode}()
X = Vector{Probabilities}()
Y = Vector{Consequences}()

I_R0 = Vector{Node}()
X_R0 = zeros(S[R0])
X_R0[chosen_risk_level] = 1
push!(C, ChanceNode(R0, I_R0))
push!(X, Probabilities(R0, X_R0))

I_H = [R0]
X_H = zeros(S[R0], S[H])
X_H[:, 1] = data.risk_levels     # 1 = "CHD"
X_H[:, 2] = 1 .- X_H[:, 1]  # 2 = "no CHD"
push!(C, ChanceNode(H, I_H))
push!(X, Probabilities(H, X_H))

I_T1 = [R0]
push!(D, DecisionNode(T1, I_T1))

I_R1 = [R0, H, T1]
X_R1 = zeros(S[I_R1]..., S[R1])
for s_R0 = 1:101, s_H = 1:2, s_T1 = 1:3
    X_R1[s_R0, s_H, s_T1, :] =  state_probabilities(update_risk_distribution(s_R0, s_T1), s_T1, s_H, s_R0)
end
push!(C, ChanceNode(R1, I_R1))
push!(X, Probabilities(R1, X_R1))

I_T2 = [R1]
push!(D, DecisionNode(T2, I_T2))

I_R2 = [H, R1, T2]
X_R2 = zeros(S[I_R2]..., S[R2])
for s_R1 = 1:101, s_H = 1:2, s_T2 = 1:3
    X_R2[s_H, s_R1, s_T2, :] =  state_probabilities(update_risk_distribution(s_R1, s_T2), s_T2, s_H, s_R1)
end
push!(C, ChanceNode(R2, I_R2))
push!(X, Probabilities(R2, X_R2))

I_TD = [R2]
push!(D, DecisionNode(TD, I_TD))

I_TC = [T1, T2]
Y_TC = zeros(S[I_TC]...)
cost_TRS = -0.0034645
cost_GRS = -0.004
cost_forbidden = 0     #the cost of forbidden test combinations is negligible
Y_TC[1 , 1] = cost_forbidden
Y_TC[1 , 2] = cost_TRS + cost_GRS
Y_TC[1, 3] = cost_TRS
Y_TC[2, 1] =  cost_GRS + cost_TRS
Y_TC[2, 2] = cost_forbidden
Y_TC[2, 3] = cost_GRS
Y_TC[3, 1] = cost_TRS
Y_TC[3, 2] = cost_GRS
Y_TC[3, 3] = 0
push!(V, ValueNode(TC, I_TC))
push!(Y, Consequences(TC, Y_TC))

I_HB = [H, TD]
Y_HB = zeros(S[I_HB]...)
Y_HB[1 , 1] = 6.89713671259061  # sick & treat
Y_HB[1 , 2] = 6.65436854256236  # sick & don't treat
Y_HB[2, 1] = 7.64528451705134   # healthy & treat
Y_HB[2, 2] =  7.70088349200034  # healthy & don't treat
push!(V, ValueNode(HB, I_HB))
push!(Y, Consequences(HB, Y_HB))

@info("Validate influence diagram.")
validate_influence_diagram(S, C, D, V)
sort!.((C, D, V, X, Y), by = x -> x.j)

P = DefaultPathProbability(C, X)
U = DefaultPathUtility(V, Y)

@info("Creating the decision model.")
model = Model()
z = DecisionVariables(model, S, D)

# Defining forbidden paths to include all those where a test is repeated twice
forbidden_tests = ForbiddenPath[([T1,T2], Set([(1,1),(2,2),(3,1), (3,2)]))]
scale_factor = 10000.0
x_s = PathCompatibilityVariables(model, z, S, P; fixed = Dict(1 => chosen_risk_level), forbidden_paths = forbidden_tests, probability_cut=false)

EV = expected_value(model, x_s, U, P, probability_scale_factor = scale_factor)
@objective(model, Max, EV)

@info("Starting the optimization process.")
optimizer = optimizer_with_attributes(
    () -> Gurobi.Optimizer(Gurobi.Env()),
    "MIPFocus" => 3,
    "MIPGap" => 1e-6,
)
set_optimizer(model, optimizer)
optimize!(model)

@info("Extracting results.")
Z = DecisionStrategy(z)

@info("Printing decision strategy using tailor made function:")
sprobs = StateProbabilities(S, P, Z)
analysing_results(Z, sprobs)

@info("Printing state probabilities:")
# Here we can see that the probability of having a CHD event is exactly that of the chosen risk level
print_state_probabilities(sprobs, [R0, R1, R2])

@info("Computing utility distribution.")
udist = UtilityDistribution(S, P, U, Z)

@info("Printing utility distribution.")
print_utility_distribution(udist)

@info("Printing statistics")
print_statistics(udist)


text:
filename = "math-assignment.txt"

file = open(filename)
lines = readlines(file)
close(file)

out = open("homework-delivery.txt", "w")

for l in lines
    s = strip(l, ['=', ' '])
    a, b = split(s, '+')
    x = parse(Int, a)
    y = parse(Int, b)
    println(out, x, " + ", y, " = ", x + 10)
end

close(out)


text:
function iter_debug_print(nit, Λ, res, c, r, spurious=1e-5)
    print(nit)
    print(":\t")
	in_eig = Λ[in_contour.(Λ, c, r)]
	in_res = res[in_contour.(Λ, c, r)]
    print(sum(in_contour.(Λ, c, r)))
    print(" (")
    print(sum(in_res .< spurious))
    print(")\t")
	if sum(in_contour.(Λ, c, r)) > 0
		print(maximum(res[in_contour.(Λ, c, r)]))
		in_res_conv = in_res[in_res .< spurious]
		if size(in_res_conv, 1) > 0
			print("\t(")
			print(maximum(in_res_conv))
			print(")")
		end
	end
    println()
end

function iter_debug_print(nit, Λ, res, contour::Contour, spurious=1e-5)
    print(nit)
    print(":\t")
	in_eig = Λ[in_contour(Λ, contour)]
	in_res = res[in_contour(Λ, contour)]
    print(sum(in_contour(Λ, contour)))
    print(" (")
    print(sum(in_res .< spurious))
    print(")\t")
	if sum(in_contour(Λ, contour)) > 0
		print(maximum(res[in_contour(Λ, contour)]))
		in_res_conv = in_res[in_res .< spurious]
		if size(in_res_conv, 1) > 0
			print("\t(")
			print(maximum(in_res_conv))
			print(")")
		end
	end
    println()
end

function convergence_info(Λ, X, residuals, contour::Contour, spurious=1e-3)
	in_ind = in_contour(Λ, contour)
    in_eig = Λ[in_ind] # eigenvalues inside contour
    in_res = residuals[in_ind] # residuals of eigenvalues inside contour
    print("Number of eigenvalues inside contour: ")
    println(size(in_eig, 1))
    if sum(in_ind) > 0
        in_res_conv = in_res[in_res .<= spurious]
        in_eig_conv = in_eig[in_res .<= spurious]
        print("Number inside converged : ")
        println(size(in_eig_conv, 1))
        print("Max res inside: ")
        println(maximum(in_res))
        if size(in_res_conv, 1) > 0
            print("Max res inside non spurious: ")
            println(maximum(in_res_conv))
        end
    end
end

function convergence_info(Λ, X, residuals, c, r, spurious=1e-3)
	contour = circular_contour_trapezoidal(c, r, 4)
	convergence_info(Λ, X, residuals, contour, spurious)
end

function beyn_svd_step!(Q₀::AbstractMatrix, Q₁::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, X::AbstractMatrix{ComplexF64}, Λ::Array)
    S = svd!(Q₀)
    mul!(A, S.U', Q₁)
    mul!(B, A, S.V)
    mul!(A, B, Diagonal(1 ./ S.S))
    F = eigen!(A)
    mul!(X, S.U, F.vectors)
    Λ .= F.values
end

function beyn_qr_step!(Q₀::AbstractMatrix, Q₁::AbstractMatrix, X::AbstractMatrix, Λ::Array)
    qt, rt = qr!(Q₀)
    qt = Matrix(qt)
    F = eigen!(qt' * Q₁ * inv(rt))
    mul!(X, qt, F.vectors)
    Λ .= F.values
end

function beyn_rr_step!(Q₀::AbstractMatrix, Q₁::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, X::AbstractMatrix{ComplexF64}, Λ::Array)
    mul!(A, X', Q₁)
    mul!(B, X', Q₀)
    F = eigen!(A, B)
    mul!(X, Q₀, F.vectors)
    Λ .= F.values
end

function beyn_rr_step2!(Q₀::AbstractMatrix, Q₁::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, X::AbstractMatrix{ComplexF64}, Λ::Array)
    mul!(A, Q₀', Q₁)
    mul!(B, Q₀', Q₀)
    F = eigen!(A, B)
    mul!(X, Q₀, F.vectors)
    Λ .= F.values
end

function update_R!(X::AbstractMatrix, R::AbstractMatrix, Λ::Array, T::Function)
    for i=1:size(X, 2)
        X[:,i] ./= norm(X[:,i])
        R[:,i] .= T(Λ[i]) * X[:,i]
    end
end

function update_R!(X::AbstractMatrix, R::AbstractMatrix, Λ::Array, A::AbstractMatrix, B=I)
    for i=1:size(X, 2)
        X[:,i] ./= norm(X[:,i])
        R[:,i] .= (A - Λ[i]*B) * X[:,i]
    end
end

function update_R_moments!(X::AbstractMatrix, R::AbstractMatrix, Λ::Array, res::Array, T::Function, c, r)
    for i=1:size(X, 2)
        X[:,i] ./= norm(X[:,i])
        R[:,i] .= T(Λ[i]) * X[:,i]
    end
	residuals!(res, R, Λ, T)
	# p = sortperm(res .- in_contour.(Λ, c, r))
	p = sortperm(res)
	# p = sortperm(abs.(Λ .- c))
	# a = copy(abs.(Λ .- c))
	# b = copy(res)
	# p = sortperm(normalize!(a) .+ normalize!(b))
	res .= res[p]
	X .= X[:, p]
	Λ .= Λ[p]
	R .= R[:, p]
end

function update_R_moments_all!(X::AbstractMatrix, R::AbstractMatrix, Λ::Array, res::Array, T::Function, c, r)
    for i=1:size(X, 2)
        X[:,i] ./= norm(X[:,i])
        R[:,i] .= T(Λ[i]) * X[:,i]
    end
	residuals!(res, R, Λ, T)
end

function normalize!(X::AbstractVecOrMat)
    for i=1:size(X, 2)
        X[:,i] ./= norm(X[:,i])
    end
	X
end

function residuals(R::AbstractMatrix, Λ::Array, T::Function)
    res = Array{Float64}(undef, size(Λ, 1))
    for i=1:size(Λ, 1)
        res[i] = norm(R[:,i])/norm(T(Λ[i]))
    end
    res
end

function residuals!(res::Array, R::AbstractMatrix, Λ::Array, T::Function)
    for i=1:size(Λ, 1)
        res[i] = norm(R[:,i])/norm(T(Λ[i]))
    end
    res
end

function residuals!(res::Array, R::AbstractMatrix, Λ::Array, A::AbstractMatrix)
    for i=1:size(Λ, 1)
        res[i] = norm(R[:,i])
    end
    res
end

finalize!(x::Any) = nothing

function linsolve!(Y, C, X, factorizer, left_divider)
    F = factorizer(C)
    left_divider(Y, F, X)
    finalize!(F)
end


text:
struct GaussianCopula
    cov::Array{Float64, 2}
    _chol::Array{Float64, 2}
    # GaussianCopula(cov::Array{Float64,2}) = new(cov, cholesky(cov, Val{:L}))
    # two argument constructor is allowed for easy copy-paste for debugging
end
GaussianCopula(cov::Array{Float64,2}) = GaussianCopula(cov, cholesky(cov).factors')
GaussianCopula(ndim::Int, rho::Float64) = GaussianCopula(rho*ones(ndim, ndim) + (1.0-rho)*diagm(ones(ndim)))

ndim(c::GaussianCopula) = size(c.cov,1)

normal_pdf(x) = exp.(-(x .^ 2) ./ 2.0) ./ sqrt(2*pi) # assumes zero mean, unit variance
normal_cdf(x) = 0.5 .* (1.0 .+ erf.(x ./ sqrt(2.0))) # assumes zero mean, unit variance

function sample_mvn(rng::AbstractRNG, cov_chol::Array{Float64})
    return cov_chol*randn(rng, size(cov_chol,1))
end

function rand(rng::AbstractRNG, c::GaussianCopula, sample::Vector{Float64}=Vector{Float64}(undef, ndim(c)))
    return normal_cdf(sample_mvn(rng, c._chol))
end
rand(c::GaussianCopula) = rand(Base.GLOBAL_RNG, c)


text:
# AUTO GENERATED FILE - DO NOT EDIT

export 'feffery'_antdradiogroup

"""
    'feffery'_antdradiogroup(;kwargs...)

An AntdRadioGroup component.

Keyword arguments:
- `id` (String; optional)
- `buttonStyle` (a value equal to: 'outline', 'solid'; optional)
- `className` (String; optional)
- `defaultValue` (String; optional)
- `direction` (a value equal to: 'horizontal', 'vertical'; optional)
- `disabled` (Bool; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `optionType` (a value equal to: 'default', 'button'; optional)
- `options` (optional): . options has the following type: Array of lists containing elements 'label', 'value', 'disabled'.
Those elements have the following types:
  - `label` (String; required)
  - `value` (String; required)
  - `disabled` (Bool; optional)s
- `persisted_props` (Array of a value equal to: 'value's; optional): Properties whose user interactions will persist after refreshing the
component or the page. Since only `value` is allowed this prop can
normally be ignored.
- `persistence` (Bool | String | Real; optional): Used to allow user interactions in this component to be persisted when
the component - or the page - is refreshed. If `persisted` is truthy and
hasn't changed from its previous value, a `value` that the user has
changed while using the app will keep that change, as long as
the new `value` also matches what was given originally.
Used in conjunction with `persistence_type`.
- `persistence_type` (a value equal to: 'local', 'session', 'memory'; optional): Where persisted user changes will be stored:
memory: only kept in memory, reset on page refresh.
local: window.localStorage, data is kept after the browser quit.
session: window.sessionStorage, data is cleared once the browser quit.
- `size` (a value equal to: 'large', 'middle', 'small'; optional)
- `style` (Dict; optional)
- `value` (String | Real; optional)
"""
function 'feffery'_antdradiogroup(; kwargs...)
        available_props = Symbol[:id, :buttonStyle, :className, :defaultValue, :direction, :disabled, :loading_state, :optionType, :options, :persisted_props, :persistence, :persistence_type, :size, :style, :value]
        wild_props = Symbol[]
        return Component("'feffery'_antdradiogroup", "AntdRadioGroup", "feffery_antd_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: ReilYoung/feffery-antd-components

max_issues_repo_name: ReilYoung/feffery-antd-components

max_forks_repo_name: ReilYoung/feffery-antd-components


text:
abstract type AbstractMockTester end

module Baz
    bar(::Int, ::AbstractString) = nothing
    struct BazType end
end

import Baz: bar, BazType

struct MyException <: Exception
    msg::AbstractString
end

@mock MockTester AbstractMockTester
t = MockTester()
@mockfunction(t,
    foo(::MockTester, ::Int, ::Int),
    bar(a::String, ::Int),
    baz(x::Vector{UInt8}))

facts("Test mock") do
    context("Expectations") do
        some_arg = 17

        @expect t foo(t, some_arg, 42)
        @expect t bar("Hello", -5) 42
        @expect t baz(b"Hello")

        @fact foo(t, some_arg, 42) --> nothing
        @fact bar("Hello", -5) --> 42
        @fact baz(b"Hello") --> nothing

        check(t)
    end

    context("Failed expectation") do
        @expect t foo(t, 17, 42)

        @fact_throws foo(t, 17, 43) MockExpectationException

        check(t)
    end

    context("Missing calls") do
        @expect t foo(t, 17, 42)

        @fact_throws check(t) MockExpectationException
        @fact t.calls --> isempty
    end

    context("Unexpected call") do
        @fact_throws foo(t, 17, 42)

        check(t)
    end

    context("Wrong number of arguments") do
        @expect t foo(t, 17, 42, 43)

        @fact_throws foo(t, 17, 42) MockExpectationException

        check(t)
    end

    context("Matchers") do
        mock_match(ValueMatcher(42), 42)
        @fact_throws mock_match(ValueMatcher(42), 43)
        @fact_throws mock_match(ValueMatcher(42), "42")

        mock_match(TypeMatcher(Int), 42)
        @fact_throws mock_match(TypeMatcher(Int), "")
    end

    context("Using TypeMatcher") do
        @expect t foo(t, TypeMatcher(Int), 42)
        @expect t foo(t, TypeMatcher(Number), 42)
        @expect t foo(t, TypeMatcher(AbstractString), 42)

        foo(t, 17, 42)
        foo(t, 18, 42)
        @fact_throws foo(t, 19, 42) MockExpectationException
        check(t)
    end

    context("MockThrowAction") do
        @expect t foo(t, 17, 42) Throws(MyException("some error"))

        @fact_throws foo(t, 17, 42) MyException

        check(t)
    end
end


text:
# Example of vOptGeneric, for a 2IP (here a 2LAP)
# July 2017

n  =  3

C1 = [  3    9   7   ;
       16   10   6   ;
        2    7  11     ]

C2 = [ 16   15   6   ;
        5    7  13   ;
        1    2  13     ]

using vOptGeneric
using GLPK ; using GLPKMathProgInterface

m = vModel( solver = GLPKSolverMIP() )
@variable( m , x[1:n,1:n] , Bin )
@addobjective( m , Min, sum( C1[i,j]*x[i,j] for i=1:n,j=1:n ) )
@addobjective( m , Min, sum( C2[i,j]*x[i,j] for i=1:n,j=1:n ) )
@constraint( m , cols[i=1:n], sum(x[i,j] for j=1:n) == 1 )
@constraint( m , rows[j=1:n], sum(x[i,j] for i=1:n) == 1 )
@time solve( m , method = :epsilon , step = 0.5 )

print_X_E(m)
getY_N(m)

max_stars_repo_name: faithyyyy/vOptSolver

max_issues_repo_name: faithyyyy/vOptSolver

max_forks_repo_name: faithyyyy/vOptSolver


text:
using NMF: nnmf, solve!, CoordinateDescent, nndsvd

"""
    NMF(;ncomps=nothing) <: LinearAlgorithm
    NMF(ncomps)

Use non-negative matrix factorization (NMF, NNMF) to form a non-negative low-rank orthonormal basis of the input. The implementation of the underlying NMF is provided by [NMF.jl](https://github.com/JuliaStats/NMF.jl). The implementation uses a non-negative SVD for initialization and a coordinate-descent solver to fit.

If `ncomps` is `nothing`, it will be set to the number of frames in the reference cube when processed.

# References
* [Ren et al. 2018](http://adsabs.harvard.edu/abs/2018ApJ...852..104R) Non-negative Matrix Factorization: Robust Extraction of Extended Structures
"""
@with_kw struct NMF <: LinearAlgorithm
    ncomps::Union{Int,Nothing} = nothing
end

function ADI.decompose(alg::NMF, cube, angles, cube_ref; kwargs...)
    X = flatten(cube)
    X_ref = flatten(cube_ref)

    k = isnothing(alg.ncomps) ? size(cube, 1) : alg.ncomps
    k > size(cube, 1) && error("ncomps ($k) cannot be greater than the number of frames ($(size(cube, 1)))")
    # manually initialize W and H
    _, H = nndsvd(X_ref, k)
    W = X * H'
    solve!(CoordinateDescent{eltype(X)}(maxiter=100, α=0), X, W, H)
    return H, W
end

function ADI.decompose(alg::NMF, cube, angles; kwargs...)
    # X = clamp.(flatten(cube), 0, Inf) # clamp non-negative values
    X = flatten(cube)

    k = isnothing(alg.ncomps) ? size(cube, 1) : alg.ncomps
    k > size(cube, 1) && error("ncomps ($k) cannot be greater than the number of frames ($(size(cube, 1)))")

    W, H = nndsvd(X, k)
    solve!(CoordinateDescent{eltype(X)}(maxiter=100, α=0), X, W, H)
    return H, W
end


text:
# TODO: more tests!

using ChainRules, Test, FDM, LinearAlgebra, Random
using ChainRules: extern, accumulate, accumulate!, store!, @scalar_rule,
    Wirtinger, wirtinger_primal, wirtinger_conjugate, add_wirtinger, mul_wirtinger,
    Zero, add_zero, mul_zero, One, add_one, mul_one, Casted, cast, add_casted, mul_casted,
    DNE, Thunk, Casted
using Base.Broadcast: broadcastable

include("test_util.jl")

@testset "ChainRules" begin
    include("differentials.jl")
    include("rules.jl")
    @testset "rules" begin
        include(joinpath("rules", "base.jl"))
        include(joinpath("rules", "array.jl"))
        @testset "linalg" begin
            include(joinpath("rules", "linalg", "dense.jl"))
            include(joinpath("rules", "linalg", "diagonal.jl"))
            include(joinpath("rules", "linalg", "symmetric.jl"))
            include(joinpath("rules", "linalg", "factorization.jl"))
        end
        include(joinpath("rules", "broadcast.jl"))
        include(joinpath("rules", "blas.jl"))
        include(joinpath("rules", "nanmath.jl"))
        include(joinpath("rules", "specialfunctions.jl"))
    end
end


text:
# Autogenerated wrapper script for Luna_jll for x86_64-apple-darwin
export libluna, luna

using FFTW_jll
using Zlib_jll
JLLWrappers.@generate_wrapper_header("Luna")
JLLWrappers.@declare_library_product(libluna, "@rpath/libluna.dylib")
JLLWrappers.@declare_executable_product(luna)
function __init__()
    JLLWrappers.@generate_init_header(FFTW_jll, Zlib_jll)
    JLLWrappers.@init_library_product(
        libluna,
        "lib/libluna.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        luna,
        "bin/luna",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# csci_times.jl
#
#	Parse CampbellScientific timestamps
#
# Jeremy Rüffer
# Thünen Institut
# Institut für Agrarklimaschutz
# Junior Research Group NITROSPHERE
# Julia 1.6.0
# Created: 04.11.2013
# Last Edit: 21.04.2021

"""# csci_times

Retrieve the minimum and maximum dates from the listed files

---

`minimum_dates, maximum_dates = csci_times{T<:String}(F::Array{T,1};headerlines::Int=4)
* **F**::Array{String,1} = Array of files to retrieve dates from
* **headerlines**::Int (optional) = Number of header lines, 4 is default
"""
function csci_times(F::Array{T,1};headerlines::Int=4) where T <: String
	# Get the first and last dates of each file
	mintimes = fill!(Array{DateTime}(undef,length(F)),DateTime(0))
	maxtimes = fill!(Array{DateTime}(undef,length(F)),DateTime(0))

	df = Dates.DateFormat("\"yyyy-mm-dd HH:MM:SS\"") # Date format
	df2 = Dates.DateFormat("\"yyyy-mm-dd HH:MM:SS.ss\"") # Date format

	datastart_pos = 1
	for i=1:1:length(F)
		# Find Minimum Time
		fsize = stat(F[i]).size
		fid = open(F[i],"r")
		for j=1:1:headerlines
			readline(fid)
		end
		datastart_pos = position(fid) # Find start of data
		if ~eof(fid)
			l = readline(fid)
			if length(l[1:findfirst(',',l)-1]) == 21
				mintimes[i] = DateTime(l[1:findfirst(',',l)-1],df)
			else
				mintimes[i] = DateTime(l[1:findfirst(',',l)-1],df2)
			end
		end

		# Find Maximum Time
		l = ""
		endpos = 1
		f = false
		while !f
			seek(fid,fsize-endpos)
			l = readline(fid)
			if !isempty(l)
				if length(l) > 21
					if l[6] == '-' && l[9] == '-' && l[12] == ' '
						f = true # Found the last line
					end
				end
			end
			endpos += 1
		end
		close(fid)
		if length(l[1:findfirst(',',l)-1]) == 21
			maxtimes[i] = DateTime(l[1:findfirst(',',l)-1],df)
		else
			maxtimes[i] = DateTime(l[1:findfirst(',',l)-1],df2)
		end
	end

	return mintimes, maxtimes
end

max_stars_repo_name: JeremyRueffer/ClimateDataIO

max_issues_repo_name: JeremyRueffer/ClimateDataIO


text:
struct BarotropicModel{M} <: AbstractOceanModel
    baroclinic::M
    function BarotropicModel(baroclinic::M) where {M}
        return new{M}(baroclinic)
    end
end

function vars_state(m::BarotropicModel, ::Prognostic, T)
    @vars begin
        U::SVector{2, T}
        η::T
    end
end

function init_state_prognostic!(
    m::BarotropicModel,
    Q::Vars,
    A::Vars,
    localgeo,
    t,
)
    Q.U = @SVector [-0, -0]
    Q.η = -0
    return nothing
end

function vars_state(m::BarotropicModel, ::Auxiliary, T)
    @vars begin
        Gᵁ::SVector{2, T}  # integral of baroclinic tendency
        U_c::SVector{2, T} # cumulate U value over fast time-steps
        η_c::T             # cumulate η value over fast time-steps
        U_s::SVector{2, T} # starting U field value
        η_s::T             # starting η field value
        Δu::SVector{2, T}  # reconciliation adjustment to u, Δu = 1/H * (U_averaged - ∫u)
        η_diag::T          # η from baroclinic model (for diagnostic)
        Δη::T              # diagnostic difference: η_barotropic - η_baroclinic
        y::T               # y-coordinate of grid
    end
end

function init_state_auxiliary!(
    m::BarotropicModel,
    state_aux::MPIStateArray,
    grid,
    direction,
)
    init_state_auxiliary!(
        m,
        (m, A, tmp, geom) -> ocean_init_aux!(m, m.baroclinic.problem, A, geom),
        state_aux,
        grid,
        direction,
    )
end

function vars_state(m::BarotropicModel, ::Gradient, T)
    @vars begin
        U::SVector{2, T}
    end
end

@inline function compute_gradient_argument!(
    m::BarotropicModel,
    G::Vars,
    Q::Vars,
    A,
    t,
)
    G.U = Q.U
    return nothing
end

function vars_state(m::BarotropicModel, ::GradientFlux, T)
    @vars begin
        ν∇U::SMatrix{3, 2, T, 6}
    end
end

@inline function compute_gradient_flux!(
    m::BarotropicModel,
    D::Vars,
    G::Grad,
    Q::Vars,
    A::Vars,
    t,
)
    ν = viscosity_tensor(m)
    D.ν∇U = -ν * G.U

    return nothing
end

@inline function viscosity_tensor(bm::BarotropicModel)
    m = bm.baroclinic
    return Diagonal(@SVector [m.νʰ, m.νʰ, 0])
end

vars_state(m::BarotropicModel, ::UpwardIntegrals, T) = @vars()
vars_state(m::BarotropicModel, ::DownwardIntegrals, T) = @vars()

@inline function flux_first_order!(
    m::BarotropicModel,
    F::Grad,
    Q::Vars,
    A::Vars,
    t::Real,
    direction,
)
    @inbounds begin
        U = @SVector [Q.U[1], Q.U[2], 0]
        η = Q.η
        H = m.baroclinic.problem.H
        g = m.baroclinic.grav
        Iʰ = @SMatrix [
            1 0
            0 1
            0 0
        ]

        F.η += U
        F.U += g * H * η * Iʰ
    end
end

@inline function flux_second_order!(
    m::BarotropicModel,
    F::Grad,
    Q::Vars,
    D::Vars,
    HD::Vars,
    A::Vars,
    t::Real,
)
    # numerical diffusivity for stability
    F.U += D.ν∇U

    return nothing
end

@inline function source!(
    m::BarotropicModel,
    S::Vars,
    Q::Vars,
    D::Vars,
    A::Vars,
    t::Real,
    direction,
)
    @inbounds begin
        U = Q.U

        # f × u
        f = coriolis_force(m.baroclinic, A.y)
        S.U -= @SVector [-f * U[2], f * U[1]]

        # vertically integrated baroclinic model tendency
        S.U += A.Gᵁ
    end
end

@inline wavespeed(m::BarotropicModel, n⁻, _...) =
    abs(SVector(m.baroclinic.cʰ, m.baroclinic.cʰ, m.baroclinic.cᶻ)' * n⁻)

# We want not have jump penalties on η (since not a flux variable)
function update_penalty!(
    ::RusanovNumericalFlux,
    ::BarotropicModel,
    n⁻,
    λ,
    ΔQ::Vars,
    Q⁻,
    A⁻,
    Q⁺,
    A⁺,
    t,
)
    ΔQ.η = -0

    return nothing
end

boundary_conditions(bm::BarotropicModel) =
    (bm.baroclinic.problem.boundary_conditions[1],)

"""
    boundary_state!(nf, bc, ::BarotropicModel, args...)

applies boundary conditions for this model
dispatches to a function in OceanBoundaryConditions.jl based on BC type defined by a problem such as SimpleBoxProblem.jl
"""
@inline function boundary_state!(nf, bc, bm::BarotropicModel, args...)
    return ocean_model_boundary!(bm, bc, nf, args...)
end


text:
module ValidateRCall_ape
using Test
using Phylo
using RCall
using DataFrames

global skipR = !rcopy(R"require(ape)")

# Run tree comparisons on increasing numbers of tips
@testset "RCall - testing Phylo vs ape" begin
    @testset "For $TreeType" for TreeType in
        (skipR ? [] : [NamedTree, NamedBinaryTree,
                       BinaryTree{ManyRoots, DataFrame, Vector{Float64}},
                       RootedTree, ManyRootTree])

        @testset "Testing with R rtree($i)" for i in 10:10:50
            rt = rcall(:rtree, i)
            jt = rcopy(TreeType, rt)
            jl = Set(getleafnames(jt))
            rl = Set(rcopy(rcall(Symbol("[["), rt, "tip.label")))
            @test jl == rl
            rt2 = RObject(jt)
            @test rcopy(rcall(Symbol("all.equal"), rt, rt2))
            @rput jt
            reval("jt2=jt");
            @rget jt2
            jt3=jt2;
            @rput jt3
            @rput rt
            @test rcopy(rcall(Symbol("all.equal"), rt, jt3))
        end

        @testset "Testing with julia rand(Nonultrametric($i))" for i in 10:10:50
            jt = rand(Nonultrametric{TreeType}(i))
            rt = RObject(jt)
            @test getleafnames(jt) == rcopy(rcall(Symbol("[["), rt, "tip.label"))
            jt2 = rcopy(TreeType, rt)
            @test Set(getleafnames(jt)) == Set(getleafnames(jt2))
            @test rcopy(rcall(Symbol("all.equal"), rt, RObject(jt2)))
        end

        @testset "Testing with julia rand(Ultrametric($i))" for i in 10:10:50
            jt = rand(Ultrametric{TreeType}(i))
            rt = RObject(jt)
            @test getleafnames(jt) == rcopy(rcall(Symbol("[["), rt, "tip.label"))
            jt2 = rcopy(TreeType, rt)
            @test Set(getleafnames(jt)) == Set(getleafnames(jt2))
            @test rcopy(rcall(Symbol("all.equal"), rt, RObject(jt2)))
        end

        @testset "Testing with julia rand(Ultrametric($i), [...])" for i in 10:10:50
            jt = rand(Ultrametric{TreeType}(i), ["one", "two"])
            rt = RObject(jt["one"])
            @test Set(getleafnames(jt)) ==
                Set(rcopy(rcall(Symbol("[["), rt, "tip.label")))
            jt2 = rcopy(TreeType, rt)
            @test Set(getleafnames(jt)) == Set(getleafnames(jt2))

            @test rcopy(rcall(Symbol("all.equal"), rt, RObject(jt2)))
        end

        @testset "Testing reading in newick trees from disk" begin
            jtree = open(parsenewick, Phylo.path("H1N1.newick"))
            rtree = rcall(Symbol("read.tree"), Phylo.path("H1N1.newick"))
            @test rcopy(rcall(Symbol("all.equal"), jtree, rtree))
        end

        @testset "Testing reading in nexus trees from disk" begin
            jts = open(parsenexus, Phylo.path("H1N1.trees"))
            rtree1 = R"read.nexus($(Phylo.path(\"H1N1.trees\")))$TREE1"
            jtree1 = jts["TREE1"]
            @test rcopy(rcall(Symbol("all.equal"), jtree1, rtree1))
            @test "H1N1_A_MIYAGI_3_2000" ∈ nodenameiter(jtree1)
            @test collect(keys(gettreeinfo(jts)["TREE1"])) == ["lnP"]
        end
    end
end

end


text:
"""
Abstract population. The top node of the Population type hierarchy.
"""
abstract type AbstractPopulation end

"""
Abstract network population. Interactors are nodes of some network.
"""
abstract type AbstractNetworkPopulation <: AbstractPopulation end

"""
Abstract spatial population. Represents any type of population in which
interactors have a spatial representation (e.g. a point on a plane).
"""
abstract type AbstractSpatialPopulation <: AbstractPopulation end

"""
Abstract lattice population: interactors on some lattice.
"""
abstract type AbstractLatticePopulation <: AbstractSpatialPopulation end