text:
using Test
using NQCDynamics
using OrdinaryDiffEq
using StaticArrays
using LinearAlgebra: Eigen
using Random: seed!
using NQCDynamics: DynamicsMethods, DynamicsUtils, Calculators
using NQCDynamics.SurfaceHoppingMethods: SurfaceHoppingMethods
using NQCDynamics.DynamicsUtils: get_positions, get_velocities

@test SurfaceHoppingMethods.FSSH{Float64}(2, :standard) isa SurfaceHoppingMethods.FSSH
atoms = Atoms(2)

@testset "FSSH" begin
    sim = Simulation{FSSH}(atoms, NQCModels.DoubleWell())
    sim.method.state = 1

    r = zeros(size(sim))
    v = randn(size(sim))
    u = DynamicsVariables(sim, v, r, SingleState(1, Adiabatic()))
    du = zero(u)

    @test DynamicsUtils.get_quantum_subsystem(u) ≈ Complex.([1 0; 0 0])

    DynamicsMethods.motion!(du, u, sim, 0.0)

    Calculators.evaluate_nonadiabatic_coupling!(sim.calculator)
    @test sim.calculator.nonadiabatic_coupling ≈ -sim.calculator.nonadiabatic_coupling'

    problem = ODEProblem(DynamicsMethods.motion!, u, (0.0, 1.0), sim)
    integrator = init(problem, Tsit5(), callback=SurfaceHoppingMethods.HoppingCallback)
    DynamicsMethods.SurfaceHoppingMethods.evaluate_hopping_probability!(sim, u, get_proposed_dt(integrator))

    σ = DynamicsUtils.get_quantum_subsystem(u)
    dσ = zero(σ)
    DynamicsUtils.set_quantum_derivative!(dσ, v, σ, sim)

    @testset "get_diabatic_population" begin
        population = Estimators.diabatic_population(sim, u)
        @test population ≈ [0.5, 0.5]
    end

    @testset "select_new_state" begin
        sim.method.hopping_probability .= [0, 1.0]
        u.state = 1
        @test 2 == SurfaceHoppingMethods.select_new_state(sim, u)
        sim.method.hopping_probability .= [0, 1.0]
        u.state = 2
        @test 2 == SurfaceHoppingMethods.select_new_state(sim, u)
        sim.method.hopping_probability .= [1.0, 0.0]
        u.state = 2
        @test 1 == SurfaceHoppingMethods.select_new_state(sim, u)
    end

    @testset "rescale_velocity!" begin
        get_velocities(integrator.u) .= 0.0 # Set momentum to zero to force frustrated hop
        integrator.p.method.new_state = 2
        cont = SurfaceHoppingMethods.rescale_velocity!(integrator.p, integrator.u)
        @test cont == false # Check hop is rejected

        get_velocities(integrator.u) .= 1e5 # Set momentum to big to force hop
        cont = SurfaceHoppingMethods.rescale_velocity!(integrator.p, integrator.u)
        @test cont == true # Check hop is accepted
    end

    @testset "calculate_potential_energy_change" begin
        integrator.p.calculator.eigen = Eigen(SVector{2}([0.9, -0.3]), SMatrix{2,2}(1, 0, 0.0, 1))
        eigs = SurfaceHoppingMethods.get_hopping_eigenvalues(integrator.p)
        @test 1.2 ≈ SurfaceHoppingMethods.calculate_potential_energy_change(eigs, 1, 2)
        @test -1.2 ≈ SurfaceHoppingMethods.calculate_potential_energy_change(eigs, 2, 1)
    end

    @testset "execute_hop!" begin
        Calculators.update_electronics!(integrator.p.calculator, get_positions(integrator.u))
        DynamicsUtils.get_velocities(integrator.u) .= 2 # Set high momentum to ensure successful hop
        integrator.u.state = 1
        integrator.p.method.new_state = 2
        eigs = SurfaceHoppingMethods.get_hopping_eigenvalues(integrator.p)
        ΔE = SurfaceHoppingMethods.calculate_potential_energy_change(eigs, 2, 1)

        KE_initial = DynamicsUtils.classical_kinetic_energy(integrator.p, get_velocities(integrator.u))
        H_initial = DynamicsUtils.classical_hamiltonian(integrator.p, integrator.u)
        @test integrator.u.state == 1

        SurfaceHoppingMethods.execute_hop!(integrator)

        KE_final = DynamicsUtils.classical_kinetic_energy(integrator.p, get_velocities(integrator.u))
        H_final = DynamicsUtils.classical_hamiltonian(integrator.p, integrator.u)
        ΔKE = KE_final - KE_initial

        @test integrator.u.state == 2 # Check state has changed
        @test H_final ≈ H_initial
        @test ΔKE ≈ -ΔE rtol=1e-3 # Test for energy conservation
    end

    @testset "run_trajectory" begin
        atoms = Atoms(2000)
        sim = Simulation{FSSH}(atoms, NQCModels.TullyModelTwo())
        v = hcat(100 / 2000)
        r = hcat(-10.0)
        u = DynamicsVariables(sim, v, r, SingleState(1, Adiabatic()))
        solution = run_trajectory(u, (0.0, 500.0), sim, output=(:hamiltonian, :state), reltol=1e-6)
        @test solution.hamiltonian[1] ≈ solution.hamiltonian[end] rtol=1e-2
    end
end

@testset "RPSH" begin
    sim = RingPolymerSimulation{FSSH}(atoms, NQCModels.DoubleWell(), 5)

    r = RingPolymerArray(zeros(size(sim)))
    v = RingPolymerArray(randn(size(sim)))
    u = DynamicsVariables(sim, v, r, SingleState(1, Adiabatic()))
    sim.method.state = u.state

    problem = ODEProblem(DynamicsMethods.motion!, u, (0.0, 1.0), sim)
    integrator = init(problem, Tsit5(), callback=SurfaceHoppingMethods.HoppingCallback)
    SurfaceHoppingMethods.evaluate_hopping_probability!(sim, u, get_proposed_dt(integrator))

    @testset "rescale_velocity!" begin
        get_velocities(integrator.u) .= 0.0 # Set momentum to zero to force frustrated hop
        integrator.p.method.new_state = 2
        cont = SurfaceHoppingMethods.rescale_velocity!(integrator.p, integrator.u)
        @test cont == false # Check hop is rejected

        get_velocities(integrator.u) .= 1e5 # Set momentum to big to force hop
        cont = SurfaceHoppingMethods.rescale_velocity!(integrator.p, integrator.u)
        @test cont == true # Check hop is accepted
    end

    @testset "calculate_potential_energy_change" begin
        integrator.p.calculator.centroid_eigen = Eigen(SVector{2}([0.9, -0.3]), SMatrix{2,2}(1, 0, 0.0, 1))
        eigs = SurfaceHoppingMethods.get_hopping_eigenvalues(integrator.p)
        @test 1.2 ≈ SurfaceHoppingMethods.calculate_potential_energy_change(eigs, 1, 2)
        @test -1.2 ≈ SurfaceHoppingMethods.calculate_potential_energy_change(eigs, 2, 1)
    end

    @testset "execute_hop!" begin
        Calculators.update_electronics!(integrator.p.calculator, get_positions(integrator.u))
        get_velocities(integrator.u) .= 5 # Set high momentum to ensure successful hop
        integrator.u.state = 1
        integrator.p.method.new_state = 2
        eigs = SurfaceHoppingMethods.get_hopping_eigenvalues(integrator.p)
        ΔE = SurfaceHoppingMethods.calculate_potential_energy_change(eigs, 2, 1)

        KE_initial = DynamicsUtils.classical_kinetic_energy(integrator.p, get_velocities(integrator.u))
        H_initial = DynamicsUtils.classical_hamiltonian(integrator.p, integrator.u)
        potential = H_initial - KE_initial
        @test integrator.u.state == 1

        SurfaceHoppingMethods.execute_hop!(integrator)

        KE_final = DynamicsUtils.classical_kinetic_energy(integrator.p, get_velocities(integrator.u))
        H_final = DynamicsUtils.classical_hamiltonian(integrator.p, integrator.u)
        potential = H_final - KE_final

        @test integrator.u.state == 2 # Check state has changed
        @test H_initial ≈ H_final # Check ring polymer hamiltonian conserved
        ΔKE = (KE_final - KE_initial) / length(sim.beads)
        @test ΔKE ≈ -sum(ΔE) rtol=1e-3 # Test for energy conservation
    end

    @testset "run_trajectory" begin
        atoms = Atoms(2000)
        sim = RingPolymerSimulation{FSSH}(atoms, TullyModelTwo(), 5; temperature=0.01)
        v = fill(100 / 2000, size(sim))
        r = fill(-10.0, size(sim)) .+ randn(1,1,5)
        u = DynamicsVariables(sim, v, r, SingleState(1, Adiabatic()))
        seed!(1)
        solution = run_trajectory(u, (0.0, 1000.0), sim, output=(:hamiltonian), dt=0.1)
        seed!(1)
        solution2 = run_trajectory(u, (0.0, 1000.0), sim, output=(:hamiltonian), algorithm=Tsit5(), abstol=1e-10, reltol=1e-10, saveat=0:0.1:1000.0)
        # Ring polymer Hamiltonian is not strictly conserved during hoppping
        @test solution.hamiltonian[1] ≈ solution.hamiltonian[end] rtol=1e-2
    end

end


text:
function initialize_avr!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::PSY.DynamicGenerator{M, S, PSY.AVRFixed, TG, P},
) where {M <: PSY.Machine, S <: PSY.Shaft, TG <: PSY.TurbineGov, P <: PSY.PSS}
    #In AVRFixed, V_ref is used as Vf
    Vf = get_inner_vars(dynamic_device)[Vf_var]
    PSY.set_V_ref!(PSY.get_avr(dynamic_device), Vf)
    PSY.set_Vf!(PSY.get_avr(dynamic_device), Vf)
    #Update Control Refs
    PSY.get_ext(dynamic_device)[CONTROL_REFS][V_ref_index] = Vf
end

function initialize_avr!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::PSY.DynamicGenerator{M, S, PSY.AVRSimple, TG, P},
) where {M <: PSY.Machine, S <: PSY.Shaft, TG <: PSY.TurbineGov, P <: PSY.PSS}
    #In AVRFixed, V_ref is used as Vf
    Vf0 = get_inner_vars(dynamic_device)[Vf_var]
    #Obtain measured terminal voltage
    Vm = sqrt(
        get_inner_vars(dynamic_device)[VR_gen_var]^2 +
        get_inner_vars(dynamic_device)[VI_gen_var]^2,
    )
    #Solve V_ref
    V_ref = Vm
    #Set Vf state equals to Vf0
    avr_ix = get_local_state_ix(dynamic_device, PSY.AVRSimple)
    avr_states = @view device_states[avr_ix]
    avr_states[1] = Vf0 #δ
    #Set V_ref
    PSY.set_V_ref!(PSY.get_avr(dynamic_device), V_ref)
    #Update Control Refs
    PSY.get_ext(dynamic_device)[CONTROL_REFS][V_ref_index] = V_ref
end

function initialize_avr!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::PSY.DynamicGenerator{M, S, PSY.AVRTypeI, TG, P},
) where {M <: PSY.Machine, S <: PSY.Shaft, TG <: PSY.TurbineGov, P <: PSY.PSS}
    #Obtain Vf0 solved from Machine
    Vf0 = get_inner_vars(dynamic_device)[Vf_var]
    #Obtain measured terminal voltage
    Vm = sqrt(
        get_inner_vars(dynamic_device)[VR_gen_var]^2 +
        get_inner_vars(dynamic_device)[VI_gen_var]^2,
    )

    #Get parameters
    avr = PSY.get_avr(dynamic_device)
    Ka = PSY.get_Ka(avr)
    Ke = PSY.get_Ke(avr)
    Kf = PSY.get_Kf(avr)
    Tf = PSY.get_Tf(avr)
    Ae = PSY.get_Ae(avr)
    Be = PSY.get_Be(avr)
    #Obtain saturated Vf
    Se_Vf = Ae * exp(Be * abs(Vf0))

    #States of AVRTypeI are Vf, Vr1, Vr2, Vm
    #To solve V_ref, Vr1, Vr2
    function f!(out, x)
        V_ref = x[1]
        Vr1 = x[2]
        Vr2 = x[3]

        out[1] = Vf0 * (Ke + Se_Vf) - Vr1 #16.12c
        out[2] = Ka * (V_ref - Vm - Vr2 - (Kf / Tf) * Vf0) - Vr1 #16.12a
        out[3] = (Kf / Tf) * Vf0 + Vr2 #16.12b
    end
    x0 = [1.0, Vf0, Vf0]
    sol = NLsolve.nlsolve(f!, x0)
    if !NLsolve.converged(sol)
        @warn("Initialization in AVR failed")
    else
        sol_x0 = sol.zero
        #Update V_ref
        PSY.set_V_ref!(avr, sol_x0[1])
        PSY.get_ext(dynamic_device)[CONTROL_REFS][V_ref_index] = sol_x0[1]
        #Update AVR states
        avr_ix = get_local_state_ix(dynamic_device, PSY.AVRTypeI)
        avr_states = @view device_states[avr_ix]
        avr_states[1] = Vf0 #δ
        avr_states[2] = sol_x0[2] #ω
        avr_states[3] = sol_x0[3]
        avr_states[4] = Vm
    end
end

