text:
function load2eqNK()
    phi_H = 0.9
    phi_L = 1.5
    rho = 0.9

    parameters = ["rho","phi"]
    parameter_values = [rho, (phi_L, phi_H)]
    regime_transmatrix = [0.3 0.7; 0.3 0.7]
    vars = ["c","inflation"]
    shocks = ["e_c"]
    expressions = ["1/exp(c) = 1/exp(c(+1))*((exp(inflation))^(phi_R))/exp(inflation(+1))",
    "c = rho*c(-1) + e_c"]

    RSDSGEModel(vars,shocks,parameters,expressions,parameter_values,regime_transmatrix)
end

function testloading2eqNK()
    load2eqNK()
    return true
end

function loadFRWZ()
    beta = 0.9976
    mu = (0.005,0.00004)
    sigma = 0.0002
    alpha = 0.33
    delta = 0.025
    transmatrix = [0.9 0.1; 0.1 0.9]

    vars = ["c","k","eps"]
    parameters = ["betta","mu","sigma","alpha","delta"]
    parameter_values = [beta,mu,sigma,alpha,delta]

    shocks = ["e_eps"]

    expressions = ["1/c = betta*1/c(+1)*exp((mu_R + sigma*eps)/(alpha-1))*(alpha*exp(mu_RP + sigma*eps(+1))*k^(alpha-1)+1-delta)",
    "c + k*exp((mu_R + sigma*eps)/(1-alpha)) = exp(mu_R+sigma*eps)*(k(-1))^(alpha) + (1-delta)*k(-1)",
    "eps = e_eps"]

    return RSDSGEModel(vars,shocks,parameters,expressions,parameter_values,transmatrix,ones(Float64,3))
end

function testloadingRBCfromFRWZ()
    loadFRWZ()
    return true
end

function testloadingfromfile()
    RSDSGEModel("../examples/RBC.txt")
    return true
end

@test testloading2eqNK()
@test testloadingRBCfromFRWZ()
@test testloadingfromfile()

twoEQNK = load2eqNK()
FRWZ = loadFRWZ()


text:
export Material

module Material

abstract type AbstractMaterial{T}                       end
abstract type AbstractElement{T} <: AbstractMaterial{T} end

include("atom.jl")

end

import .Material

ext: jl


text:
function box_count(x,y,z)
    #= BOX COUNT algorithm, assumes a long time series of x,y,z that are normalized in the (0,1) cube. =#

    # for n in nvec: n = 1/ϵ where n is the number of boxes in one direction, ϵ is the scale
    nvec = 2.^(3:8)
    Nmax = zeros(nvec)

    # largest distance between consequtive points - to compare with box size ϵ
    #Δx = maximum(sqrt.((x[1:end-1]-x[2:end]).^2 + (y[1:end-1]-y[2:end]).^2 + (z[1:end-1]-z[2:end]).^2))

    Nt = length(x)

    for (ni,n) in enumerate(nvec)
        # preallocate
        B = zeros(Int,n,n,n)    # could also be set up as boolean, indicates whether a box has been visited or not
        N = zeros(Int,Nt+1)     # counter array, set up as array to also check for convergence behaviour if necessary

        for i = 1:Nt
            xi = Int(ceil(x[i]*n))  # allocate a point into its box via rounding
            yi = Int(ceil(y[i]*n))
            zi = Int(ceil(z[i]*n))
            N[i+1] = N[i]+(1-B[xi,yi,zi])   # only add to counter if box was previously not visited
            B[xi,yi,zi] = 1     # mark box as visited
        end

        Nmax[ni] = N[end]
    end

    # estimate slope via linear regression
    interscept,slope = linreg(-log10.(1./nvec),log10.(Nmax))

    return slope
end

max_stars_repo_name: milankl/misc

max_issues_repo_name: milankl/misc

max_forks_repo_name: milankl/misc


text:
@testset "collect" begin
  ZZ7 = ResidueRing(ZZ, 7)
  @test collect(ZZ7) == ZZ7.(0:6)
  ZZ7_fmpz = ResidueRing(ZZ, ZZ(7))
  @test collect(ZZ7_fmpz) == ZZ7_fmpz.(0:6)
  F7 = GF(7)
  @test collect(F7) == F7.(0:6)
end
@testset "#465" begin
  G, mG = unit_group(quo(ZZ, 8)[1]);
  H = collect(G)
  X = map(mG, H)
  @test H == map(x->preimage(mG, x), X)
end


text:
using Base
using MatlabCompat
using MatlabCompat.Io
using MatlabCompat.ImageTools
import MatlabCompat.MathTools
import MatlabCompat.MathTools: max
using Base.Test
using Images
using FixedPointNumbers
using Tk

const TEST_FOLDER = dirname(@__FILE__)
const JANUS_PATH = joinpath(TEST_FOLDER, "janus.m");
const TESTMAT_PATH = joinpath(TEST_FOLDER, "test.mat");

# basic test
println("foo test")
x = function ()
    return 1 == 1
end
@test x()

################################
#     ImageTools Module Tests
################################

# make a test image
array = [0.1 0.001 0 0; 0.1 0.401 0 0; 0.1 0.01 0 0; 0.1 0.901 0 0.6];
img = grayim(array);
properties(img)["IMcs"] = "Gray";

# test the fuctions
println("setting true tests results");
trueThreshold = 0.11100752480522075;
trueBWSum= 3;

# println("testing graythresh()");
# @test graythresh(img) == trueThreshold;

println("testing im2bw()");
@test sum(reinterpret(Float32, float32(im2bw(img, trueThreshold)))[:]) == trueBWSum;

println("testing imread()");
img2 = MatlabCompat.imread("http://matlabcompat.github.io/img/example.tif");
expected = nothing
if VERSION < v"0.4-"
    expected = "Image{Gray{UfixedBase{Uint8,8}},2,Array{Gray{UfixedBase{Uint8,8}},2}}";
else
    expected = "Images.Image{ColorTypes.Gray{FixedPointNumbers.UFixed{UInt8,8}},2,Array{ColorTypes.Gray{FixedPointNumbers.UFixed{UInt8,8}},2}}"
	end
@test "$(typeof(img2))" == expected;

################################
#     Support Module Tests
################################

println("testing im2mat()");
mat = im2mat(img)
@test "$(typeof(mat))" == "Array{Float32,2}";

println("testing mat2im()");
expected = nothing
if (VERSION < v"0.4-")
    expected = "Image{Float64,2,Array{Float64,2}}"
else
    expected = "Images.Image{Float64,2,Array{Float64,2}}"
end
@test "$(typeof(mat2im(img)))" == expected;

println("testing rosetta()");
@test "$(typeof(rosetta(JANUS_PATH)))" == "Array{UTF8String,1}";

################################
#     Io Module Tests
################################

println("testing load()");
println("$(typeof(MatlabCompat.load(TESTMAT_PATH)))")
@test "$(typeof(MatlabCompat.load(TESTMAT_PATH)))" == "Dict{ASCIIString,Any}";
# "$(typeof(A))" == "Dict{Any,Any}"

################################
#     StringTools Module Tests
################################
println("testing disp()");
@test disp("test") == true;

println("testing num2str()");
@test num2str(3) == "3";

println("testing strcat()");
@test strcat("1", "2", "3") == "123";

################################
#     MathTools Module Tests
################################
matrix = [1 2 3];

println("testing numel()");
@test numel(matrix[:]) == 3;

println("testing max()");

@test max(matrix) == maximum(matrix);


text:
"Utility functions"

function AR_precisionmat(γ, order)
    mW = 1e8*Matrix{Float64}(I, order, order)
    mW[1,1] = γ
    return mW
end

function shift_mat(dim)
    S = Matrix{Float64}(I, dim, dim)
    for i in dim:-1:2
           S[i,:] = S[i-1, :]
    end
    S[1, :] = zeros(dim)
    return S
end

function shift_vec(dim, pos=1)
    s = zeros(dim)
    s[pos] = 1
    return dim == 1 ? s[pos] : s
end

max_stars_repo_name: wmkouw/LARX

max_issues_repo_name: wmkouw/LARX

max_forks_repo_name: wmkouw/LARX


text:
using NMRData
using Test

@testset "NMRData.jl" begin
    # Write your tests here.
end


text:
struct AssemblerSparsityPattern{Tv,Ti}
    K::SparseMatrixCSC{Tv,Ti}
    f::Vector{Tv}
    permutation::Vector{Int}
    sorteddofs::Vector{Int}
end

"""
    start_assemble(K::SparseMatrixCSC, [f::Vector])

Create an `AssemblerSparsityPattern` from the sparsity pattern in `K`
and optionally a "force" vector.
"""
function start_assemble(K::SparseMatrixCSC, f = eltype(K)[])
    if !isempty(f)
        @assert size(K, 1) == length(f)
    end
    AssemblerSparsityPattern(K, f, Int[], Int[])
end

"""
    assemble_local_vector!(f, dofs, fe)
Assembles the element residual `fe` into the global residual vector `f`.
"""
Base.@propagate_inbounds function assemble_local_vector!(f::AbstractVector, dofs::AbstractVector{Int}, fe::AbstractVector)
    @boundscheck checkbounds(f, dofs)
    @inbounds for i in 1:length(dofs)
        f[dofs[i]] += fe[i]
    end
end

function assemble_local!(A::AssemblerSparsityPattern, dofs::AbstractVector, Ke::AbstractMatrix, fe::AbstractVector)
    @assert !isempty(A.f)
    @boundscheck checkbounds(A.K, dofs, dofs)
    @boundscheck checkbounds(A.f, dofs)
    @inbounds assemble_local_vector!(A.f, dofs, fe)
    @inbounds assemble_local_matrix!(A, dofs, Ke)
end

"""
    assemble!(A::AssemblerSparsityPattern, dofs2, Ke)

Assemble a local dense element matrix `Ke` into the sparse matrix
wrapped by the assembler `A`.
location given by lists of indices `dofs1` and `dofs2`.

# Example

```julia
using SparseArrays
sparsity_pattern = sparse([1. 0 1; 1 0 1; 1 1 1])
fill!(sparsity_pattern, 0)
assembler = FEMSparse.start_assemble(sparsity_pattern)
dofs = [1, 3]
Ke = [1.0 2.0; 3.0 4.0]
FEMSparse.FEMSparse.assemble_local_matrix!(assembler, dofs, Ke)
Matrix(sparsity_pattern)

# output
julia> Matrix(sparsity_pattern)
3×3 Array{Float64,2}:
 1.0  0.0  2.0
 0.0  0.0  0.0
 3.0  0.0  4.0
```
"""
Base.@propagate_inbounds function assemble_local_matrix!(A::AssemblerSparsityPattern, dofs::AbstractVector{Int}, Ke::AbstractMatrix)
    permutation = A.permutation
    sorteddofs = A.sorteddofs
    K = A.K

    @boundscheck checkbounds(K, dofs, dofs)
    resize!(permutation, length(dofs))
    resize!(sorteddofs, length(dofs))
    copyto!(sorteddofs, dofs)
    sortperm2!(sorteddofs, permutation)

    current_col = 1
    @inbounds for Kcol in sorteddofs
        maxlookups = length(dofs)
        current_idx = 1
        for r in nzrange(K, Kcol)
            Kerow = permutation[current_idx]
            if K.rowval[r] == dofs[Kerow]
                Kecol = permutation[current_col]
                K.nzval[r] += Ke[Kerow, Kecol]
                current_idx += 1
            end
            current_idx > maxlookups && break
        end
        if current_idx <= maxlookups
            error("some row indices were not found")
        end
        current_col += 1
    end
