text:
# Saving the LogL landscape data for S4 Fig
using CSV
using SurNoR_2020
using DataFrames
using Query

# A function to read the parameters' name and their optimization range given
# their index
function Func_if_VarName(num_param)
    if num_param==1
        # P_change              # The name in the codes
        VarName = "m/(1+m)"     # The name in the paper
        x = 0.00:0.0025:0.3
    elseif num_param==2
        # Eps
        VarName = "eps"
        x = -12:0.2:0.2
        x = 10 .^ x
    elseif num_param==3
        # β_R
        VarName = "beta_1"
        x = 0.1:0.2:15
    elseif num_param==4
        # β_RN1
        VarName = "beta_N1"
        x = 0.00:0.005:0.6
    elseif num_param==5
        # β_RN2
        VarName = "beta_N2"
        x = 0.00:0.005:0.6
    elseif num_param==6
        # γ_R
        VarName = "λ_R"
        x = 0.85:0.01:0.99
    elseif num_param==7
        # γ_N
        VarName = "λ_N"
        x = 0.20:0.05:0.99
    elseif num_param==8
        # λ_R
        VarName = "μ_R"
        x = 0.85:0.01:1
    elseif num_param==9
        # λ_N
        VarName = "μ_N"
        x = 0.20:0.05:1
    elseif num_param==10
        # α
        VarName = "ρ_b"
        x = 0.0:0.02:0.2
    elseif num_param==11
        # Δα
        VarName = "delta ρ"
        x = 0.0:0.05:(1-Param1.α)
    elseif num_param==12
        # Q_N_0
        VarName = "Q_N0"
        x = 0.0:0.25:12
    elseif num_param==13
        # N_PS
        VarName = "T_PS"
        x = cat(0:14,15:10:60,dims=1)
    elseif num_param==14
        # w_2
        VarName = "w_0"
        x = 0:0.05:1
    elseif num_param==15
        # w_11
        VarName = "w_11"
        x = 0:0.05:1
    elseif num_param==16
        # w_12
        VarName = "w_12"
        x = 0:0.05:1
    elseif num_param==17
        # w_scaler
        VarName = "w_scale"
        x = 0:0.1:6
    elseif num_param==18
        # β_R_2
        VarName = "beta_2"
        x = 0.1:0.2:15
    end
    return (VarName,x)
end

# ------------------------------------------------------------------------------
# Reading Data and making input
BehavData = CSV.read("./data/BehavData_SwitchState.csv")

# --------------------------------------------------------------------------
# Defining Hyper Parameters
update_form = "VarSMiLe"

if update_form == "SMiLe"
        surprise_form="CC"
else
        surprise_form="BF"
end

novelty_form="log-count"
update_form_model_free = "Q_learning"

# ----------------------------------------------------------------------
# Defining Hyper Parameters
HyperParam = SurNoR_2020.Str_HyperParam(; gamma_form="satur",
                                surprise_form=surprise_form,
                                novelty_form=novelty_form,
                                update_form=update_form,
                                update_form_model_free=update_form_model_free)

# ----------------------------------------------------------------------
# Subject Information
Block_Set = 1:2
Sub_Set = [1,2,3,4,5,6,7,8,9,11,12,13]
Epi_Set = 1:5

# ------------------------------------------------------------------------------
# Loading Parameters
Path_Params = "data/Fitted_parameters/Overall/OverAll_SurNoR_BestParams_diff5.csv"
Param_Data = CSV.read(Path_Params)
Row = 3

# ------------------------------------------------------------------------------
# Agents
Param1 = SurNoR_2020.Str_Param(; num_act=4, num_obs=11,
                                 count_decay=1, Q_R_0 = 0)

Param1.p_change = Param_Data[Row,1+1]
Param1.ϵ = Param_Data[Row,1+2]
Param1.β_R = Param_Data[Row,1+3]
Param1.β_ratio_N2R = Param_Data[Row,1+4]
Param1.γ_RL_R = Param_Data[Row,1+6]
Param1.γ_RL_N = Param_Data[Row,1+7]
Param1.λ_elig_R = Param_Data[Row,1+8]
Param1.λ_elig_N = Param_Data[Row,1+9]
Param1.α = Param_Data[Row,1+10]
Param1.Δα = Param_Data[Row,1+11]
Param1.Q_N_0 = Param_Data[Row,1+12]
Param1.N_prior_sweep = Param_Data[Row,1+13]
Param1.w_MF_to_MB_2 = Param_Data[Row,1+14]
Param1.w_MF_to_MB_1 = Param_Data[Row,1+15]
Param1.w_MF_to_MB_Scaler = Param_Data[Row,1+17]

Param2 = deepcopy(Param1)
Param2.β_ratio_N2R = Param_Data[Row,1+5]
Param2.w_MF_to_MB_1 = Param_Data[Row,1+16]
Param2.β_R = Param_Data[Row,1+18]

Param = [Param1, Param2]

Param_perm_order = 1:18
for i = 11
    num_param_i = Param_perm_order[i]

    Param_temp = zeros(18,)
    for i_dummy = 1:18
        Param_temp[i_dummy] = Param_Data[Row, 1 .+ i_dummy]
    end

    (VarName_i, x_i) = Func_if_VarName(num_param_i)
    @show VarName_i

    LogL_range = zeros(length(x_i))
    for x_ind_i = 1:length(x_i)
        Param_temp[num_param_i] = x_i[x_ind_i]

        Param1.p_change = Param_temp[1]
        Param1.ϵ = Param_temp[2]
        Param1.β_R = Param_temp[3]
        Param1.β_ratio_N2R = Param_temp[4]
        Param1.γ_RL_R = Param_temp[6]
        Param1.γ_RL_N = Param_temp[7]
        Param1.λ_elig_R = Param_temp[8]
        Param1.λ_elig_N = Param_temp[9]
        Param1.α = Param_temp[10]
        Param1.Δα = Param_temp[11]
        Param1.Q_N_0 = Param_temp[12]
        Param1.N_prior_sweep = Param_temp[13]
        Param1.w_MF_to_MB_2 = Param_temp[14]
        Param1.w_MF_to_MB_1 = Param_temp[15]
        Param1.w_MF_to_MB_Scaler = Param_temp[17]

        Param2 = deepcopy(Param1)
        Param2.β_ratio_N2R = Param_temp[5]
        Param2.w_MF_to_MB_1 = Param_temp[16]
        Param2.β_R = Param_temp[18]

        Param = [Param1, Param2]

        (Agents, LogL) = SurNoR_2020.Func_train_set_of_Agents(BehavData;
                                                               Sub_Set = Sub_Set,
                                                               Epi_Set=Epi_Set,
                                                               Block_Set=Block_Set,
                                                               Param=deepcopy(Param),
                                                               HyperParam=HyperParam)
        LogL_range[x_ind_i] = LogL
    end

    A = DataFrame(x=Array(x_i),y=LogL_range)
    Path_Save = "data/LogL_Landscape/"
    CSV.write(string(Path_Save, "LogL_range_i", string(num_param_i),"_Row", string(Row), ".CSV" ),A)
end


text:
#### Pressure model kernels

"""
    perturbation_pressure(
        m::AtmosModel{FT},
        press::PressureModel,
        state::Vars,
        diffusive::Vars,
        aux::Vars,
        t::Real,
        env,
        i,
    ) where {FT}
Returns the value of perturbation pressure gradient
for updraft i following He et al. (JAMES, 2020), given:
 - `m`, an `AtmosModel`
 - `press`, a `PressureModel`
 - `state`, state variables
 - `diffusive`, additional variables
 - `aux`, auxiliary variables
 - `t`, the time
 - `env`, NamedTuple of environment variables
 - `i`, index of the updraft
"""
function perturbation_pressure(
    m::AtmosModel{FT},
    press::PressureModel,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    env,
    i,
) where {FT}

    # Alias convention:
    gm = state
    en = state.turbconv.environment
    up = state.turbconv.updraft
    up_aux = aux.turbconv.updraft
    up_dif = diffusive.turbconv.updraft

    N_up = n_updrafts(m.turbconv)
    w_up_i = up[i].ρaw / up[i].ρa

    nh_press_buoy = press.α_b * up_aux[i].buoyancy
    nh_pressure_adv = -press.α_a * w_up_i * up_dif[i].∇w[3]
    nh_pressure_drag =
        press.α_d * (w_up_i - env.w) * abs(w_up_i - env.w) / press.H_up

    dpdz = nh_press_buoy + nh_pressure_adv + nh_pressure_drag

    return dpdz
end;

ext: jl

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
module NLSolvers

import Base: show

# We use often use the LinearAlgebra functions dot and norm for operations rela-
# ted to assessing angles between vectors, size of vectors and so on.
using LinearAlgebra: dot, I, norm,
                     mul!, rmul!,
                     cholesky,
                     UniformScaling, Symmetric,
                     ldiv!, factorize,
                     eigen, Hermitian,
                     diag, # for trust region diagonal manipulation
                     opnorm # for NWI safe guards

using RandomNumbers # for better random number generators and also rand!

using StaticArrays

include("Manifolds.jl")
include("objectives.jl")
export NonDiffed, OnceDiffed, TwiceDiffed

# make this struct that has scheme and approx
abstract type QuasiNewton{T1} end

abstract type HessianApproximation end
struct Inverse <: HessianApproximation end
struct Direct <: HessianApproximation end
export Inverse, Direct

struct OptOptions{T1, T2}
    g_tol::T1
    maxiter::T2
    show_trace::Bool
end
export OptOptions

OptOptions(; g_tol=1e-8, maxiter=10000, show_trace=false) =
OptOptions(g_tol, maxiter, show_trace)

# Globalization strategies
abstract type LineSearch end
include("globalization/linesearches/root.jl")
export backtracking, Backtracking, TwoPointQuadratic
include("globalization/trs_solvers/root.jl")
export NWI, Dogleg#, TRSolver

# Quasi-Newton (including Newton and gradient descent) functionality
include("quasinewton/quasinewton.jl")
export DBFGS, BFGS, SR1, DFP, GradientDescent, Newton

# Include the actual functions that expose the functionality in this package.
include("optimize/linesearch/linesearch.jl")
include("optimize/randomsearch/randomsearch.jl")
include("optimize/directsearch/directsearch.jl")
export NelderMead
include("optimize/trustregions/trustregions.jl")
export minimize, minimize!, OptProblem

include("nlsolve/root.jl")
export nlsolve, nlsolve!, NLEProblem
export ResidualKrylov, ResidualKrylovProblem
# Forcing Terms
export FixedForceTerm, DemboSteihaug, EisenstatWalkerA, EisenstatWalkerB
end # module