function initialize_avr!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::PSY.DynamicGenerator{M, S, PSY.AVRTypeII, TG, P},
) where {M <: PSY.Machine, S <: PSY.Shaft, TG <: PSY.TurbineGov, P <: PSY.PSS}
    #Obtain Vf0 solved from Machine
    Vf0 = get_inner_vars(dynamic_device)[Vf_var]
    #Obtain measured terminal voltage
    Vm = sqrt(
        get_inner_vars(dynamic_device)[VR_gen_var]^2 +
        get_inner_vars(dynamic_device)[VI_gen_var]^2,
    )

    #Get parameters
    avr = PSY.get_avr(dynamic_device)
    K0 = PSY.get_K0(avr)
    T1 = PSY.get_T1(avr)
    T2 = PSY.get_T2(avr)
    T3 = PSY.get_T3(avr)
    T4 = PSY.get_T4(avr)
    Tr = PSY.get_Tr(avr)
    Ae = PSY.get_Ae(avr)
    Be = PSY.get_Be(avr)
    #Obtain saturated Vf
    Se_Vf = Ae * (exp(Be * abs(Vf0)))

    #States of AVRTypeII are Vf, Vr1, Vr2, Vm
    #To solve V_ref, Vr1, Vr2
    function f!(out, x)
        V_ref = x[1]
        Vr1 = x[2]
        Vr2 = x[3]

        Vr = K0 * Vr2 + (T4 / T3) * (Vr1 + K0 * (T2 / T1) * (V_ref - Vm)) #16.21

        out[1] = (Vf0 * (1.0 + Se_Vf) - Vr) #16.18
        out[2] = K0 * (1.0 - (T2 / T1)) * (V_ref - Vm) - Vr1 #16.14
        out[3] = (1.0 - (T4 / T3)) * (Vr1 + K0 * (T2 / T1) * (V_ref - Vm)) - K0 * Vr2  #16.20
    end
    x0 = [1.0, Vf0, Vf0]
    sol = NLsolve.nlsolve(f!, x0)
    if !NLsolve.converged(sol)
        @warn("Initialization in AVR failed")
    else
        sol_x0 = sol.zero
        #Update V_ref
        PSY.set_V_ref!(avr, sol_x0[1])
        PSY.get_ext(dynamic_device)[CONTROL_REFS][V_ref_index] = sol_x0[1]
        #Update AVR states
        avr_ix = get_local_state_ix(dynamic_device, PSY.AVRTypeII)
        avr_states = @view device_states[avr_ix]
        avr_states[1] = Vf0 #Vf
        avr_states[2] = sol_x0[2] #Vr1
        avr_states[3] = sol_x0[3] #Vr2
        avr_states[4] = Vm #Vm
    end
end

function initialize_avr!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::PSY.DynamicGenerator{M, S, PSY.ESAC1A, TG, P},
) where {M <: PSY.Machine, S <: PSY.Shaft, TG <: PSY.TurbineGov, P <: PSY.PSS}
    #Obtain Vf0 solved from Machine
    Vf0 = get_inner_vars(dynamic_device)[Vf_var]
    #Obtain I_fd obtained from Machine:
    Xad_Ifd0 = get_inner_vars(dynamic_device)[Xad_Ifd_var]
    #Obtain measured terminal voltage
    Vm0 = sqrt(
        get_inner_vars(dynamic_device)[VR_gen_var]^2 +
        get_inner_vars(dynamic_device)[VI_gen_var]^2,
    )

    #Get parameters
    avr = PSY.get_avr(dynamic_device)
    Tr = PSY.get_Tr(avr)
    Tb = PSY.get_Tb(avr)
    Tc = PSY.get_Tc(avr)
    Ka = PSY.get_Ka(avr)
    Ta = PSY.get_Ta(avr)
    Va_min, Va_max = PSY.get_Va_lim(avr)
    Te = PSY.get_Te(avr)
    Kf = PSY.get_Kf(avr)
    Tf = PSY.get_Tf(avr)
    Kc = PSY.get_Kc(avr)
    Kd = PSY.get_Kd(avr)
    Ke = PSY.get_Ke(avr)
    E1, E2 = PSY.get_E_sat(avr)
    SE1, SE2 = PSY.get_Se(avr)
    Vr_min, Vr_max = PSY.get_Vr_lim(avr)
    #Obtain saturation
    #Se_Vf = saturation_function(Vm)

    #Solve Ve from rectifier function
    function f_Ve!(out, x)
        V_e0 = x[1]
        I_N0 = Kc * Xad_Ifd0 / V_e0

        out[1] = Vf0 - V_e0 * rectifier_function(I_N0)
    end
    x0 = [1.0]
    sol = NLsolve.nlsolve(f_Ve!, x0)
    if !NLsolve.converged(sol)
        @error("Initialization in V_e0 failed")
    else
        sol_x0 = sol.zero
        V_e0 = sol_x0[1]
        I_N0 = Kc * Xad_Ifd0 / V_e0
    end
    Se0 = saturation_function(avr, V_e0) #To do saturation function
    V_FE0 = Kd * Xad_Ifd0 + Ke * V_e0 + Se0 * V_e0
    V_r20 = V_FE0
    if (V_r20 > Vr_max) || (V_r20 < Vr_min)
        @error("Regulator Voltage V_R = $(V_r20) outside the limits")
    end
    Tc_Tb_ratio = Tb <= eps() ? 0.0 : Tc / Tb
    V_r30 = -(Kf / Tf) * V_FE0
    V_r10 = (V_r20 / Ka) - (Tc_Tb_ratio) * (V_r20 / Ka)
    V_ref0 = Vm0 + V_r20 / Ka

    #States of ESAC1A are Vm, Vr1, Vr2, Ve, Vr3
    #To solve V_ref, Vr1, Vr2, Ve, Vr3
    function f!(out, x)
        V_ref = x[1]
        Vr1 = x[2]
        Vr2 = x[3]
        Ve = x[4]
        Vr3 = x[5]

        #Compute auxiliary variables

        I_N = Kc * Xad_Ifd0 / Ve
        Se = saturation_function(avr, Ve)
        V_FE = Kd * Xad_Ifd0 + Ke * Ve + Ve * Se
        V_F = Vr3 + (Kf / Tf) * V_FE
        V_in = V_ref - Vm0 - V_F
        V_out = Vr1 + (Tc_Tb_ratio) * V_in

        # Time Constants eliminated because don't appear in SteadyState
        out[1] = (V_in * (1 - Tc_Tb_ratio) - Vr1) #dVr1/dt
        out[2] = (Ka * V_out - Vr2) #dVr2/dt
        out[3] = (-(Kf / Tf) * V_FE - Vr3) #dVr3/dt
        out[4] = (Vr2 - V_FE) #dVe/dt
        out[5] = Vf0 - Ve * rectifier_function(I_N)
    end
    x0 = [V_ref0, V_r10, V_r20, V_e0, V_r30]
    sol = NLsolve.nlsolve(f!, x0)
    if !NLsolve.converged(sol)
        @warn("Initialization in AVR failed")
    else
        sol_x0 = sol.zero
        #Update V_ref
        PSY.set_V_ref!(avr, sol_x0[1])
        PSY.get_ext(dynamic_device)[CONTROL_REFS][V_ref_index] = sol_x0[1]
        #Update AVR states
        avr_ix = get_local_state_ix(dynamic_device, typeof(avr))
        avr_states = @view device_states[avr_ix]
        avr_states[1] = Vm0
        avr_states[2] = sol_x0[2] #Vr1
        avr_states[3] = sol_x0[3] #Vr2
        avr_states[4] = sol_x0[4] #Ve
        avr_states[5] = sol_x0[5] #Vr3
    end
end

function initialize_avr!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::PSY.DynamicGenerator{M, S, PSY.SEXS, TG, P},
) where {M <: PSY.Machine, S <: PSY.Shaft, TG <: PSY.TurbineGov, P <: PSY.PSS}
    #Obtain Vf0 solved from Machine
    Vf0 = get_inner_vars(dynamic_device)[Vf_var]
    #Obtain measured terminal voltage
    Vm = sqrt(
        get_inner_vars(dynamic_device)[VR_gen_var]^2 +
        get_inner_vars(dynamic_device)[VI_gen_var]^2,
    )

    #Get parameters
    avr = PSY.get_avr(dynamic_device)
    Ta_Tb = PSY.get_Ta_Tb(avr)
    K = PSY.get_K(avr)

    #States of AVRTypeI are Vf, Vr1, Vr2, Vm
    #To solve V_ref, Vr
    function f!(out, x)
        V_ref = x[1]
        Vr = x[2]

        V_in = V_ref - Vm
        V_LL = Vr + Ta_Tb * V_in

        out[1] = K * V_LL - Vf0
        out[2] = V_in * (1 - Ta_Tb) - Vr
    end
    x0 = [1.0, Vf0]
    sol = NLsolve.nlsolve(f!, x0)
    if !NLsolve.converged(sol)
        @warn("Initialization in AVR failed")
    else
        sol_x0 = sol.zero
        #Update V_ref
        PSY.set_V_ref!(avr, sol_x0[1])
        PSY.get_ext(dynamic_device)[CONTROL_REFS][V_ref_index] = sol_x0[1]
        #Update AVR states
        avr_ix = get_local_state_ix(dynamic_device, PSY.SEXS)
        avr_states = @view device_states[avr_ix]
        avr_states[1] = Vf0 #Vf
        avr_states[2] = sol_x0[2] #Vr
    end
end

ext: jl


text:
module xla
  const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
  include("xla_data_pb.jl")
  include("hlo_pb.jl")
  include("xla_pb.jl")
end


text:
using IterTools, ProgressBars
using JLD
using LaTeXStrings, PyPlot

include("synthetic.jl")

RUNS    = 1
EPOCHS  = 1000
T       = n * EPOCHS

SEED    = 42
Random.seed!(SEED)

function gen_syn_data_all()
    gaussian_samples = Dict()
    # for each distribution, pre-generate the stream of examples
    # to be used in T iterations of online SGD
    # for a specified number of runs
    iter_qrngs = ProgressBar(qrngs)
    for (qrng_name, ~) in iter_qrngs
        println(iter_qrngs, "generating data for $qrng_name...")

        data = []
        iter_runs = ProgressBar(1:RUNS)
        for run in iter_runs
            println(iter_runs, "working on run $run")
            data_run = gen_syn_data(T, qrng_name, d+1, SEED+run)
            push!(data, data_run)
        end
        gaussian_samples[qrng_name] = data
    end
    return gaussian_samples
end

println("generating data...")
samples_online = gen_syn_data_all()

println("begin training...")
result = Dict()
for qrng_name in ProgressBar(qrng_names)
    subopt_list = zeros(RUNS, T)

    for run in 1:RUNS
        # run sgd for online mode
        samples = samples_online[qrng_name][run]
        subopt_list[run, :] = sgd(T, w0, wgen, samples, collect(1:T))
    end
    result[qrng_name] = subopt_list
end

println("training completed")
println("plotting...")

# final plotting
figure(figsize = (6,6))
for qrng_name in ["IID Uniform", "Sobol"]
    subopt_list = result[qrng_name]
    subopt_mean = mean(subopt_list, dims=1)'[100:100:T]

    loglog(collect(100:100:T), subopt_mean, label=qrng_name)
end

loglog(collect(100:100:T), d * (collect(100:100:T)).^(-1), label=L"1/t", color="#d62728", linewidth=3, linestyle="dashed")
loglog(collect(100:100:T), 16 * d^2 *(collect(100:100:T)).^(-2), label=L"1/t^2", color="#2ca02c", linewidth=3, linestyle="dashed")
legend(loc="lower left", fontsize=18)

xlabel(string("iterations ", L"t"), fontsize=20)
ylabel(L"\Vert w_t-w^*\Vert^2", fontsize=20)
title("Online", fontsize=24)
xticks(fontsize=20)
yticks(fontsize=20)

ax = gca()
locmin = matplotlib.ticker.LogLocator(base=10.0, subs=(0.25,0.5,0.75), numticks=50)
ax.yaxis.set_minor_locator(locmin)
ax.yaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())

tight_layout()
plot_fn = joinpath(pwd(),"figs", "qmc-online.pdf")
savefig(plot_fn)
println("plot saved at $plot_fn")

max_stars_repo_name: EugeneLYC/qmc-ordering

max_issues_repo_name: EugeneLYC/qmc-ordering

max_forks_repo_name: EugeneLYC/qmc-ordering


text:
function discretizeheart!(system::CVSystem)
    append!(system.heart.lv.V,zeros(system.solverparams.numsteps+1))
    append!(system.heart.lv.P,zeros(system.solverparams.numsteps+1))
    append!(system.heart.lv.E,zeros(system.solverparams.numsteps+1))

    append!(system.heart.rv.V,zeros(system.solverparams.numsteps+1))
    append!(system.heart.rv.P,zeros(system.solverparams.numsteps+1))
    append!(system.heart.rv.Q,zeros(system.solverparams.numsteps+1))
    append!(system.heart.rv.E,zeros(system.solverparams.numsteps+1))

    append!(system.heart.la.V,zeros(system.solverparams.numsteps+1))
    append!(system.heart.la.P,zeros(system.solverparams.numsteps+1))
    append!(system.heart.la.Q,zeros(system.solverparams.numsteps+1))

    append!(system.heart.ra.V,zeros(system.solverparams.numsteps+1))
    append!(system.heart.ra.P,zeros(system.solverparams.numsteps+1))
    append!(system.heart.ra.Q,zeros(system.solverparams.numsteps+1))

    append!(system.heart.av.zeta,zeros(system.solverparams.numsteps+1))
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
#
# This file is part of the SlowActors.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

using SlowActors, Test, .Threads

a = Ref{Int}()
a[] = 1

inca(a, b) = a[] = a[] + b

slow = spawn(Func(inca, a))

send!(slow, 1)
sleep(0.1)
@test a[] == 2
@test request!(slow, 1) == 3
@test a[] == 3
become!(slow, threadid)
@test request!(slow) > 1

slow1 = spawn(Func(threadid), sticky=true)
@test request!(slow1) == 1


text:
using Gadfly

plot(x=rand(20), y=rand(20),
     color=vcat(fill("A", 5), fill("B", 5), fill("C", 5), fill("D", 5)),
     ygroup=vcat(fill("U", 10), fill("V", 10)),
     Geom.subplot_grid(Scale.color_discrete, Geom.line))


text:
using Random
using CUDA
using Plots