end

##################
# Sort utilities #
##################
# We bundle a few sorting utilities here because the ones
# in have some unacceptable overhead.

# Sorts B and stores the permutation in `ii`
function sortperm2!(B, ii)
   @inbounds for i = 1:length(B)
      ii[i] = i
   end
   quicksort!(B, ii)
   return
end

function quicksort!(A, order, i=1,j=length(A))
    @inbounds if j > i
        if  j - i <= 12
           # Insertion sort for small groups is faster than Quicksort
           insertionsort!(A, order, i, j)
           return A
        end

        pivot = A[div(i+j,2)]
        left, right = i, j
        while left <= right
            while A[left] < pivot
                left += 1
            end
            while A[right] > pivot
                right -= 1
            end
            if left <= right
                A[left], A[right] = A[right], A[left]
                order[left], order[right] = order[right], order[left]

                left += 1
                right -= 1
            end
        end  # left <= right

        quicksort!(A,order, i,   right)
        quicksort!(A,order, left,j)
    end  # j > i

    return A
end

function insertionsort!(A, order, ii=1, jj=length(A))
    @inbounds for i = ii+1 : jj
        j = i - 1
        temp  = A[i]
        itemp = order[i]

        while true
            if j == ii-1
                break
            end
            if A[j] <= temp
                break
            end
            A[j+1] = A[j]
            order[j+1] = order[j]
            j -= 1
        end

        A[j+1] = temp
        order[j+1] = itemp
    end  # i
    return
end


text:
using FHIRClientXML
using Test

import BSON
import BrokenRecord
import HTTP
import XMLDict

@testset "FHIRClientXML.jl" begin
    include("unit-tests.jl")
    include("recorded-tests.jl")
end


text:
include("../../examples/matern.jl")
include("../plotting/gnuplot_utils.jl")

const VGRID = range(0.25, 10.0, length=101) # to hit "near-integer" v, which is the hardest.
const XGRID = range(0.0,  30.0, length=201)[2:end]

# First and second derivatives with FD:
fdbesselkdv(v, x) = (besselk(v+1e-6, x) - besselk(v+1e-6, x))/1e-6
fdbesselkdvdv(v, x) = (besselk(v+2e-6, x) - 2*besselk(v+1e-6, x) + besselk(v, x))/1e-12

# First and second derivatives with AD:
adbesselkdv(v, x) = ForwardDiff.derivative(_v->BesselK._besselk(_v, x), v)
adbesselkdvdv(v, x) = ForwardDiff.derivative(_v->adbesselkdv(_v, x), v)

function assemble_matrix(fn, pts, p)
  out = Array{Float64}(undef, length(pts), length(pts))
  Threads.@threads for k in 1:length(pts)
    out[k,k] = p[1]
    @inbounds for j in 1:(k-1)
      out[j,k] = fn(pts[j], pts[k], p)
    end
  end
  Symmetric(out)
end

atolfun(tru, est) = isnan(est) ? NaN : (isinf(tru) ? 0.0 : abs(tru-est))
atolfun(tru_est)  = atolfun(tru_est[1], tru_est[2])
rtolfun(tru, est) = atolfun(tru, est)/abs(tru)
rtolfun(tru_est)  = rtolfun(tru_est[1], tru_est[2])


text:
function test_single_stage_sequential(in_memory)
    template_ed = get_template_nomin_ed_simulation()
    c_sys = PSB.build_system(PSITestSystems, "c_sys5_uc")
    models = SimulationModels([
        DecisionModel(template_ed, c_sys, name="ED", optimizer=ipopt_optimizer),
    ])
    test_sequence =
        SimulationSequence(models=models, ini_cond_chronology=InterProblemChronology())
    sim_single = Simulation(
        name="consecutive",
        steps=2,
        models=models,
        sequence=test_sequence,
        simulation_folder=mktempdir(cleanup=true),
    )
    build_out = build!(sim_single)
    @test build_out == PSI.BuildStatus.BUILT
    execute_out = execute!(sim_single, in_memory=in_memory)
    @test execute_out == PSI.RunStatus.SUCCESSFUL
end

@testset "Single stage sequential tests" begin
    for in_memory in (true, false)
        test_single_stage_sequential(in_memory)
    end
end

function test_2_stage_decision_models_with_feedforwards(in_memory)
    template_uc = get_template_basic_uc_simulation()
    template_ed = get_template_nomin_ed_simulation()
    set_device_model!(template_ed, InterruptibleLoad, StaticPowerLoad)
    set_device_model!(template_ed, HydroEnergyReservoir, HydroDispatchReservoirBudget)
    set_network_model!(template_uc, NetworkModel(
        CopperPlatePowerModel,
        # MILP "duals" not supported with free solvers
        # duals = [CopperPlateBalanceConstraint],
    ))
    set_network_model!(
        template_ed,
        NetworkModel(
            CopperPlatePowerModel,
            duals=[CopperPlateBalanceConstraint],
            use_slacks=true,
        ),
    )
    c_sys5_hy_uc = PSB.build_system(PSITestSystems, "c_sys5_hy_uc")
    c_sys5_hy_ed = PSB.build_system(PSITestSystems, "c_sys5_hy_ed")
    models = SimulationModels(
        decision_models=[
            DecisionModel(template_uc, c_sys5_hy_uc; name="UC", optimizer=GLPK_optimizer),
            DecisionModel(template_ed, c_sys5_hy_ed; name="ED", optimizer=ipopt_optimizer),
        ],
    )

    sequence = SimulationSequence(
        models=models,
        feedforwards=Dict(
            "ED" => [
                SemiContinuousFeedforward(
                    component_type=ThermalStandard,
                    source=OnVariable,
                    affected_values=[ActivePowerVariable],
                ),
                EnergyLimitFeedforward(
                    component_type=HydroEnergyReservoir,
                    source=ActivePowerVariable,
                    affected_values=[ActivePowerVariable],
                    number_of_periods=12,
                ),
            ],
        ),
        ini_cond_chronology=InterProblemChronology(),
    )
    sim = Simulation(
        name="no_cache",
        steps=2,
        models=models,
        sequence=sequence,
        simulation_folder=mktempdir(cleanup=true),
    )

    build_out = build!(sim; console_level=Logging.Error)
    @test build_out == PSI.BuildStatus.BUILT
    execute_out = execute!(sim, in_memory=in_memory)
    @test execute_out == PSI.RunStatus.SUCCESSFUL
end

@testset "2-Stage Decision Models with FeedForwards" begin
    for in_memory in (true, false)
        test_2_stage_decision_models_with_feedforwards(in_memory)
    end
end

function test_2_stages_with_storage_ems(in_memory)
    template_uc =
        get_template_hydro_st_uc(NetworkModel(CopperPlatePowerModel, use_slacks=true))
    template_ed =
        get_template_hydro_st_ed(NetworkModel(CopperPlatePowerModel, use_slacks=true))
    set_device_model!(template_ed, InterruptibleLoad, StaticPowerLoad)
    c_sys5_hy_uc = PSB.build_system(PSITestSystems, "c_sys5_hy_ems_uc")
    c_sys5_hy_ed = PSB.build_system(PSITestSystems, "c_sys5_hy_ems_ed")
    models = SimulationModels(
        decision_models=[
            DecisionModel(template_uc, c_sys5_hy_uc; name="UC", optimizer=GLPK_optimizer),
            DecisionModel(template_ed, c_sys5_hy_ed; name="ED", optimizer=GLPK_optimizer),
        ],
    )

    sequence_cache = SimulationSequence(
        models=models,
        feedforwards=Dict(
            "ED" => [
                SemiContinuousFeedforward(
                    component_type=ThermalStandard,
                    source=OnVariable,
                    affected_values=[ActivePowerVariable],
                ),
                EnergyLimitFeedforward(
                    component_type=HydroEnergyReservoir,
                    source=ActivePowerVariable,
                    affected_values=[ActivePowerVariable],
                    number_of_periods=12,
                ),
            ],
        ),
        ini_cond_chronology=InterProblemChronology(),
    )
    sim_cache = Simulation(
        name="cache",
        steps=2,
        models=models,
        sequence=sequence_cache,
        simulation_folder=mktempdir(cleanup=true),
    )
    build_out = build!(sim_cache)
    @test build_out == PSI.BuildStatus.BUILT
    execute_out = execute!(sim_cache, in_memory=in_memory)
    @test execute_out == PSI.RunStatus.SUCCESSFUL
end

@testset "Simulation with 2-Stages with Storage EMS" begin
    for in_memory in (true, false)
        test_2_stages_with_storage_ems(in_memory)
    end
end