text:
@testset "single_knockout" begin
    m = StandardModel()
    add_metabolite!(m, Metabolite("A"))
    add_metabolite!(m, Metabolite("B"))

    add_gene!(m, Gene("g1"))
    add_gene!(m, Gene("g2"))
    add_reaction!(
        m,
        Reaction("v1", metabolites = Dict("A" => -1.0, "B" => 1.0), grr = [["g1"]]),
    )
    add_reaction!(
        m,
        Reaction("v2", metabolites = Dict("A" => -1.0, "B" => 1.0), grr = [["g1", "g2"]]),
    )
    add_reaction!(
        m,
        Reaction("v3", metabolites = Dict("A" => -1.0, "B" => 1.0), grr = [["g1"], ["g2"]]),
    )
    add_reaction!(
        m,
        Reaction(
            "v4",
            metabolites = Dict("A" => -1.0, "B" => 1.0),
            grr = [["g1", "g2"], ["g2"]],
        ),
    )

    opt_model = make_optimization_model(m, Tulip.Optimizer)
    knockout("g1")(m, opt_model)

    # Knockout should remove v1
    @test normalized_rhs(opt_model[:lbs][1]) == 0
    @test normalized_rhs(opt_model[:ubs][1]) == 0

    # Knockout should remove [g1, g2] (AND) and thus remove reaction
    @test normalized_rhs(opt_model[:lbs][2]) == 0
    @test normalized_rhs(opt_model[:ubs][2]) == 0

    # Knockout should remove [g1], but keep reaction (OR)
    @test normalized_rhs(opt_model[:lbs][3]) == 1000
    @test normalized_rhs(opt_model[:ubs][3]) == 1000

    # Knockout should remove [g1, g2] (AND), but keep reaction (OR)
    @test normalized_rhs(opt_model[:lbs][4]) == 1000
    @test normalized_rhs(opt_model[:ubs][4]) == 1000
end

@testset "multiple_knockouts" begin
    m = StandardModel()
    add_metabolite!(m, Metabolite("A"))
    add_metabolite!(m, Metabolite("B"))
    add_gene!(m, Gene("g1"))
    add_gene!(m, Gene("g2"))
    add_gene!(m, Gene("g3"))
    add_reaction!(
        m,
        Reaction("v1", metabolites = Dict("A" => -1.0, "B" => 1.0), grr = [["g1"], ["g3"]]),
    )
    add_reaction!(
        m,
        Reaction(
            "v2",
            metabolites = Dict("A" => -1.0, "B" => 1.0),
            grr = [["g1", "g2"], ["g3"]],
        ),
    )
    add_reaction!(
        m,
        Reaction(
            "v3",
            metabolites = Dict("A" => -1.0, "B" => 1.0),
            grr = [["g1"], ["g2"], ["g3"]],
        ),
    )

    opt_model = make_optimization_model(m, Tulip.Optimizer)
    knockout(["g1", "g3"])(m, opt_model)

    # Reaction 1 should be knocked out, because both
    # gene1 and gene 3 are knocked out
    @test normalized_rhs(opt_model[:lbs][1]) == 0
    @test normalized_rhs(opt_model[:ubs][1]) == 0

    # Reaction 2 should be knocked out, because both
    # [g1, g2] is an AND relationship
    @test normalized_rhs(opt_model[:lbs][1]) == 0
    @test normalized_rhs(opt_model[:ubs][1]) == 0

    # Reaction 3 should stay, because gene2 is still
    # available (the arrays have an OR relationship)
    @test normalized_rhs(opt_model[:lbs][3]) == 1000
    @test normalized_rhs(opt_model[:ubs][3]) == 1000
end

@testset "Knockouts on realistic models" begin
    for model in [
        load_model(StandardModel, model_paths["e_coli_core.json"]), #test on standardModel
        load_model(model_paths["e_coli_core.json"]), #then on JSONModel with the same contents
    ]

        sol = flux_balance_analysis_dict(
            model,
            Tulip.Optimizer;
            modifications = [
                change_objective("BIOMASS_Ecoli_core_w_GAM"),
                change_constraint("EX_glc__D_e"; lb = -12, ub = -12),
                change_sense(MAX_SENSE),
                change_optimizer_attribute("IPM_IterationsLimit", 110),
                knockout(["b0978", "b0734"]), # knockouts out cytbd
            ],
        )
        @test isapprox(
            sol["BIOMASS_Ecoli_core_w_GAM"],
            0.2725811189335953,
            atol = TEST_TOLERANCE,
        )

        sol = flux_balance_analysis_dict(
            model,
            Tulip.Optimizer;
            modifications = [
                change_objective("BIOMASS_Ecoli_core_w_GAM"),
                change_constraint("EX_glc__D_e"; lb = -12, ub = -12),
                change_sense(MAX_SENSE),
                change_optimizer_attribute("IPM_IterationsLimit", 110),
                knockout("b2779"), # knockouts out enolase
            ],
        )
        @test isapprox(sol["BIOMASS_Ecoli_core_w_GAM"], 0.0, atol = TEST_TOLERANCE)
    end
end


text:
"""
    Simply cuts the elements up into n_procs chunks
    and assigns them in increasing order based on the
    element IDs

    Note: This is generally a REALLY bad idea!
    (but it will get us going)
"""
abstract SimplePartitioner <: Partitioner

function partition!(mesh::Mesh, ::Type{SimplePartitioner})
    n_procs = MPI.Comm_size(MPI.COMM_WORLD)
    proc_id = MPI.Comm_rank(MPI.COMM_WORLD)

    n_elems = length(mesh.elements)

    n_elems_per_proc = div(n_elems,n_procs)

    # We'll just give the extra elements to the last processor
    extra_elems = rem(n_elems, n_procs)

    elems_on_proc = 0
    current_proc = 0

    # Loop over all of the elements and assign their processor IDs
    for elem in mesh.elements
        elem.processor_id = current_proc

        elems_on_proc += 1

        if elems_on_proc >= n_elems_per_proc && current_proc != n_procs-1
            elems_on_proc = 0
            current_proc += 1
        end
    end

    assignNodeProcessorIDs(mesh)
end


text:
using RDatasets, Gadfly

set_default_plot_size(6inch, 3inch)

plot(dataset("plm", "Cigar"), x=:Year, y=:Sales, Scale.x_discrete, Geom.bar)


text:
using BinaryBuilder

# Collection of pre-build pandoc binaries
name = "pandoc"
version = v"2.9.2"
pandoc_ver = "2.9.2.1"