function generate_random_walks!(positions, cache, T)
    positions .= zero(eltype(positions))

    # Make the code typesafe!
    limit = convert(eltype(cache), 0.5);
    el_two = convert(eltype(positions), 2);
    el_one = one(eltype(positions))

    for _ in 1:T
        rand!(cache) # Put random numbers in the cache
        positions .+= el_two .* (cache .< limit) .- el_one
    end
    nothing
end

# positions = zeros(10^6);
# cache = similar(positions);
# generate_random_walks!(positions, cache, 20)

# using CUDA # supplies cu function
# positions_gpu = cu(zeros(10^6));
# cache_gpu = similar(positions_gpu);
# generate_random_walks!(positions_gpu, cache_gpu, 20)

# histogram(positions)

# histogram(Array(positions_gpu))

function time_function(fn, params...; repeats = 5)
    min_time = Inf64
    for _ in 1:repeats
        new_time = @elapsed fn(params...)
        min_time = min(new_time, min_time)
    end
    return min_time
end

function measure_array_times(samples, T)
    function measure_cpu(n)
        dist = zeros(Int32, n)
        cache = zeros(Float32, n)
        return time_function(generate_random_walks!, dist, cache, T)
    end
    function measure_gpu(n)
        dist = cu(zeros(Int32, n))
        cache = cu(zeros(Float32, n))
        return time_function(generate_random_walks!, dist, cache, T)
    end
    cpu_times = zeros(Float64, length(samples))
    gpu_times = zeros(Float64, length(samples))
    for (i, n) in enumerate(samples)
        cpu_times[i] = measure_cpu(n)
        gpu_times[i] = measure_gpu(n)
        # Save memory
        GC.gc()
        CUDA.reclaim()
    end
    return cpu_times, gpu_times
end

function plot_threaded_gpu(T, samples, cpu_threaded_times, gpu_times)
    plt = plot(samples, cpu_threaded_times, label="Threaded"; markershape=:utriangle)
    plot!(samples, gpu_times, label="GPU"; markershape=:square, xscale=:log10, yscale=:log10)
    xlabel!("n")
    ylabel!("Time (s) for T=$T")
    return plt
end

samples = 10 .^ collect(0: 9);
cpu_times_array, gpu_times_array = measure_array_times(samples, 20);

plt_threaded_vs_gpu_array = plot_threaded_gpu(20, samples, cpu_times_array, gpu_times_array);
display(plt_threaded_vs_gpu_array)

max_stars_repo_name: JamieMair/intro-to-julia-talk

max_issues_repo_name: JamieMair/intro-to-julia-talk

max_forks_repo_name: JamieMair/intro-to-julia-talk


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

## The LAPACK module of interfaces to LAPACK subroutines
module LAPACK

const liblapack = Base.liblapack_name

import Base.@blasfunc

import Base.LinAlg: BlasFloat, Char, BlasInt, LAPACKException,
    DimensionMismatch, SingularException, PosDefException, chkstride1, checksquare

#Generic LAPACK error handlers
"""
Handle only negative LAPACK error codes

*NOTE* use only if the positive error code is useful.
"""
function chkargsok(ret::BlasInt)
    if ret < 0
        throw(ArgumentError("invalid argument #$(-ret) to LAPACK call"))
    end
end

"Handle all nonzero info codes"
function chklapackerror(ret::BlasInt)
    if ret == 0
        return
    elseif ret < 0
        throw(ArgumentError("invalid argument #$(-ret) to LAPACK call"))
    else # ret > 0
        throw(LAPACKException(ret))
    end
end

function chknonsingular(ret::BlasInt)
    if ret > 0
        throw(SingularException(ret))
    end
end

function chkposdef(ret::BlasInt)
    if ret > 0
        throw(PosDefException(ret))
    end
end

"Check that upper/lower (for special matrices) is correctly specified"
function chkuplo(uplo::Char)
    if !(uplo == 'U' || uplo == 'L')
        throw(ArgumentError("uplo argument must be 'U' (upper) or 'L' (lower), got $uplo"))
    end
    uplo
end

"Check that {c}transpose is correctly specified"
function chktrans(trans::Char)
    if !(trans == 'N' || trans == 'C' || trans == 'T')
        throw(ArgumentError("trans argument must be 'N' (no transpose), 'T' (transpose), or 'C' (conjugate transpose), got $trans"))
    end
    trans
end

"Check that left/right hand side multiply is correctly specified"
function chkside(side::Char)
    if !(side == 'L' || side == 'R')
        throw(ArgumentError("side argument must be 'L' (left hand multiply) or 'R' (right hand multiply), got $side"))
    end
    side
end

"Check that unit diagonal flag is correctly specified"
function chkdiag(diag::Char)
    if !(diag == 'U' || diag =='N')
        throw(ArgumentError("diag argument must be 'U' (unit diagonal) or 'N' (non-unit diagonal), got $diag"))
    end
    diag
end

subsetrows(X::AbstractVector, Y::AbstractArray, k) = Y[1:k]
subsetrows(X::AbstractMatrix, Y::AbstractArray, k) = Y[1:k, :]

function chkfinite(A::StridedMatrix)
    for a in A
        if !isfinite(a)
            throw(ArgumentError("matrix contains Infs or NaNs"))
        end
    end
    return true
end

# LAPACK version number
function laver()
    major = Ref{BlasInt}(0)
    minor = Ref{BlasInt}(0)
    patch = Ref{BlasInt}(0)
    ccall((@blasfunc(ilaver_), liblapack), Void,
        (Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}),
        major, minor, patch)
    return major[], minor[], patch[]
end

## (ST) Symmetric tridiagonal - eigendecomposition
for (stev, stebz, stegr, stein, elty) in
    ((:dstev_,:dstebz_,:dstegr_,:dstein_,:Float64),
     (:sstev_,:sstebz_,:sstegr_,:sstein_,:Float32)
#     , (:zstev_,:Complex128)  Need to rewrite for ZHEEV, rwork, etc.
#     , (:cstev_,:Complex64)
     )
    @eval begin
        function stev!(job::Char, dv::AbstractVector{$elty}, ev::AbstractVector{$elty},
                       Zmat::StridedMatrix{$elty}, work::Array{$elty}, info::Ref{BlasInt})
            n = length(dv)

            ldz = stride(Zmat, 2)

            ccall((@blasfunc($stev), liblapack), Void,
                  (Ptr{UInt8}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty},
                   Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}),
                  &job, &n, dv, ev, Zmat, &ldz, work, info)
            chklapackerror(info[])
        end

        function stev!(job::Char, dv::Vector{$elty}, ev::Vector{$elty}, query::Bool = false)
            n = length(dv)
            if length(ev) != n - 1
                throw(DimensionMismatch("ev has length $(length(ev)) but needs one less than dv's length, $n)"))
            end

            Zmat = similar(dv, $elty, (n, job != 'N' ? n : 0))
            work = Array($elty, max(1, 2n-2))

            info = Ref{BlasInt}()

            query && return Zmat, work, info
            stev!(job, dv, ev, Zmat, work, info)
            dv, Zmat
        end

        function stegr!(jobz::Char, range::Char,
                        dv::AbstractVector{$elty}, eev::AbstractVector{$elty},
                        vl::Real, vu::Real, il::Integer, iu::Integer,
                        abstol::Array{$elty}, m::Array{BlasInt}, w::Vector{$elty}, Z::StridedMatrix{$elty},
                        isuppz::Vector{BlasInt},
                        work::Array{$elty}, lwork::BlasInt,
                        iwork::Array{BlasInt}, liwork::BlasInt,
                        info::Ref{BlasInt})
            n = length(dv)

            ldz = stride(Z, 2)

            ccall((@blasfunc($stegr), liblapack), Void,
                  (Ptr{UInt8}, Ptr{UInt8}, Ptr{BlasInt}, Ptr{$elty},
                   Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt},
                   Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty},
                   Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty},
                   Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}),
                  &jobz, &range, &n, dv,
                  eev, &vl, &vu, &il,
                  &iu, abstol, m, w,
                  Z, &ldz, isuppz, work,
                  &lwork, iwork, &liwork, info)
            chklapackerror(info[])
        end

        function stegr!(jobz::Char, range::Char,
                        dv::Vector{$elty}, ev::Vector{$elty},
                        vl::Real, vu::Real, il::Integer, iu::Integer,
                        query::Bool = false)
            n = length(dv)
            if length(ev) != n - 1
                throw(DimensionMismatch("ev has length $(length(ev)) but needs one less than dv's length, $n)"))
            end
            eev = [ev; zero($elty)]

            abstol = Array($elty, 1)
            m = Array(BlasInt, 1)
            w = similar(dv, $elty, n)
            ldz = jobz == 'N' ? 1 : n
            Z = similar(dv, $elty, ldz, n)
            isuppz = similar(dv, BlasInt, 2n)
            work = Array($elty, 1)
            lwork = BlasInt(-1)
            iwork = Array(BlasInt, 1)
            liwork = BlasInt(-1)

            info = Ref{BlasInt}()
            stegr!(jobz, range, dv, eev, vl, vu, il, iu, abstol, m, w, Z, isuppz, work, lwork, iwork, liwork, info)

            lwork = BlasInt(work[1])
            work = Array($elty, lwork)
            liwork = iwork[1]
            iwork = Array(BlasInt, liwork)

            query && return abstol, m, w, Z, isuppz, work, lwork, iwork, liwork, info
            stegr!(jobz, range, dv, eev, vl, vu, il, iu, abstol, m, w, Z, isuppz, work, lwork, iwork, liwork, info)

            w[1:m[1]], Z[:,1:m[1]]
        end
    end
end
stegr!(jobz::Char, dv::Vector, ev::Vector) = stegr!(jobz, 'A', dv, ev, 0.0, 0.0, 0, 0)

"""
    stev!(job, dv, ev) -> (dv, Zmat)

Computes the eigensystem for a symmetric tridiagonal matrix with `dv` as
diagonal and `ev` as off-diagonal. If `job = N` only the eigenvalues are
found and returned in `dv`. If `job = V` then the eigenvectors are also found
and returned in `Zmat`.
"""
stev!(job::Char, dv::Vector, ev::Vector)

"""
    stegr!(jobz, range, dv, ev, vl, vu, il, iu) -> (w, Z)

Computes the eigenvalues (`jobz = N`) or eigenvalues and eigenvectors
(`jobz = V`) for a symmetric tridiagonal matrix with `dv` as diagonal
and `ev` as off-diagonal. If `range = A`, all the eigenvalues
are found. If `range = V`, the eigenvalues in the half-open interval
`(vl, vu]` are found. If `range = I`, the eigenvalues with indices between
`il` and `iu` are found. The eigenvalues are returned in `w` and the eigenvectors
in `Z`.
"""
stegr!(jobz::Char, range::Char, dv::Vector, ev::Vector, vl::Real, vu::Real, il::Integer, iu::Integer)

end # module


text:
# Functions for creating a Cayley digraph based on S_n and permutation
# generators.

using Permutations, Combinatorics
using SimpleGraphs

function Cayley(gens::Array{Permutation,1})
    ng = length(gens)
    if ng == 0
        return SimpleDigraph{Permutation}()
    end

    n = length(gens[1])

    for k = 2:ng
        if length(gens[k]) != length(gens[k-1])
            error("All generators must be same length")
        end
    end

    verts = [Permutation(p) for p in permutations(1:n)]

    G = SimpleDigraph{Permutation}()

    for v in verts
        add!(G, v)
        for p in gens
            add!(G, v, v * p)
        end
    end

    return G
end

Cayley(perms...) = Cayley(collect(perms))


text:
using Statistics, Printf

using Oceananigans
using Oceananigans.Fields
using Oceananigans.TurbulenceClosures
using Oceananigans.OutputWriters
using Oceananigans.Diagnostics
using Oceananigans.Utils

""" Friction velocity. See equation (16) of Vreugdenhil & Taylor (2018). """
function uτ(model, Uavg, U_wall)
    Nz, Hz, Δz = model.grid.Nz, model.grid.Hz, model.grid.Δz
    ν = model.closure.ν

    compute!(Uavg)
    U = Array(interior(Uavg))  # Exclude average of halo region.

    # Use a finite difference to calculate dU/dz at the top and bottomtom walls.
    # The distance between the center of the cell adjacent to the wall and the
    # wall itself is Δz/2.
    uτ²_top    = ν * abs(U[1] - U_wall)    / (Δz/2)  # Top wall    where u = +U_wall
    uτ²_bottom = ν * abs(-U_wall  - U[Nz]) / (Δz/2)  # Bottom wall where u = -U_wall

    uτ_top, uτ_bottom = √uτ²_top, √uτ²_bottom

    return uτ_top, uτ_bottom
end

""" Heat flux at the wall. See equation (16) of Vreugdenhil & Taylor (2018). """
function q_wall(model, Tavg, Θ_wall)
    Nz, Hz, Δz = model.grid.Nz, model.grid.Hz, model.grid.Δz
    κ = model.closure.κ.T

    compute!(Tavg)
    Θ = Array(interior(Tavg)) # Exclude average of halo region.

    # Use a finite difference to calculate dθ/dz at the top and bottomtom walls.
    # The distance between the center of the cell adjacent to the wall and the
    # wall itself is Δz/2.
    q_wall_top    = κ * abs(Θ[1] - Θ_wall)   / (Δz/2)  # Top wall    where Θ = +Θ_wall
    q_wall_bottom = κ * abs(-Θ_wall - Θ[Nz]) / (Δz/2)  # Bottom wall where Θ = -Θ_wall

    return q_wall_top, q_wall_bottom
end

struct FrictionReynoldsNumber{H, U}
    Uavg :: H
    U_wall :: U
end

struct NusseltNumber{H, T}
    Tavg :: H
    Θ_wall :: T
end

""" Friction Reynolds number. See equation (20) of Vreugdenhil & Taylor (2018). """
function (Reτ::FrictionReynoldsNumber)(model)
    ν = model.closure.ν
    h = model.grid.Lz / 2
    uτ_top, uτ_bottom = uτ(model, Reτ.Uavg, Reτ.U_wall)

    return h * uτ_top / ν, h * uτ_bottom / ν