@testset "Test Simulation Utils" begin
    template_uc = get_template_basic_uc_simulation()
    set_device_model!(template_uc, ThermalStandard, ThermalStandardUnitCommitment)
    set_network_model!(template_uc, NetworkModel(
        CopperPlatePowerModel,
        # MILP "duals" not supported with free solvers
        # duals = [CopperPlateBalanceConstraint],
    ))

    template_ed = get_template_nomin_ed_simulation(
        NetworkModel(
            CopperPlatePowerModel;
            # Added because of data issues
            use_slacks=true,
            duals=[CopperPlateBalanceConstraint],
        ),
    )
    set_device_model!(template_ed, InterruptibleLoad, StaticPowerLoad)
    set_device_model!(template_ed, HydroEnergyReservoir, HydroDispatchReservoirBudget)
    c_sys5_hy_uc = PSB.build_system(PSITestSystems, "c_sys5_hy_uc")
    c_sys5_hy_ed = PSB.build_system(PSITestSystems, "c_sys5_hy_ed")
    models = SimulationModels(
        decision_models=[
            DecisionModel(template_uc, c_sys5_hy_uc; name="UC", optimizer=GLPK_optimizer),
            DecisionModel(template_ed, c_sys5_hy_ed; name="ED", optimizer=ipopt_optimizer),
        ],
    )

    sequence = SimulationSequence(
        models=models,
        feedforwards=Dict(
            "ED" => [
                SemiContinuousFeedforward(
                    component_type=ThermalStandard,
                    source=OnVariable,
                    affected_values=[ActivePowerVariable],
                ),
                EnergyLimitFeedforward(
                    component_type=HydroEnergyReservoir,
                    source=ActivePowerVariable,
                    affected_values=[ActivePowerVariable],
                    number_of_periods=12,
                ),
            ],
        ),
        ini_cond_chronology=InterProblemChronology(),
    )
    sim = Simulation(
        name="aggregation",
        steps=2,
        models=models,
        sequence=sequence,
        simulation_folder=mktempdir(cleanup=false),
    )

    build_out = build!(sim; console_level=Logging.Error)
    @test build_out == PSI.BuildStatus.BUILT
    execute_out = execute!(sim)
    @test execute_out == PSI.RunStatus.SUCCESSFUL

    @testset "Verify simulation events" begin
        file = joinpath(PSI.get_simulation_dir(sim), "recorder", "simulation_status.log")
        @test isfile(file)
        events = PSI.list_simulation_events(
            PSI.InitialConditionUpdateEvent,
            PSI.get_simulation_dir(sim);
            step=1,
        )
        @test length(events) == 0
        events = PSI.list_simulation_events(
            PSI.InitialConditionUpdateEvent,
            PSI.get_simulation_dir(sim);
            step=2,
        )
        # This value needs to be checked
        @test length(events) == 20

        PSI.show_simulation_events(
            devnull,
            PSI.InitialConditionUpdateEvent,
            PSI.get_simulation_dir(sim),
            ;
            step=2,
        )
        events = PSI.list_simulation_events(
            PSI.InitialConditionUpdateEvent,
            PSI.get_simulation_dir(sim);
            step=1,
            model_name="UC",
        )
        @test length(events) == 0
        events = PSI.list_simulation_events(
            PSI.InitialConditionUpdateEvent,
            PSI.get_simulation_dir(sim),
            ;
            step=2,
            model_name="UC",
        )
        @test length(events) == 20
        PSI.show_simulation_events(
            devnull,
            PSI.InitialConditionUpdateEvent,
            PSI.get_simulation_dir(sim),
            ;
            step=2,
            model_name="UC",
        )
    end

    # @testset "Check Serialization - Deserialization of Sim" begin
    #     path = mktempdir()
    #     files_path = PSI.serialize_simulation(sim; path = path)
    #     deserialized_sim = Simulation(files_path, stage_info)
    #     build_out = build!(deserialized_sim)
    #     @test build_out == PSI.BuildStatus.BUILT
    #     for stage in values(PSI.get_stages(deserialized_sim))
    #         @test PSI.is_stage_built(stage)
    #     end
    # end

end

function test_3_stage_simulation_with_feedforwards(in_memory)
    sys_rts_da = PSB.build_system(PSITestSystems, "modified_RTS_GMLC_DA_sys")
    sys_rts_rt = PSB.build_system(PSITestSystems, "modified_RTS_GMLC_RT_sys")
    sys_rts_ha = deepcopy(sys_rts_rt)

    PSY.transform_single_time_series!(sys_rts_da, 36, Hour(24))
    PSY.transform_single_time_series!(sys_rts_ha, 24, Hour(1))
    PSY.transform_single_time_series!(sys_rts_rt, 12, Hour(1))

    template_uc = get_template_standard_uc_simulation()
    set_network_model!(template_uc, NetworkModel(CopperPlatePowerModel))
    template_ha = deepcopy(template_uc)
    # network slacks added because of data issues
    template_ed = get_thermal_dispatch_template_network(
        NetworkModel(CopperPlatePowerModel, use_slacks=true),
    )

    models = SimulationModels(
        decision_models=[
            DecisionModel(
                template_uc,
                sys_rts_da;
                name="UC",
                optimizer=HiGHS_optimizer,
                initialize_model=false,
            ),
            DecisionModel(
                template_ha,
                sys_rts_ha;
                name="HA",
                optimizer=HiGHS_optimizer,
                initialize_model=false,
            ),
            DecisionModel(
                template_ed,
                sys_rts_rt;
                name="ED",
                optimizer=HiGHS_optimizer,
                initialize_model=false,
            ),
        ],
    )

    sequence = SimulationSequence(
        models=models,
        feedforwards=Dict(
            "ED" => [
                SemiContinuousFeedforward(
                    component_type=ThermalStandard,
                    source=OnVariable,
                    affected_values=[ActivePowerVariable],
                ),
            ],
        ),
        ini_cond_chronology=InterProblemChronology(),
    )

    sim = Simulation(
        name="3stage_feedforward",
        steps=1,
        models=models,
        sequence=sequence,
        simulation_folder=mktempdir(cleanup=true),
    )
    build_out = build!(sim)
    @test build_out == PSI.BuildStatus.BUILT
    # execute_out = execute!(sim, in_memory = in_memory)
    # @test execute_out == PSI.RunStatus.SUCCESSFUL
end

@testset "Test 3 stage simulation with FeedForwards" begin
    for in_memory in (true, false)
        test_3_stage_simulation_with_feedforwards(in_memory)
    end
end

@testset "Test HybridSystem simulations" begin
    sys_uc = PSB.build_system(PSITestSystems, "c_sys5_hybrid_uc")

    template_uc = get_template_standard_uc_simulation()
    set_device_model!(template_uc, HybridSystem, BasicHybridDispatch)
    set_network_model!(template_uc, NetworkModel(CopperPlatePowerModel, use_slacks=true))

    models = SimulationModels(
        decision_models=[
            DecisionModel(
                template_uc,
                sys_uc;
                name="UC",
                optimizer=HiGHS_optimizer,
                initialize_model=false,
            ),
        ],
    )

    sequence =
        SimulationSequence(models=models, ini_cond_chronology=InterProblemChronology())

    sim = Simulation(
        name="hybrid_test",
        steps=2,
        models=models,
        sequence=sequence,
        simulation_folder=mktempdir(cleanup=true),
    )
    build_out = build!(sim)
    @test build_out == PSI.BuildStatus.BUILT
    @test execute!(sim) == PSI.RunStatus.SUCCESSFUL
end

@testset "Test HybridSystem simulations" begin
    sys_uc = PSB.build_system(PSITestSystems, "c_sys5_hybrid_uc")
    sys_ed = PSB.build_system(PSITestSystems, "c_sys5_hybrid_ed")

    template_uc = get_template_standard_uc_simulation()
    set_device_model!(template_uc, HybridSystem, BasicHybridDispatch)
    set_network_model!(template_uc, NetworkModel(CopperPlatePowerModel, use_slacks=true))
    template_ed = get_thermal_dispatch_template_network(
        NetworkModel(CopperPlatePowerModel, use_slacks=true),
    )
    set_device_model!(template_ed, HybridSystem, BasicHybridDispatch)

    models = SimulationModels(
        decision_models=[
            DecisionModel(
                template_uc,
                sys_uc;
                name="UC",
                optimizer=HiGHS_optimizer,
                initialize_model=false,
            ),
            DecisionModel(
                template_ed,
                sys_ed;
                name="ED",
                optimizer=HiGHS_optimizer,
                initialize_model=false,
            ),
        ],
    )

    sequence =
        SimulationSequence(models=models, ini_cond_chronology=InterProblemChronology())

    sim = Simulation(
        name="hybrid_test",
        steps=2,
        models=models,
        sequence=sequence,
        simulation_folder=mktempdir(cleanup=true),
    )
    build_out = build!(sim)
    @test build_out == PSI.BuildStatus.BUILT
    execute_out = execute!(sim)
    @test execute_out == PSI.RunStatus.SUCCESSFUL
end

@testset "UC with MarketBid Cost in ThermalGenerators simulations" begin
    template = get_thermal_dispatch_template_network(
        NetworkModel(CopperPlatePowerModel, use_slacks=true),
    )
    set_device_model!(template, DeviceModel(ThermalStandard, ThermalDispatchNoMin))
    set_device_model!(template, DeviceModel(ThermalMultiStart, ThermalBasicUnitCommitment))

    models = SimulationModels(
        decision_models=[
            DecisionModel(
                UnitCommitmentProblem,
                template,
                PSB.build_system(PSITestSystems, "c_market_bid_cost");
                optimizer=cbc_optimizer,
                initialize_model=false,
            ),
        ],
    )

    sequence =
        SimulationSequence(models=models, ini_cond_chronology=InterProblemChronology())

    sim = Simulation(
        name="pwl_cost_test",
        steps=2,
        models=models,
        sequence=sequence,
        simulation_folder=mktempdir(cleanup=true),
    )

    @test build!(sim) == PSI.BuildStatus.BUILT
    @test execute!(sim) == PSI.RunStatus.SUCCESSFUL
    # TODO: Add more testing of resulting values
end


text:
"""
    GeneralFieldModule

Module for general fields.
"""
module GeneralFieldModule

__precompile__(true)

using ..FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import ..FieldModule.AbstractField
import ..FieldModule.@add_Field_fields

"""
    GeneralField{T<:Number} <: AbstractField

General field, meaning the entities can be anything.
"""
mutable struct GeneralField{T<:Number} <: AbstractField
    @add_Field_fields()
end

"""
    GeneralField(data::FMat{T}=[]) where {T<:Number}

Constructor of general field.  The values of the field are given by the array
on input, `data`. This array needs to have as many rows as there are entities,
and as many columns as there are degrees of freedom per entities.
"""
function GeneralField(data::FMat{T}=[]) where {T<:Number}
    values = deepcopy(data)
    dofnums = 0*similar(values,FInt)
    is_fixed = similar(values,Bool)
    fill!(is_fixed, 0)
    fixed_values = zeros(T,size(values))
    nfreedofs = 0
    return GeneralField(values, dofnums, is_fixed, fixed_values, nfreedofs)
end

"""
    GeneralField(data::FVec{T}) where {T<:Number}

Constructor of general field.  The values of the field are given by the vector
on input, `data`. This vector needs to have as many rows as there are entities.
"""
function GeneralField(data::FVec{T}) where {T<:Number}
    return GeneralField(reshape(data, length(data), 1))
end

end


text:
# CodePoints iterator
#
# Copyright 2017-2020 Gandalf Software, Inc., Scott P. Jones
# Licensed under MIT License, see LICENSE.md

struct CodePoints{T<:AbstractString}
    xs::T
end

"""
    codepoints(str)

An iterator that generates the code points of a string

# Examples
```jldoctest
julia> a = str("abc\U1f596")

julia> collect(a)

julia> collect(codepoints(a))
```
"""
codepoints(xs) = CodePoints(xs)
eltype(::Type{<:CodePoints{S}}) where {S} = eltype(S)
length(it::CodePoints) = length(it.xs)
iterate(it::CodePoints, state=1) = iterate(it.xs, state)


text:
import HTTP, JSON2

"""
    Agent{}

Agent is an abstract type that represents a solution running on the cluster.
"""
abstract type Agent{} end

function deploy(agent::Agent)
end

"""
    Runtime(agent::Agent, openapispec::Dict{String, Any})

Runtime holds the common runtime configuration for agents.
"""
mutable struct Runtime
    agent::Agent
    openapispec::Dict{String, Any}
end

function run(runtime::Runtime)
    # should setup a basic server based on the openapi definition

    # the different types of agents should be responsible for their own update functions.
end