sources = [
    ArchiveSource("https://github.com/jgm/pandoc/releases/download/$pandoc_ver/pandoc-$pandoc_ver-linux-amd64.tar.gz", "5b61a981bd2b7d48c1b4ba5788f1386631f97e2b46d0d1f1a08787091b4b0cf8"; unpack_target = "x86_64-linux-gnu"),
    ArchiveSource("https://github.com/jgm/pandoc/releases/download/$pandoc_ver/pandoc-$pandoc_ver-macOS.zip", "c4847f7a6e6a02a7d1b8dc17505896d8a6e4c2ee9e8b325e47a0468036675307"; unpack_target = "x86_64-apple-darwin14"),
    ArchiveSource("https://github.com/jgm/pandoc/releases/download/$pandoc_ver/pandoc-$pandoc_ver-windows-i386.zip", "db5a8533b7e2ef38114e9788e56530bb6be330c326731692f236218682017d4d"; unpack_target = "i686-w64-mingw32"),
    ArchiveSource("https://github.com/jgm/pandoc/releases/download/$pandoc_ver/pandoc-$pandoc_ver-windows-x86_64.zip", "223f7ef1dd926394ee57b6b5893484e51100be8527bd96eec26e284774863084"; unpack_target = "x86_64-w64-mingw32"),
    FileSource("https://raw.githubusercontent.com/jgm/pandoc/master/COPYING.md", "2b0d4dda4bf8034e1506507a67f80f982131137afe62bf144d248f9faea31da4"),
    FileSource("https://raw.githubusercontent.com/jgm/pandoc/master/COPYRIGHT", "39be98cc4d2906dd44abf8573ab91557e0b6d51f503d3a889dab0e8bcca1c43f"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/
mkdir -p "${bindir}"
if [[ "${target}" != *-mingw* ]]; then
    subdir="bin/"
fi
cp ${target}/pandoc-*/${subdir}pandoc${exeext} ${bindir}
cp ${target}/pandoc-*/${subdir}pandoc-citeproc${exeext} ${bindir}
chmod +x ${bindir}/*
install_license COPYRIGHT
install_license COPYING.md
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64),
    MacOS(),
    Windows(:x86_64),
    Windows(:i686),
]

# The products that we will ensure are always built
products = [
    ExecutableProduct("pandoc", :pandoc),
    ExecutableProduct("pandoc-citeproc", :pandoc_citeproc),
]

# Dependencies that must be installed before this package can be built
dependencies = [
	Dependency("Zlib_jll"),
	Dependency("Libiconv_jll")
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: tbeason/Yggdrasil

max_issues_repo_name: tbeason/Yggdrasil

max_forks_repo_name: tbeason/Yggdrasil


text:
using MyProject
using Test

@testset "MyProject.jl" begin
    # Write your tests here.
    @test foo(0) <1E-16
end


text:
function logsumexp_dim1(x_km)
  max_m = maximum(x_km, dims=1)
  exp_km = exp.(x_km .- max_m)
  s_m = log.(sum(exp_km, dims=1)) .+ max_m
end

function logsumexp_dim2(x_km)
  max_k = maximum(x_km, dims=2)
  exp_km = exp.(x_km .- max_k)
  s_k = log.(sum(exp_km, dims=2)) .+ max_k
end

function wham_equation!(F, x, K, M, bias_km, log_Neff_km, log_heff_m)
    log_prob_m = reshape(x[1:M], 1, M)
    f_k = reshape(x[(M+1):(M+K)], K, 1)

    # 1st equation
    log_denominator_km = f_k .- bias_km
    log_denominator_km = log_denominator_km .+ log_Neff_km
    log_denominator_m  = logsumexp_dim1(log_denominator_km)
    log_numerator_m = log_heff_m
    F[1:M] = - log_prob_m .+ log_numerator_m .- log_denominator_m

    # 2nd equation
    F[(M+1):(M+K)] = f_k .+ logsumexp_dim2(log_prob_m .- bias_km)
    F[(M+1):(M+K)] = F[(M+1):(M+K)] .- F[M+1]
end

"""
wham

estimates (reduced) free energies of umbrella-windows and potential of mean force in data-bins by using the WHAM equations.

References
[1] S. Kumar, D. Bouzida, R. H. Swendsen, P. A. Kollman, and
    J. M. Rosenberg, J. Comput. Chem. 13, 1011 (1992).
"""
function wham(h_km, bias_km; ftol=1e-10, iterations=10^5)
    # K: number of umbrella-windows
    K = size(h_km, 1)

    # M: number of bins
    M = size(h_km, 2)

    # conversion from Int64 to Float64
    h_km = map(Float64, h_km)

    # calculate counts (N_k)
    # N_k: number of samples counted in the histogram from umbrella-window (k)
    N_k = sum(h_km, dims=2)
    log_N_k = zeros(Float64, K, 1)
    log_N_k[N_k .> 0.0] .= log.(N_k[N_k .> 0.0])

    # calculate effective histogram (heff_m)
    # heff_m: effective number of independent samples in data-bin (m)
    g_km = ones(Float64, K, M)
    heff_m = sum(h_km./g_km, dims=1)
    log_heff_m = zeros(Float64, 1, M)
    log_heff_m[:, :] .= -10.0 #double(log(eps(Float32)))
    log_heff_m[heff_m .> 0.0] .= log.(heff_m[heff_m .> 0.0])

    # calculate effective counts (Neff_km)
    # Neff_km: effective number of independent samples for umbrella-window (k) in data-bin (m)
    Neff_km = N_k ./ g_km
    log_Neff_km = zeros(Float64, K, M)
    log_Neff_km[:, :] .= -10.0 #double(log(eps(Float32)))
    log_Neff_km[Neff_km .> 0] .= log.(Neff_km[Neff_km .> 0])

    #F_tmp = zeros(Float64, M+K)
    #initial_x = zeros(Float64, M+K)
    #wham_equation!(F_tmp, initial_x, K, M, bias_km, log_Neff_km, log_heff_m)

    f!(F, x) = wham_equation!(F, x, K, M, bias_km, log_Neff_km, log_heff_m)
    initial_x = zeros(Float64, M+K)
    sol = nlsolve(f!, initial_x, autodiff=:forward, ftol=ftol, iterations=iterations)

    log_prob_m = reshape(sol.zero[1:M], 1, M)
    f_k = reshape(sol.zero[(M+1):(M+K)], K, 1)
    pmf_m = - log_prob_m
    f_k .- f_k[1, 1], pmf_m .- pmf_m[1, 1]
end

"""
wham_iteration

Old version of wham(). this function iteratively solve wham equations.

References
[1] S. Kumar, D. Bouzida, R. H. Swendsen, P. A. Kollman, and
    J. M. Rosenberg, J. Comput. Chem. 13, 1011 (1992).
"""
function wham_iteration(h_km, bias_km)
    TOLERANCE = 10^(-8)

    # K: number of umbrella-windows
    K = size(h_km, 1)

    # M: number of bins
    M = size(h_km, 2)

    # conversion from Int64 to Float64
    h_km = map(Float64, h_km)

    # calculate counts (N_k)
    # N_k: number of samples counted in the histogram from umbrella-window (k)
    N_k = sum(h_km, dims=2)
    log_N_k = zeros(Float64, K, 1)
    log_N_k[N_k .> 0.0] .= log.(N_k[N_k .> 0.0])

    # calculate effective histogram (heff_m)
    # heff_m: effective number of independent samples in data-bin (m)
    g_km = ones(Float64, K, M)
    heff_m = sum(h_km./g_km, dims=1)
    log_heff_m = zeros(Float64, 1, M)
    log_heff_m[:, :] .= -10.0 #double(log(eps(Float32)))
    log_heff_m[heff_m .> 0.0] .= log.(heff_m[heff_m .> 0.0])

    # calculate effective counts (Neff_km)
    # Neff_km: effective number of independent samples for umbrella-window (k) in data-bin (m)
    Neff_km = N_k ./ g_km
    log_Neff_km = zeros(Float64, K, M)
    log_Neff_km[:, :] .= -10.0 #double(log(eps(Float32)))
    log_Neff_km[Neff_km .> 0] .= log.(Neff_km[Neff_km .> 0])

    # solve the WHAM equations by self-consistent iteration
    f_k = zeros(Float64, K, 1)
    check_convergence = Inf

    log_prob_m = zeros(Float64, 1, M)
    count_iteration = 0
    while check_convergence > TOLERANCE

      # 1st equation
      log_denominator_km = f_k .- bias_km
      log_denominator_km = log_denominator_km .+ log_Neff_km
      log_denominator_m  = logsumexp_dim1(log_denominator_km)

      log_numerator_m = log_heff_m

      log_prob_m = log_numerator_m .- log_denominator_m

      # 2nd equation
      f_k_new = - logsumexp_dim2(log_prob_m .- bias_km)
      f_k_new = f_k_new .- f_k_new[1]

      # check convergence
      check_convergence = maximum(abs.(f_k_new - f_k))./std(f_k_new)
      f_k = f_k_new;

      count_iteration = count_iteration + 1
      if mod(count_iteration, 100) == 0
        @printf("%dth iteration  delta = %e  tolerance = %e\n", count_iteration, check_convergence, TOLERANCE)
        @printf("free energies = ")
        for k = 1:K
          @printf("%f ", f_k[k])
        end
        @printf("\n")
        @printf("\n")
      end
    end

    pmf_m = - log_prob_m
    f_k .- f_k[1, 1], pmf_m .- pmf_m[1, 1]
end


text:
@testset "Testing DScribe SOAP Descriptor with dc Si." begin
@info("Testing DScribe SOAP Descriptor for dc Si.")
using DescriptorZoo, JuLIP, Test

at = bulk(:Si, cubic=true)
desc = dscribe_soap(at, 6.5, n_max=4, l_max=4)
soap_ref = [0.47285277, 1.7203233, -1.6327636, 4.0070906, 6.258845, -5.940287, 14.578514, 5.637943, -13.836509, 33.957237, 3.5743687e-33, -1.0012961e-32, 1.5516055e-32, -4.7689123e-33, 2.8058653e-32, -4.353206e-32, 1.3339438e-32, 6.878232e-32, -2.1905732e-32, 8.337165e-33, 7.119684e-33, -2.2345982e-32, 2.8977972e-32, -2.1937631e-32, 7.015937e-32, -9.1644326e-32, 6.978438e-32, 1.4435643e-31, -1.2865728e-31, 1.2816743e-31, 0.025995005, -0.113665685, 0.20420578, 0.19071767, 0.49701422, -0.8929096, -0.8339315, 1.6041543, 1.4981973, 1.3992391, 0.00081991457, -0.0046709054, 0.009903433, 0.037733737, 0.026609303, -0.056418065, -0.21496229, 0.11961975, 0.4557713, 1.7365651]
soap_now = vcat(desc[1,:]...)
tst = [isapprox(soap_now[x], soap_ref[x]; atol=0.00001) for x=1:length(soap_ref)]
println(@test all(tst))
end


text:
export KarateClub

"""
    Zachary's Karate Club

The Karate Club dataset originally appeared in Ref [1].

The network contains 34 nodes (members of the karate club).
The nodes are connected by 78 undirected and unweighted edges.
The edges indicate if the two members interacted outside the club.

The node labels indicate which community or the karate club the member belongs to.
The club based labels are as per the original dataset in Ref [1].
The community labels are obtained by modularity-based clustering following Ref [2].
The data is retrieved from Ref [3] and [4].
One node per unique label is used as training data.

# References

[1]: [An Information Flow Model for Conflict and Fission in Small Groups](http://www1.ind.ku.dk/complexLearning/zachary1977.pdf)
[2]: [Semi-supervised Classification with Graph Convolutional Networks](https://arxiv.org/abs/1609.02907)
[3]: [PyTorch Geometric Karate Club Dataset](https://pytorch-geometric.readthedocs.io/en/latest/_modules/torch_geometric/datasets/karate.html#KarateClub)
[4]: [NetworkX Zachary's Karate Club Dataset](https://networkx.org/documentation/stable/_modules/networkx/generators/social.html#karate_club_graph)
"""
struct KarateClub
    metadata::Dict{String, Any}
    graphs::Vector{Graph}
end

function KarateClub()
    src = [
        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2,
        2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5,
        6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 11, 11,
        11, 12, 13, 13, 14, 14, 14, 14, 14, 15, 15, 16, 16, 17, 17, 18, 18,
        19, 19, 20, 20, 20, 21, 21, 22, 22, 23, 23, 24, 24, 24, 24, 24, 25,
        25, 25, 26, 26, 26, 27, 27, 28, 28, 28, 28, 29, 29, 29, 30, 30, 30,
        30, 31, 31, 31, 31, 32, 32, 32, 32, 32, 32, 33, 33, 33, 33, 33, 33,
        33, 33, 33, 33, 33, 33, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34,
        34, 34, 34, 34, 34, 34
    ]
    target = [
        2, 3, 4, 5, 6, 7, 8, 9, 11, 12, 13, 14, 18, 20, 22, 32, 1, 3, 4, 8,
        14, 18, 20, 22, 31, 1, 2, 4, 8, 9, 10, 14, 28, 29, 33, 1, 2, 3, 8,
        13, 14, 1, 7, 11, 1, 7, 11, 17, 1, 5, 6, 17, 1, 2, 3, 4, 1, 3, 31,
        33, 34, 3, 34, 1, 5, 6, 1, 1, 4, 1, 2, 3, 4, 34, 33, 34, 33, 34, 6,
        7, 1, 2, 33, 34, 1, 2, 34, 33, 34, 1, 2, 33, 34, 26, 28, 30, 33,
        34, 26, 28, 32, 24, 25, 32, 30, 34, 3, 24, 25, 34, 3, 32, 34, 24,
        27, 33, 34, 2, 9, 33, 34, 1, 25, 26, 29, 33, 34, 3, 9, 15, 16, 19,
        21, 23, 24, 30, 31, 32, 34, 9, 10, 14, 15, 16, 19, 20, 21, 23, 24,
        27, 28, 29, 30, 31, 32, 33
    ]

    labels_clubs = [
            1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1,
            0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    labels_comm = [
                1, 1, 1, 1, 3, 3, 3, 1, 0, 1, 3, 1, 1, 1, 0, 0, 3, 1, 0, 1, 0, 1,
                0, 0, 2, 2, 0, 0, 2, 0, 0, 2, 0, 0]

    node_data = (; labels_clubs, labels_comm)
    g = Graph(; num_nodes=34, num_edges=156, edge_index=(src, target), node_data)

    metadata = Dict{String, Any}()
    return KarateClub(metadata, [g])
end

Base.length(d::KarateClub) = length(d.graphs)
Base.getindex(d::KarateClub) = d.graphs
Base.getindex(d::KarateClub, i) = getindex(d.graphs, i)


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__()

module CUDex

# package code goes here
const libcudex = Libdl.find_library(["libcudex"], [dirname(@__FILE__)])

dir(path...) = joinpath(dirname(dirname(@__FILE__)),path...)

export gpu, DexArray, relu, invx, sigm, logp, conv_f, pool_f, rnn_ft

include("gpu.jl")
include("dexptr.jl")
include("dexarray.jl")
include("basic_bi.jl")
include("basic_uni.jl")
include("cublas.jl")
include("cudnn.jl")

# See if we have a gpu at initialization:
function __init__()
    try
        r = gpu(true)
        info(r >= 0 ? "Using GPU $r" : "No GPU found, Using CPU")
    catch e
        warn("Using CPU: $e")
        gpu(false)
    end
end

end # module


text:
# Copyright 2018 IBM Corp.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

module WiringRDF
export wiring_diagram_to_rdf, semantic_graph_to_rdf

using Parameters
using Serd
using Catlab
using Catlab.WiringDiagrams, Catlab.WiringDiagrams.WiringDiagramSerialization

using ...Doctrine
using ..RDFUtils

const R = RDF.Resource

# Configuration
###############

default_box_rdf_node(box::String) = RDF.Blank(box)
default_port_rdf_node(box::String, port::String) = RDF.Blank("$box:$port")
default_wire_rdf_node(wire::String) = RDF.Blank(wire)

function default_value_to_rdf(node::RDF.Node, value)
  if isnothing(value)
    RDF.Statement[]
  else
    [ RDF.Triple(node, R("monocl","value"), RDF.Literal(string(value))) ]
  end
end

@with_kw struct RDFConfig
  box_rdf_node::Function = default_box_rdf_node
  port_rdf_node::Function = default_port_rdf_node
  wire_rdf_node::Function = default_wire_rdf_node
  box_value_to_rdf::Function = default_value_to_rdf
  port_value_to_rdf::Function = default_value_to_rdf
  wire_value_to_rdf::Function = default_value_to_rdf
  include_provenance::Bool = false
end

# Wiring diagrams
#################

""" Convert a wiring diagram to RDF.
"""
function wiring_diagram_to_rdf(diagram::WiringDiagram; kw...)
  config = RDFConfig(; kw...)
  box_to_rdf(config, diagram, Int[])
end

function box_to_rdf(config::RDFConfig, diagram::WiringDiagram, path::Vector{Int})
  node = config.box_rdf_node(box_id(path))
  stmts = RDF.Statement[
    [ RDF.Triple(node, R("rdf","type"), R("monocl","WiringDiagram")) ];
    ports_to_rdf(config, diagram, path);
  ]

  # Add RDF for boxes.
  for v in box_ids(diagram)
    box_node, box_stmts = box_to_rdf(config, box(diagram, v), [path; v])
    push!(stmts, RDF.Triple(node, R("monocl","hasBox"), box_node))
    append!(stmts, box_stmts)
  end

  # Add RDF for wires.
  box_rdf_node = port::Port -> config.box_rdf_node(
    box_id(diagram, [path; port.box]))
  port_rdf_node = port::Port -> config.port_rdf_node(
    box_id(diagram, [path; port.box]), port_name(diagram, port))

  for (i, wire) in enumerate(wires(diagram))
    wire_node = config.wire_rdf_node(wire_id(path, i))
    src_port_node = port_rdf_node(wire.source)
    tgt_port_node = port_rdf_node(wire.target)
    append!(stmts, [
      RDF.Triple(node, R("monocl","hasWire"), wire_node),
      RDF.Triple(wire_node, R("rdf","type"), R("monocl","Wire")),
      RDF.Triple(wire_node, R("monocl","source"), src_port_node),
      RDF.Triple(wire_node, R("monocl","target"), tgt_port_node),
      RDF.Triple(src_port_node, R("monocl","wire"), tgt_port_node),
    ])
    append!(stmts, config.wire_value_to_rdf(wire_node, wire.value))
    if config.include_provenance
      src_box_node = box_rdf_node(wire.source)
      tgt_box_node = box_rdf_node(wire.target)
      append!(stmts, [
        RDF.Triple(tgt_box_node, R("prov","wasInformedBy"), src_box_node),
        RDF.Triple(tgt_port_node, R("prov","wasDerivedFrom"), src_port_node),
      ])
    end
  end

  (node, stmts)
end

function box_to_rdf(config::RDFConfig, box::Box, path::Vector{Int})
  node = config.box_rdf_node(box_id(path))
  stmts = RDF.Statement[
    [ RDF.Triple(node, R("rdf","type"), R("monocl","Box")) ];
    config.box_value_to_rdf(node, box.value);
    ports_to_rdf(config, box, path);
  ]
  (node, stmts)
end

function ports_to_rdf(config::RDFConfig, box::AbstractBox, path::Vector{Int})
  name = box_id(path)
  node = config.box_rdf_node(name)
  port_node = port -> config.port_rdf_node(name, port)
  inputs, outputs = input_ports(box), output_ports(box)
  nin, nout = length(inputs), length(outputs)
  owl_inputs_outputs(node, port_node, nin, nout, index=true) do cell, is_input, i
    port_value = (is_input ? inputs : outputs)[i]
    RDF.Statement[
      [ RDF.Triple(cell, R("rdf","type"), R("monocl","Port")) ];
      config.port_value_to_rdf(cell, port_value);
      if config.include_provenance
        if is_input
          [ RDF.Triple(node, R("prov","used"), cell) ]
        else
          [ RDF.Triple(cell, R("prov","wasGeneratedBy"), node) ]
        end
      else [] end;
    ]
  end
end

# Semantic flow graphs
######################

""" Convert a semantic flow graph to RDF.
"""
function semantic_graph_to_rdf(diagram::WiringDiagram, concept_rdf_node::Function;
    include_provenance::Bool=false, kw...)

  function box_expr_to_rdf(node::RDF.Node, expr::Monocl.Hom)
    concept = if head(expr) == :generator
      concept_rdf_node(expr)
    elseif head(expr) == :construct
      concept_rdf_node(codom(expr))
    else
      error("Cannot serialize Monocl morphism of type: ", head(expr))
    end
    RDF.Statement[
      [ RDF.Triple(node, R("monocl","isConcept"), concept) ];
      include_provenance && head(expr) == :generator ?
        [ RDF.Triple(node, R("rdf","type"), concept) ] : [];
    ]
  end

  function port_expr_to_rdf(node::RDF.Node, expr::Monocl.Ob)
    concept = concept_rdf_node(expr)
    RDF.Statement[
      [ RDF.Triple(node, R("monocl","isConcept"), concept) ];
      include_provenance ?
        [ RDF.Triple(node, R("rdf","type"), concept) ] : [];
    ]
  end

  wiring_diagram_to_rdf(diagram;
    box_value_to_rdf = box_expr_to_rdf,
    port_value_to_rdf = port_expr_to_rdf,
    include_provenance = include_provenance,
    kw...)
end

end

max_stars_repo_name: cclauss/semanticflowgraph

max_issues_repo_name: cclauss/semanticflowgraph

max_forks_repo_name: cclauss/semanticflowgraph


text:
using .Distributions # Requires.jl replace usual syntax with this

# Define recipe for 2d MvNormal

@recipe function f(dist::MvNormal, alpha=0.0001)
    @assert length(dist) == 2 "MvNormal Plot supports only 2d distribution"

    seriestype --> :contours
    x_mdist = Normal(dist.μ[1], sqrt(dist.Σ[1,1]))
    y_mdist = Normal(dist.μ[2], sqrt(dist.Σ[2,2]))
    xmin, xmax = quantile.(x_mdist, [alpha, 1-alpha])
    ymin, ymax = quantile.(y_mdist, [alpha, 1-alpha])
    (range(xmin, xmax, length=100), range(ymin, ymax, length=100), (x,y)->pdf(dist, [x,y]))
end

function find_limit(mdist_list, prior_p, alpha)
    q_l = [quantile(mdist, alpha) for mdist in mdist_list]
    l = minimum(q_l)
    r = maximum(q_l)

    for i in 1:100
        s = (l+r)/2
        sum_cdf = sum([cdf(mdist, s) * p for (mdist, p) in zip(mdist_list, prior_p)])
        if sum_cdf < alpha
            l = s
        else
            r = s
        end
    end

    left = (l + r) / 2
end

function find_limit2(mdist_list, prior_p, alpha)
    find_limit(mdist_list, prior_p, alpha), find_limit(mdist_list, prior_p, 1-alpha)
end

@recipe function f(mix_dist::MixtureModel{Multivariate, Continuous, <:MvNormal}, alpha=0.0001)
    @assert length(mix_dist.components[1]) == 2 "MixtureModel Plot supports only 2d component distribution"

    seriestype --> :contours

    x_mdist_list = [Normal(comp.μ[1], sqrt(comp.Σ[1,1])) for comp in mix_dist.components]
    y_mdist_list = [Normal(comp.μ[2], sqrt(comp.Σ[2,2])) for comp in mix_dist.components]

    prior_p = mix_dist.prior.p

    xmin, xmax = find_limit2(x_mdist_list, prior_p, alpha)
    ymin, ymax = find_limit2(y_mdist_list, prior_p, alpha)

    (range(xmin, xmax, length=100), range(ymin, ymax, length=100), (x,y)->pdf(mix_dist, [x,y]))
end


text:
function mw_lrstat_vec(x, cov, cov_a)
  nsim = size(x, 2)
  q = size(x, 1)
  # Compute likelihood under ho @
  lf_h0 = mw_marg_lf(x, cov)
  lf_ha = NaN*zeros(size(x,2), size(cov_a, 3))
  for i ∈ 1:size(cov_a, 3)
    lf_ha[:, i] = mw_marg_lf(x, cov_a[:,:,i])
  end
  lf_ha_avg = mean(lf_ha, 2)
  lrstat = lf_ha_avg./lf_h0
  return lrstat
end


text:
config = MOI.Test.TestConfig()
optimize!(mock) = MOIU.mock_optimize!(mock, MOI.INFEASIBLE, tuple(),
                                      MOI.INFEASIBILITY_CERTIFICATE)
for mock in mocks(optimize!)
    Tests.choi_test(mock, config)
end


text:
using GibbsSampling
using Test
using CSV

## Tips de RAz para ser bien pro
# Em;iezas un Pkg ocn PkgTemplates.jl
# generas un paquete y le haces `dev MyPkg` para que genere el path correcto
# haces todo tu codigo en runtests.jl e iteras con tests

@testset "GibbsSampling.jl" begin
    # Write your own tests here.
    @test 1 == 1
end

abstract type AbstractVector end
struct PVec{T,S} <: AbstractVector

    x :: T
    y :: S
end

pvecx(v :: PVec{T,S}) where {T,S} = v.x
pvecy(v :: PVec{T,S}) where {T,S} = v.y
# TODO Setters en vez de getters

p1 = PVec(1,1)
p2 = PVec(1,0)
p3 = PVec(true,false)

@test p1.x==1
@test p1.y==1
@test p2.x==1
@test p2.y==0
@test p3.x==true
@test p3.y==false

@test typeof(p3.x) <: Bool && typeof(p3.y) <: Bool

path = "/home/pepefv97/.julia/dev/GibbsSampling/src/"
cd(path)
xa = CSV.read("centroids.csv",header=false)
xa = convert(Matrix, xa)
# Si entradas estan al reves correr metodo distinto para escribir matriz desde numpy
# probalemente cambiar el a.flatten()
mat = [PVec(xa[i],xa[i+1]) for i in 1:2:length(xa)-1]
mat = reshape(mat,192,202)

kronecker(i,j, x) = i == j ? one(x) : zero(x)
Base.one(v :: PVec{T,S}) where {T,S} = one(pvecx(v))
Base.zero(v :: PVec{T,S}) where {T,S} = zero(pvecx(v))

@test kronecker(1,1,1) == 1
@test kronecker(1,1,1.0) == 1.0
@test kronecker(1,2,1) == 0
@test kronecker(1,2,1.0) == 0.0
@test kronecker(1,1,true) == true
@test kronecker(1,2,true) == false

p4 = PVec(4.0,Bool)
@test one(1) == 1.0
@test zero(p4) == 0.0

choice() = rand(1:10)
@test 1 <= choice() <= 10

function energia(i,j,mat, r)
    @assert 2 == length(size(mat))
    res1 :: = interactuar_vecinos(i, j, mat, r)
    res2 = res1 - delta_vecinos(i, j, mat, r)
end

delta_vecinos(i,j,mat, r) = kronecker(pvecx(mat[i,j]), r, mat[i,j])
@test delta_vecinos(6,4,mat, choice()) == 0.0

#@inline function interactuar_vecinos(i,j,x, r) :: typeof(one(x[1,1]))
#    sum(skipmissing((vecino(findvecinoup(i,j,x),r) ,vecino(findvecinodown(i,j,x),r),vecino(findvecinoleft(i,j,x),r),vecino(findvecinoright(i,j,x),r))))
#end
@inline function interactuar_vecinos(i,j,x,r)
    m=vecino(findvecinoup(i,j,x),r)+vecino(findvecinodown(i,j,x),r)+vecino(findvecinoleft(i,j,x),r)+vecino(findvecinoright(i,j,x),r)
end
# Para cambiar a forntera toroidal usar mod1
@inbounds findvecinoup(i,j,x) =  i > 1 ? x[i-1,j] : missing
@inbounds findvecinodown(i,j,x) = i < size(x)[1] ? x[i+1,j] : missing
@inbounds findvecinoleft(i,j,x) = j > 1 ? x[i,j-1] : missing
@inbounds findvecinoright(i,j,x) = j < size(x)[2] ? x[i,j+1] : missing

vecino(v :: PVec{T,S}, r) where {T,S} = (pvecx(v) - r)^2
#vecino(x :: Missing, r) = missing
vecino(x :: Missing, r) = 0

# TODO

@test findvecinoup(1,1,mat) === missing
@test findvecinodown(192,1,mat) === missing
@test findvecinoright(1,202,mat) === missing
@test findvecinoleft(1,1,mat) === missing

@test findvecinoup(1,10,mat) === missing

abstract type AbstractSimulation end
struct CircularSimulation <: AbstractSimulation
    N :: Int
    mat :: Matrix
end
struct MCSimulation <: AbstractSimulation
    N :: Int
    mat :: Matrix
end

# energias(s :: CircularSimulation, mat) = println("tada")
# energias(s :: MCSimulation, mat) = println("asa")

"""
exp(- energia(i,j,mat,r) * log(i + j + k))/ exp( - energia(i,j,mat,pvecx(mat[i,j])) * log(i + j + k))

counter incrementa con cada recorrido entero de la matriz
"""
@inbounds function cociente(i,j,mat,counter,r)
    set_zero_subnormals(true)
    exp(- energia(i,j,mat,r) * log(i + j + counter) + energia(i,j,mat,pvecx(mat[i,j])) * log(i + j + counter))
end

@inbounds function update!(i,j,mat, α, r)
    if α >= 1.0
        mat[i,j] = PVec(convert(typeof(one(mat[i,j])),r),pvecy(mat[i,j])) # TODO Setters
    else
        ρ = rand() # esto vive en [0,1] uniforme
        if ρ <= α
            return
        else
            mat[i,j] = PVec(convert(typeof(one(mat[i,j])),r), pvecy(mat[i,j]))
        end
    end
end

using UnicodePlots
miniplot(mat) = spy(sim.mat .|> pvecx)
@inbounds function simulate(s :: AbstractSimulation)
    # TODO usar Random.seed!(7)
    I,J = [size(s.mat)...]
    N = s.N

    @inbounds for counter in 1:N
        for idxj in 1:J
            for idxi in 1:I

                # Desechar algunos casos
                if pvecy(mat[idxi,idxj]) == one(mat[1,1])
                    continue
                end
                r = choice()
                # En casos de interes, checar propension energetica
                α = cociente(idxi,idxj,mat, counter, r)

                # Reemplazar valores
                update!(idxi,idxj, mat, α, r)
            end
        end
    end

    return mat
    # TODO guardar las ultimas mil
end

### Benchmarking

### Threads
mattest = deepcopy(mat2)
Random.seed!(7)

@btime PVec(1,2)
@btime cociente(1,1,$mattest,1,4)

@test update!(1,1,mattest, )


text:
workspace()
Base.Test.@testset "controller.jl" begin
    include("controller.jl")
end

workspace()
Base.Test.@testset "params.jl" begin
    include("params.jl")
end

workspace()
Base.Test.@testset "octo/assoc.jl" begin
    include("octo/assoc.jl")
end

workspace()
Base.Test.@testset "octo/changeset.jl" begin
    include("octo/changeset.jl")
end

workspace()
Base.Test.@testset "octo/natural.jl" begin
    include("octo/natural.jl")
end

workspace()
Base.Test.@testset "view.jl" begin
    include("view.jl")
end

workspace()
Base.Test.@testset "renderers/text.jl" begin
    include("renderers/text.jl")
end

workspace()
Base.Test.@testset "renderers/json.jl" begin
    include("renderers/json.jl")
end

workspace()
Base.Test.@testset "renderers/html.jl" begin
    include("renderers/html.jl")
end

workspace()
Base.Test.@testset "renderers/html/tag.jl" begin
    include("renderers/html.jl")
end

workspace()
Base.Test.@testset "renderers/html/form.jl" begin
    include("renderers/html.jl")
end

workspace()
Base.Test.@testset "renderers/markdown.jl" begin
    include("renderers/markdown.jl")
end

workspace()
Base.Test.@testset "renderers/mustache.jl" begin
    include("renderers/mustache.jl")
end

workspace()
Base.Test.@testset "logger.jl" begin
    include("logger.jl")
end

workspace()
Base.Test.@testset "router.jl" begin
    include("router.jl")
end

workspace()
Base.Test.@testset "scope.jl" begin
    include("scope.jl")
end

workspace()
Base.Test.@testset "endpoint.jl" begin
    include("endpoint.jl")
end

workspace()
Base.Test.@testset "config.jl" begin
    include("config.jl")
end

workspace()
Base.Test.@testset "plugins.jl" begin
    include("plugins.jl")
end

workspace()
Base.Test.@testset "plug/static.jl" begin
    include("plug/static.jl")
end

workspace()
Base.Test.@testset "plug/logger.jl" begin
    include("plug/logger.jl")
end

workspace()
Base.Test.@testset "server.jl" begin
    include("server.jl")
end


text:
function ADScalarObjective(f,x0::AbstractArray)
    Hres = DiffResults.HessianResult(x0)
    function _g(df,x,Hresult)
        ForwardDiff.gradient!(Hresult,f,x)
        df .= DiffResults.gradient(Hresult)
        df
    end
    function _fg(df,x,Hresult)
        ForwardDiff.gradient!(Hresult,f,x)
        df .= DiffResults.gradient(Hresult)
        fx = DiffResults.value(Hresult)
        return fx,df
    end

    function _fgh(df,d2f,x,Hresult)
        ForwardDiff.hessian!(Hresult,f,x)
        d2f .= DiffResults.hessian(Hresult)
        df .= DiffResults.gradient(Hresult)
        fx = DiffResults.value(Hresult)
        return fx,df,d2f
    end

    function h(d2f,x)
        ForwardDiff.hessian!(d2f,f,x)
        d2f
    end
    g(df,x) = _g(df,x,Hres)
    fg(df,x) = _fg(df,x,Hres)
    fgh(df,d2f,x) =  _fgh(df,d2f,x,Hres)
    return ScalarObjective(f=f,
    g=g,
    fg=fg,
    fgh=fgh,
    h=h)
end

function ADScalarObjective(f,x0::Number)
    function g(x)
        return derivative(f,x)
    end

    function fg(∂fx,x)
        ∂fx,x = f∂f(f,x)
        return ∂fx,x
    end

    function fgh(∂fx,∂2fx,x)
        fx,∂fx,∂2fx = f∂f∂2f(f,x)
        return fx,∂fx,∂2fx
    end

    function h(∂2fx,x)
        ∂2fx = derivative(g,x)
        return ∂2fx
    end

    return ScalarObjective(f=f,
    g=g,
    fg=fg,
    fgh=fgh,
    h=h)
end
"""
    function optimize(f,x0,method=LineSearch(Newton()), options=OptimizationOptions())
"""
function optimize(f,x0::Number,method=LineSearch(Newton()),options=OptimizationOptions())
    scalarobj = ADScalarObjective(f,x0)
    optprob = OptimizationProblem(scalarobj; inplace=false)
    return NLSolvers.solve(optprob, x0, method,options)
end


text:
using Flux, DiffEqFlux, OrdinaryDiffEq,  StatsBase, RecursiveArrayTools
using Flux: onehotbatch
using MLDatasets:MNIST
using Base.Iterators: repeated,partition

batch_size=10

train_x, train_y = MNIST.traindata();
test_x, test_y = MNIST.testdata();
train_y_hot = onehotbatch(train_y,0:9);

train_data = [(reshape(train_x[:,:,i],(28,28,1,batch_size)),train_y_hot[:,i]) for i in partition(1:60_000,batch_size)];

downsample = Chain(
                   Conv((3,3),1=>32,stride=1),
                   BatchNorm(32,relu),
                   Conv((4,4),32=>32,stride=2,pad=1),
                   BatchNorm(32,relu),
                   Conv((4,4),32=>32,stride=2,pad=1)
                  )
dudt = Chain(
          BatchNorm(32,relu),
          Conv((3,3),32=>32,stride =1, pad=1),
          BatchNorm(32),
          Conv((3,3),32=>32,stride =1, pad=1),
          BatchNorm(32)
          )
classify = Chain(
                 BatchNorm(32,relu),
                 MeanPool((6,6)),
                 x->view(x,1,1,:,:), #Need More General Flatten Here
                 Dense(32,10)
                )

ps = Flux.params(dudt)

function n_ode(batch_u0, batch_t)
    neural_ode(dudt,batch_u0,(0.,25.),Tsit5(),
               save_start=false,
               saveat=batch_t,      # Ugly way to only get sol[end]
               reltol=1e-3,abstol=1e-3)
end

model = Chain(downsample,u->n_ode(u,[25.])[:,:,:,:,end],classify) # Further ugliness getting sol[end]

loss(x,y) = Flux.mse(model(x),y)

loss(train_data[1]...) # Works and is tracked
# This fails with neural_ode_rd with T undefined error

opt = ADAM(0.1)
Ps = params(model,ps)

Tracker.gradient(()->loss(xt,yt),Ps) # Fails with BoundsError
#ERROR: BoundsError: attempt to access 6×6×32×10×1 Array{Float32,5} at index [Base.Slice(Base.OneTo(6)), 1]

Flux.train!(loss, Ps, train_data, opt, cb = ()) # Fails with BoundsError

function ODEfunc(du,u,p,t)
    du.= dudt(u)
end
prob = ODEProblem(ODEfunc,downsample(xt),(0.,25.))

solve(prob,Tsit5(),saveat=25.)

max_stars_repo_name: jessebett/model-zoo

max_issues_repo_name: jessebett/model-zoo

max_forks_repo_name: jessebett/model-zoo


text:
abstract Cut{R}
export Cut, addcuts, getcuts

(::Type{T}){T<:Cut{1}}() = T(Cut{1}, CutCount{1}())
(::Type{T}){T<:Cut{2}}() = T(Cut{2}, CutCount{2}())
(::Type{T}){T<:Cut{3}}() = T(Cut{3}, CutCount{3}())

# inelegant code to extract the "native" order for a cut;
# for example, transitivity is enforced by constraining every
# triplet of boxes, and so extract_cut_order(Transitivity) = 3
extract_cut_order{C<:Cut}(::C) = _cut_order(super(C))
_cut_order{R}(::Type{Cut{R}}) = R

function get_groups{T}(prob::Problem, cut::Cut, cutcount::CutCount{T})
    N = prob.N
    R = extract_cut_order(cut)
    if R < T
        error("Cannot add 𝒪^$R cuts when only 𝒪^$T are available")
    elseif R == T
        return combinations(1:N, R)
    else
        function _lt(p1,p2)
            v1 = mapreduce(+, combinations(p1,2)) do x
                get_metrics(cutcount.metric, prob, x...)
            end
            v2 = mapreduce(+, combinations(p2,2)) do x
                get_metrics(cutcount.metric, prob, x...)
            end
            v1 < v2
        end
        ordered_combos = sort!(collect(combinations(1:N,R)), lt=_lt, rev=true)
        numcuts = min(cutcount.num_rounds*N, length(ordered_combos))
        return ordered_combos[1:numcuts]
    end
end

get_groups(prob::Problem, cut::Cut, v::Vector{Tuple}) = v

function addcuts(model::Model, prob::Problem, vars::Variables, cut::Cut)
    groups = get_groups(prob, cut, cut.groups)
    for group in groups
        cuttingplanes = getcut(prob, vars, cut, group...)
        for c in cuttingplanes
            addConstraint(model, c)
        end
    end
end

function getcuts(model::Model, prob::Problem, vars::Variables, cut::Cut)
    groups = get_groups(prob, cut, cut.groups)
    cuts = LinearConstraint[]
    sizehint!(cuts, length(groups))
    for group in groups
        append!(cuts, getcut(prob, vars, cut, group...))
    end
    return cuts
end

include("integer_variable_inequalities.jl")
include("literature_cuts.jl")
include("objective_cuts.jl")
include("three_box_cuts.jl")
include("tightened_inequalities.jl")
include("upperbound_inequalities.jl")


text:
using OffsetArrays

@testset "Iteration Bound Tests" begin
  function masktest_incr1_none1start!(y, x)
    @turbo for i ∈ 0:20
      y[i] = x[i] + 2
    end
  end
  function masktest_incr2_none1start!(y, x)
    @turbo for i ∈ 0:2:20
      y[i] = x[i] + 2
    end
  end

  x = OffsetVector(rand(24), 0:23)
  y = copy(x)
  masktest_incr1_none1start!(y, x)
  @test y == x .+ ifelse.(axes(x, 1) .≤ 20, 2, 0)
  @turbo y .= x
  @test y == x
  masktest_incr2_none1start!(y, x)
  @test y == x .+ ifelse.((axes(x, 1) .≤ 20) .& iseven.(axes(x, 1)), 2, 0)

  # issue #290
  function my_gemm_noturbo!(out, s::Matrix{UInt8}, V, srows, scols, Vcols, μ)
    k = srows >> 2
    rem = srows & 3
    fill!(out, 0)

    for c = 1:Vcols
      for i = 1:scols
        for l = 1:k
          block = s[l, i]
          for p = 1:4
            Aij = (block >> (2 * (p - 1))) & 3
            out[i, c] += (((Aij >= 2) + (Aij == 3) + (Aij == 1) * μ[i]) * V[4*(l-1)+p, c])
          end
        end
      end
    end
    nothing
  end
  function my_gemm!(out, s::Matrix{UInt8}, V, srows, scols, Vcols, μ)
    k = srows >> 2
    rem = srows & 3
    fill!(out, 0)

    @avx for c = 1:Vcols
      for i = 1:scols
        for l = 1:k
          block = s[l, i]
          for p = 1:4
            Aij = (block >> (2 * (p - 1))) & 3
            out[i, c] += (((Aij >= 2) + (Aij == 3) + (Aij == 1) * μ[i]) * V[4*(l-1)+p, c])
          end
        end
      end
    end
    nothing
  end

  out_true = Matrix{Float64}(undef, 100, 100)
  out_test1 = similar(out_true)
  # out_test2 = zeros(100, 100)
  μ = rand(100)
  s = rand(UInt8, 100, 100)
  V = rand(400, 100)

  my_gemm_noturbo!(out_true, s, V, 400, 100, 100, μ)
  my_gemm!(out_test1, s, V, 400, 100, 100, μ)
  @test out_true ≈ out_test1
end


text:
include("euler/euler.jl")
using .Numbers: get_digits
using Combinatorics: with_replacement_combinations
using StatsBase: countmap

function compute(n::Int)::Int
    squares = (0:9) .* (0:9)
    result = 0
    for i ∈ with_replacement_combinations(0:9, n)
        combination = parse(Int, join(i))
        while true
            combination = sum(squares[get_digits(combination) .+ 1])
            if combination == 89
                result += factorial(n) ÷ prod(factorial.(values(countmap(i))))
                break
            elseif combination ≤ 1
                break
            end
        end
    end
    return result
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
module TestWDIJSON

using Test
using WorldBankData
using DataFrames
using JSON

# http://api.worldbank.org/countries/US/indicators/NY.GNP.PCAP.CD?date=1990:1997&format=json
testjsonstr = "[{\"page\":1,\"pages\":1,\"per_page\":\"50\",\"total\":8},[{\"indicator\":{\"id\":\"NY.GNP.PCAP.CD\",\"value\":\"GNI per capita, Atlas method (current US\$)\"},\"country\":{\"id\":\"US\",\"value\":\"United States\"},\"value\":\"31390\",\"decimal\":\"0\",\"date\":\"1997\"},{\"indicator\":{\"id\":\"NY.GNP.PCAP.CD\",\"value\":\"GNI per capita, Atlas method (current US\$)\"},\"country\":{\"id\":\"US\",\"value\":\"United States\"},\"value\":\"30380\",\"decimal\":\"0\",\"date\":\"1996\"},{\"indicator\":{\"id\":\"NY.GNP.PCAP.CD\",\"value\":\"GNI per capita, Atlas method (current US\$)\"},\"country\":{\"id\":\"US\",\"value\":\"United States\"},\"value\":\"29150\",\"decimal\":\"0\",\"date\":\"1995\"},{\"indicator\":{\"id\":\"NY.GNP.PCAP.CD\",\"value\":\"GNI per capita, Atlas method (current US\$)\"},\"country\":{\"id\":\"US\",\"value\":\"United States\"},\"value\":\"27750\",\"decimal\":\"0\",\"date\":\"1994\"},{\"indicator\":{\"id\":\"NY.GNP.PCAP.CD\",\"value\":\"GNI per capita, Atlas method (current US\$)\"},\"country\":{\"id\":\"US\",\"value\":\"United States\"},\"value\":\"26480\",\"decimal\":\"0\",\"date\":\"1993\"},{\"indicator\":{\"id\":\"NY.GNP.PCAP.CD\",\"value\":\"GNI per capita, Atlas method (current US\$)\"},\"country\":{\"id\":\"US\",\"value\":\"United States\"},\"value\":\"25780\",\"decimal\":\"0\",\"date\":\"1992\"},{\"indicator\":{\"id\":\"NY.GNP.PCAP.CD\",\"value\":\"GNI per capita, Atlas method (current US\$)\"},\"country\":{\"id\":\"US\",\"value\":\"United States\"},\"value\":\"24370\",\"decimal\":\"0\",\"date\":\"1991\"},{\"indicator\":{\"id\":\"NY.GNP.PCAP.CD\",\"value\":\"GNI per capita, Atlas method (current US\$)\"},\"country\":{\"id\":\"US\",\"value\":\"United States\"},\"value\":\"24150\",\"decimal\":\"0\",\"date\":\"1990\"}]]"

dfref = DataFrame(iso2c=["US" for _ in range(1, length=8)],
  country=["United States" for _ in range(1, length=8)],
  NY_GNP_PCAP_CD=[31390.0, 30380.0, 29150.0, 27750.0, 26480.0, 25780.0, 24370.0, 24150.0],
  year=[ 1997.0, 1996.0, 1995.0, 1994.0, 1993.0, 1992.0, 1991.0, 1990.0])

#│ Row │ iso2c │ country       │ NY_GNP_PCAP_CD │ year   │
#├─────┼───────┼───────────────┼────────────────┼────────┤
#│ 1   │ US    │ United States │ 31390.0        │ 1997.0 │
#│ 2   │ US    │ United States │ 30380.0        │ 1996.0 │
#│ 3   │ US    │ United States │ 29150.0        │ 1995.0 │
#│ 4   │ US    │ United States │ 27750.0        │ 1994.0 │
#│ 5   │ US    │ United States │ 26480.0        │ 1993.0 │
#│ 6   │ US    │ United States │ 25780.0        │ 1992.0 │
#│ 7   │ US    │ United States │ 24370.0        │ 1991.0 │
#│ 8   │ US    │ United States │ 24150.0        │ 1990.0

testjson = JSON.parse(testjsonstr)[2]

df = WorldBankData.parse_wdi("NY.GNP.PCAP.CD", testjson, 1990, 1997)

@test df == dfref

end


text:
# journal.jl
# ===============
#
# A type representing a journal of an array, string, or sequence-like container,
# against a reference of some other array, string, or sequence-like container.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

type Journal{T}
    "The reference object for the journal."
    reference::T

    "A buffer containing elements required for insertions."
    buffer::T

    "A sorted vector of the journal entries."
    entries::Vector{JE}

    "The length of the journaled object."
    len::Int
end

function Journal{T}(reference::T)
    l = length(reference)
    return Journal{T}(reference, T(), init_entries(l), l)
end

@inline function init_entries(l::Int)
    v = Vector{JE}(1)
    v[1] = JournalEntry(SRC_REFERENCE, 1, 1, 1, l)
    return v
end

function Base.show{T}(io::IO, ::MIME"text/plain", j::Journal{T})
    println(io, "A ", T, " of length ", length(j), " in a journalled format.")
    println(io, "Reference:")
    println(io, reference(j))
end

Base.length(j::Journal) = j.len

@inline function reinit_entries!(j::Journal, l::Int)
    j.entries = init_entries(l)
end

reference(j::Journal) = j.reference
function reference!{T}(j::Journal{T}, ref::T)
    j.reference = ref
    l = length(ref)
    j.len = l
    reinit_entries!(j, l)
end

function clear{T}(j::Journal{T})
    l = length(reference(j))
    reinit_entries!(j, l)
    j.buffer = T()
    j.len = l
end

function insert_record!(tree::Vector{JE}, vpos, bufferpos, len)
    # We find the last entry in the journal with a virtual
    # position that is lower or equal to vpos.
    refEntry = JournalEntry(SRC_NULL, 0, vpos, 0, 0)
    insertpoint = searchsortedlast(tree, refEntry)
    ie = tree[insertpoint]
    ievpos = vposition(ie)

    # Now we need to create new journal entries.

    # If the journal entry found above, contains the virtual position we
    # want to insert into.
    if (ievpos + length(ie)) > vpos
        if ievpos == vpos

        end

    end

end

"""
    insert!(j, ins)

Insert a piece of data into the reference.

This adds an insertion entry into the journal, and updates other entries.
"""
function insert!{T}(j::Journal{T}, ins::T, idx::Int)
    # The length of the journalled data increases.
    l = length(ins)
    j.len += l
    # Note the length of the buffer before insertion.
    bufferstart = length(j.buffer)
    # Add the insertion onto the buffer.
    append!(j.buffer, ins)
    # Insert a record into the tree.
    insert_record!(j.entries, idx, bufferstart, l)
end


text:
module ExampleAside

["a aside note"]

end


text:
mutable struct DDSolution
    solvetime::Float64
    objval::Float64
    best_lower_bound::Float64
    gap::Float64
    numiterations::Int64
    termination::Symbol

    # Iteration Data
    iterval::Array{Float64,1}
    iterbound::Array{Float64,1}
    itertime::Array{Float64,1} # in seconds
    clocktime::Array{Float64,1}

    # DD data
    alpha::Vector{Float64}
    step::Vector{Float64}
end

DDSolution() = DDSolution(0.0,NaN,NaN,NaN,0,:NotExecuted,Float64[],Float64[],Float64[],Float64[],Float64[],Float64[])

#TODO
# function saveiteration(s::Solution,tstamp::Float64,arr::Array{Float64,1},n=1)
#     push!(s.iterval,arr[1])
#     push!(s.iterbound,arr[2])
#     push!(s.itertime,arr[3])
#     push!(s.clocktime,arr[4])
#     if length(arr) == 6 && s.method  == :dual_decomposition
#         push!(s.α,arr[5])
#         push!(s.step,arr[6])
#     end
#     if n == 1
#         s.bestbound = maximum(s.iterbound)
#         s.objval = minimum(s.iterval)
#     else
#         s.bestbound = minimum(s.iterbound)
#         s.objval = maximum(s.iterval)
#     end
#     s.numiterations += 1
#     s.gap = abs(s.objval - s.bestbound)/s.objval
#     s.solvetime = tstamp
# end
#
# function Base.show(io::IO,s::Solution)
#     println("-----------------")
#     println(" SOLUTION SUMMARY")
#     println("-----------------")
#     println("Method: $(s.method)")
#     println("Objective Value : $(s.objval)")
#     println("Best Bound : $(s.bestbound)")
#     println("Gap : $(round(s.gap*100,digits = 2)) %")
#     println("Iterations : $(s.numiterations)")
#     println("Solution Time : $(round(s.solvetime,digits = 2)) s")
#     println("Termination : $(s.termination)")
# end
#
# function printiterationsummary(s::Solution; singleline=false)
#     if singleline
#         print("$(s.numiterations)  ")
#         print("| $(round(s.iterval[end],digits = 4)) \t")
#         print("$(round(s.iterbound[end],digits = 4)) \t")
#         print("| $(round(s.objval,digits = 4)) \t")
#         print("$(round(s.bestbound,digits = 4)) \t")
#         print("| $(round(s.gap*100,digits = 2)) %\t")
#         print("| $(round(s.clocktime[end],digits = 2)) s\n")
#     else
#         println("-----------------")
#         println(" Iteration $(s.numiterations)")
#         println("-----------------")
#         println("Iteration Value : $(s.iterval[end])")
#         println("Iteration Bound : $(s.iterbound[end])")
#         println("Objective Value : $(s.objval)")
#         println("Best Bound : $(s.bestbound)")
#         println("Gap : $(round(s.gap*100,digits = 2)) %")
#         println("Iteration Time : $(round(s.itertime[end],digits = 2)) s")
#         println("Elapsed Time : $(round(s.clocktime[end],digits = 2)) s")
#     end
# end
#
# function iterationheader()
#     print("    |    Iteration\t")
#     print("|     Best Bounds\t|\n")
#     print("It  ")
#     print("| Value \t")
#     print("Bound \t")
#     print("| Value \t")
#     print("Bound \t")
#     print("| Gap\t")
#     print("| Time\n")
# end


text:
# This file is a part of RadiationDetectorSignals.jl, licensed under the MIT License (MIT).

# TODO: Replace by true custom types to avoid type piracy (e.g. with plotting
# recipes) on NamedTuple and TypedTables.Table (even though the types aliased
# here are quite specific, so the risk of a conflict is low).

"""
    DetectorHit = NamedTuple{(:evtno, :detno, :thit, :edep, :pos), ...)

Representation of an localized energy deposition in a detector.

Use [`DetectorHitEvents`](@ref) for arrays of `DetectorHit` that have a
compact memory layout.
"""
const DetectorHit = NamedTuple{
    (:evtno, :detno, :thit, :edep, :pos),
    <:Tuple{Integer, Integer, RealQuantity, RealQuantity, AbstractVector{<:RealQuantity}}
}
export DetectorHit

const DetectorHits = TypedTables.Table{<:DetectorHit,1}
export DetectorHits
DetectorHits(x::NamedTuple{(:evtno, :detno, :thit, :edep, :pos)}) = TypedTables.Table(x)::DetectorHits
DetectorHits(x::Any) = TypedTables.Table(detno = x.detno, thit = x.thit, edep = x.edep, pos = x.pos)::DetectorHits

"""
    DetectorHitEvents = TypedTables.Table{...}

A `TypedTables.Table` of [`DetectorHit`](@ref) events.
"""
const DetectorHitEvents = TypedTables.Table{
    <:NamedTuple{
        (:evtno, :detno, :thit, :edep, :pos),
        <:Tuple{
            Union{Integer, AbstractVector{<:Integer}},
            Union{Integer, AbstractVector{<:Integer}},
            AbstractVector{<:RealQuantity},
            AbstractVector{<:RealQuantity},
            AbstractVector{<:AbstractVector{<:RealQuantity}}
        }
    }
}

export DetectorHitEvents

DetectorHitEvents(x::NamedTuple{(:evtno, :detno, :thit, :edep, :pos)}) = TypedTables.Table(x)::DetectorHitEvents
DetectorHitEvents(x::Any) = TypedTables.Table(evtno = x.evtno, detno = x.detno, thit = x.thit, edep = x.edep, pos = x.pos)::DetectorHitEvents

function _to_scalar(A::AbstractArray)
    x = first(A)
    !all(isequal(x), A) && throw(ArgumentError("Array values are not all equal, can't convert to scalar"))
    x
end

function group_by_evtno(hits::DetectorHits)
    issorted(hits.evtno) || throw(ArgumentError("Hits are not sorted by evtno"))
    grouped_cols = consgroupedview(hits.evtno, Tables.columns(hits))
    Table(merge(grouped_cols, (evtno = _to_scalar.(grouped_cols.evtno),)))
end

function ungroup_by_evtno(events::DetectorHitEvents)
    expanded_evtno = ((evtno, edep) -> Fill(evtno, size(edep))).(events.evtno, events.edep)
    Table(merge(
        map(flatview, Tables.columns(events)),
        (evtno = collect(flatview(expanded_evtno)),)
    ))
end

function group_by_evtno_and_detno(hits::DetectorHits)
    k = ((a,b) -> (b,a)).(hits.evtno, hits.detno)
    perm = sortperm(k)
    k_sorted = k[perm]
    sorted_hits = hits[perm]
    grouped_cols = consgroupedview(k_sorted, Tables.columns(sorted_hits))
    Table(merge(
        grouped_cols,
        (
            evtno = _to_scalar.(grouped_cols.evtno),
            detno = _to_scalar.(grouped_cols.detno),
        )
    ))
end


text:
module brexit_investment

using Revise

using Optim, Parameters, QuantEcon, BenchmarkTools,Dates,JLD,Interpolations
using .Threads #so we don't have to write Threads.@threads every time

#Type definitions and their functions (brexDefs.jl)
export pars,stat_equil,check_par,fillV

#Functions from brexComp.jl
export SE_compute!,ED

#Export functions from brexTools.jl
export genFolderName,prep0,saveAll,saveInFile,loadFromFile

### Include source files ###
include("brexDefs.jl")
include("brexComp.jl")
include("brexTools.jl")

end # module


text:
struct Resource
    relative_package_path::Union{Nothing, Vector{String}}
    dev_package_path::Union{Nothing, Vector{String}}
    external_url::Union{Nothing, Vector{String}}
    type::Symbol
    async::Symbol # :none, :eager, :lazy May be we should use enum
    function Resource(;relative_package_path, dev_package_path = nothing, external_url = nothing, type = :js, dynamic = nothing, async=nothing)
        (!isnothing(dynamic) && !isnothing(async)) && throw(ArgumentError("Can't have both 'dynamic' and 'async'"))
        !in(type, [:js, :css]) &&  throw(ArgumentError("type must be `:js` or `:css`"))
        async_symbol = :none
        if !isnothing(dynamic)
            dynamic == true && (async_symbol = :lazy)
        elseif !isnothing(async) && async != :false
            async_symbol = async == :lazy ? :lazy : :eager
        end
        return new(_path_to_vector(relative_package_path), _path_to_vector(dev_package_path), _path_to_vector(external_url), type, async_symbol)
    end
end

_path_to_vector(s::Nothing) = nothing
_path_to_vector(s::String) = [s]
_path_to_vector(s::Vector{String}) = s

has_relative_path(r::Resource) = !isnothing(r.relative_package_path)
has_dev_path(r::Resource) = !isnothing(r.dev_package_path)
has_external_url(r::Resource) = !isnothing(r.external_url)

get_type(r::Resource) = r.type
get_external_url(r::Resource) = r.external_url
get_dev_path(r::Resource) = r.dev_package_path
get_relative_path(r::Resource) = r.relative_package_path

isdynamic(resource::Resource, eager_loading::Bool) = resource.async == :lazy || (resource.async == :eager && !eager_loading)

struct ResourcePkg
    namespace ::String
    path ::String
    resources ::Vector{Resource}
    version ::String
    ResourcePkg(namespace, path, resources = Resource[]; version = "")  = new(namespace, path, resources, version)
end

struct ResourcesRegistry
    dash_dependency ::NamedTuple{(:dev, :prod), Tuple{ResourcePkg,ResourcePkg}}
    dash_renderer ::ResourcePkg
    components ::Dict{String, ResourcePkg}
    ResourcesRegistry(;dash_dependency, dash_renderer) = new(dash_dependency,  dash_renderer, Dict{String, ResourcePkg}())
end

function register_package!(registry::ResourcesRegistry, pkg::ResourcePkg)
    registry.components[pkg.namespace] = pkg
end

get_dash_dependencies(registry::ResourcesRegistry, prop_check::Bool) = prop_check ?
                                                                    registry.dash_dependency[:dev] :
                                                                    registry.dash_dependency[:prod]

get_componens_pkgs(registry::ResourcesRegistry) = values(registry.components)
get_dash_renderer_pkg(registry::ResourcesRegistry) = registry.dash_renderer

const RESOURCE_PATH = realpath(joinpath(ROOT_PATH, "resources"))

const resources_registry = ResourcesRegistry(
    dash_dependency = (
        dev = ResourcePkg(
            "dash_renderer",
            RESOURCE_PATH, version = "1.5.0",
            [
                Resource(
                relative_package_path = "react@16.13.0.js",
                external_url = "https://unpkg.com/react@16.13.0/umd/react.development.js"
                ),
                Resource(
                    relative_package_path = "react-dom@16.13.0.js",
                    external_url = "https://unpkg.com/react-dom@16.13.0/umd/react-dom.development.js"
                ),
                Resource(
                    relative_package_path = "polyfill@7.8.7.min.js",
                    external_url = "https://unpkg.com/@babel/polyfill@7.8.7/dist/polyfill.min.js"
                ),
                Resource(
                    relative_package_path = "prop-types@15.7.2.js",
                    external_url = "https://unpkg.com/prop-types@15.7.2/prop-types.js",
                ),
            ]
        ),
        prod = ResourcePkg(
            "dash_renderer",
            RESOURCE_PATH, version = "1.2.2",
            [
                Resource(
                relative_package_path = "react@16.13.0.min.js",
                external_url = "https://unpkg.com/react@16.13.0/umd/react.production.min.js"
                ),
                Resource(
                    relative_package_path = "react-dom@16.13.0.min.js",
                    external_url = "https://unpkg.com/react-dom@16.13.0/umd/react-dom.production.min.js"
                ),
                Resource(
                    relative_package_path = "polyfill@7.8.7.min.js",
                    external_url = "https://unpkg.com/@babel/polyfill@7.8.7/dist/polyfill.min.js"
                ),
                Resource(
                    relative_package_path = "prop-types@15.7.2.min.js",
                    external_url = "https://unpkg.com/prop-types@15.7.2/prop-types.min.js"
                ),
            ]
        )
    ),
    dash_renderer = ResourcePkg(
        "dash_renderer",
        RESOURCE_PATH, version = "1.5.0",
        [
            Resource(
                relative_package_path = "dash_renderer.min.js",
                dev_package_path = "dash_renderer.dev.js",
                external_url = "https://unpkg.com/dash-renderer@1.5.0/dash_renderer/dash_renderer.min.js"
            ),
            Resource(
                relative_package_path = "dash_renderer.min.js.map",
                dev_package_path = "dash_renderer.dev.js.map",
                dynamic = true,
            ),
        ]
    )
)

register_package(pkg::ResourcePkg) = register_package!(resources_registry, pkg)

main_registry() = resources_registry


text:
__precompile__()
module Compounds

export Compound,getConstPropList,newCompound,
GetCompoundList,getConstPropList,property

StringPropertiesBaseDict = begin
    Dict{String,Union{Missing,String}}(
    "ID"=>missing,
    "commonName"=>missing,
    "iupacName"=>missing,
    "casRegistryNumber"=>missing,
    "chemicalFormula"=>missing,
    "SMILESformula"=>missing,
    "compoundType"=>"Standard") #For future solid, ionic, bacterial and other interactions
end
ValuePropertiesBaseDict = begin
    Dict{String,Float64}(
    "molecularWeight"=>NaN,
    "criticalTemperature"=>NaN,
    "criticalPressure"=>NaN,
    "criticalVolume"=>NaN,
    "criticalCompressibilityFactor"=>NaN,
    "criticalDensity"=>NaN,
    "acentricFactor"=>NaN,
    "dipoleMoment"=>NaN,
    "parachor"=>NaN,
    "gyrationRadius"=>NaN,
    "associationParameter"=>NaN,
    "diffusionVolume"=>NaN,
    "vanderwaalsArea"=>NaN,
    "energyLennardJones"=>NaN,
    "lengthLennardJones"=>NaN,
    "normalBoilingPoint"=>NaN,
    "heatOfVaporizationAtNormalBoilingPoint"=>NaN,
    "normalFreezingPoint"=>NaN,
    "heatOfFusionAtNormalFreezingPoint"=>NaN,
    "liquidDensityAt25C"=>NaN,
    "liquidVolumeAt25C"=>NaN,
    "idealGasGibbsFreeEnergyOfFormationAt25C"=>NaN,
    "idealGasEnthalpyOfFormationAt25C"=>NaN,
    "standardFormationEnthalpySolid"=>NaN,
    "standardFormationEnthalpyLiquid"=>NaN,
    "standardFormationEnthalpyGas"=>NaN,
    "standardFormationGibbsEnergySolid"=>NaN,
    "standardFormationGibbsEnergyLiquid"=>NaN,
    "standardFormationGibbsEnergyGas"=>NaN,
    "standardEntropySolid"=>NaN,
    "standardEntropyLiquid"=>NaN,
    "standardEntropyGas"=>NaN,
    "triplePointTemperature"=>NaN,
    "triplePointPressure"=>NaN,
    "BornRadius"=>NaN,
    "charge"=>NaN,
    "StandardEnthalpyAqueousDilution"=>NaN,
    "StandardGibbsAqueousDilution"=>NaN)
end

ArrayPropertiesBaseDict = begin
    Dict{String,Dict{String,Float64}}(
    "DebugConstants"=>Dict{String,Float64}(
    "empty"=>NaN)
)
end
###########################################
struct Compound
    #stores an individual compound, and its constant properties
    StringProperties::Dict{String,Union{Missing,String}}
    ValueProperties::Dict{String,Float64}
    ArrayProperties::Dict{String,Dict{String,Float64}}
    function Compound(str=StringPropertiesBaseDict,val=ValuePropertiesBaseDict,arr=ArrayPropertiesBaseDict;ID="")
    x1 = merge(StringPropertiesBaseDict,str)
    x2 = merge(ValuePropertiesBaseDict,val)
    x3= merge(ArrayPropertiesBaseDict,arr)
    if ID!=""
        x1["ID"]=ID
    end
    return new(x1,x2,x3)
end
end

struct CompoundList
    #stores a set of compounds, and its constant properties
    Compounds::Array{Compound,1}
    CompoundsID::Array{String,1}
    function CompoundList(Compoundsz::Array{Compound})
        IDs = map(x->x.StringProperties["ID"],Compoundsz)
        if allunique(IDs) || throw("one or more IDs of the compound list are diferent.")
            return new(Compoundsz,IDs)
        else
       end
    end
    #strArr:: Array{String,1}
end

########################
#Indexing to CompoundList
Base.length(S::CompoundList) = length(S.CompoundsID)
Base.iterate(S::CompoundList, state=1) = state > length(S) ? nothing : (S.compounds[state], state+1)

function Base.getindex(S::CompoundList, i::Number)
    1 <= i <= length(S) || throw(BoundsError(S, i))
    return S.compounds[i]
end

#function Base.getindex(S::CompoundList, I)
 #   return [S[i] for i in I]
#end

function Base.getindex(list::CompoundList, prop::String)
    if haskey(ValuePropertiesBaseDict, prop)
        xxx = map(x->x.ValueProperties[prop],list.Compounds)
    elseif haskey(StringPropertiesBaseDict, prop)
        xxx= map(x->x.StringProperties[prop],list.Compounds)
    elseif haskey(ArrayPropertiesBaseDict, prop)
        xxx= map(x->x.ArrayProperties[prop],list.Compounds)
    else throw("Property doesnt exist.")
    end
    return xxx
end

Base.firstindex(S::CompoundList)=1
Base.lastindex(S::CompoundList)= length(S)

Base.iterate(S::CompoundList, state=1) = state > length(S) ? nothing : (S.compounds[state], state+1)

function getCompound(list::CompoundList)
return(list.Compounds)
end

function getCompound(list::CompoundList,ID::String="")
    if ID==""
    return(list.Compounds)
    else
        list.Compounds(ID)
    end
#SetCompoundProperty()
end

#function garr(Compoundsz::Array{Compound}) #transform compound list to Array

function newCompoundList(Compoundsz::Array{Compound})
    IDs = map(x->x.StringProperties["ID"],Compoundsz)
    if allunique(IDs) || throw("one or more IDs of the compound list are diferent.")
        return CompoundList(Compoundsz,IDs)
    else

   end
end

#function GetCompoundConstant(;props::Array{String};compIds::Array{String}=[])

function property(list::CompoundList,prop::String)
if haskey(ValuePropertiesBaseDict, prop)
    xxx = map(x->x.ValueProperties[prop],list.Compounds)
elseif haskey(StringPropertiesBaseDict, prop)
    xxx= map(x->x.StringProperties[prop],list.Compounds)
elseif haskey(ArrayPropertiesBaseDict, prop)
    xxx= map(x->x.ArrayProperties[prop],list.Compounds)
else throw("Property doesnt exist.")
end
return xxx
end

#Cape OPEN functions
function GetCompoundList(list::CompoundList)
return[
    property(list,"ID"),
    property(list,"chemicalFormula"),
    property(list,"commonName"),
    property(list,"normalBoilingPoint"),
    property(list,"molecularWeight"),
    property(list,"casRegistryNumber"),
    ]
end

function GetNumCompounds(list::CompoundList)
    return length(list)
end

getConstPropList() = vcat(collect(keys(StringPropertiesBaseDict)),collect(keys(ValuePropertiesBaseDict)))

## IMPLEMENT THE FOLLOWING ON ALL OBJECTS (examples from julialang docs)
#Base.iterate(S::Squares, state=1) = state > S.count ? nothing : (state*state, state+1)
#Base.eltype(::Type{Squares}) = Int # Note that this is defined for the type

#END OF MODULE
end


text:
"""
**Basic gradient descent with fixed learning rate**
```julia
VanillaGradDescent(; η::Real=0.01)
```

Algorithm
```math
\\Delta x_t = \\eta g_t
```

"""
mutable struct VanillaGradDescent <: Optimizer
    opt_type::String
    t::Int64
    η::Float64
end

function VanillaGradDescent(; η::Real=0.01)
    @assert η > 0.0 "η must be greater than 0"

    VanillaGradDescent("Vanilla Gradient Descent", 0, η)
end

params(opt::VanillaGradDescent) = "η=$(opt.η)"

function update(opt::VanillaGradDescent, g_t::AbstractArray{T}) where {T<:Real}
    # update timestep
    opt.t += 1
    return opt.η * g_t
end


text:
"""
Reptile: A Scalable Meta-Learning Algorithm [1]

[1] https://blog.openai.com/reptile/
"""
function reptile(model; meta_opt=Descent(0.1), inner_opt=Descent(0.02), epochs=30_000,
                 train_batch_size=10, eval_batch_size=10, eval_interval=1000)

    weights = params(model)
    dist = Uniform(-5, 5)
    testx = Float32.(range(-5, stop=5, length=50))
    x = testx

    for i in 1:epochs
        prev_weights = deepcopy(Flux.data.(weights))
        task = SineWave()

        # Train on task for k steps on the dataset
        y = task(x)
        for idx in partition(randperm(length(x)), train_batch_size)
            l = Flux.mse(model(x[idx]'), y[idx]')
            Flux.back!(l)
            Flux.Optimise._update_params!(inner_opt, weights)
        end

        # Reptile update
        for (w1, w2) in zip(weights, prev_weights)
            gw = Flux.Optimise.apply!(meta_opt, w2, w1.data - w2)
            @. w1.data = w2 + gw
        end

        if i % eval_interval == 0
            @printf("Iteration %d, evaluating model on random task...\n", i)
            evalx = Float32.(rand(dist, eval_batch_size))
            eval_model(model, evalx, testx, SineWave())
        end

    end
end

max_stars_repo_name: JuliaTagBot/MetaLearning

max_issues_repo_name: JuliaTagBot/MetaLearning

max_forks_repo_name: JuliaTagBot/MetaLearning