end

""" Nusselt number. See equation (20) of Vreugdenhil & Taylor (2018). """
function (Nu::NusseltNumber)(model)
    κ = model.closure.κ.T
    h = model.grid.Lz / 2

    q_wall_top, q_wall_bottom = q_wall(model, Nu.Tavg, Nu.Θ_wall)

    return (q_wall_top * h)/(κ * Nu.Θ_wall), (q_wall_bottom * h)/(κ * Nu.Θ_wall)
end

"""
    simulate_stratified_couette_flow(; Nxy, Nz, h=1, U_wall=1, Re=4250, Pr=0.7,
                                     Ri, Ni=10, end_time=1000)

Simulate stratified plane Couette flow with `Nxy` grid cells in each horizontal
direction, `Nz` grid cells in the vertical, in a domain of size (4πh, 2πh, 2h),
with wall velocities of `U_wall` at the top and -`U_wall` at the bottom, at a Reynolds
number `Re, Prandtl number `Pr`, and Richardson number `Ri`.

`Ni` is the number of "intermediate" time steps taken at a time before printing a progress
statement and updating the time step.
"""
function simulate_stratified_couette_flow(; Nxy, Nz, arch=GPU(), h=1, U_wall=1,
                                          Re=4250, Pr=0.7, Ri, Ni=10, end_time=1000)
    #####
    ##### Computed parameters
    #####

         ν = U_wall * h / Re    # From Re = U_wall h / ν
    Θ_wall = Ri * U_wall^2 / h  # From Ri = L Θ_wall / U_wall²
         κ = ν / Pr             # From Pr = ν / κ

    #####
    ##### Impose boundary conditions
    #####

    grid = RegularRectilinearGrid(size = (Nxy, Nxy, Nz), extent = (4π*h, 2π*h, 2h))

    Tbcs = FieldBoundaryConditions(top = BoundaryCondition(Value,  Θ_wall),
                                   bottom = BoundaryCondition(Value, -Θ_wall))

    ubcs = FieldBoundaryConditions(top = BoundaryCondition(Value,  U_wall),
                                   bottom = BoundaryCondition(Value, -U_wall))

    vbcs = FieldBoundaryConditions(top = BoundaryCondition(Value, 0),
                                   bottom = BoundaryCondition(Value, 0))

    #####
    ##### Non-dimensional model setup
    #####

    model = NonhydrostaticModel(
               architecture = arch,
                       grid = grid,
                   buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(α=1.0, β=0.0)),
                    closure = AnisotropicMinimumDissipation(ν=ν, κ=κ),
        boundary_conditions = (u=ubcs, v=vbcs, T=Tbcs)
    )

    #####
    ##### Set initial conditions
    #####

    # Add a bit of surface-concentrated noise to the initial condition
    ε(σ, z) = σ * randn() * z/model.grid.Lz * (1 + z/model.grid.Lz)

    # We add a sinusoidal initial condition to u to encourage instability.
    T₀(x, y, z) = 2Θ_wall * (1/2 + z/model.grid.Lz) * (1 + ε(5e-1, z))
    u₀(x, y, z) = 2U_wall * (1/2 + z/model.grid.Lz) * (1 + ε(5e-1, z)) * (1 + 0.5*sin(4π/model.grid.Lx * x))
    v₀(x, y, z) = ε(5e-1, z)
    w₀(x, y, z) = ε(5e-1, z)

    set!(model, u=u₀, v=v₀, w=w₀, T=T₀)

    #####
    ##### Print simulation banner
    #####

    @printf(
        """
        Simulating stratified plane Couette flow

                N : %d, %d, %d
                L : %.3g, %.3g, %.3g
               Re : %.3f
               Ri : %.3f
               Pr : %.3f
                ν : %.3g
                κ : %.3g
           U_wall : %.3f
           Θ_wall : %.3f

        """, model.grid.Nx, model.grid.Ny, model.grid.Nz,
             model.grid.Lx, model.grid.Ly, model.grid.Lz,
             Re, Ri, Pr, ν, κ, U_wall, Θ_wall)

    #####
    ##### Set up field output writer
    #####

    base_dir = @sprintf("stratified_couette_flow_data_Nxy%d_Nz%d_Ri%.2f", Nxy, Nz, Ri)
    prefix = @sprintf("stratified_couette_flow_Nxy%d_Nz%d_Ri%.2f", Nxy, Nz, Ri)

    function init_save_parameters_and_bcs(file, model)
        file["parameters/reynolds_number"] = Re
        file["parameters/richardson_number"] = Ri
        file["parameters/prandtl_number"] = Pr
        file["parameters/viscosity"] = ν
        file["parameters/diffusivity"] = κ
        file["parameters/wall_velocity"] = U_wall
        file["parameters/wall_temperature"] = Θ_wall
    end

    fields = Dict(
        :u => model -> Array(model.velocities.u.data.parent),
        :v => model -> Array(model.velocities.v.data.parent),
        :w => model -> Array(model.velocities.w.data.parent),
        :T => model -> Array(model.tracers.T.data.parent),
   :kappaT => model -> Array(model.diffusivity_fields.κₑ.T.data.parent),
       :nu => model -> Array(model.diffusivity_fields.νₑ.data.parent))

    field_writer =
        JLD2OutputWriter(model, fields, dir=base_dir, prefix=prefix * "_fields",
                         init=init_save_parameters_and_bcs, schedule=TimeInterval(10),
                         force=true, verbose=true)

    #####
    ##### Set up profile output writer
    #####

    Uavg = AveragedField(model.velocities.u,            dims=(1, 2))
    Vavg = AveragedField(model.velocities.v,            dims=(1, 2))
    Wavg = AveragedField(model.velocities.w,            dims=(1, 2))
    Tavg = AveragedField(model.tracers.T,               dims=(1, 2))
    νavg = AveragedField(model.diffusivity_fields.νₑ,   dims=(1, 2))
    κavg = AveragedField(model.diffusivity_fields.κₑ.T, dims=(1, 2))

    profiles = Dict(
         :u => Uavg,
         :v => Vavg,
         :w => Wavg,
         :T => Tavg,
        :nu => νavg,
    :kappaT => κavg)

    profile_writer =
        JLD2OutputWriter(model, profiles, dir=base_dir, prefix=prefix * "_profiles",
                         init=init_save_parameters_and_bcs, schedule=TimeInterval(1),
                         force=true, verbose=true)

    #####
    ##### Set up statistic output writer
    #####

    Reτ = FrictionReynoldsNumber(Uavg, U_wall)
     Nu = NusseltNumber(Tavg, Θ_wall)

    statistics = Dict(
        :Re_tau => model -> Reτ(model),
        :Nu     => model -> Nu(model))

    statistics_writer =
        JLD2OutputWriter(model, statistics, dir=base_dir, prefix=prefix * "_statistics",
                         init=init_save_parameters_and_bcs, schedule=TimeInterval(1/2),
                         force=true, verbose=true)

    #####
    ##### Time stepping
    #####

    wizard = TimeStepWizard(cfl=0.02, Δt=0.0001, max_change=1.1, max_Δt=0.02)

    # We will ramp up the CFL used by the adaptive time step wizard during spin up.
    cfl(t) = min(0.01t, 0.1)

    function print_progress(simulation)
        model = simulation.model
        clock = model.clock

        wizard.cfl = cfl(model.clock.time)

        progress = 100 * (clock.time / end_time)

        umax = maximum(abs, model.velocities.u.data.parent)
        vmax = maximum(abs, model.velocities.v.data.parent)
        wmax = maximum(abs, model.velocities.w.data.parent)
        CFL = wizard.Δt / cell_advection_timescale(model)

        Δ = min(model.grid.Δx, model.grid.Δy, model.grid.Δz)
        νmax = maximum(model.diffusivity_fields.νₑ.data.parent)
        κmax = maximum(model.diffusivity_fields.κₑ.T.data.parent)
        νCFL = wizard.Δt / (Δ^2 / νmax)
        κCFL = wizard.Δt / (Δ^2 / κmax)

        @printf("[%06.2f%%] i: %d, t: %.2e, umax: (%.2e, %.2e, %.2e), CFL: %.2e, νκmax: (%.2e, %.2e), νκCFL: (%.2e, %.2e), next Δt: %.2e, wall time: %s\n",
                progress, model.clock.iteration, model.clock.time, umax, vmax, wmax,
                CFL, νmax, κmax, νCFL, κCFL, wizard.Δt, prettytime(simulation.run_time))
    end

    simulation = Simulation(model, Δt=wizard, stop_time=end_time,
                            progress=print_progress, iteration_interval=Ni)
    push!(simulation.output_writers, field_writer, profile_writer, statistics_writer)
    run!(simulation)

    return simulation
end


text:
using DatasetsCF
using Persa
using DataFrames

datasets = [DatasetsCF.MovieLens,
            DatasetsCF.MovieTweeting,
            DatasetsCF.CiaoDVD,
            DatasetsCF.FilmTrust,
            DatasetsCF.YahooMusic,
            DatasetsCF.LastFM]
######
ds = datasets[5]()

hist = Persa.histogram(ds)
df = DataFrame()

df[:rating] = Array{eltype(ds.preferences)}(size(ds.preferences))
df[:a] = Array{Int}(size(ds.preferences))
df[:b] = Array{Float64}(size(ds.preferences))

i = 1

for rating in sort(ds.preferences.possibles)

    df[:rating][i] = rating
    df[:a][i] = hist[rating]
    i = i + 1
end

df[:b] = df[:a] ./ sum(df[:a])

df[:rating]
df[:a]
df[:b]

##############
datasets
ds = datasets[5]()

hist_user = zeros(Int, ds.users)
hist_item = zeros(Int, ds.items)

for (u, v, r) in ds
    hist_user[u] = hist_user[u] + 1
    hist_item[v] = hist_item[v] + 1
end

sort!(hist_user, rev = true)
sort!(hist_item, rev = true)
###
factor = convert(Int, round(length(hist_user) / 200))

open("t.txt", "w") do f
    write(f, "n\tcount\n")
    for i = 1:length(hist_user)
        if i % factor == 1
            write(f, "$i\t$(hist_user[i])\n")
        end
    end
end

###
factor = convert(Int, round(length(hist_item) / 200))

open("t.txt", "w") do f
    write(f, "n\tcount\n")
    for i = 1:length(hist_item)
        if i % factor == 1
            write(f, "$i\t$(hist_item[i])\n")
        end
    end
end

#####
ds = datasets[3]()
sort!(ds.file, cols = :timestamp)
ds.file[:timestamp]

moments = Dict{Int,Int}()

times = unique(ds.file[:timestamp])

for i = 1:length(times)
    moments[times[i]] = 0
end

for i = 1:length(ds.file[:timestamp])
    moments[ds.file[:timestamp][i]] += 1
end

moments2 = hcat([keys(moments)...], [values(moments)...])

moments2_sorted = sortrows(moments2)

moments2_sorted[:,1] = moments2_sorted[:,1] .- moments2_sorted[1,1]

factor = convert(Int, round(size(moments2_sorted)[1] / 200))

open("t.txt", "w") do f
    write(f, "n\tcount\n")
    for i = 1:size(moments2_sorted)[1]
        if i % factor == 1
            write(f, "$(moments2_sorted[i,1])\t$(sum(moments2_sorted[1:i,2]) ./ length(ds))\n")
        end
    end
end

###############
sort!(ds.file, cols = :timestamp)
ds.file[:timestamp]

times = unique(ds.file[:timestamp])

qnt = collect(1:length(times))

for i = 1:length(times)
    repeats = length(find(r->r == times[i], ds.file[:timestamp]))
    if repeats > 1
        qnt[i] = qnt[i] + repeats - 1
    end
end

times = times .- times[1]

################
using PyPlot
ds = datasets[1]()
sort!(ds.file, cols = :timestamp)
ds.file[:timestamp] = ds.file[:timestamp] .- ds.file[:timestamp][1]

moments = Dict{Int,Array{Int}}()

times = unique(ds.file[:timestamp])

for i = 1:length(times)
    moments[times[i]] = Array{Int}(0)
end

for (u, v, r, t) in ds
    push!(moments[t], u)
end

moments2 = Array{Tuple{Int,Int,Int}}(0)

for i = 1:length(times)
    users = moments[times[i]]
    for j = 1:length(unique(users))
        push!(moments2, (times[i], users[j], length(find(r->r == users[j], users))))
        println(length(find(r->r == users[j], users)))
    end
end

x = Array{Int}(length(moments2))
y = Array{Int}(length(moments2))
z = Array{Int}(length(moments2))

for i = 1:length(moments2)
    x[i] = moments2[i][1]
    y[i] = moments2[i][2]
    z[i] = moments2[i][3]
end

###
moments = hcat(ds.file[:user], ds.file[:timestamp])
moments[:,2] = moments[:,2] .- moments[1,2]

select = find(r->r in [1:10...], y)
select = (length(x) - 1000):length(x)
###
scatter(x[select], y[select], s = z[select] .* 25, alpha = 0.5)

for i = 1:10
    plot(x[select], repeat([i], inner = length(x[select])))
end

grid("on")
title("MovieLens")
xlabel("Tempo")
ylabel("Usuário")


text:
"""

$(TYPEDEF)

Unit conversions.

$(TYPEDFIELDS)

"""
@with_kw struct Units
  mole = 6.022140857e23
  Angs3tocm3 = 1e24
  Angs3toL = 1e27
  Angs3tocm3permol = mole/Angs3tocm3
  Angs3toLpermol = mole/Angs3toL
  SitesperAngs3tomolperL = Angs3toL/mole
end
units = Units()


text:
using UnitJuMP
using Test
import HiGHS