# use a plain `Dict` as a "data store"
const ANIMALS = Dict{Int, Animal}()

function processhandler(req::HTTP.Request)
    input = Dict{String, Any}()
    input = JSON2.read(IOBuffer(HTTP.payload(req)), input)
    # TODO: validate input off openapi spec

    # process

    output = Dict{String, Any}()

    # TODO: validate output off openapi spec
    return HTTP.Response(200, JSON2.write(output))
end

function errorhandler(req::HTTP.Request)
    error = Dict{String, Any}()
    error = JSON2.read(IOBuffer(HTTP.payload(req)), error)

    # process error

    return HTTP.Response(200)
end

mutable struct Stats
    kv::Dict{string, Any}
end

function statshandler(req::HTTP.Request)
    # get stats

    return HTTP.Response(200, JSON2.write(stats))
end

# define REST endpoints to dispatch to "service" functions
const AGENT_ROUTER = HTTP.Router()
HTTP.@register(AGENT_ROUTER, "POST", "/process", processhandler)
HTTP.@register(AGENT_ROUTER, "POST", "/error", errorhandler)
HTTP.@register(AGENT_ROUTER, "GET", "/stats", statshandler)

HTTP.serve(AGENT_ROUTER, Sockets.localhost, 8080)

max_stars_repo_name: pbarker/kubemeta

max_issues_repo_name: pbarker/kubemeta

max_forks_repo_name: pbarker/kubemeta


text:
import ACTRSimulators: draw_object!
"""
** Tracking **

- `width`: screen width
- `height`: screen height
- `scheduler`: event scheduler
- `screen`: visual objects on screen
- `canvas`: GTK canvas
- `window`: GTK window
- `visible`: GUI visible
- `speed`: real time speed

Function Signature

````julia
PVT(;n_trials=10, trial=1, lb=2.0, ub=10.0, width=600.0, height=600.0, scheduler=nothing,
    screen=Vector{VisualObject}(), window=nothing, canvas=nothing, visible=false, speed=1.0)
````
"""
@concrete mutable struct Tracking <: AbstractTask
    width::Float64
    height::Float64
    scheduler
    screen::Vector{VisualObject}
    canvas
    window
    dot
    cursor
    visible::Bool
    realtime::Bool
    speed::Float64
    Δθ::Float64
    λθ::Float64
end

Broadcast.broadcastable(x::Tracking) = Ref(x)

function Tracking(;
    width = 600.0,
    height = 600.0,
    scheduler = nothing,
    screen = Vector{VisualObject}(),
    window = nothing,
    canvas = nothing,
    dot = Dot(;width, height),
    cursor = Cursor(),
    visible = false,
    realtime = false,
    speed = 1.0,
    λθ = .01,
    revs_per_sec = 10,
    Δθ = (2 * π) * λθ * (1 / revs_per_sec)
    )
    visible ? ((canvas,window) = setup_window(width)) : nothing
    visible ? Gtk.showall(window) : nothing
    return Tracking(
        width,
        height,
        scheduler,
        screen,
        canvas,
        window,
        dot,
        cursor,
        visible,
        realtime,
        speed,
        Δθ,
        λθ
    )
end

mutable struct Dot
    x::Float64
    y::Float64
    c_x::Float64
    c_y::Float64
    r::Float64
    θ::Float64
end

function Dot(;
    width,
    height,
    r = (width + height) / 8,
    θ = rand(Uniform(0, 2 * π))
    )
    c_x = width / 2
    c_y = height / 2
    x,y = compute_position(c_x, c_y, r, θ)
    Dot(x, y, c_x, c_y, r, θ)
end

function compute_position!(dot)
    x,y = compute_position(dot.c_x, dot.c_y, dot.r, dot.θ)
    dot.x = x
    dot.y = y
end

compute_position(dot) = compute_position(dot.c_x, dot.c_y, dot.r, dot.θ)

function compute_position(c_x, c_y, r, θ)
    x = c_x + r * cos(θ)
    y = c_y + r * sin(θ)
    return x, y
end

function get_position(dot)
    return dot.x, dot.y
end

function update_angle!(dot, Δθ)
    dot.θ = mod(dot.θ + Δθ, 2 * π)
end

function update_dot!(task)
    update_angle!(task.dot, task.Δθ)
    compute_position!(task.dot)
end

function update_vo!(models, task)
    x,y = get_position(task.dot)
    move_vo!.(models, x, y)
end

function move_dot!(models, task)
    update_dot!(task)
    update_vo!(models, task)
    repaint!(task, models)
end

mutable struct Cursor
    x::Float64
    y::Float64
end

Cursor(;x=0.0, y=0.0) = Cursor(x, y)

function draw_object!(task, dot::Dot)
    draw_object!(task, "O", dot.x, dot.y)
end

function start!(task::Tracking, models)
    x,y = compute_position(task.dot)
    present_stimulus(task, models, task.dot)
    run_trial!(task, models)
end

function present_stimulus(task, models, dot)
    vo = VisualObject(;x = dot.x, y = dot.y)
    add_to_visicon!(models, vo; stuff=true)
    push!(task.screen, vo)
    task.visible ? repaint!(task, models) : nothing
end

function run_trial!(task, models)
    λ = task.λθ
    register!(task, move_dot!, every, λ, models, task)
end

function press_key!(task::Tracking, actr, key)
    println("press_key not implimented")
end

function repaint!(task::Tracking, models)
    clear!(task)
    draw_object!(task, task.dot)
    draw_attention!.(task, models)
end