function test_simple_example()
    model = Model(HiGHS.Optimizer)
    set_silent(model)
    max_speed = 2u"m/s"
    a = 9.8u"ft/s^2"
    M = 0.5u"m/s"
    no_speed = 0.0u"m/s"
    data = 4
    @variables(model, begin
        x >= 0, u"m/s"
        y >= 0, u"s"
        z >= 0, u"km/s"
        v[1:4], u"m/s"
        w, Bin
    end)
    @constraints(model, begin
        x + a * y <= max_speed
        z + x <= max_speed
        x + 0.5 * a * y <= max_speed
        x <= M * w
        -3x + 2z <= no_speed
        x == 0.15 * max_speed
        con7[i = 1:4], v[i] <= max_speed
        con8, sum(v[i] for i in 1:3) <= max_speed, u"ft/s"
        -3 * data * x + 2z <= no_speed
    end)
    @objective(model, Max, x + 0.2 * z)
    optimize!(model)
    @test value(x) == 0.3u"m/s"
    @test value(y) == 0.0u"s"
    @test value(z) ≈ 0.00045u"km/s"
    @test value(x + 0.2 * z) ≈ 0.39u"m/s"
    return
end

test_simple_example()


text:
module RandomisationModels

using DocStringExtensions
using Random: AbstractRNG, GLOBAL_RNG
using StatsBase: sample, Weights, FrequencyWeights
using Distances: euclidean
using LinearAlgebra: normalize

import Base.deepcopy
import Base.show
import Base.==

export MultiArmRandomisationModel
export CompleteRandomisation
export MassWeightedUrn
export PermutedBlock
export BlockUrn

export target, narms, sequence, dist, imbalance, predictability, randomise, randomize, randomise!, randomize!
export prob, update!
export mass, blocksize, weights

# Internal helper functions
include("util.jl")

include("types.jl")
include("models/complete-randomisation.jl")
include("models/mass-weighted-urn.jl")
include("models/permuted-block.jl")
include("models/block-urn.jl")

end # module


text:
module ArtifactTests
import ..Pkg # ensure we are using the correct Pkg

using Test, Random, Pkg.Artifacts, Base.BinaryPlatforms, Pkg.PlatformEngines
import Pkg.Artifacts: pack_platform!, unpack_platform, with_artifacts_directory, ensure_all_artifacts_installed, extract_all_hashes
using TOML, Dates
import Base: SHA1

using ..Utils

# Helper function to create an artifact, then chmod() the whole thing to 0o644.  This is
# important to keep hashes stable across platforms that have different umasks, changing
# the permissions within a tree hash, breaking our tests.
function create_artifact_chmod(f::Function)
    create_artifact() do path
        f(path)

        # Change all files to have 644 permissions, leave directories alone
        for (root, dirs, files) in walkdir(path)
            for f in files
                f = joinpath(root, f)
                islink(f) || chmod(f, 0o644)
            end
        end
    end
end

@testset "Artifact Creation" begin
    # We're going to ensure that our artifact creation does in fact give git-tree-sha1's.
    creators = [
        # First test the empty artifact
        (path -> begin
            # add no contents
        end, "4b825dc642cb6eb9a060e54bf8d69288fbee4904"),

        # Next test creating a single file
        (path -> begin
            open(joinpath(path, "foo"), "w") do io
                print(io, "Hello, world!")
            end
        end, "339aad93c0f854604248ea3b7c5b7edea20625a9"),

        # Next we will test creating multiple files
        (path -> begin
            open(joinpath(path, "foo1"), "w") do io
                print(io, "Hello")
            end
            open(joinpath(path, "foo2"), "w") do io
                print(io, "world!")
            end
        end, "98cda294312216b19e2a973e9c291c0f5181c98c"),

        # Finally, we will have nested directories and all that good stuff
        (path -> begin
            mkpath(joinpath(path, "bar", "bar"))
            open(joinpath(path, "bar", "bar", "foo1"), "w") do io
                print(io, "Hello")
            end
            open(joinpath(path, "bar", "foo2"), "w") do io
                print(io, "world!")
            end
            open(joinpath(path, "foo3"), "w") do io
                print(io, "baz!")
            end

            # Empty directories do nothing to effect the hash, so we create one with a
            # random name to prove that it does not get hashed into the rest.  Also, it
            # turns out that life is cxomplex enough that we need to test the nested
            # empty directories case as well.
            rand_dir = joinpath(path, Random.randstring(8), "inner")
            mkpath(rand_dir)

            # Symlinks are not followed, even if they point to directories
            symlink("foo3", joinpath(path, "foo3_link"))
            symlink("../bar", joinpath(path, "bar", "infinite_link"))
        end, "86a1ce580587d5851fdfa841aeb3c8d55663f6f9"),
    ]

    # Enable the following code snippet to figure out the correct gitsha's:
    if false
        for (creator, blah) in creators
            mktempdir() do path
                creator(path)
                for (root, dirs, files) in walkdir(path)
                    for f in files
                        chmod(joinpath(root, f), 0o644)
                    end
                end
                cd(path) do
                    read(`git init .`)
                    read(`git add . `)
                    read(`git commit --allow-empty -m "foo"`)
                    hash = chomp(String(read(`git log -1 --pretty='%T' HEAD`)))
                    println(hash)
                end
            end
        end
    end

    for (creator, known_hash) in creators
        # Create artifact
        hash = create_artifact_chmod(creator)

        # Ensure it hashes to the correct gitsha:
        @test all(hash.bytes .== hex2bytes(known_hash))

        # Test that we can look it up and that it sits in the right place
        @test basename(dirname(artifact_path(hash))) == "artifacts"
        @test basename(artifact_path(hash)) == known_hash
        @test artifact_exists(hash)

        # Test that the artifact verifies
        if !Sys.iswindows()
            @test verify_artifact(hash)
        end
    end

    @testset "File permissions" begin
        mktempdir() do artifacts_dir
            with_artifacts_directory(artifacts_dir) do
                subdir = "subdir"
                file1 = "file1"
                file2 = "file2"
                dir_link = "dir_link"
                file_link = "file_link"
                hash = create_artifact() do dir
                    # Create files, links, and directories
                    mkpath(joinpath(dir, subdir))
                    touch(joinpath(dir, subdir, file1))
                    touch(joinpath(dir, subdir, file2))
                    symlink(basename(subdir), joinpath(dir, dir_link))
                    symlink(basename(file1), joinpath(dir, subdir, file_link))
                end
                artifact_dir = artifact_path(hash)
                # Make sure only files are read-only
                @test iszero(filemode(joinpath(artifact_dir, file1)) & 0o222)
                @test iszero(filemode(joinpath(artifact_dir, file2)) & 0o222)
                @test iszero(filemode(joinpath(artifact_dir, file_link)) & 0o222)
                @test !iszero(filemode(joinpath(artifact_dir, subdir)) & 0o222)
                @test !iszero(filemode(joinpath(artifact_dir, dir_link)) & 0o222)
                # Make sure we can delete the artifact directory without having
                # to manually change permissions
                rm(artifact_dir; recursive=true)
            end
        end
    end
end

@testset "with_artifacts_directory()" begin
    mktempdir() do art_dir
        with_artifacts_directory(art_dir) do
            hash = create_artifact() do path
                touch(joinpath(path, "foo"))
            end
            @test startswith(artifact_path(hash), art_dir)
        end
    end
end

@testset "Artifacts.toml Utilities" begin
    # First, let's test our ability to find Artifacts.toml files;
    ATS = joinpath(@__DIR__, "test_packages", "ArtifactTOMLSearch")
    test_modules = [
        joinpath(ATS, "pkg.jl") =>  joinpath(ATS, "Artifacts.toml"),
        joinpath(ATS, "sub_module", "pkg.jl") =>  joinpath(ATS, "Artifacts.toml"),
        joinpath(ATS, "sub_package", "pkg.jl") =>  joinpath(ATS, "sub_package", "Artifacts.toml"),
        joinpath(ATS, "julia_artifacts_test", "pkg.jl") =>  joinpath(ATS, "julia_artifacts_test", "JuliaArtifacts.toml"),
        joinpath(@__DIR__, "test_packages", "BasicSandbox", "src", "Foo.jl") => nothing,
    ]
    for (test_src, artifacts_toml) in test_modules
        temp_pkg_dir() do pkg_dir
            # Test that the Artifacts.toml that was found is what we expected
            @test find_artifacts_toml(test_src) == artifacts_toml

            # Load `arty` and check its gitsha
            if artifacts_toml !== nothing
                arty_hash = SHA1("43563e7631a7eafae1f9f8d9d332e3de44ad7239")
                @test artifact_hash("arty", artifacts_toml) == arty_hash

                @test arty_hash in extract_all_hashes(artifacts_toml)

                # Ensure it's installable (we uninstall first, to make sure)
                @test !artifact_exists(arty_hash)

                @test ensure_artifact_installed("arty", artifacts_toml) == artifact_path(arty_hash)
                if !Sys.iswindows()
                    @test verify_artifact(arty_hash)
                end

                # Make sure doing it twice "just works"
                @test ensure_artifact_installed("arty", artifacts_toml) == artifact_path(arty_hash)

                # clean up after thyself
                remove_artifact(arty_hash)
                @test !verify_artifact(arty_hash)
            end
        end
    end

    # Test binding/unbinding
    temp_pkg_dir() do path
        hash = create_artifact() do path
            open(joinpath(path, "foo.txt"), "w") do io
                print(io, "hello, world!")
            end
        end

        # Bind this artifact to something
        artifacts_toml = joinpath(path, "Artifacts.toml")
        @test artifact_hash("foo_txt", artifacts_toml) == nothing
        bind_artifact!(artifacts_toml, "foo_txt", hash)

        # Test that this binding worked
        @test artifact_hash("foo_txt", artifacts_toml) == hash
        @test ensure_artifact_installed("foo_txt", artifacts_toml) == artifact_path(hash)

        # Test that binding caused an entry in the manifest_usage.toml
        usage = TOML.parsefile(joinpath(Pkg.logdir(), "artifact_usage.toml"))
        @test any(x -> startswith(x, artifacts_toml), keys(usage))

        # Test that we can overwrite bindings
        hash2 = create_artifact() do path
            open(joinpath(path, "foo.txt"), "w") do io
                print(io, "goodbye, world!")
            end
        end
        @test_throws ErrorException bind_artifact!(artifacts_toml, "foo_txt", hash2)
        @test artifact_hash("foo_txt", artifacts_toml) == hash
        bind_artifact!(artifacts_toml, "foo_txt", hash2; force=true)
        @test artifact_hash("foo_txt", artifacts_toml) == hash2

        # Test that we can un-bind
        unbind_artifact!(artifacts_toml, "foo_txt")
        @test artifact_hash("foo_txt", artifacts_toml) == nothing

        # Test platform-specific binding and providing download_info
        download_info = [
            ("http://google.com/hello_world", "0"^64),
            ("http://microsoft.com/hello_world", "a"^64),
        ]

        # First, test the binding of things with various platforms and overwriting and such works properly
        linux64 = Platform("x86_64", "linux")
        win32 = Platform("i686", "windows")
        bind_artifact!(artifacts_toml, "foo_txt", hash; download_info=download_info, platform=linux64)
        @test artifact_hash("foo_txt", artifacts_toml; platform=linux64) == hash
        @test artifact_hash("foo_txt", artifacts_toml; platform=Platform("x86_64", "macos")) == nothing
        @test_throws ErrorException bind_artifact!(artifacts_toml, "foo_txt", hash2; download_info=download_info, platform=linux64)
        bind_artifact!(artifacts_toml, "foo_txt", hash2; download_info=download_info, platform=linux64, force=true)
        bind_artifact!(artifacts_toml, "foo_txt", hash; download_info=download_info, platform=win32)
        @test artifact_hash("foo_txt", artifacts_toml; platform=linux64) == hash2
        @test artifact_hash("foo_txt", artifacts_toml; platform=win32) == hash
        @test ensure_artifact_installed("foo_txt", artifacts_toml; platform=linux64) == artifact_path(hash2)
        @test ensure_artifact_installed("foo_txt", artifacts_toml; platform=win32) == artifact_path(hash)

        # Next, check that we can get the download_info properly:
        meta = artifact_meta("foo_txt", artifacts_toml; platform=win32)
        @test meta["download"][1]["url"] == "http://google.com/hello_world"
        @test meta["download"][2]["sha256"] == "a"^64
    end

    # Let's test some known-bad Artifacts.toml files
    badifact_dir = joinpath(@__DIR__, "artifacts", "bad")

    # First, parsing errors
    @test_logs (:error, r"contains no `git-tree-sha1`") artifact_meta("broken_artifact", joinpath(badifact_dir, "no_gitsha.toml"))
    @test_logs (:error, r"malformed, must be array or dict!") artifact_meta("broken_artifact", joinpath(badifact_dir, "not_a_table.toml"))

    # Next, test incorrect download errors
    if !Sys.iswindows()
        for ignore_hash in (false, true); withenv("JULIA_PKG_IGNORE_HASHES" => ignore_hash ? "1" : nothing) do; mktempdir() do dir
            with_artifacts_directory(dir) do
                @test artifact_meta("broken_artifact", joinpath(badifact_dir, "incorrect_gitsha.toml")) != nothing
                @test_logs (:error, r"Tree Hash Mismatch!") match_mode=:any begin
                    if !ignore_hash
                        @test_throws ErrorException ensure_artifact_installed("broken_artifact", joinpath(badifact_dir, "incorrect_gitsha.toml"))
                    else
                        path = ensure_artifact_installed("broken_artifact", joinpath(badifact_dir, "incorrect_gitsha.toml"))
                        @test endswith(path, "0000000000000000000000000000000000000000")
                        @test isdir(path)
                    end
                end
            end
        end end end
    end

    mktempdir() do dir
        with_artifacts_directory(dir) do
            @test artifact_meta("broken_artifact", joinpath(badifact_dir, "incorrect_sha256.toml")) != nothing
            @test_logs (:error, r"Hash Mismatch!") match_mode=:any begin
                @test_throws ErrorException ensure_artifact_installed("broken_artifact", joinpath(badifact_dir, "incorrect_sha256.toml"))
            end
        end
    end
end

@testset "Artifact archival" begin
    mktempdir() do art_dir
        with_artifacts_directory(art_dir) do
            hash = create_artifact(p -> touch(joinpath(p, "foo")))
            tarball_path = joinpath(art_dir, "foo.tar.gz")
            archive_artifact(hash, tarball_path)
            @test "foo" in PlatformEngines.list_tarball_files(tarball_path)

            # Test archiving something that doesn't exist fails
            remove_artifact(hash)
            @test_throws ErrorException archive_artifact(hash, tarball_path)
        end
    end
end

@testset "Artifact Usage" begin
    # Don't use `temp_pkg_dir()` here because we need `Pkg.test()` to run in the
    # same package context as the one we're running in right now.  Yes, this pollutes
    # the global artifact namespace and package list, but it should be harmless.
    mktempdir() do project_path
        with_pkg_env(project_path) do
            path = git_init_package(project_path, joinpath(@__DIR__, "test_packages", "ArtifactInstallation"))
            add_this_pkg()
            Pkg.add(Pkg.Types.PackageSpec(
                name="ArtifactInstallation",
                uuid=Base.UUID("02111abe-2050-1119-117e-b30112b5bdc4"),
                path=path,
            ))

            # Run test harness
            Pkg.test("ArtifactInstallation")

            # Also manually do it
            Core.eval(Module(:__anon__), quote
                using ArtifactInstallation
                do_test()
            end)
        end
    end

    # Ensure that `instantiate()` installs artifacts:
    temp_pkg_dir() do project_path
        copy_test_package(project_path, "ArtifactInstallation")
        Pkg.activate(joinpath(project_path, "ArtifactInstallation"))
        add_this_pkg()
        Pkg.instantiate(; verbose=true)

        # Manual test that artifact is installed by instantiate()
        artifacts_toml = joinpath(project_path, "ArtifactInstallation", "Artifacts.toml")
        c_simple_hash = artifact_hash("c_simple", artifacts_toml)
        @test artifact_exists(c_simple_hash)
    end

    # Ensure that porous platform coverage works with ensure_all_installed()
    temp_pkg_dir() do project_path
        copy_test_package(project_path, "ArtifactInstallation")
        artifacts_toml = joinpath(project_path, "ArtifactInstallation", "Artifacts.toml")

        # Install artifacts such that `c_simple` is not installed properly
        # because of the platform we requested, but `socrates` is.
        missing_platform = Platform("powerpc64le", "linux")
        artifacts = select_downloadable_artifacts(artifacts_toml; platform=missing_platform)
        for name in keys(artifacts)
            ensure_artifact_installed(name, artifacts[name], artifacts_toml; platform=missing_platform)
        end

        # Test that c_simple doesn't even show up
        c_simple_hash = artifact_hash("c_simple", artifacts_toml; platform=missing_platform)
        @test c_simple_hash == nothing

        # Test that socrates shows up, but is not installed
        socrates_hash = artifact_hash("socrates", artifacts_toml; platform=missing_platform)
        @test !artifact_exists(socrates_hash)

        # Test that collapse_the_symlink is installed
        cts_hash = artifact_hash("collapse_the_symlink", artifacts_toml; platform=missing_platform)
        @test artifact_exists(cts_hash)
    end

    # Ensure that platform augmentation hooks work.  We will switch between two arbitrary artifacts for this,
    # by inspecting an environment variable in our package hook.
    engaged_hash = SHA1("a5f8161ca1ab2e94fedd3578586fe06d7906177c")
    engaged_url = "https://github.com/JuliaBinaryWrappers/HelloWorldGo_jll.jl/releases/download/HelloWorldGo-v1.0.4+0/HelloWorldGo.v1.0.4.aarch64-linux-musl.tar.gz"
    engaged_sha256 = "9b66d6b02a370d0170a8c217a872cd1f3d53de267d4e63c22a40b49f04367f8a"
    disengaged_hash = SHA1("ea8ea92ecd57aa602d254ca6c637309642202768")
    disengaged_url = "https://github.com/JuliaBinaryWrappers/HelloWorldGo_jll.jl/releases/download/HelloWorldGo-v1.0.4+0/HelloWorldGo.v1.0.4.i686-w64-mingw32.tar.gz"
    disengaged_sha256 = "5c96a327fc6f0dc71d533373bc6cc6719a1e477c72319b800f29abf1b1e7d812"

    function generate_flooblegrank_artifacts(ap_path)
        # Bind both "engaged" and "disengaged" variants of our `gooblebox` artifact to generate an Artifacts.toml file
        artifacts_toml = joinpath(ap_path, "Artifacts.toml")
        engaged_platform = HostPlatform()
        engaged_platform["flooblecrank"] = "engaged"
        Pkg.Artifacts.bind_artifact!(
            artifacts_toml,
            "gooblebox",
            engaged_hash;
            download_info = [(engaged_url, engaged_sha256)],
            platform = engaged_platform,
        )
        disengaged_platform = HostPlatform()
        disengaged_platform["flooblecrank"] = "disengaged"
        Pkg.Artifacts.bind_artifact!(
            artifacts_toml,
            "gooblebox",
            disengaged_hash;
            download_info = [(disengaged_url, disengaged_sha256)],
            platform = disengaged_platform,
        )
    end

    for flooblecrank_status in ("engaged", "disengaged")
        # Ensure that they're both missing at first so tests can fail properly
        temp_pkg_dir() do project_path
            copy_test_package(project_path, "AugmentedPlatform")
            ap_path = joinpath(project_path, "AugmentedPlatform")
            generate_flooblegrank_artifacts(ap_path)

            Pkg.activate(ap_path)

            @test !isdir(artifact_path(engaged_hash))
            @test !isdir(artifact_path(disengaged_hash))

            if flooblecrank_status == "engaged"
                right_hash = engaged_hash
                wrong_hash = disengaged_hash
            else
                right_hash = disengaged_hash
                wrong_hash = engaged_hash
            end

            withenv("FLOOBLECRANK" => flooblecrank_status) do
                add_this_pkg()
                @test isdir(artifact_path(right_hash))
                @test !isdir(artifact_path(wrong_hash))
            end

            # Manual test that artifact is installed by instantiate()
            artifacts_toml = joinpath(ap_path, "Artifacts.toml")
            p = HostPlatform()
            p["flooblecrank"] = flooblecrank_status
            flooblecrank_hash = artifact_hash("gooblebox", artifacts_toml; platform=p)
            @test flooblecrank_hash == right_hash
            @test artifact_exists(flooblecrank_hash)

            # Test that if we load the package, it knows how to find its own artifact,
            # because it feeds the right `Platform` object through to `@artifact_str()`
            cmd = setenv(`$(Base.julia_cmd()) --color=yes --project=$(ap_path) -e 'using AugmentedPlatform; print(get_artifact_dir("gooblebox"))'`,
                         "JULIA_DEPOT_PATH" => join(Base.DEPOT_PATH, Sys.iswindows() ? ";" : ":"),
                         "FLOOBLECRANK" => flooblecrank_status)
            using_output = chomp(String(read(cmd)))
            @test success(cmd)
            @test artifact_path(right_hash) == using_output

            mkpath("/tmp/foo/$(flooblecrank_status)")
            rm("/tmp/foo/$(flooblecrank_status)"; recursive=true, force=true)
            cp(project_path, "/tmp/foo/$(flooblecrank_status)")
            cp(Base.DEPOT_PATH[1], "/tmp/foo/$(flooblecrank_status)/depot")
        end
    end

    # Also run a test of "cross-installation"
    temp_pkg_dir() do project_path
        copy_test_package(project_path, "AugmentedPlatform")
        ap_path = joinpath(project_path, "AugmentedPlatform")
        generate_flooblegrank_artifacts(ap_path)

        Pkg.activate(ap_path)
        @test !isdir(artifact_path(engaged_hash))
        @test !isdir(artifact_path(disengaged_hash))

        # Instantiate with the environment variable set, but with an explicit
        # tag set in the platform object, which overrides.
        withenv("FLOOBLECRANK" => "disengaged") do
            p = HostPlatform()
            p["flooblecrank"] = "engaged"
            add_this_pkg(; platform=p)
            @test isdir(artifact_path(engaged_hash))
            @test !isdir(artifact_path(disengaged_hash))
        end
    end
end

@testset "Artifact GC collect delay" begin
    temp_pkg_dir() do tmpdir
        live_hash = create_artifact_chmod() do path
            open(joinpath(path, "README.md"), "w") do io
                print(io, "I will not go quietly into that dark night.")
            end
            open(joinpath(path, "binary.data"), "w") do io
                write(io, rand(UInt8, 1024))
            end
        end
        die_hash = create_artifact_chmod() do path
            open(joinpath(path, "README.md"), "w") do io
                print(io, "Let me sleep!")
            end
            open(joinpath(path, "binary.data"), "w") do io
                write(io, rand(UInt8, 1024))
            end
        end

        # We have created two separate artifacts
        @test live_hash != die_hash

        # Test that artifact_usage.toml does not exist
        usage_path = joinpath(Pkg.logdir(), "artifact_usage.toml")
        @test !isfile(usage_path)

        # We bind them here and now
        artifacts_toml = joinpath(tmpdir, "Artifacts.toml")
        bind_artifact!(artifacts_toml, "live", live_hash)
        bind_artifact!(artifacts_toml, "die", die_hash)

        # Now test that the usage file exists, and contains our Artifacts.toml
        usage = TOML.parsefile(usage_path)
        @test any(x -> startswith(x, artifacts_toml), keys(usage))

        # Test that a gc() doesn't remove anything
        @test artifact_exists(live_hash)
        @test artifact_exists(die_hash)
        Pkg.gc()
        @test artifact_exists(live_hash)
        @test artifact_exists(die_hash)

        # Test that unbinding the `die_hash` and running `gc()` again still doesn't
        # remove it, but it does add it to the orphan list
        unbind_artifact!(artifacts_toml, "die")
        Pkg.gc()
        @test artifact_exists(live_hash)
        @test artifact_exists(die_hash)

        orphaned_path = joinpath(Pkg.logdir(), "orphaned.toml")
        orphanage = TOML.parsefile(orphaned_path)
        @test any(x -> startswith(x, artifact_path(die_hash)), keys(orphanage))

        # Now, sleep for 0.2 seconds, then gc with a collect delay of 0.1 seconds
        # This should reap the `die_hash` immediately, as it has already been moved to
        # the orphaned list.
        sleep(0.2)
        Pkg.gc(;collect_delay=Millisecond(100))
        @test artifact_exists(live_hash)
        @test !artifact_exists(die_hash)

        # die_hash should still be listed within the orphan list, but one more gc() will
        # remove it; this is intentional and allows for robust removal scheduling.
        orphanage = TOML.parsefile(orphaned_path)
        @test any(x -> startswith(x, artifact_path(die_hash)), keys(orphanage))
        Pkg.gc()
        orphanage = TOML.parsefile(orphaned_path)
        @test !any(x -> startswith(x, artifact_path(die_hash)), keys(orphanage))

        # Next, unbind the live_hash, then run with collect_delay=0, and ensure that
        # things are cleaned up immediately.
        unbind_artifact!(artifacts_toml, "live")
        Pkg.gc(;collect_delay=Second(0))
        @test !artifact_exists(live_hash)
        @test !artifact_exists(die_hash)
    end
end