text:
using Test, LinearAlgebra
using DiffEqSensitivity, SteadyStateDiffEq, DiffEqBase, NLsolve
using OrdinaryDiffEq
using ForwardDiff, Calculus
@testset "Adjoint sensitivities of steady state solver" begin
  function f!(du,u,p,t)
    du[1] = p[1] + p[2]*u[1]
    du[2] = p[3]*u[1] + p[4]*u[2]
  end

  function jac!(J,u,p,t) #df/dx
    J[1,1] = p[2]
    J[2,1] = p[3]
    J[1,2] = 0
    J[2,2] = p[4]
    nothing
  end

  function paramjac!(fp,u,p,t) #df/dp
    fp[1,1] = 1
    fp[2,1] = 0
    fp[1,2] = u[1]
    fp[2,2] = 0
    fp[1,3] = 0
    fp[2,3] = u[1]
    fp[1,4] = 0
    fp[2,4] = u[2]
    nothing
  end

  function dg!(out,u,p,t,i)
    (out.=-2.0.+u)
  end

  function g(u,p,t)
    sum((2.0.-u).^2)/2 + sum(p.^2)/2
  end

  u0 = zeros(2)
  p = [2.0,-2.0,1.0,-4.0]
  prob = SteadyStateProblem(f!,u0,p)
  abstol = 1e-10
  @testset "for p" begin
    println("Calculate adjoint sensitivities from Jacobians")

    sol_analytical = [-p[1]/p[2], p[1]*p[3]/(p[2]*p[4])]

    J = zeros(2,2)
    fp = zeros(2,4)
    gp = zeros(4)
    gx = zeros(1,2)
    delg_delp = copy(p)

    jac!(J,sol_analytical,p,nothing)
    dg!(vec(gx),sol_analytical,p,nothing,nothing)
    paramjac!(fp,sol_analytical,p,nothing)

    lambda = J' \gx'
    res_analytical = delg_delp'-lambda' * fp # = -gx*inv(J)*fp

    @info "Expected result" sol_analytical, res_analytical, delg_delp'-gx*inv(J)*fp

    @info "Calculate adjoint sensitivities from autodiff & numerical diff"
    function G(p)
      tmp_prob = remake(prob,u0=convert.(eltype(p),prob.u0),p=p)
      sol = solve(tmp_prob,
        SSRootfind(nlsolve = (f!,u0,abstol) -> (res=NLsolve.nlsolve(f!,u0,autodiff=:forward,method=:newton,iterations=Int(1e6),ftol=1e-14);res.zero))
        )
      A = convert(Array,sol)
      g(A,p,nothing)
    end
    res1 = ForwardDiff.gradient(G,p)
    res2 = Calculus.gradient(G,p)
    #@info res1, res2, res_analytical

    @test res1 ≈ res_analytical' rtol = 1e-7
    @test res2 ≈ res_analytical' rtol = 1e-7
    @test res1 ≈ res2 rtol = 1e-7

    @info "Adjoint sensitivities"

    # with jac, param_jac
    f1 = ODEFunction(f!;jac=jac!, paramjac=paramjac!)
    prob1 = SteadyStateProblem(f1,u0,p)
    sol1 = solve(prob1,DynamicSS(Rodas5()))
    res1a = adjoint_sensitivities(sol1,DynamicSS(Rodas5()),sensealg=SteadyStateAdjoint(),g,dg!)
    res1b = adjoint_sensitivities(sol1,DynamicSS(Rodas5()),sensealg=SteadyStateAdjoint(),g,nothing)
    res1c = adjoint_sensitivities(sol1,DynamicSS(Rodas5()),sensealg=SteadyStateAdjoint(autodiff=false),g,nothing)
    res1d = adjoint_sensitivities(sol1,DynamicSS(Rodas5()),sensealg=SteadyStateAdjoint(autojacvec=TrackerVJP()),g,nothing)
    res1e = adjoint_sensitivities(sol1,DynamicSS(Rodas5()),sensealg=SteadyStateAdjoint(autojacvec=ReverseDiffVJP()),g,nothing)
    res1f = adjoint_sensitivities(sol1,DynamicSS(Rodas5()),sensealg=SteadyStateAdjoint(autojacvec=ZygoteVJP()),g,nothing)
    res1g = adjoint_sensitivities(sol1,DynamicSS(Rodas5()),sensealg=SteadyStateAdjoint(autodiff=false,autojacvec=false),g,nothing)

    # with jac, without param_jac
    f2 = ODEFunction(f!;jac=jac!)
    prob2 = SteadyStateProblem(f2,u0,p)
    sol2 = solve(prob2,DynamicSS(Rodas5()))
    res2a = adjoint_sensitivities(sol2,DynamicSS(Rodas5()),sensealg=SteadyStateAdjoint(),g,dg!)
    res2b = adjoint_sensitivities(sol2,DynamicSS(Rodas5()),sensealg=SteadyStateAdjoint(),g,nothing)
    res2c = adjoint_sensitivities(sol2,DynamicSS(Rodas5()),sensealg=SteadyStateAdjoint(autodiff=false),g,nothing)
    res2d = adjoint_sensitivities(sol2,DynamicSS(Rodas5()),sensealg=SteadyStateAdjoint(autojacvec=TrackerVJP()),g,nothing)
    res2e = adjoint_sensitivities(sol2,DynamicSS(Rodas5()),sensealg=SteadyStateAdjoint(autojacvec=ReverseDiffVJP()),g,nothing)
    res2f = adjoint_sensitivities(sol2,DynamicSS(Rodas5()),sensealg=SteadyStateAdjoint(autojacvec=ZygoteVJP()),g,nothing)
    res2g = adjoint_sensitivities(sol2,DynamicSS(Rodas5()),sensealg=SteadyStateAdjoint(autodiff=false,autojacvec=false),g,nothing)

    # without jac, without param_jac
    f3 = ODEFunction(f!)
    prob3 = SteadyStateProblem(f3,u0,p)
    sol3 = solve(prob3,DynamicSS(Rodas5()))
    res3a = adjoint_sensitivities(sol3,DynamicSS(Rodas5()),sensealg=SteadyStateAdjoint(),g,dg!)
    res3b = adjoint_sensitivities(sol3,DynamicSS(Rodas5()),sensealg=SteadyStateAdjoint(),g,nothing)
    res3c = adjoint_sensitivities(sol3,DynamicSS(Rodas5()),sensealg=SteadyStateAdjoint(autodiff=false),g,nothing)
    res3d = adjoint_sensitivities(sol3,DynamicSS(Rodas5()),sensealg=SteadyStateAdjoint(autojacvec=TrackerVJP()),g,nothing)
    res3e = adjoint_sensitivities(sol3,DynamicSS(Rodas5()),sensealg=SteadyStateAdjoint(autojacvec=ReverseDiffVJP()),g,nothing)
    res3f = adjoint_sensitivities(sol3,DynamicSS(Rodas5()),sensealg=SteadyStateAdjoint(autojacvec=ZygoteVJP()),g,nothing)
    res3g = adjoint_sensitivities(sol3,DynamicSS(Rodas5()),sensealg=SteadyStateAdjoint(autodiff=false,autojacvec=false),g,nothing)

    @test norm(res_analytical' .- res1a) < 1e-7
    @test norm(res_analytical' .- res1b) < 1e-7
    @test norm(res_analytical' .- res1c) < 1e-7
    @test norm(res_analytical' .- res1d) < 1e-7
    @test norm(res_analytical' .- res1e) < 1e-7
    @test_broken norm(res_analytical' .- res1f) < 1e-7
    @test norm(res_analytical' .- res1g) < 1e-7
    @test norm(res_analytical' .- res2a) < 1e-7
    @test norm(res_analytical' .- res2b) < 1e-7
    @test norm(res_analytical' .- res2c) < 1e-7
    @test norm(res_analytical' .- res2d) < 1e-7
    @test norm(res_analytical' .- res2e) < 1e-7
    @test_broken norm(res_analytical' .- res2f) < 1e-7
    @test norm(res_analytical' .- res2g) < 1e-7
    @test norm(res_analytical' .- res3a) < 1e-7
    @test norm(res_analytical' .- res3b) < 1e-7
    @test norm(res_analytical' .- res3c) < 1e-7
    @test norm(res_analytical' .- res3d) < 1e-7
    @test norm(res_analytical' .- res3e) < 1e-7
    @test_broken norm(res_analytical' .- res3f) < 1e-7
    @test norm(res_analytical' .- res3g) < 1e-7

    @info "oop checks"
    function foop(u,p,t)
      dx = p[1] + p[2]*u[1]
      dy = p[3]*u[1] + p[4]*u[2]
      [dx,dy]
    end
    proboop = SteadyStateProblem(foop,u0,p)
    soloop = solve(proboop,DynamicSS(Rodas5()))

    res4a = adjoint_sensitivities(soloop,DynamicSS(Rodas5()),sensealg=SteadyStateAdjoint(),g,dg!)
    res4b = adjoint_sensitivities(soloop,DynamicSS(Rodas5()),sensealg=SteadyStateAdjoint(),g,nothing)
    res4c = adjoint_sensitivities(soloop,DynamicSS(Rodas5()),sensealg=SteadyStateAdjoint(autodiff=false),g,nothing)
    res4d = adjoint_sensitivities(soloop,DynamicSS(Rodas5()),sensealg=SteadyStateAdjoint(autojacvec=TrackerVJP()),g,nothing)
    res4e = adjoint_sensitivities(soloop,DynamicSS(Rodas5()),sensealg=SteadyStateAdjoint(autojacvec=ReverseDiffVJP()),g,nothing)
    res4f = adjoint_sensitivities(soloop,DynamicSS(Rodas5()),sensealg=SteadyStateAdjoint(autojacvec=ZygoteVJP()),g,nothing)
    res4g = adjoint_sensitivities(soloop,DynamicSS(Rodas5()),sensealg=SteadyStateAdjoint(autodiff=false,autojacvec=false),g,nothing)
    res4h = adjoint_sensitivities(soloop,DynamicSS(Rodas5()),sensealg=SteadyStateAdjoint(autodiff=true,autojacvec=false),g,nothing)

    @test norm(res_analytical' .- res4a) < 1e-7
    @test norm(res_analytical' .- res4b) < 1e-7
    @test norm(res_analytical' .- res4c) < 1e-7
    @test norm(res_analytical' .- res4d) < 1e-7
    @test norm(res_analytical' .- res4e) < 1e-7
    @test norm(res_analytical' .- res4f) < 1e-7
    @test norm(res_analytical' .- res4g) < 1e-7
    @test norm(res_analytical' .- res4h) < 1e-7
  end

  @testset "for u0: (should be zero, steady state does not depend on initial condition)" begin
    res5 = ForwardDiff.gradient(prob.u0) do u0
      tmp_prob = remake(prob,u0=u0)
      sol = solve(tmp_prob,
        SSRootfind(nlsolve = (f!,u0,abstol) -> (res=NLsolve.nlsolve(f!,u0,autodiff=:forward,method=:newton,iterations=Int(1e6),ftol=1e-14);res.zero) )
        )
      A = convert(Array,sol)
      g(A,p,nothing)
    end
    @test abs(dot(res5,res5)) < 1e-7
  end
end

using Zygote
@testset "concrete_solve derivatives steady state solver" begin

  function g1(u,p,t)
    sum(u)
  end

  function g2(u,p,t)
    sum((2.0.-u).^2)/2
  end

  u0 = zeros(2)
  p = [2.0,-2.0,1.0,-4.0]

  @testset "iip" begin
    function f!(du,u,p,t)
      du[1] = p[1] + p[2]*u[1]
      du[2] = p[3]*u[1] + p[4]*u[2]
    end
    prob = SteadyStateProblem(f!,u0,p)

    sol = solve(prob,DynamicSS(Rodas5()))
    res1 = adjoint_sensitivities(sol,DynamicSS(Rodas5()),sensealg=SteadyStateAdjoint(),g1,nothing)
    res2 = adjoint_sensitivities(sol,DynamicSS(Rodas5()),sensealg=SteadyStateAdjoint(),g2,nothing)

    dp1 = Zygote.gradient(p->sum(solve(prob,DynamicSS(Rodas5()),u0=u0,p=p,sensealg=SteadyStateAdjoint())),p)
    dp2 = Zygote.gradient(p->sum((2.0.-solve(prob,DynamicSS(Rodas5()),u0=u0,p=p,sensealg=SteadyStateAdjoint())).^2)/2.0,p)

    @test res1 ≈ dp1[1] rtol=1e-12
    @test res2 ≈ dp2[1] rtol=1e-12

    res1 = Zygote.gradient(p->sum(Array(solve(prob,DynamicSS(Rodas5()),u0=u0,p=p,sensealg=SteadyStateAdjoint()))[1]),p)
    dp1 = Zygote.gradient(p->sum(solve(prob,DynamicSS(Rodas5()),u0=u0,p=p,save_idxs=1:1,sensealg=SteadyStateAdjoint())),p)
    dp2 = Zygote.gradient(p->solve(prob,DynamicSS(Rodas5()),u0=u0,p=p,save_idxs=1,sensealg=SteadyStateAdjoint())[1],p)
    @test res1[1] ≈ dp1[1] rtol=1e-10
    @test res1[1] ≈ dp2[1] rtol=1e-10
  end

  @testset "oop" begin
    function f(u,p,t)
      dx = p[1] + p[2]*u[1]
      dy = p[3]*u[1] + p[4]*u[2]
      [dx,dy]
    end
    proboop = SteadyStateProblem(f,u0,p)

    soloop = solve(proboop,DynamicSS(Rodas5()))
    res1oop = adjoint_sensitivities(soloop,DynamicSS(Rodas5()),sensealg=SteadyStateAdjoint(),g1,nothing)
    res2oop = adjoint_sensitivities(soloop,DynamicSS(Rodas5()),sensealg=SteadyStateAdjoint(),g2,nothing)

    dp1oop = Zygote.gradient(p->sum(solve(proboop,DynamicSS(Rodas5()),u0=u0,p=p,sensealg=SteadyStateAdjoint())),p)
    dp2oop = Zygote.gradient(p->sum((2.0.-solve(proboop,DynamicSS(Rodas5()),u0=u0,p=p,sensealg=SteadyStateAdjoint())).^2)/2.0,p)

    @test res1oop ≈ dp1oop[1] rtol=1e-12
    @test res2oop ≈ dp2oop[1] rtol=1e-12

    res1oop = Zygote.gradient(p->sum(Array(solve(proboop,DynamicSS(Rodas5()),u0=u0,p=p,sensealg=SteadyStateAdjoint()))[1]),p)
    dp1oop = Zygote.gradient(p->sum(solve(proboop,DynamicSS(Rodas5()),u0=u0,p=p,save_idxs=1:1,sensealg=SteadyStateAdjoint())),p)
    dp2oop = Zygote.gradient(p->solve(proboop,DynamicSS(Rodas5()),u0=u0,p=p,save_idxs=1,sensealg=SteadyStateAdjoint())[1],p)
    @test res1oop[1] ≈ dp1oop[1] rtol=1e-10
    @test res1oop[1] ≈ dp2oop[1] rtol=1e-10
  end
end


text:
#include("ar/test.jl")
#include("spline/test.jl")

# TODO: Reexport Distribution
using Filtering, Distributions, Test, Random, StatsBase

include("interface.jl")

include("particle/ar.jl")
include("particle/score.jl")

include("regressions/ar.jl")


text:
export OfflinePolicy, JuliaRLTransition, gen_JuliaRL_dataset

export calculate_CQL_loss, maximum_mean_discrepancy_loss

struct JuliaRLTransition
    state
    action
    reward
    terminal
    next_state
end

Base.@kwdef struct OfflinePolicy{L,T} <: AbstractPolicy
    learner::L
    dataset::T
    continuous::Bool
    batch_size::Int
end

(π::OfflinePolicy)(env) = π(env, ActionStyle(env), action_space(env))

function (π::OfflinePolicy)(env, ::MinimalActionSet, ::Base.OneTo)
    if π.continuous
        π.learner(env)
    else
        findmax(π.learner(env))[2]
    end
end
(π::OfflinePolicy)(env, ::FullActionSet, ::Base.OneTo) = findmax(π.learner(env), legal_action_space_mask(env))[2]

function (π::OfflinePolicy)(env, ::MinimalActionSet, A)
    if π.continuous
        π.learner(env)
    else
        A[findmax(π.learner(env))[2]]
    end
end
(π::OfflinePolicy)(env, ::FullActionSet, A) = A[findmax(π.learner(env), legal_action_space_mask(env))[2]]

function RLBase.update!(
    p::OfflinePolicy,
    traj::AbstractTrajectory,
    ::AbstractEnv,
    ::PreExperimentStage,
)
    l = p.learner
    if in(:pretrain_step, fieldnames(typeof(l)))
        println("Pretrain...")
        for _ in 1:l.pretrain_step
            inds, batch = sample(l.rng, p.dataset, p.batch_size)
            update!(l, batch)
        end
    end
end

function RLBase.update!(
    p::OfflinePolicy,
    traj::AbstractTrajectory,
    ::AbstractEnv,
    ::PreActStage,
)
    l = p.learner
    l.update_step += 1

    if in(:target_update_freq, fieldnames(typeof(l))) && l.update_step % l.target_update_freq == 0
        copyto!(l.target_approximator, l.approximator)
    end

    l.update_step % l.update_freq == 0 || return

    inds, batch = sample(l.rng, p.dataset, p.batch_size)

    update!(l, batch)
end

"""
    gen_JuliaRL_dataset(alg::Symbol, env::Symbol, type::AbstractString; dataset_size)

Generate the dataset by trajectory from the trajectory obtained from the experiment (`alg` + `env`). `type` represents the method of collecting data. Possible values: random/medium/expert. `dataset_size` is the size of the generated dataset.
"""
function gen_JuliaRL_dataset(alg::Symbol, env::Symbol, type::AbstractString; dataset_size::Int)
    dataset_ex = Experiment(
            Val(:GenDataset),
            Val(alg),
            Val(env),
            type;
            dataset_size = dataset_size)

    run(dataset_ex)

    dataset = []
    s, a, r, t = dataset_ex.policy.trajectory.traces
    for i in 1:dataset_size
        push!(dataset, JuliaRLTransition(s[:, i], a[i], r[i], t[i], s[:, i+1]))
    end
    dataset
end

function StatsBase.sample(rng::AbstractRNG, dataset::Vector{T}, batch_size::Int) where {T}
    valid_range = 1:length(dataset)
    inds = rand(rng, valid_range, batch_size)
    batch_data = dataset[inds]
    s_length = size(batch_data[1].state)[1]
    a_type = typeof(batch_data[1].action)

    s = Array{Float32}(undef, s_length, batch_size)
    s′ = Array{Float32}(undef, s_length, batch_size)
    a = Array{a_type}(undef, batch_size)
    r = Array{Float32}(undef, batch_size)
    t = Array{Float32}(undef, batch_size)
    for (i, data) in enumerate(batch_data)
        s[:, i] = data.state
        a[i] = data.action
        s′[:, i] = data.next_state
        r[i] = data.reward
        t[i] = data.terminal
    end
    batch = NamedTuple{SARTS}((s, a, r, t, s′))
    inds, batch
end

"""
    calculate_CQL_loss(q_value, action; method)
See paper: [Conservative Q-Learning for Offline Reinforcement Learning](https://arxiv.org/abs/2006.04779)
"""
function calculate_CQL_loss(q_value::Matrix{T}, action::Vector{R}; method = "CQL(H)") where {T, R}
    if method == "CQL(H)"
        cql_loss = mean(log.(sum(exp.(q_value), dims=1)) .- q_value[action])
    else
        @error Wrong method parameter
    end
    return cql_loss
end

function maximum_mean_discrepancy_loss(raw_sample_action, raw_actor_action, type::Symbol, mmd_σ::Float32=10.0f0)
    A, B, N = size(raw_sample_action)
    diff_xx = reshape(raw_sample_action, A, B, N, 1) .- reshape(raw_sample_action, A, B, 1, N)
    diff_xy = reshape(raw_sample_action, A, B, N, 1) .- reshape(raw_actor_action, A, B, 1, N)
    diff_yy = reshape(raw_actor_action, A, B, N, 1) .- reshape(raw_actor_action, A, B, 1, N)
    diff_xx = calculate_sample_distance(diff_xx, type, mmd_σ)
    diff_xy = calculate_sample_distance(diff_xy, type, mmd_σ)
    diff_yy = calculate_sample_distance(diff_yy, type, mmd_σ)
    mmd_loss = sqrt.(diff_xx .+ diff_yy .- 2.0f0 .* diff_xy .+ 1.0f-6)
end

function calculate_sample_distance(diff, type::Symbol, mmd_σ::Float32)
    if type == :gaussian
        diff = diff .^ 2
    elseif type == :laplacian
        diff = abs.(diff)
    else
        error("Wrong parameter.")
    end
    return vec(mean(exp.(-sum(diff, dims=1) ./ (2.0f0 * mmd_σ)), dims=(3, 4)))
end


text:
import LinearAlgebra: diagind
"""
    euclidean_distance(out, A, B; tmp_m=nothing, tmp_n=nothing)

Computes all pairwise distances between two sets of data, A(p x n) and B(p x m).
temp memory: length n, length m
Output: n x m.
"""
function euclidean_distance!(out::AbstractArray, A::AbstractMatrix, B::AbstractMatrix; tmp_n::Union{AbstractArray, Nothing}=nothing, tmp_m::Union{AbstractArray, Nothing}=nothing)
    AT = typeof(B) # this choice is intentional due to the method below.
    ATN =  hasproperty(AT, :name) ? AT.name.wrapper : A

    p, n = size(A)
    p2, m = size(B)
    @assert p == p2
    @assert size(out) == (n, m)
    T = eltype(A)

    tmp_n = (tmp_n == nothing) ? ATN{T}(undef, n) : tmp_n
    @assert length(tmp_n) == n
    tmp_m = (tmp_m == nothing) ? ATN{T}(undef, m) : tmp_m
    @assert length(tmp_m) == m

    tmp_n .= reshape(sum(A.^2; dims=1), n)
    tmp_m .= reshape(sum(B.^2; dims=1), m)

    LinearAlgebra.mul!(out, transpose(A), B)
    out .= sqrt.(max.((-2out .+ transpose(tmp_m)) .+ tmp_n, zero(T)))

    if A === B
        out[diagind(out)] .= zero(T)
    end
    out
end

"""
    euclidean_distance!(out, A)

Computes all pairwise distances between data in A (p x [n]). Output: n x [n].
"""
function euclidean_distance!(out::MPIMatrix{T,A}, data::MPIMatrix{T,A};
                             tmp_data::Union{A, Nothing}=nothing, tmp_dist::Union{A, Nothing}=nothing,
                             tmp_vec1::Union{A,Nothing}=nothing, tmp_vec2::Union{A,Nothing}=nothing) where {T,A}
    p, n = size(data)
    @assert size(out) == (n, n)
    local_len = n ÷ Size()
    remainder = n % Size()

    local_len_p = remainder > 0 ? local_len + 1 : local_len

    tmp_data = (tmp_data == nothing) ? A{T}(undef, p*(local_len_p)) : tmp_data
    @assert length(tmp_data) >= p * local_len_p
    tmp_dist = (tmp_dist == nothing) ? A{T}(undef, local_len_p^2) : tmp_dist
    @assert length(tmp_dist) >= local_len_p^2
    tmp_vec1 = (tmp_vec1 == nothing) ? A{T}(undef, local_len_p) : tmp_vec1
    @assert length(tmp_vec1) >= local_len_p
    tmp_vec2 = (tmp_vec2 == nothing) ? A{T}(undef, local_len_p) : tmp_vec2
    @assert length(tmp_vec2) >= local_len_p

    for r in 0:Size()-1
        this = data.localarray
        if r == Rank()
            other = data.localarray
            sync()
            Bcast!(reshape(data.localarray, :); root=r)
        else
            other = r < remainder ? @view(tmp_data[1:(p*(local_len+1))]) : @view(tmp_data[1:p*local_len])
            sync()
            Bcast!(other; root=r)
            other = reshape(other, p, :)
        end

        tmp_dist_cols = (Rank() < remainder) ? local_len + 1 : local_len
        tmp_dist_rows = (r < remainder) ? local_len + 1 : local_len
        tmp_dist_view = reshape(@view(tmp_dist[1:(tmp_dist_rows * tmp_dist_cols)]),
                                      tmp_dist_rows, tmp_dist_cols)
        euclidean_distance!(tmp_dist_view, other, this;
                            tmp_n = @view(tmp_vec1[1:tmp_dist_rows]),
                            tmp_m = @view(tmp_vec2[1:tmp_dist_cols]))
        out.localarray[out.partitioning[r+1][2], :] .= tmp_dist_view
    end
    out
end


text:
using Pkg; Pkg.activate("."); Pkg.instantiate()
using Revise

# include(".\\examples\\nonlin_statics\\twisting_circle_examples.jl");
# pl = twisting_circle_examples.twisting_circle_1()
# savejson(pl, "plot.json")
# open("plot.json", "r") do f; plot(JSON.parse(Plot, String(read(f)))) end

# include(".\\examples\\nonlin_statics\\twisting_circle_examples.jl");
# twisting_circle_examples.twisting_circle()

# snapnum(n) = begin
#     s = ""
#     if n <= 9
#         s = "00$(n)"
#     elseif n <= 99
#         s = "0$(n)"
#     else
#         s = "$(n)"
#     end
#     s
# end
# run(`cp ../twi*.png .`)
# K = 1
# for i in 1:5:501
#     if mod(i, 10) == 1
#         run(`cp "twisting_circle-$(i).png" "snap$(snapnum(K)).png"`)
#         global K += 1
#     end
# end
# run(`magick snap*.png anim.gif`)

# include(".\\examples\\nonlin_transient\\twisting_circle_examples.jl");
# twisting_circle_examples.allrun()

include(".\\examples\\nonlin_statics\\curved_cantilever_examples.jl");
curved_cantilever_examples.curved_cantilever_thin()

include(".\\examples\\nonlin_statics\\curved_cantilever_examples.jl");
curved_cantilever_examples.curved_cantilever()

include(".\\examples\\meshing\\mesh_examples.jl");
mesh_examples.allrun()

include(".\\examples\\modal\\beam_modal_examples.jl");
beam_modal_examples.allrun()

include(".\\examples\\modal\\argyr_l_frame_modal_examples.jl");
# argyr_l_frame_modal_examples.allrun()
argyr_l_frame_modal_examples.argyr_l_frame_modal_anim()

include(".\\examples\\linear_buckling\\tippling_examples.jl");
tippling_examples.allrun()

include(".\\examples\\nonlin_transient\\fast_top_examples.jl");
@time fast_top_examples.fasttop2()

include(".\\examples\\nonlin_transient\\slow_top_examples.jl");
@time slow_top_examples.slowtop3()

include(".\\examples\\nonlin_transient\\argyr_swing_examples.jl");
argyr_swing_examples.argyr_swing_animated()

include(".\\examples\\nonlin_transient\\argyr_swing_examples.jl");
argyr_swing_examples.argyr_swing_compare()


text:
"""
    box dynamics
    - 3D particle with orientation and 8 corners subject to contact forces

    - configuration: q = (px, py, pz, rx, ry, rz) ∈ R⁶
        - orientation : modified Rodrigues angles

    - contacts (8x)
        - impact force (magnitude): n ∈ R₊
        - friction force: b ∈ R²
            - contact force: λ = (b, n) ∈ R² × R₊
            - friction coefficient: μ ∈ R₊

    Discrete Mechanics and Variational Integrators
        pg. 363
"""

struct Box
	n::Int
	m::Int
	d::Int

    mass # mass
    J # inertia
    μ # friction coefficient
    g # gravity

    r             # corner length
	n_corners     # number of corners
    corner_offset # precomputed corner offsets

    nq # configurations dimension
end

# Methods
function mass_matrix(model::Box)
	Diagonal(@SVector [model.mass, model.mass, model.mass,
		model.J, model.J, model.J])
end

function gravity(model::Box)
	@SVector [0., 0., model.mass * model.g, 0., 0., 0.]
end

function kinematics(model::Box, q)
    p = view(q, 1:3)
    r = view(q, 4:6)

    R = MRP(r...)

    SVector{24}([(p + R * model.corner_offset[1])...,
                 (p + R * model.corner_offset[2])...,
                 (p + R * model.corner_offset[3])...,
                 (p + R * model.corner_offset[4])...,
                 (p + R * model.corner_offset[5])...,
                 (p + R * model.corner_offset[6])...,
                 (p + R * model.corner_offset[7])...,
                 (p + R * model.corner_offset[8])...])
end

function jacobian(model::Box, q)
	k(z) = kinematics(model, z)
	ForwardDiff.jacobian(k, q)
end

function signed_distance(model::Box, q)
	idx = collect([3, 6, 9, 12, 15, 18, 21, 24])
	kinematics(model, q)[idx]
end

function P_func(model::Box, q)
	idx = collect([1, 2, 4, 5, 7, 8, 10, 11, 13, 14, 16, 17, 19, 20, 22, 23])
	k(z) = kinematics(model, z)[idx]
	ForwardDiff.jacobian(k, q)
end

function B_func(model::Box, q)
	# p = view(q, 1:3)
    r = view(q, 4:6)

    R = MRP(r...)
    # @SMatrix [0. 0. 0. R[1,1] R[2,1] R[3,1];
    #           0. 0. 0. R[1,2] R[2,2] R[3,2];
    #           0. 0. 0. R[1,3] R[2,3] R[3,3]]
	@SMatrix [1. 0. 0. 0.0 0.0 0.0;
			  0. 1. 0. 0.0 0.0 0.0;
			  0. 0. 1. 0.0 0.0 0.0;
			  0. 0. 0. R[1,1] R[2,1] R[3,1];
		      0. 0. 0. R[1,2] R[2,2] R[3,2];
		      0. 0. 0. R[1,3] R[2,3] R[3,3]]
end

# dynamics
function dynamics(model, q1, q2, q3, u1, λ, h)
      nq = model.nq
      SVector{nq}(mass_matrix(model) * (2.0 * q2 - q1 - q3) / h
            - h * gravity(model)
            + h * jacobian(model, q3)' * λ
			+ h * transpose(B_func(model, q3)) * u1)
end

# Kinematics
num_contacts = 8
d = 0.5
c8 = @SVector [d, d, d]
c2 = @SVector [d, d, -d]
c3 = @SVector [d, -d, d]
c4 = @SVector [d, -d, -d]
c5 = @SVector [-d, d, d]
c6 = @SVector [-d, d, -d]
c7 = @SVector [-d, -d, d]
c1 = @SVector [-d, -d, -d]

corner_offset = @SVector [c1, c2, c3, c4, c5, c6, c7, c8]

# Model
nq = 6
nu = 6
model = Box(2 * nq, nu, 0,
			1.0,
			1.0 / 12.0 * 1.0 * ((2.0 * d)^2 + (2.0 * d)^2),
 			0.5,
			9.81,
            d, num_contacts, corner_offset,
            nq)

# qq = rand(nq)
# kinematics(model, qq)
# ϕ_func(model, qq)
# jacobian(model, qq)
# P_func(model, qq)
num_contacts
# var
num_var = nq + num_contacts + num_contacts + 3 * num_contacts + 3 * num_contacts
function unpack(z)
	q = view(z, 1:nq)
	n = view(z, nq .+ (1:num_contacts))
	sϕ = view(z, nq + num_contacts .+ (1:num_contacts))
	b = view(z, nq + 2 * num_contacts .+ (1:3 * num_contacts))
	sb = view(z, nq + 2 * num_contacts + 3 * num_contacts .+ (1:3 * num_contacts))

	b_traj = [view(b, (i - 1) * 3 .+ (1:3)) for i = 1:num_contacts]
	sb_traj = [view(sb, (i - 1) * 3 .+ (1:3)) for i = 1:num_contacts]

	return q, n, sϕ, b_traj, sb_traj
end

function initialize(q2, num_var; z_init = 1.0)
	z = z_init * ones(num_var)

	z[1:nq] = copy(q2)
	z[nq + num_contacts .+ (1:num_contacts)] = signed_distance(model, q2)
	z[nq + 2 * num_contacts .+ (1:3 * num_contacts)] = vcat([[model.μ * view(z, nq .+ (1:num_contacts))[i]; 0.1 * model.μ * view(z, nq .+ (1:num_contacts))[i]; 0.1 * model.μ * view(z, nq .+ (1:num_contacts))[i]] for i = 1:num_contacts]...)
	z[nq + 2 * num_contacts + 3 * num_contacts .+ (1:3 * num_contacts)] = vcat([[10.0 * z_init; 1.0 * z_init; 1.0 * z_init] for i = 1:num_contacts]...)
	return z
end

function visualize!(vis, model::Box, q;
        Δt = 0.1)

	default_background!(vis)

	setobject!(vis["box"], GeometryBasics.Rect(Vec(-1.0 * model.r,
		-1.0 * model.r,
		-1.0 * model.r),
		Vec(2.0 * model.r, 2.0 * model.r, 2.0 * model.r)),
		MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))

    for i = 1:model.n_corners
        setobject!(vis["corner$i"], GeometryBasics.Sphere(Point3f0(0),
            convert(Float32, 0.05)),
            MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0.0, 1.0)))
    end

    anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

    for t = 1:length(q)
        MeshCat.atframe(anim, t) do

            settransform!(vis["box"],
				compose(Translation(q[t][1:3]...), LinearMap(MRP(q[t][4:6]...))))

            for i = 1:model.n_corners
                settransform!(vis["corner$i"],
                    Translation((q[t][1:3] + MRP(q[t][4:6]...) * (corner_offset[i]))...))
            end
        end
    end
    # settransform!(vis["/Cameras/default"], compose(Translation(-1, -1, 0),LinearMap(RotZ(pi/2))))
    MeshCat.setanimation!(vis, anim)
end


text:
module Builtins

# Raw pointer accesses
typealias DeviceArray Ptr

import Base.getindex
@inline getindex{T}(a::DeviceArray{T}, i0::Real) =
    Base.unsafe_load(a, Base.to_index(i0))::T
import Base.setindex!
@inline setindex!{T}(a::DeviceArray{T}, x::T, i0::Real) =
    Base.unsafe_store!(a, x, Base.to_index(i0))

# intrinsic functions
using Core.Intrinsics.llvmcall

const BUILTINS = Symbol[]

is_builtin(b::Symbol) = in(b, BUILTINS)
is_builtin(b) = false

const int_intrinsics = [
    ("@_Z13get_global_idj", Int64),
    ("@_Z12get_local_idj", Int64),
    ("@_Z17get_global_offsetj", Int64),
    ("@_Z12get_group_idj", Int64),
    ("@_Z15get_global_sizej", Int64),
    ("@_Z14get_local_sizej", Int64),
    ("@_Z14get_num_groupsj", Int64),
    ("@_Z12get_work_dimv", Int64),
    ("@_Z23get_enqueued_local_sizej", Int64),
    ("@_Z20get_global_linear_idv", Int64),
    ("@_Z19get_local_linear_idv", Int64),
    ]

type IntrinsicDef
    name
    ret_type
    arg_types

    function IntrinsicDef()
        new("", Void, Tuple{})
    end
end

const typeMap = Dict(
    "j" => Int32,
    "v" => Void
    )

const llvmTypeMap = Dict(
    Int32 => "i32",
    Int64 => "i64"
    )

function parse_mangled_name(name)
    def = IntrinsicDef()

    len_str = match(r"([0-9]+)(.+)", name)

    assert(len_str != nothing)

    len = parse(Int, len_str.captures[1])

    name_tail = len_str.captures[2];

    def.name = name_tail[1:len]

    arg_typecodes = name_tail[len+1:end]

    arg_types = filter(
    (x) -> x != Void,
        [ typeMap[string(x)] for x in arg_typecodes ]
    )

    def.arg_types = arg_types

    return def
end

function intrinsic_impl()
    for intr in int_intrinsics
        mangled_name = intr[1]

        def = parse_mangled_name(mangled_name)
        def.ret_type = intr[2]

        llvm_argtypes = map(
            (x) -> llvmTypeMap[x],
            def.arg_types
        )

        arg_names = Array(Symbol, length(def.arg_types))
        arg_exprs = Array(Expr, length(def.arg_types))

        arg = 'a'
        for i = 1:length(arg_names)
            arg_name = symbol(arg)
            arg_type = def.arg_types[i]

            arg_exprs[i] = :($arg_name :: $arg_type)

            arg_names[i] = arg_name
            arg += 1
        end

        intr_name = symbol(def.name)

        call_expr = Expr(
            :call,
            intr_name,
            arg_exprs...)

        llvm_argtype_string = join(llvm_argtypes, ",")
        llvm_argcall_string = join([ "$(llvm_argtypes[i]) %$(i-1)" for i = 1:length(llvm_argtypes)], ",")

        tuple_expr = Expr(
            :tuple,
            "declare spir_func i64 $mangled_name($llvm_argtype_string) readnone nounwind",
            """%res = tail call spir_func i64 $mangled_name($llvm_argcall_string)
               ret i64 %res""")
        argtype_expr = Expr(
            :curly,
            :Tuple,
            def.arg_types...
            )

        llvmcall_expr = Expr(
            :call,
            :llvmcall,
            tuple_expr,
            def.ret_type,
            argtype_expr,
            arg_names...
            )

        impl_expr = Expr(
            :(=),
            call_expr,
            llvmcall_expr
        )

        final_expr = quote
            $impl_expr
        end

        # Add final code to the module
        eval(final_expr)

        push!(BUILTINS, intr_name)
    end
end

intrinsic_impl()
end # module Builtins


text:
# Beg et al. (2007): https://doi.org/10.1073/pnas.0609845104.
module BegData

    using DataFrames
    import ..Chemostat_Heerden2013
    const ChH = Chemostat_Heerden2013
    import CSV

    import UtilsJL
    const UJL = UtilsJL
    UJL.gen_sub_proj(@__MODULE__)

    # The average crowding coefficient a was fitted to obtain the minimum square
    # deviation between the measured and model predicted growth rates,
    # resulting in a = 0.0040 ± 0.0005 h•g/mmol, in which g is grams dry weight.
    # However, the maximum growth rates on glucose and glycerol are more consistent
    # with a = 0.0031 ± 0.0001 h•g/mmol and a = 0.0053 ± 0.0001 h•g/mmol, respectively.
    const ave_a  = 0.0031 # I will use this to set the unknown costs

    # enzymatic costs from Beg et al. (2007): https://doi.org/10.1073/pnas.0609845104.
    function load_enz_data()
        datfile = rawdir("beg2007___enzymatic_data.tsv")
        beg_enz_data = CSV.read(datfile, DataFrame)
    end

    function __init__()
        UJL.create_proj_dirs(@__MODULE__)
    end

end


text:
using Gtk.ShortNames

function clickwindow()
    clicks = 0
    win = Window("Click Counter", 300, 100) |> (Frame() |> (vbox = Box(:v)))
    lab = Label("There have been no clicks yet.")
    but = Button("click me")
    push!(vbox, lab)
    push!(vbox, but)
    setproperty!(vbox, :expand, lab, true)
    setproperty!(vbox, :spacing, 20)
    callback(w) = (clicks += 1; setproperty!(lab, :label, "There have been $clicks button clicks."))
    id = signal_connect(callback, but, "clicked")
    c = Condition()
    endit(w) = notify(c)
    signal_connect(endit, win, :destroy)
    showall(win)
    wait(c)
end

clickwindow()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Routes
using Base.Test

r1 = Routes.route("./RESULT/",:day=>r"(\d\d\d\d-\d\d-\d\d)",".gz")
@test r1|>typeof == Routes.Route
@test r1.regex == r"./RESULT/(?P<day>(\d\d\d\d-\d\d-\d\d)).gz"

m1 = match( r1.regex, "./RESULT/2016-06-11.gz" )
@test m1 != nothing
@test m1.match == "./RESULT/2016-06-11.gz"
@test Routes.file( r1, m1.match) == "./RESULT/./RESULT/2016-06-11.gz.gz"

max_stars_repo_name: closescreen/Routes

max_issues_repo_name: closescreen/Routes

max_forks_repo_name: closescreen/Routes


text:
# Arrays

initn(dims...) = randn(dims...)/100

flatten(xs) = reshape(xs, size(xs, 1), :)

unsqueeze(xs, dim) = reshape(xs, (size(xs)[1:dim-1]..., 1, size(xs)[dim:end]...))

stack(xs, dim) = cat(dim, unsqueeze.(xs, dim)...)
unstack(xs, dim) = [slicedim(xs, dim, i) for i = 1:size(xs, dim)]

"""
    chunk(xs, n)

Split `xs` into `n` parts.

```julia
julia> chunk(1:10, 3)
3-element Array{Array{Int64,1},1}:
 [1, 2, 3, 4]
 [5, 6, 7, 8]
 [9, 10]
```
"""
chunk(xs, n) = collect(Iterators.partition(xs, ceil(Int, length(xs)/n)))

batchindex(xs, i) = (reverse(Base.tail(reverse(indices(xs))))..., i)

"""
  batch(xs)

Batch the arrays in `xs` into a single array.

```julia
julia> batch([[1,2,3],[4,5,6]])
3×2 Array{Int64,2}:
 1  4
 2  5
 3  6
```
"""
function batch(xs)
  data = first(xs) isa AbstractArray ?
    similar(first(xs), size(first(xs))..., length(xs)) :
    Vector{eltype(xs)}(length(xs))
  for (i, x) in enumerate(xs)
    data[batchindex(data, i)...] = x
  end
  return data
end

Base.rpad(v::AbstractVector, n::Integer, p) = [v; fill(p, max(n - length(v), 0))]

"""
    batchseq(seqs, pad)

Take a list of `N` sequences, and turn them into a single sequence where each
item is a batch of `N`. Short sequences will be padded by `pad`.

```julia
julia> batchseq([[1, 2, 3], [4, 5]], 0)
3-element Array{Array{Int64,1},1}:
 [1, 4]
 [2, 5]
 [3, 0]
```
"""
function batchseq(xs, pad = nothing, n = maximum(length(x) for x in xs))
  xs_ = [rpad(x, n, pad) for x in xs]
  [batch([xs_[j][i] for j = 1:length(xs_)]) for i = 1:n]
end

# Other

function accuracy(m, data)
  n = 0
  correct = 0
  for (x, y) in data
    x, y = tobatch.((x, y))
    n += size(x, 1)
    correct += sum(argmax(m(x)) .== argmax(y))
  end
  return correct/n
end

"""
Returns a function that when invoked, will only be triggered at most once
during `timeout` seconds. Normally, the throttled function will run
as much as it can, without ever going more than once per `wait` duration;
but if you'd like to disable the execution on the leading edge, pass
`leading=false`. To enable execution on the trailing edge, ditto.
"""
function throttle(f, timeout; leading=true, trailing=false)
  cooldown = true
  later = nothing

  function throttled(args...; kwargs...)
    yield()

    if cooldown
      if leading
        f(args...; kwargs...)
      else
        later = () -> f(args...; kwargs...)
      end

      cooldown = false
      @schedule try
        while (sleep(timeout); later != nothing)
          later()
          later = nothing
        end
      finally
        cooldown = true
      end
    elseif trailing
      later = () -> f(args...; kwargs...)
    end

    nothing
  end
end


text:
using Test
using TravisTest

println("Starting Test 1")
##Test 1, check the local func works (defined in the module script)
#Function doubles input arg
Y=MrFunc(2);
@test Y==4
println("Test 1 passed")

println("Starting Test 2")
##Test 2, check the external one works (imported into the module script)
#Function halves input arg
Y=TravisTest.MrsFunc(8);
@test Y==4
println("Test 2 passed")

max_stars_repo_name: Timmmdavis/JuliaTravisTest

max_issues_repo_name: Timmmdavis/JuliaTravisTest

max_forks_repo_name: Timmmdavis/JuliaTravisTest


text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mb03bd(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    KMAX = 6
    NMAX = 50
    LDA1 = NMAX
    LDA2 = NMAX
    LDQ1 = NMAX
    LDQ2 = NMAX
    LDWORK = KMAX + max( 2*NMAX, 8*KMAX )
    LIWORK = 2*KMAX + NMAX
    QIND = Array{BlasInt,1}(undef, KMAX)
    S = Array{BlasInt,1}(undef, KMAX)
    A = Array{Float64,3}(undef, LDA1,LDA2,KMAX)
    Q = Array{Float64,3}(undef, LDQ1,LDQ2,KMAX)
    ALPHAR = Array{Float64,1}(undef, NMAX)
    ALPHAI = Array{Float64,1}(undef, NMAX)
    BETA = Array{Float64,1}(undef, NMAX)
    SCAL = Array{BlasInt,1}(undef, NMAX)
    IWORK = Array{BlasInt,1}(undef, LIWORK)
    DWORK = Array{Float64,1}(undef, LDWORK)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    JOB =  vs[1][1]
    DEFL =  vs[2][1]
    COMPQ =  vs[3][1]
    K = parse(BlasInt, vs[4])
    N = parse(BlasInt, vs[5])
    H = parse(BlasInt, vs[6])
    ILO = parse(BlasInt, vs[7])
    IHI = parse(BlasInt, vs[8])
    if ( N<0 || N>NMAX )
        @error "Illegal N=$N"
    end

    vs = String[]
    _isz = K
    while length(vs) < _isz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    S[1:_isz] .= parsex.(BlasInt, vs)

    vs = String[]
    _isz,_jsz,_ksz = (N,N,K)
    while length(vs) < _isz*_jsz*_ksz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for k in 1:_ksz
      for i in 1:_isz
        _i0 = (i-1)*_jsz + (k-1)*_jsz*_isz
        A[i,1:_jsz,k] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
      end
    end
    if ( LSAME( COMPQ, 'P' ) )

    vs = String[]
    _isz = K
    while length(vs) < _isz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    QIND[1:_isz] .= parsex.(BlasInt, vs)

    end # if
    close(f)
# interp call 1

    INFO, IWARN = SLICOT.mb03bd!(JOB, DEFL, COMPQ, QIND, K, N, H, ILO, IHI, S, A, Q, ALPHAR, ALPHAI, BETA, SCAL, LIWORK, LDWORK)
    @test INFO == 0
    println(io, "IWARN = $IWARN")

    if ( LSAME( JOB, 'S' ) || LSAME( JOB, 'T' ) )
# interp output 1
    println(io, "A:")
    _nc = N
    _nr = N
    _nk = K
    show(io, "text/plain", A[1:_nr,1:_nc,1:_nk])
    println(io)

    end # if
    if ( LSAME( COMPQ, 'U' ) || LSAME( COMPQ, 'I' ) )
# interp output 2
    println(io, "Q:")
    _nc = N
    _nr = N
    _nk = K
    show(io, "text/plain", Q[1:_nr,1:_nc,1:_nk])
    println(io)

    elseif ( LSAME( COMPQ, 'P' ) )
        for L in 1:K
            if ( QIND[L]>0 )
                println(io, "factor ",QIND[L])
#   write QIND( L )
# unable to translate write loop:
#  write ( I, J, QIND( L ) ), J = 1, N
# interp output 3
                show(io, "text/plain", Q[1:N,1:N,QIND[L]])
                println(io)

            end # if
        end # for
    end # if
# interp output 4
    println(io, "ALPHAR:")
    _nr = N
    show(io, "text/plain", ALPHAR[1:_nr])
    println(io)

# interp output 5
    println(io, "ALPHAI:")
    _nr = N
    show(io, "text/plain", ALPHAI[1:_nr])
    println(io)

# interp output 6
    println(io, "BETA:")
    _nr = N
    show(io, "text/plain", BETA[1:_nr])
    println(io)

# interp output 7
    println(io, "SCAL:")
    _nr = N
    show(io, "text/plain", SCAL[1:_nr])
    println(io)
end # run_mb03bd()


text:
using Revise
using JLD2
using ViewerGL
using LarSurf

# created by surface_extraction_parallel_ircad01.jl
@load "liver01.jld2" V FV

FVtri = LarSurf.triangulate_quads(FV)
EV = LarSurf.Smoothing.get_EV_quads(FV)
# Plasm.view(val)
ViewerGL.VIEW([
    # ViewerGL.GLGrid(V,FVtri,ViewerGL.Point4d(1,0,1,0.1))
    ViewerGL.GLGrid(V,EV,ViewerGL.Point4d(0.5,0,1,0.5))
	ViewerGL.GLAxis(ViewerGL.Point3d(-1,-1,-1),ViewerGL.Point3d(1,1,1))
])