@testset "Override.toml" begin
    # We are going to test artifact overrides by creating an overlapping set of depots,
    # each with some artifacts installed within, then checking via things like
    # `artifact_path()` to ensure that our overrides are actually working.

    mktempdir() do depot_container
        depot1 = joinpath(depot_container, "depot1")
        depot2 = joinpath(depot_container, "depot2")
        depot3 = joinpath(depot_container, "depot3")

        make_foo(dir) = open(io -> print(io, "foo"), joinpath(dir, "foo"), "w")
        make_bar(dir) = open(io -> print(io, "bar"), joinpath(dir, "bar"), "w")
        make_baz(dir) = open(io -> print(io, "baz"), joinpath(dir, "baz"), "w")

        foo_hash = SHA1("2f42e2c1c1afd4ef8c66a2aaba5d5e1baddcab33")
        bar_hash = SHA1("64d0b4f8d9c004b862b38c4acfbd74988226995c")
        baz_hash = SHA1("087d8c93bff2f2b05f016bcd6ec653c8def76568")

        # First, create artifacts in each depot, with some overlap
        with_artifacts_directory(joinpath(depot3, "artifacts")) do
            @test create_artifact_chmod(make_foo) == foo_hash
            @test create_artifact_chmod(make_bar) == bar_hash
            @test create_artifact_chmod(make_baz) == baz_hash
        end
        with_artifacts_directory(joinpath(depot2, "artifacts")) do
            @test create_artifact_chmod(make_bar) == bar_hash
            @test create_artifact_chmod(make_baz) == baz_hash
        end
        with_artifacts_directory(joinpath(depot1, "artifacts")) do
            @test create_artifact_chmod(make_baz) == baz_hash
        end

        # Next, set up our depot path, with `depot1` as the "innermost" depot.
        old_depot_path = DEPOT_PATH
        empty!(DEPOT_PATH)
        append!(DEPOT_PATH, [depot1, depot2, depot3])

        # First sanity check; does our depot path searching code actually work properly?
        @test startswith(artifact_path(foo_hash), depot3)
        @test startswith(artifact_path(bar_hash), depot2)
        @test startswith(artifact_path(baz_hash), depot1)

        # Our `test/test_packages/ArtifactOverrideLoading` package contains some artifacts
        # that will not load unless they are properly overridden
        aol_uuid = Base.UUID("7b879065-7f74-5fa4-bdd5-9b7a15df8941")

        # Create an arbitrary absolute path for `barty`
        barty_override_path = abspath(joinpath(depot_container, "a_wild_barty_appears"))
        mkpath(barty_override_path)

        # Next, let's start spitting out some Overrides.toml files.  We'll make
        # one in `depot2`, then eventually create one in `depot1` to override these
        # overrides!
        open(joinpath(depot2, "artifacts", "Overrides.toml"), "w") do io
            overrides = Dict(
                # Override baz_hash to point to `bar_hash`
                bytes2hex(baz_hash.bytes) => bytes2hex(bar_hash.bytes),

                # Override "ArtifactOverrideLoading.arty" to point to `bar_hash` as well.
                # Override "ArtifactOverrideLoading.barty" to point to a location on disk
                string(aol_uuid) => Dict(
                    "arty" => bytes2hex(bar_hash.bytes),
                    "barty" => barty_override_path,
                )
            )
            TOML.print(io, overrides)
        end

        # Force Pkg to reload what it knows about artifact overrides
        @inferred Union{Nothing,Dict{Symbol,Any}} Pkg.Artifacts.load_overrides(;force=true)

        # Verify that the hash-based override worked
        @test artifact_path(baz_hash) == artifact_path(bar_hash)
        @test !endswith(artifact_path(baz_hash), bytes2hex(baz_hash.bytes))

        # Verify that the name-based override worked; extract paths from module that
        # loads overridden package artifacts.
        Pkg.activate(depot_container) do
            copy_test_package(depot_container, "ArtifactOverrideLoading")
            git_init_and_commit(joinpath(depot_container, "ArtifactOverrideLoading"))
            add_this_pkg()
            Pkg.add(Pkg.Types.PackageSpec(
                name="ArtifactOverrideLoading",
                uuid=aol_uuid,
                path=joinpath(depot_container, "ArtifactOverrideLoading"),
            ))

            (arty_path, barty_path) = Core.eval(Module(:__anon__), quote
                using ArtifactOverrideLoading
                arty_path, barty_path
            end)

            @test arty_path == artifact_path(bar_hash)
            @test barty_path == barty_override_path
        end

        # Excellent.  Let's add another Overrides.toml, this time in `depot1`, to muck
        # with the two overrides we put in previously, as well as `foo_hash`.
        open(joinpath(depot1, "artifacts", "Overrides.toml"), "w") do io
            overrides = Dict(
                # Override `foo` to an absolute path, then remove all overrides on `baz`
                bytes2hex(foo_hash.bytes) => barty_override_path,
                bytes2hex(baz_hash.bytes) => "",

                # Override "ArtifactOverrideLoading.arty" to point to `barty_override_path` as well.
                string(aol_uuid) => Dict(
                    "arty" => barty_override_path,
                )
            )
            TOML.print(io, overrides)
        end

        # Force Pkg to reload what it knows about artifact overrides
        Pkg.Artifacts.load_overrides(;force=true)

        # Force Julia to re-load ArtifactOverrideLoading from scratch
        pkgid = Base.PkgId(aol_uuid, "ArtifactOverrideLoading")
        delete!(Base.loaded_modules, pkgid)

        # Verify that the hash-based overrides (and clears) worked
        @test artifact_path(foo_hash) == barty_override_path
        @test endswith(artifact_path(baz_hash), bytes2hex(baz_hash.bytes))

        # Verify that the name-based override worked; extract paths from module that
        # loads overridden package artifacts.
        Pkg.activate(depot_container) do
            (arty_path, barty_path) = Core.eval(Module(:__anon__), quote
                using ArtifactOverrideLoading
                arty_path, barty_path
            end)

            @test arty_path == barty_override_path
            @test barty_path == barty_override_path
        end

        # Finally, let's test some invalid overrides:
        function test_invalid_override(overrides::Dict, msg)
            open(joinpath(depot1, "artifacts", "Overrides.toml"), "w") do io
                TOML.print(io, overrides)
            end
            @test_logs (:error, msg) match_mode=:any Pkg.Artifacts.load_overrides(;force=true)
        end

        # Mapping to a non-absolute path or SHA1 hash
        test_invalid_override(
            Dict("0"^40 => "invalid override target"),
            r"must map to an absolute path or SHA1 hash!",
        )
        test_invalid_override(
            Dict("0"^41 => "1"^40),
            r"Invalid SHA1 hash",
        )
        test_invalid_override(
            Dict("invalid UUID" => Dict("0"^40 => "1"^40)),
            r"Invalid UUID",
        )
        test_invalid_override(
            Dict("0"^40 => ["not", "a", "string", "or", "dict"]),
            r"failed to parse entry",
        )
        empty!(DEPOT_PATH)
        append!(DEPOT_PATH, old_depot_path)
    end
end

end # module


text:
abstract type LinearProcess{T} <: ContinuousTimeProcess{T}
end

b(s, x, P::LinearProcess) = B(s, P)*x + β(s, P)
σ(s, x, P::LinearProcess) = σ(s, P)
a(s, P::LinearProcess) = outer(σ(s, P))

symmetrize(A) = (A + A')/2
#####################

mutable struct Ptilde{T} <: ContinuousTimeProcess{T}
    cs::CSpline{T}
    σ
    a
    Γ
    Ptilde{T}(cs, σ) where T = new(cs, σ, σ*σ', inv(σ*σ'))
end
b(t::Float64, x, P::Ptilde) = P.cs(t)

B(t::Float64, P::Ptilde) = 0.0
β(t::Float64, P::Ptilde) = P.cs(t)

mu(s, x, t, P::Ptilde) = x + integrate(P.cs, s, t)
σ(t, x, P::Ptilde) = P.σ
a(t, x, P::Ptilde) = P.a
a(t, P::Ptilde) = P.a
gamma(t, x, P::Ptilde) = P.Γ
constdiff(::Ptilde) = true

"""
    Ptilde(cs::CSpline, σ)

Affine diffusion ``dX = cs(t) dt + σdW``
with cs a cubic spline `::CSpline`.

"""
Ptilde(cs::CSpline{T}, σ) where {T} = Ptilde{T}(cs, σ)

function lp(s, x, t, y, P::Ptilde{T}) where T
    logpdfnormal(y - mu(s,x,t,P), (t-s)*P.a)
end

transitionprob(s,x,t,P::Ptilde) = Gaussian(mu(s,x,t,P), (t-s)*P.a)

function V(t, T, v, P::Ptilde)
    v - integrate(P.cs, t, T)
end

function dotV(t, T, v, P::Ptilde)
     P.cs(t)
end

function H(t, T, P::Ptilde)
    P.Γ/(T-t)
end
function H(t, T, P::Ptilde, x)
    P.Γ*x/(T-t)
end

mutable struct LinPro{S,T,U} <: ContinuousTimeProcess{T}
    B::S
    μ::T
    σ::U
    a::U
    Γ::U
    lambda::S # stationary covariance
    function LinPro(B::S, μ::T, σ::U) where {S,T,U}
        a = σ*σ'
        return new{S,T,U}(B, μ, σ, a, inv(a), symmetrize(lyap(B, a)))
    end
end

B(t, P::LinPro) = P.B
β(t, P::LinPro) = -P.B*P.μ
b(t::Float64, x, P::LinPro) = P.B*(x .- P.μ)
σ(t, x, P::LinPro) = P.σ
bderiv(t, x, P::LinPro) = P.B
σderiv(t, x, P::LinPro) = 0*(P.σ)

a(t, x, P::LinPro) = P.a
a(t, P::LinPro) = P.a
constdiff(::LinPro) = true

"""
    LinPro(B, μ::T, σ)

Linear diffusion ``dX = B(X - μ)dt + σdW``.

"""
LinPro

function lp(s, x, t, y, P::LinPro{T}) where T
    logpdfnormal(y - mu(s,x,t,P), K(s, t, P::LinPro))
end

transitionprob(s,x,t,P::LinPro) = Gaussian(mu(s,x,t,P), K(s, t, P::LinPro))

function Phi(t, T, P::LinPro)
    exp((T-t)*P.B)
end

function mu(t, x, T, P::LinPro)
    phi = exp((T-t)*P.B)
    phi*(x - P.μ) + P.μ
end

function K(t, T, P::LinPro)
    phi = exp((T-t)*P.B)
    P.lambda - phi*P.lambda*phi'
end

function H(t, T, P::LinPro, x)
     phim = exp(-(T-t)*P.B)
     (phim*P.lambda*phim'-P.lambda)\x
end
function H(t, T, P::LinPro)
     phim = exp(-(T-t)*P.B)
     inv(phim*P.lambda*phim'-P.lambda)
end

function V(t, T, v, P::LinPro)
    phim = exp(-(T-t)*P.B)
    phim*(v - P.μ) + P.μ
end

function dotV(t, T, v, P::LinPro)
    exp(-(T-t)*P.B)*P.B*(v - P.μ)
end

"""
    LinProBridge

Bridge process of `P::LinPro` with `μ == 0` conditional on ending in `v` at time `t`.
"""
struct LinProBridge{T,S<:LinPro} <: ContinuousTimeProcess{T}
    t::Float64  # end time
    v::T        # end point
    P::S
    LinProBridge(t, v::T, P::S) where {T,S<:LinPro} = !iszero(P.μ) ? throw(ArgumentError("μ ≠ 0")) : new{T,S}(t,v,P)
end

b(s, x, P::LinProBridge) = P.P.B * (x - P.P.μ) + P.P.a * H(s, P.t, P.P, V(s, P.t, P.v, P.P) - x)
B(s, P::LinProBridge) = P.P.B - P.P.a * H(s, P.t, P.P)

σ(t, x, P::LinProBridge) = P.P.σ
bderiv(t, x, P::LinProBridge) = B(t, P)
σderiv(t, x, P::LinProBridge) = 0*(P.P.σ)

a(t, x, P::LinProBridge) = P.P.a
a(t, P::LinProBridge) = P.P.a
constdiff(::LinProBridge) = true

Phi(t, P::LinProBridge) = P.P.lambda*exp((P.t - t)*P.P.B)' - exp(-(P.t - t)*P.P.B)*P.P.lambda

Phi(t, T, Ps::LinProBridge) = Phi(T, Ps)*inv(Phi(t, Ps))

mu(s, x, t, P::LinProBridge{Float64}) = inv(sinh(P.P.B*(P.t - s)))*( P.v*sinh(P.P.B*(t - s)) + x*sinh(P.P.B*(P.t - t)) )
K(s, t, P::LinProBridge{Float64}) = 2 * P.P.lambda * csch(-P.P.B*(P.t - s)) * sinh(-P.P.B*(t-s))*sinh(-P.P.B*(P.t - t))
# Wolfram alpha: integral_(1/e)^ gamma (sinh^2(-(π - gamma )))/(sinh^2(-(π - s))) ds = sinh(1/e - gamma ) sinh( gamma - π) (-csch(1/e - π))≈0.170692

#function mu(s, x, t, Ps::LinProBridge)
#    T = Ps.t
#    P = Ps.P
#    inv(P.a - expm(P.B*(T-s))*P.a*expm(P.B*(T-s))')*(P.a - expm(P.B*(T-t))*P.a*expm(P.B*(T-t))')*x
#end

transitionprob(s, x, t, P::LinProBridge{Float64}) = Gaussian(mu(s, x, t, P), K(s, t, P))

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

"""
    LinearAppr(tt, B, β, a)
"""
struct LinearAppr{R,S,T} <: ContinuousTimeProcess{T}
    tt::Vector{Float64}
    xx::Vector{T}
    B::Vector{S}
    b::Vector{T}
    Σ::Vector{R}
end

_b((i,s)::IndexedTime, x, P::LinearAppr) = P.B[i]*(x - P.xx[i]) + P.b[i]
B((i,s)::IndexedTime, P::LinearAppr) = P.B[i]
β((i,s)::IndexedTime, P::LinearAppr) = P.b[i] - P.B[i]*P.xx[i]
a((i,s)::IndexedTime, x, P::LinearAppr) = outer(P.Σ[i])
a((i,s)::IndexedTime, P::LinearAppr) = outer(P.Σ[i])
constdiff(::LinearAppr) = false

linearappr(Y, P) = LinearAppr(Y.tt, Y.yy, map((t,x) -> Bridge.bderiv(t, x, P), Y.tt, Y.yy), map((t,x) -> Bridge.b(t, x, P), Y.tt, Y.yy), map((t,x) -> Bridge.σ(t, x, P), Y.tt, Y.yy))
function linearappr!(Pt::LinearAppr, Y, P)
    Pt.tt[:] = Y.tt
    Pt.xx[:] = Y.yy
    Pt.B[:] = map((t,x) -> Bridge.bderiv(t, x, P), Y.tt, Y.yy)
    Pt.b[:] = map((t,x) -> Bridge.b(t, x, P), Y.tt, Y.yy)
    Pt.Σ[:] = map((t,x) -> Bridge.σ(t, x, P), Y.tt, Y.yy)
    Pt
end


text:
# Julia wrapper for header: /usr/include/libavcodec/vorbis_parser.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_vorbis_parse_init,
    av_vorbis_parse_free,
    av_vorbis_parse_frame_flags,
    av_vorbis_parse_frame,
    av_vorbis_parse_reset

function av_vorbis_parse_init(extradata, extradata_size::Integer)
    ccall((:av_vorbis_parse_init, libavcodec), Ptr{AVVorbisParseContext}, (Ptr{UInt8}, Cint), extradata, extradata_size)
end

function av_vorbis_parse_free(s)
    ccall((:av_vorbis_parse_free, libavcodec), Cvoid, (Ptr{Ptr{AVVorbisParseContext}},), s)
end

function av_vorbis_parse_frame_flags(s, buf, buf_size::Integer, flags)
    ccall((:av_vorbis_parse_frame_flags, libavcodec), Cint, (Ptr{AVVorbisParseContext}, Ptr{UInt8}, Cint, Ptr{Cint}), s, buf, buf_size, flags)
end

function av_vorbis_parse_frame(s, buf, buf_size::Integer)
    ccall((:av_vorbis_parse_frame, libavcodec), Cint, (Ptr{AVVorbisParseContext}, Ptr{UInt8}, Cint), s, buf, buf_size)
end

function av_vorbis_parse_reset(s)
    ccall((:av_vorbis_parse_reset, libavcodec), Cvoid, (Ptr{AVVorbisParseContext},), s)
end


text:
import Pkg; Pkg.add(Pkg.PackageSpec(url="https://github.com/JuliaComputing/JuliaAcademyData.jl"))
using JuliaAcademyData; activate("World of machine learning")

# # Character based RNN language model
# Based on http://karpathy.github.io/2015/05/21/rnn-effectiveness

#-

# ## Setup
# 1-Adds required packages to Julia.
# 2-Loads the data and a pretrained model.

cd(datapath("charlm"))
include("charlm.jl")

# ## Sample Data-1
# A random subset of the Shakespeare training data

LEN = 500
r = rand(1:length(shake_text)-LEN)
println(shake_text[r:r+LEN])
flush(stdout)

# ## Sample Generation-1
# Random Shakespeare style text generated by the model

generate(shake_model, shake_chars, LEN)

# ## Sample Data-2
# A random subset of the code in Julia base

r = rand(1:length(julia_text)-LEN)
println(julia_text[r:r+LEN])
flush(stdout)

# ## Sample Generation-2
# Random Julia code generated by the model

generate(julia_model, julia_chars, LEN)

#-

max_stars_repo_path: Courses/World of machine learning/charlm.jl

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
## This makes implementing operators simpler
# but overrided BandedMatrix directly is likely to be faster

abstract type TridiagonalOperator{T} <: Operator{T} end
abstract type SymTridiagonalOperator{T} <: TridiagonalOperator{T} end
abstract type BidiagonalOperator{T} <: TridiagonalOperator{T} end
abstract type DiagonalOperator{T} <: BidiagonalOperator{T} end
# override getindex

bandwidths(::TridiagonalOperator)=1,1
bandwidths(::BidiagonalOperator)=0,1
bandwidths(::DiagonalOperator)=0,0

Tridiagonal(J::TridiagonalOperator{T},n::Integer) where {T}=Tridiagonal(T[J.γ(k) for k=2:n],
T[J.α(k) for k=1:n],
T[J.β(k) for k=1:n-1])

function symmetrize(J::TridiagonalOperator{T},n::Integer) where T
    d=Array{T}(undef, n)
    d[1]=1
    for k=2:n
        d[k]=sqrt(J[k,k-1]/J[k-1,k])*d[k-1]
    end

   SymTridiagonal(
    T[J[k,k] for k=1:n],
    T[J[k,k+1]*d[k+1]/d[k] for k=1:n-1])
end

struct DiagIteratorOperator{IT,T} <: DiagonalOperator{T}
    iterator::IT
end

DiagIteratorOperator(it) = DiagIteratorOperator{typeof(it),eltype(it)}(it)

getindex(D::DiagIteratorOperator,k::Integer,j::Integer) =
    k==j ? D.iterator[k] : zero(eltype(D))

domainspace(D::DiagIteratorOperator) = ℓ⁰
rangespace(D::DiagIteratorOperator) = ℓ⁰


text:
using Missings

#=
General notes about handling missing data, zero MAD etc:
The idea is that bicor should switch to cor whenever it is feasible, it helps, and it is requested:
(1) if median is missing, the mean would be missing as well, so there's no point in switching to Pearson
(2) In the results, columns and rows corresponding to input with missing means/medians are missing'd out.
(3) The convention is that if zeroMAD is set to non-zero, it is the index of the column in which MAD is
    zero (plus one for C indexing)
=#

const RefUX = 0.5

rescale(x, ::Type{Val{true}}, lowQ, hiQ) = x * RefUX / lowQ
rescale(x, ::Type{Val{false}}, lowQ, hiQ) = x * RefUX / hiQ
rescale(x, ::Type{Val{missing}}, lowQ, hiQ) = missing
rescale(x, y, lowQ, hiQ) = rescale(x, Val{y}, lowQ, hiQ)
rescale(x, lowQ, hiQ) = rescale(x, x < 0, lowQ, hiQ)

function prepareColBicor(col::Vector, maxPOutliers::Float64, fallback::Int64, cosine::Bool)
    #=
    asymptCorr = 1.4826, qnorm75 = 0.6744898
    Note to self: asymptCorr * qnorm75 is very close to 1 and should equal 1 theoretically. Should
    probably leave them out completely.
    =#
    nr = length(col)
    clean_col = collect(AbstractFloat, skipmissing(col))
    nNAentries::UInt64 = nr - length(clean_col)
    NAmed::Bool = false
    zeroMAD::Bool = false
    res = copy(col)
    aux = copy(col)

    if fallback == 4
        res, nNAentries, NAmed = prepareColCor(col, cosine)
        return (res, nNAentries, NAmed, zeroMAD)
    end

    # Calculate the median of col
    med = median(clean_col)  # original one implicitly arrange NA to the end

    # Create a conditional copy of the median
    medX = (cosine)? 0.0 : med

    # calculate absolute deviations from the median
    zeroMAD = false
    NAmed = ismissing(med)
    if NAmed
        res = zeros(nr)
        zeroMAD = true
    else
        res = col .- medX
        aux = abs.(col .- med)

        # calculate mad, i.e. median absolute deviation
        mad = median(collect(AbstractFloat, skipmissing(aux)))  # original one implicitly arrange NA to the end

        # If mad is zero, value of fallback decides what is it we will do.
        if mad == 0.0
            zeroMAD = true
            if fallback == 1
                # Return after zeoring out results and setting the NAmed flag
                res = zeros(nr)
                NAmed = true
            elseif fallback == 2
                # Switch to Pearson correlation and return
                # Rprintf("mad is zero in a column. Switching to Pearson for this column.\n");
                res, nNAentries, NAmed = prepareColCor(col, cosine)
            elseif fallback == 3
                # Do nothing: the setting of *zeroMAD above is enough.
            end
            return res, nNAentries, NAmed, zeroMAD
        end

        # We now re-use aux to store a copy of the weights ux. To calculate them, first get (x-med)/(9*mad).
        # Rprintf("median: %6.4f, mad: %6.4f, cosine: %d\n", med, mad, cosine);
        aux = (col - med) / (9.0 * mad)

        # Get the low and high quantiles  of ux
        clean_aux = collect(AbstractFloat, skipmissing(aux))
        lowQ = quantile(clean_aux, maxPOutliers)
        hiQ = quantile(clean_aux, 1.0 - maxPOutliers)

        # Rprintf("prepareColBicor: lowQ=%f, hiQ = %f\n", lowQ, hiQ);
        # If the low quantile is below -1, rescale the aux (that serve as ux below)
        # such that the low quantile will fall at -1; similarly for the high quantile

        (lowQ > -RefUX) && (lowQ = -RefUX)
        (hiQ < RefUX) && (hiQ = RefUX)
        lowQ = abs(lowQ)
        aux = rescale.(aux, lowQ, hiQ)

        # Calculate the (1-ux^2)^2 * (x-median(x))
        sum::BigFloat = BigFloat(0.0)
        for k=1:nr
            if !ismissing(res[k])
                ux = (abs(aux[k]) > 1)? 1.0 : aux[k]  # sign of ux doesn't matter.
                sum += (res[k] * (1 - ux^2)^2)^2
            else
                res[k] = 0
            end
        end
        sum = √(sum)
        if sum == 0.0
            res = zeros(nr)
            NAmed = true
        else
            res = Float64.(res ./ sum)
        end
    end
    return res, nNAentries, NAmed, zeroMAD
end

function prepareColCor(x::Vector, cosine::Bool)
    nr = length(x)
    # res
    nNAentries::UInt64 = 0
    NAmean::Bool

    count::UInt64 = 0
    mean::BigFloat = BigFloat(0.0)
    sum::BigFloat = BigFloat(0.0)
    for k=1:nr
        if !ismissing(x[k])
            mean += x[k]
            sum += x[k]^2
            count += 1
        end
    end

    if count > 0
        NAmean = false
        nNAentries = nr - count
        mean = (cosine)? 0 : mean/count
        sum = √(sum - count * mean*mean)
        if sum > 0
            # Rprintf("sum: %Le\n", sum)
            for k=1:nr
                res[k] = (!ismissing(x[k]))? (x[k] - mean)/sum : 0.0
            end
        else
            # Rprintf("prepareColCor: have zero variance.\n")
            NAmean = true
            res = zeros(nr)
        end
    else
        NAmean = true
        nNAentries = nr
        res = zeros(nr)
    end
    return res, nNAentries, NAmean
end


text:
using SymmetrySnap
using Documenter

DocMeta.setdocmeta!(SymmetrySnap, :DocTestSetup, :(using SymmetrySnap); recursive=true)

makedocs(;
    modules=[SymmetrySnap],
    authors="Jay Spendlove, Brayden Becker, Jason Meziere, Lydia Serafin",
    repo="https://github.com/jmeziere/SymmetrySnap.jl/blob/{commit}{path}#{line}",
    sitename="SymmetrySnap.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://jmeziere.github.io/SymmetrySnap.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/jmeziere/SymmetrySnap.jl",
)


text:
@defcomp Discontinuity begin

    region = Index()
    y_year = Parameter(index=[time], unit="year")
    y_year_0 = Parameter(unit="year")

    rand_discontinuity = Parameter(unit="unitless", default=.5)

    irefeqdis_eqdiscimpact = Variable(index=[region], unit="%")
    wincf_weightsfactor_sea = Parameter(index=[region], unit="")
    wdis_gdplostdisc = Parameter(unit="%", default=3.)

    igdpeqdis_eqdiscimpact = Variable(index=[time,region], unit="%")
    rgdp_per_cap_NonMarketRemainGDP = Parameter(index=[time,region], unit="\$/person")
    GDP_per_cap_focus_0_FocusRegionEU = Parameter(unit="\$/person", default=34298.93698672955)
    ipow_incomeexponent = Parameter(unit="unitless", default=-0.13333333333333333)

    igdp_realizeddiscimpact = Variable(index=[time,region], unit="%")
    occurdis_occurrencedummy = Variable(index=[time], unit="unitless")
    expfdis_discdecay = Variable(index=[time], unit="unitless")

    distau_discontinuityexponent = Parameter(unit="unitless", default=20.)

    idis_lossfromdisc = Variable(index=[time], unit="degreeC")
    tdis_tolerabilitydisc = Parameter(unit="degreeC", default=1.5)
    rt_g_globaltemperature = Parameter(index=[time], unit="degreeC")
    pdis_probability = Parameter(unit="%/degreeC", default=20.)

    isatg_saturationmodification = Parameter(unit="unitless")
    isat_satdiscimpact = Variable(index=[time,region], unit="%")

    isat_per_cap_DiscImpactperCapinclSaturation = Variable(index=[time,region], unit="\$/person")
    rcons_per_cap_DiscRemainConsumption = Variable(index=[time, region], unit="\$/person")
    rcons_per_cap_NonMarketRemainConsumption = Parameter(index=[time, region], unit="\$/person")

    function run_timestep(p, v, d, t)

        v.idis_lossfromdisc[t] = max(0, p.rt_g_globaltemperature[t] - p.tdis_tolerabilitydisc)

        if is_first(t)
            if v.idis_lossfromdisc[t] * (p.pdis_probability / 100) > p.rand_discontinuity
                v.occurdis_occurrencedummy[t] = 1
            else
                v.occurdis_occurrencedummy[t] = 0
            end
            v.expfdis_discdecay[t] = exp(-(p.y_year[t] - p.y_year_0) / p.distau_discontinuityexponent)
        else
            if v.idis_lossfromdisc[t] * (p.pdis_probability / 100) > p.rand_discontinuity
                v.occurdis_occurrencedummy[t] = 1
            elseif v.occurdis_occurrencedummy[t - 1] == 1
                v.occurdis_occurrencedummy[t] = 1
            else
                v.occurdis_occurrencedummy[t] = 0
            end
            v.expfdis_discdecay[t] = exp(-(p.y_year[t] - p.y_year[t - 1]) / p.distau_discontinuityexponent)
        end

        for r in d.region
            v.irefeqdis_eqdiscimpact[r] = p.wincf_weightsfactor_sea[r] * p.wdis_gdplostdisc

            v.igdpeqdis_eqdiscimpact[t,r] = v.irefeqdis_eqdiscimpact[r] * (p.rgdp_per_cap_NonMarketRemainGDP[t,r] / p.GDP_per_cap_focus_0_FocusRegionEU)^p.ipow_incomeexponent

            if is_first(t)
                v.igdp_realizeddiscimpact[t,r] = v.occurdis_occurrencedummy[t] * (1 - v.expfdis_discdecay[t]) * v.igdpeqdis_eqdiscimpact[t,r]
            else
                v.igdp_realizeddiscimpact[t,r] = v.igdp_realizeddiscimpact[t - 1,r] + v.occurdis_occurrencedummy[t] * (1 - v.expfdis_discdecay[t]) * (v.igdpeqdis_eqdiscimpact[t,r] - v.igdp_realizeddiscimpact[t - 1,r])
            end

            if v.igdp_realizeddiscimpact[t,r] < p.isatg_saturationmodification
                v.isat_satdiscimpact[t,r] = v.igdp_realizeddiscimpact[t,r]
            else
                v.isat_satdiscimpact[t,r] = p.isatg_saturationmodification + (100 - p.isatg_saturationmodification) * ((v.igdp_realizeddiscimpact[t,r] - p.isatg_saturationmodification) / ((100 - p.isatg_saturationmodification) + (v.igdp_realizeddiscimpact[t,r] - p.isatg_saturationmodification)))
            end
            v.isat_per_cap_DiscImpactperCapinclSaturation[t,r] = (v.isat_satdiscimpact[t,r] / 100) * p.rgdp_per_cap_NonMarketRemainGDP[t,r]
            v.rcons_per_cap_DiscRemainConsumption[t,r] = p.rcons_per_cap_NonMarketRemainConsumption[t,r] - v.isat_per_cap_DiscImpactperCapinclSaturation[t,r]
        end
    end
end