text:
export SA_AMGsetup;

"""
Based on the paper
Eran Treister and Irad Yavneh, Non-Galerkin Multigrid based on Sparsified Smoothed Aggregation.
SIAM Journal on Scientific Computing, 37 (1), A30-A54, 2015.
"""
function SA_AMGsetup(AT::SparseMatrixCSC{VAL,IND},param::MGparam{VAL,IND},symm::Bool = true,nrhs::Int64 = 1,verbose::Bool=false) where {VAL,IND}
Ps = Array{SparseMatrixCSC{real(VAL),IND}}(undef,param.levels-1);
Rs = Array{SparseMatrixCSC{real(VAL),IND}}(undef,param.levels-1);
As = Array{SparseMatrixCSC{VAL,IND}}(undef,param.levels);

if symm == false
	error("not supported yet...");
end

T_time = 0;
N = size(AT,2);
As[1] = AT;
relaxPrecs = Array{Array{VAL}}(undef,param.levels-1);
Cop = nnz(AT);
for l = 1:(param.levels-1)
	if verbose
		T_time = time_ns();
	end
    AT = As[l];
    if param.relaxType=="Jac" || param.relaxType=="Jac-GMRES" || param.relaxType=="SPAI"
		relaxPrecs[l] = getRelaxPrec(AT,param.relaxType,param.relaxParam,[],false);
	else
		error("Unknown relaxation type !!!!");
	end
	P0 = getAggregation(AT,param.strongConnParam);
	Nc = size(P0,2);
	P0 = sparse(P0');
	if (size(P0,1)==size(P0,2))
		if verbose; println(string("Stopped Coarsening at level ",l)); end
		param.levels = l;
		As = As[1:l];
		Ps = Ps[1:l-1];
		Rs = Rs[1:l-1];
		relaxPrecs = relaxPrecs[1:l-1];
		break;
	else
		DAT = AT*spdiagm(0=>relaxPrecs[l]); # Julia's .* doesn't work for this multiplication. Only from the right.
		rhoDAT = min(norm(DAT,1),norm(DAT,Inf));
		PT = P0 - (1.33/rhoDAT)*P0*DAT;
		Rs[l] = sparse(PT'); # this is becasue we hold the transpose of the matrices and P = R' anyway here....
		Ps[l] = PT;

		Act = Ps[l]*AT*Rs[l];
		As[l+1] = Act;
		Cop = Cop + nnz(Act);
		#fprintf('setup level %3d, dim:%3dx%3d, nnz:%3d\n',nlevels+1,nc(1),nc(2),nnz(Ac))

		if verbose; println("MG setup: ",N," took:",(time_ns() - T_time)/1e+9); end;
		N = Nc;
	end
end
if verbose
	T_time = time_ns();
	println("MG Setup: Operator complexity = ",Cop/nnz(As[1]));
end
As[end] = As[end] + 1e-8*norm(As[end],1)*sparse(1.0I,size(As[end],2),size(As[end],2));

defineCoarsestAinv(param,As[end]);

if verbose
	println("MG setup coarsest: ",N,", done LU in ",(time_ns() - T_time)/1e+9);
end
param.As = As;
param.Ps = Ps;
param.Rs = Rs;
param.relaxPrecs = relaxPrecs;
param = adjustMemoryForNumRHS(param,nrhs,verbose);
return;
end

function getAggregation(AT::SparseMatrixCSC,strengthConnParam::Float64)
	if size(AT,2) <= 100
		return sparse(1.0I,size(AT,2),size(AT,2));
	end
	S = getStrengthMatrix(AT,strengthConnParam);
	aggr = neighborhoodAggregationNew(S);
	# println(aggr)
	return aggrArray2P(aggr);
end

function getStrengthMatrix(AT::SparseMatrixCSC,strengthConnParam::Float64)
S = -AT;
mm = 1e-16*maximum(S.nzval);
n = size(S,2);
for j = 1:n
	maxVal_j = mm;
	for gidx = S.colptr[j] : S.colptr[j+1]-1
		if S.nzval[gidx] > maxVal_j
			maxVal_j = S.nzval[gidx];
		end
	end
	scal_k = 1.0 ./maxVal_j;
	for gidx = S.colptr[j] : S.colptr[j+1]-1
		S.nzval[gidx]*=scal_k;
	end
	for gidx = S.colptr[j] : S.colptr[j+1]-1
		if S.rowval[gidx] == j
			# We have a diagonal element.
			S.nzval[gidx] = 1.0;
		end
	end
	for gidx = S.colptr[j] : S.colptr[j+1]-1
		if S.nzval[gidx] <  strengthConnParam
			S.nzval[gidx] = 0.0;
		end
	end
end
S = S + S';
end

function neighborhoodAggregationNew(S::SparseMatrixCSC)
## S is symmetric. We look only at non-zero values here...
tau = 3.0;
n = size(S,2);
aggr = zeros(Int64,n);

avg_sparsity = 0.0;
aux = zeros(n);
aux_count = zeros(Int64,n);
for k=1:n
	avg_sparsity += S.colptr[k+1] - S.colptr[k];
end
avg_sparsity /= n;

for k = 1:n
    if S.colptr[k+1] - S.colptr[k] > tau*avg_sparsity
		aux_count[k] = -1;
	end
end
#########################################################################
for k = 1:n
    Neighbors_Aggregated_flag = false;
	if aux_count[k] == -1
		continue; # ignore...
	end
    for global_idx = S.colptr[k] : S.colptr[k+1]-1
        if aggr[S.rowval[global_idx]] != 0
            Neighbors_Aggregated_flag = true;
            break;
		end
    end
	if !Neighbors_Aggregated_flag
        for global_idx = S.colptr[k] : S.colptr[k+1]-1
			if aux_count[S.rowval[global_idx]] != -1
				aggr[S.rowval[global_idx]] = k; #conversion to Matlab's indices;
				aux_count[k]+=1;
			end
        end
    end
end
#########################################################################
for k = 1:n
    Neighbors_Aggregated_flag = false;
	if aux_count[k]!=-1
		continue;
	end
	aux_count[k] = 0;
    for global_idx = S.colptr[k] : S.colptr[k+1]-1
        if aggr[S.rowval[global_idx]] != 0
            Neighbors_Aggregated_flag = true;
            break;
		end
    end
	if !Neighbors_Aggregated_flag
        for global_idx = S.colptr[k] : S.colptr[k+1]-1
			aggr[S.rowval[global_idx]] = k; #conversion to Matlab's indices;
			aux_count[k]+=1;
        end
    end
end
#########################################################################
for k = 1:n
	chosen_score = 0.0;
	chosen = 0;
	if aggr[k] == 0
		for global_idx = S.colptr[k] : S.colptr[k+1]-1
			if aggr[S.rowval[global_idx]] > 0 # aggr can be negative here....
                agg_of_neighbor = aggr[S.rowval[global_idx]];
				aux[agg_of_neighbor] += S.nzval[global_idx];
			end
			for global_idx = S.colptr[k] : S.colptr[k+1]-1
				if aggr[S.rowval[global_idx]]>0 # aggr can be negative here....
                    agg_of_neighbor = aggr[S.rowval[global_idx]];
                    if chosen_score < aux[agg_of_neighbor]/aux_count[agg_of_neighbor]
						chosen_score = aux[agg_of_neighbor]/aux_count[agg_of_neighbor];
						chosen = agg_of_neighbor;
                        aux[agg_of_neighbor] = 0;
					end
                end
			end
            aggr[k] = -chosen; # we don't want the new aggregate to spread - use only original aggregate
		end
	end
end
#########################################################################
for k = 1:n
	if aggr[k] < 0.0
		aggr[k] = -aggr[k];
	end
end

return aggr;
end

function aggrArray2P(aggr)
	fine2coarse = zeros(Int64,size(aggr));
	n = length(aggr);
	I = findall(aggr .== 1:n);
	fine2coarse[I] = 1:length(I);
	aggr = fine2coarse[aggr];
	if sum(aggr.==0)>0
		error("nodes without aggregates");
	end
	P = sparse(collect(1:n),aggr,ones(n));
	return P;
end


text:
VcovFormula(::Type{Val{:robust}}) = VcovRobustFormula()

struct VcovRobustFormula <: AbstractVcovFormula  end

struct VcovRobustMethod <: AbstractVcovMethod end
VcovMethod(::AbstractDataFrame, ::VcovRobustFormula) = VcovRobustMethod()

function vcov!(v::VcovRobustMethod, x::VcovData)
    S = shat!(v, x)
    return sandwich(x.crossmatrix, S)
end

# S_{(l-1) * K + k, (l'-1)*K + k'} = \sum_i X[i, k] res[i, l] X[i, k'] res[i, l']
function shat!(::VcovRobustMethod, x::VcovData{T, N}) where {T, N}
    dim = size(x.regressors, 2) * size(x.residuals, 2)
    X2 = fill(zero(Float64), size(x.regressors, 1), dim)
    index = 0
    for k in 1:size(x.residuals, 2)
        for j in 1:size(x.regressors, 2)
            index += 1
            @inbounds @simd for i in 1:size(x.regressors, 1)
                X2[i, index] = x.regressors[i, j] * x.residuals[i, k]
            end
        end
    end
    S2 = X2' * X2
    rmul!(S2, size(x.regressors, 1) / x.dof_residual)
    return S2
end

ext: jl


text:
linear_index(combinations, grade::Integer) = findfirst(==(grade), length.(combinations))
linear_index(combinations, indices) = findfirst(==(indices), combinations)

indices_from_linear_index(combinations, index::Integer) = combinations[index]


text:
export TabularApproximator, TabularVApproximator, TabularQApproximator

"""
    TabularApproximator(table<:AbstractArray, opt)

For `table` of 1-d, it will serve as a state value approximator. See [`TabularVApproximator`](@ref).
For `table` of 2-d, it will serve as a state-action value approximator. See [`TabularQApproximator`](@ref).

Note that actions and states should be presented to `TabularApproximator` as integers starting from
1 to be used as the index of the table. That is, e.g., [`RLBase.state_space`](@ref) is expected to
return `Base.OneTo(n_state)`, where `n_state` is the number of states.

!!! warning
    For `table` of 2-d, the first dimension is action and the second dimension is state.
"""
struct TabularApproximator{N,T<:AbstractArray{<:AbstractFloat, N},O} <: AbstractApproximator
    table::T
    optimizer::O
    function TabularApproximator(table::T, opt::O) where {T<:AbstractArray,O}
        n = ndims(table)
        n <= 2 || throw(ArgumentError("the dimension of table must be <= 2"))
        new{n,T,O}(table, opt)
    end
end

const TabularVApproximator = TabularApproximator{1}
const TabularQApproximator = TabularApproximator{2}

"""
    TabularVApproximator(; n_state, init = 0.0, opt = InvDecay(1.0))

A state value approximator represented by a 1-d table. `init` is the initial value of each state.
"""
TabularVApproximator(; n_state, init = 0.0, opt = InvDecay(1.0)) =
    TabularApproximator(fill(init, n_state), opt)
"""
    TabularQApproximator(; n_state, n_action, init = 0.0, opt = InvDecay(1.0))

An action-state value approximator represented by a 2-d table. `init` is the initial value of each
pair of action-state.
"""
TabularQApproximator(; n_state, n_action, init = 0.0, opt = InvDecay(1.0)) =
    TabularApproximator(fill(init, n_action, n_state), opt)

(app::TabularVApproximator)(s::Int) = @views app.table[s]

(app::TabularQApproximator)(s::Int) = @views app.table[:, s]
(app::TabularQApproximator)(s::Int, a::Int) = app.table[a, s]

function RLBase.update!(app::TabularVApproximator, correction::Pair{Int,Float64})
    s, e = correction
    x = @view app.table[s]
    x̄ = @view [e][1]
    Flux.Optimise.update!(app.optimizer, x, x̄)
end

function RLBase.update!(app::TabularQApproximator, correction::Pair{Tuple{Int,Int},Float64})
    (s, a), e = correction
    x = @view app.table[a, s]
    x̄ = @view [e][1]
    Flux.Optimise.update!(app.optimizer, x, x̄)
end

function RLBase.update!(app::TabularQApproximator, correction::Pair{Int,Vector{Float64}})
    s, errors = correction
    x = @view app.table[:, s]
    Flux.Optimise.update!(app.optimizer, x, errors)
end


text:
module NReco
include("types.jl")
include("recof.jl")
include("recohits.jl")
include("nemareco.jl")
include("classify_events.jl")
end

max_stars_repo_name: andLaing/NReco

max_issues_repo_name: andLaing/NReco

max_forks_repo_name: andLaing/NReco


text:
# SPDX-License-Identifier: X11
# 2020-11-20
# Popular Vote (200pt)

N,M = parse.(Int, (split ∘ readline)())
A = parse.(Int, (split ∘ readline)())
println((length ∘ filter)(x -> x >= ceil(sum(A)/(4M)), A) >= M ?
          "Yes" : "No")

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg
using BinaryBuilderBase: get_addable_spec

name = "Pfapack"
version = v"0.1.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/xrq-phys/Pfapack.git", "ed885d81e31d88c86016f6664f7df08e50c32d6a")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd Pfapack

cmake fortran \
    -DCMAKE_POSITION_INDEPENDENT_CODE=ON \
    -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release
make -j${nproc} VERBOSE=1
make install

# Manual conversion from static to dynamic lib.
# Avoid linking libgfortran (these .f files make no reference to Fortran libs.)
cc -shared $(flagon -Wl,--whole-archive) ${prefix}/lib/libpfapack.a $(flagon -Wl,--no-whole-archive) \
    -o ${libdir}/libpfapack.${dlext} \
    -L${libdir} -lblastrampoline -lm

# Copy license file
install_license LapackLicence
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
# No need to expand Fortran version since there code make no call to
# standard Fortran libraries
platforms = supported_platforms()

# Since we need to link to libblastrampoline which has seen multiple
# ABI-incompatible versions, we need to expand the julia versions we target
julia_versions = [v"1.7.0", v"1.8.0", v"1.9.0"]
function set_julia_version(platforms::Vector{Platform}, julia_version::VersionNumber)
    _platforms = deepcopy(platforms)
    for p in _platforms
        p["julia_version"] = string(julia_version)
    end
    return _platforms
end
expand_julia_versions(platforms::Vector{Platform}, julia_versions::Vector{VersionNumber}) =
    vcat(set_julia_version.(Ref(platforms), julia_versions)...)
platforms = expand_julia_versions(platforms, julia_versions)

# The products that we will ensure are always built
products = [
    LibraryProduct("libpfapack", :libpfapack)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(get_addable_spec("libblastrampoline_jll", v"3.0.4+0"); platforms=filter(p -> VersionNumber(p["julia_version"]) == v"1.7.0", platforms)),
    Dependency(get_addable_spec("libblastrampoline_jll", v"5.1.0+0"); platforms=filter(p -> VersionNumber(p["julia_version"]) >= v"1.8.0", platforms)),
]

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

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
__precompile__()
module KernelRidgeRegression

export KRR,
    FastKRR,
    RandomFourierFeatures,
    TruncatedNewtonKRR,
    NystromKRR,
    SubsetRegressorsKRR,
    fitPar

import Base: show, showcompact, display
import MLKernels
import StatsBase
import StatsBase: fit, fitted, predict, nobs, predict!, RegressionModel, sample, weights

abstract AbstractKRR{T} <: RegressionModel

function fit(::Type{AbstractKRR}) error("not implemented") end

StatsBase.fitted(KRR::AbstractKRR) = predict(KRR, KRR.X)

include("krr.jl")
include("util.jl")

end # module KernelRidgeRegression


text:
using CairoMakie
using Colors
using Makie.GeometryBasics
using JLD2

training_data = [
    (-5e-4, -1e-8),
    (-5e-4, -2e-8),
    (-5e-4, -3e-8),

    (-3.5e-4, -1e-8),
    (-3.5e-4, -2e-8),
    (-3.5e-4, -3e-8),

    (-2e-4, -1e-8),
    (-2e-4, -2e-8),
    (-2e-4, -3e-8),

    (-5e-4, 1e-8),
    (-5e-4, 2e-8),
    (-5e-4, 3e-8),

    (-3.5e-4, 1e-8),
    (-3.5e-4, 2e-8),
    (-3.5e-4, 3e-8),

    (-2e-4, 1e-8),
    (-2e-4, 2e-8),
    (-2e-4, 3e-8),
]

train_files = [
    "wind_-5e-4_heating_-1e-8_new"
    "wind_-5e-4_heating_-2e-8_new"
    "wind_-5e-4_heating_-3e-8_new"

    "wind_-3.5e-4_heating_-1e-8_new"
    "wind_-3.5e-4_heating_-2e-8_new"
    "wind_-3.5e-4_heating_-3e-8_new"

    "wind_-2e-4_heating_-1e-8_new"
    "wind_-2e-4_heating_-2e-8_new"
    "wind_-2e-4_heating_-3e-8_new"

    "wind_-5e-4_cooling_1e-8_new"
    "wind_-5e-4_cooling_2e-8_new"
    "wind_-5e-4_cooling_3e-8_new"

    "wind_-3.5e-4_cooling_1e-8_new"
    "wind_-3.5e-4_cooling_2e-8_new"
    "wind_-3.5e-4_cooling_3e-8_new"

    "wind_-2e-4_cooling_1e-8_new"
    "wind_-2e-4_cooling_2e-8_new"
    "wind_-2e-4_cooling_3e-8_new"
]

interpolating_data = [
    (-4.5e-4, 2.5e-8),
    (-4.5e-4, 1.5e-8),
    (-2.5e-4, 2.5e-8),
    (-2.5e-4, 1.5e-8),

    (-4.5e-4, -2.5e-8),
    (-4.5e-4, -1.5e-8),
    (-2.5e-4, -2.5e-8),
    (-2.5e-4, -1.5e-8),
]

interpolating_files = [
    "wind_-4.5e-4_cooling_2.5e-8"
    "wind_-4.5e-4_cooling_1.5e-8"
    "wind_-2.5e-4_cooling_2.5e-8"
    "wind_-2.5e-4_cooling_1.5e-8"

    "wind_-4.5e-4_heating_-2.5e-8"
    "wind_-4.5e-4_heating_-1.5e-8"
    "wind_-2.5e-4_heating_-2.5e-8"
    "wind_-2.5e-4_heating_-1.5e-8"
]

extrapolating_data = [
    (-1.5e-4, 3.5e-8),
    (-5.5e-4, 3.5e-8),
    (-5.5e-4, 0),
    (-5.5e-4, -3.5e-8),
    (-1.5e-4, -3.5e-8),
]

extrapolating_files = [
    "wind_-1.5e-4_cooling_3.5e-8"
    "wind_-5.5e-4_cooling_3.5e-8"

    "wind_-5.5e-4_new"

    "wind_-5.5e-4_heating_-3.5e-8"
    "wind_-1.5e-4_heating_-3.5e-8"
]

diurnal_data = [
    (-5.5e-4, 5.5e-8),
    (-1.5e-4, 5.5e-8),
]

diurnal_files = [
    "wind_-5.5e-4_diurnal_5.5e-8"
    "wind_-1.5e-4_diurnal_5.5e-8"
]

momentum_fluxes_training = [data[1] for data in training_data]
buoyancy_fluxes_training = [data[2] for data in training_data]

momentum_fluxes_interpolating = [data[1] for data in interpolating_data]
buoyancy_fluxes_interpolating = [data[2] for data in interpolating_data]

momentum_fluxes_extrapolating = [data[1] for data in extrapolating_data]
buoyancy_fluxes_extrapolating = [data[2] for data in extrapolating_data]

momentum_fluxes_diurnal = [data[1] for data in diurnal_data]
buoyancy_fluxes_diurnal = [data[2] for data in diurnal_data]

fig = CairoMakie.Figure(resolution=(1000, 650))

ax =fig[1,1] = CairoMakie.Axis(fig, xlabel="Buoyancy Flux / m² s⁻³", ylabel="Momentum Flux / m² s⁻²")
color_palette = distinguishable_colors(4, [RGB(1,1,1), RGB(0,0,0)], dropseed=true)

rectangle = CairoMakie.poly!(ax, Point2f0[(-3e-8, -2e-4), (-3e-8, -5e-4), (3e-8, -5e-4), (3e-8, -2e-4)], color=("paleturquoise3", 0.5))

training_points = CairoMakie.scatter!(ax, buoyancy_fluxes_training, momentum_fluxes_training, color=color_palette[1])
interpolating_points = CairoMakie.scatter!(ax, buoyancy_fluxes_interpolating, momentum_fluxes_interpolating, color=color_palette[2])
extrapolating_points = CairoMakie.scatter!(ax, buoyancy_fluxes_extrapolating, momentum_fluxes_extrapolating, color=color_palette[3])

diurnal_line_1 = CairoMakie.lines!(ax, [-5.5e-8, 5.5e-8], [-1.5e-4, -1.5e-4], color=color_palette[4])
diurnal_line_2 = CairoMakie.lines!(ax, [-5.5e-8, 5.5e-8], [-5.5e-4, -5.5e-4], color=color_palette[4])

# diurnal_points = CairoMakie.scatter!(ax, buoyancy_fluxes_diurnal, momentum_fluxes_diurnal, color=color_palette[4])

legend = fig[2,1] = CairoMakie.Legend(fig, [training_points, interpolating_points, extrapolating_points, diurnal_line_1, rectangle],
                                            ["Training", "Interpolating", "Extrapolating", "Diurnal Fluxes", "Interpolation Region"], orientation=:horizontal)

rowsize!(fig.layout, 1, CairoMakie.Relative(0.95))
trim!(fig.layout)
fig
save("final_results/data_WCWH.png", fig, px_per_unit = 4)

T_loss(data) = sum(data) / 1153

losses_training = []

for train_file in train_files
    file = jldopen("final_results/18sim_old/train_$(train_file)/profiles_fluxes_oceananigans.jld2")
    data = file["NDE_profile"]
    close(file)

    loss = T_loss(data["T_losses"])
    loss_mpp = T_loss(data["T_losses_modified_pacanowski_philander"])
    loss_kpp = T_loss(data["T_losses_kpp"])

    loss_min = argmin([loss, loss_mpp, loss_kpp])

    loss_str = ["NDE", "mpp", "kpp"]
    push!(losses_training, loss_str[loss_min])
end

losses_training

losses_interpolating = []

for interpolating_file in interpolating_files
    file = jldopen("final_results/18sim_old/test_$(interpolating_file)/profiles_fluxes_oceananigans.jld2")
    data = file["NDE_profile"]
    close(file)

    loss = T_loss(data["T_losses"])
    loss_mpp = T_loss(data["T_losses_modified_pacanowski_philander"])
    loss_kpp = T_loss(data["T_losses_kpp"])

    @show loss, loss_mpp, loss_kpp

    loss_min = argmin([loss, loss_mpp, loss_kpp])

    loss_str = ["NDE", "mpp", "kpp"]
    push!(losses_interpolating, loss_str[loss_min])
end

losses_interpolating

losses_extrapolating = []

for extrapolating_file in extrapolating_files
    file = jldopen("final_results/18sim_old/test_$(extrapolating_file)/profiles_fluxes_oceananigans.jld2")
    data = file["NDE_profile"]
    close(file)

    loss = T_loss(data["T_losses"])
    loss_mpp = T_loss(data["T_losses_modified_pacanowski_philander"])
    loss_kpp = T_loss(data["T_losses_kpp"])

    @show loss, loss_mpp, loss_kpp

    loss_min = argmin([loss, loss_mpp, loss_kpp])

    loss_str = ["NDE", "mpp", "kpp"]
    push!(losses_extrapolating, loss_str[loss_min])
end

losses_extrapolating

losses_diurnal = []

for diurnal_file in diurnal_files
    file = jldopen("final_results/18sim_old/test_$(diurnal_file)/profiles_fluxes_oceananigans.jld2")
    data = file["NDE_profile"]
    close(file)

    loss = T_loss(data["T_losses"])
    loss_mpp = T_loss(data["T_losses_modified_pacanowski_philander"])
    loss_kpp = T_loss(data["T_losses_kpp"])

    @show loss, loss_mpp, loss_kpp

    loss_min = argmin([loss, loss_mpp, loss_kpp])

    loss_str = ["NDE", "mpp", "kpp"]
    push!(losses_diurnal, loss_str[loss_min])
end

losses_diurnal

fig = CairoMakie.Figure(resolution=(1000, 650))

ax =fig[1,1] = CairoMakie.Axis(fig, xlabel="Buoyancy Flux / m² s⁻³", ylabel="Momentum Flux / m² s⁻²")
color_palette = distinguishable_colors(4, [RGB(1,1,1), RGB(0,0,0)], dropseed=true)

loss_colors = Dict(
    "NDE" => color_palette[1],
    "mpp" => color_palette[2],
    "kpp" => color_palette[3],
)

rectangle = CairoMakie.poly!(ax, Point2f0[(-3e-8, -2e-4), (-3e-8, -5e-4), (3e-8, -5e-4), (3e-8, -2e-4)], color=("paleturquoise3", 0.5))

NDE_point = CairoMakie.scatter!(ax, [buoyancy_fluxes_training[1]], [momentum_fluxes_training[1]], color=loss_colors["NDE"])
mpp_point = CairoMakie.scatter!(ax, [buoyancy_fluxes_training[1]], [momentum_fluxes_training[1]], color=loss_colors["mpp"])
kpp_point = CairoMakie.scatter!(ax, [buoyancy_fluxes_training[1]], [momentum_fluxes_training[1]], color=loss_colors["kpp"])

training_points = [CairoMakie.scatter!(ax, [buoyancy_fluxes_training[i]], [momentum_fluxes_training[i]], color=loss_colors[losses_training[i]]) for i in 1:length(losses_training)]
interpolating_points = [CairoMakie.scatter!(ax, [buoyancy_fluxes_interpolating[i]], [momentum_fluxes_interpolating[i]], color=loss_colors[losses_interpolating[i]]) for i in 1:length(losses_interpolating)]
extrapolating_points = [CairoMakie.scatter!(ax, [buoyancy_fluxes_extrapolating[i]], [momentum_fluxes_extrapolating[i]], color=loss_colors[losses_extrapolating[i]]) for i in 1:length(losses_extrapolating)]

# extrapolating_points = CairoMakie.scatter!(ax, buoyancy_fluxes_extrapolating, momentum_fluxes_extrapolating, color=color_palette[3])

diurnal_line_1 = CairoMakie.lines!(ax, [-5.5e-8, 5.5e-8], [-5.5e-4, -5.5e-4], color=loss_colors[losses_diurnal[1]])
diurnal_line_2 = CairoMakie.lines!(ax, [-5.5e-8, 5.5e-8], [-1.5e-4, -1.5e-4], color=loss_colors[losses_diurnal[2]])

["Ri-based Diffusivity Only", "K-Profile Parameterisation", "NN Embedded in Oceananigans.jl"]
legend = fig[2,1] = CairoMakie.Legend(fig, [mpp_point, kpp_point, NDE_point],
                                ["Ri-based Diffusivity Only", "K-Profile Parameterisation", "NN Embedded in Oceananigans.jl"], orientation=:horizontal)

rowsize!(fig.layout, 1, CairoMakie.Relative(0.95))
trim!(fig.layout)
fig
save("final_results/data_WCWH_loss_results.png", fig, px_per_unit = 4)


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

"""
Unit tests runner for TestPackage.jl package.
"""

# --- Set up Julia environment

test_tools_dir = joinpath(@__DIR__, "..", "..", "..", "..", "..")
push!(LOAD_PATH, test_tools_dir)

# --- Imports

# Standard library
using Test

# External packages
using TestTools: jltest

# Local package
include(joinpath(@__DIR__, "..", "src", "TestPackage.jl"))
using .TestPackage

# --- Run tests

jltest.run_tests(@__DIR__)

# --- Clean up

# Restore LOAD_PATH
filter!(x -> x != test_tools_dir, LOAD_PATH)


text:
using RecyclePermanent
using RNGPool
using JLD2
using DataFrames
using BenchmarkTools
include("../models/gkModel.jl")
include("compare.jl")

@load "data/gkData.jld" ys θ
model = makeGKModel(θ, ys)

@load "compare/gktau.jld" ns trials outR outS outKM1 outKM2 outKM3

PR = Int.(round.((x->x[1]).(outR), sigdigits=2))
NS = Int.(round.((x->x[1]).(outS), sigdigits=2))
PS = Int.(round.(ns .*(x->x[1]).(outS), sigdigits=2))
PKM1 = Int.(round.((x->x[1]).(outKM1), sigdigits=2))
PKM2 = Int.(round.((x->x[1]).(outKM2), sigdigits=2))
PKM3 = Int.(round.((x->x[1]).(outKM3), sigdigits=2))

df = DataFrame(n = ns, PR = PR, PS = PS, PKM1 = PKM1, PKM2 = PKM2, PKM3 = PKM3)
show(df)
println()

m = length(ns)

tR = (i -> @belapsed rpEstimate(model, $(RPIO{model.Tx}(PR[i],ns[i],1)))).(1:m)
tS = (i -> @belapsed simpleEstimate(model, $(RPIO{model.Tx}(NS[i],ns[i],1)))).(1:m)
tKM1 = (i -> @belapsed kmEstimate(model, $(KMIO{model.Tx}(PKM1[i],ns[i])),1.0)).(1:m)
tKM2 = (i -> @belapsed kmEstimate(model, $(KMIO{model.Tx}(PKM2[i],ns[i])),0.001)).(1:m)
tKM3 = (i -> @belapsed kmEstimate(model, $(KMIO{model.Tx}(PKM3[i],ns[i])),0.0001)).(1:m)

TR = round.(tR*1000, digits=1)
TS = round.(tS*1000, digits=1)
TKM1 = round.(tKM1*1000, digits=1)
TKM2 = round.(tKM2*1000, digits=1)
TKM3 = round.(tKM3*1000, digits=1)

df2 = DataFrame(n = ns, TR = TR, TS = TS, TKM1 = TKM1, TKM2 = TKM2, TKM3 = TKM3)
show(df2)
println()


text:
#' @exported
#' @description
#'
#' Render a DataTableRow to an IO system. Each column of the DataTableRow
#' is printed on a separate line.
#'
#' @param io::IO The IO system where rendering will take place.
#' @param r::DataTableRow The DataTableRow to be rendered to `io`.
#'
#' @returns o::Void A `nothing` value.
#'
#' @examples
#'
#' dt = DataTable(A = 1:3, B = ["x", "y", "z"])
#' for r in eachrow(dt)
#'     show(STDOUT, r)
#' end
function Base.show(io::IO, r::DataTableRow)
    labelwidth = mapreduce(n -> length(string(n)), max, _names(r)) + 2
    @printf(io, "DataTableRow (row %d)\n", r.row)
    for (label, value) in r
        println(io, rpad(label, labelwidth, ' '), value)
    end
end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct CustomModel_CustomModelParamValue <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function CustomModel_CustomModelParamValue(; kwargs...)
        obj = new(meta(CustomModel_CustomModelParamValue), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct CustomModel_CustomModelParamValue
const __meta_CustomModel_CustomModelParamValue = Ref{ProtoMeta}()
function meta(::Type{CustomModel_CustomModelParamValue})
    ProtoBuf.metalock() do
        if !isassigned(__meta_CustomModel_CustomModelParamValue)
            __meta_CustomModel_CustomModelParamValue[] = target = ProtoMeta(CustomModel_CustomModelParamValue)
            fnum = Int[10,20,30,40,50,60]
            allflds = Pair{Symbol,Union{Type,String}}[:doubleValue => Float64, :stringValue => AbstractString, :intValue => Int32, :longValue => Int64, :boolValue => Bool, :bytesValue => Vector{UInt8}]
            oneofs = Int[1,1,1,1,1,1]
            oneof_names = Symbol[Symbol("value")]
            meta(target, CustomModel_CustomModelParamValue, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, oneofs, oneof_names)
        end
        __meta_CustomModel_CustomModelParamValue[]
    end
end
function Base.getproperty(obj::CustomModel_CustomModelParamValue, name::Symbol)
    if name === :doubleValue
        return (obj.__protobuf_jl_internal_values[name])::Float64
    elseif name === :stringValue
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :intValue
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :longValue
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :boolValue
        return (obj.__protobuf_jl_internal_values[name])::Bool
    elseif name === :bytesValue
        return (obj.__protobuf_jl_internal_values[name])::Vector{UInt8}
    else
        getfield(obj, name)
    end
end

mutable struct CustomModel_ParametersEntry <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function CustomModel_ParametersEntry(; kwargs...)
        obj = new(meta(CustomModel_ParametersEntry), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct CustomModel_ParametersEntry (mapentry)
const __meta_CustomModel_ParametersEntry = Ref{ProtoMeta}()
function meta(::Type{CustomModel_ParametersEntry})
    ProtoBuf.metalock() do
        if !isassigned(__meta_CustomModel_ParametersEntry)
            __meta_CustomModel_ParametersEntry[] = target = ProtoMeta(CustomModel_ParametersEntry)
            allflds = Pair{Symbol,Union{Type,String}}[:key => AbstractString, :value => CustomModel_CustomModelParamValue]
            meta(target, CustomModel_ParametersEntry, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_CustomModel_ParametersEntry[]
    end
end
function Base.getproperty(obj::CustomModel_ParametersEntry, name::Symbol)
    if name === :key
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :value
        return (obj.__protobuf_jl_internal_values[name])::CustomModel_CustomModelParamValue
    else
        getfield(obj, name)
    end
end

mutable struct CustomModel <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function CustomModel(; kwargs...)
        obj = new(meta(CustomModel), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct CustomModel
const __meta_CustomModel = Ref{ProtoMeta}()
function meta(::Type{CustomModel})
    ProtoBuf.metalock() do
        if !isassigned(__meta_CustomModel)
            __meta_CustomModel[] = target = ProtoMeta(CustomModel)
            fnum = Int[10,30,40]
            allflds = Pair{Symbol,Union{Type,String}}[:className => AbstractString, :parameters => Base.Dict{AbstractString,CustomModel_CustomModelParamValue}, :description => AbstractString]
            meta(target, CustomModel, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_CustomModel[]
    end
end
function Base.getproperty(obj::CustomModel, name::Symbol)
    if name === :className
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :parameters
        return (obj.__protobuf_jl_internal_values[name])::Base.Dict{AbstractString,CustomModel_CustomModelParamValue}
    elseif name === :description
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

export CustomModel_CustomModelParamValue, CustomModel_ParametersEntry, CustomModel
# mapentries: "CustomModel_ParametersEntry" => ("AbstractString", "CustomModel_CustomModelParamValue")


text:
function select_optimal_demand_curve(t_mask, p)
    local row_opt

    tau = 1                      # start time of the current interval
    hist_d = []
    hist_df = DataFrame()
    row_opt = 0
    for t in range(1, stop=T-1)  # simulation loop
        old_price = p
        realized_d = sample_actual_demand(p)
        append!(hist_d, realized_d)

        if t_mask[t] !== t_mask[t + 1]  # end of the interval
            interval_mean_d = mean(hist_d[tau : t])

            min_dist = Inf
            for i in 1:nrow(h_vec)  # search for the best hypothesis
                dist = abs(interval_mean_d - h_vec[i, :d](p))
                if dist < min_dist
                    min_dist = dist
                    row_opt = i
                end
            end

            p = h_vec[row_opt, :p_opt] # set price for the next interval
            tau = t + 1                # switch to the next interval
        end
        append!(hist_df,
            DataFrame(:price => old_price, :realized_d => realized_d,
                :opt_row => row_opt))
    end
    realized_d = sample_actual_demand(p)
    append!(hist_df,
        DataFrame(:price => p, :realized_d => realized_d,
            :opt_row => row_opt))
    (p, row_opt, hist_df)
end


text:
module TestTranspose
using ParallelAccelerator

@acc transpose_t(A) = A'

function test()
    A = [1. 2. 3.; 4. 5. 6.]
    B = transpose_t(A)
    return B
end

@acc function test2()
    A = [1 2 3; 4 5 6]
    B = transpose(A)
    return B
end

end

using Base.Test
println("testing transpose...")
@test_approx_eq TestTranspose.test() [1. 4.; 2. 5.; 3. 6.]
@test TestTranspose.test2()==[1 4; 2 5; 3 6]
println("Done testing transpose.")


text:
# Important information about oldTraj:
# ======================================
# oldTraj.oldTraj contains all state information of one lap. It is structured like follows:
# The first <prebuf> values are the end of one lap before the next lap starts (necessary for inter-lap-system-ID)
# The last value of <prebuf> ends with the last value *before* the finish line (s < s_target)
# The s-values within <prebuf> need to be below zero (-> before the finish line). Otherwise they can be mistaken as the end of the current (not previous) lap.
# After <prebuf> follow the recorded states of the trajectory (as many as there were during one lap)
# The number of recorded states during one trajectory (0 <= s < s_target) is equal to <oldCost>.
# After the recorded trajectory, the rest of the vector (until <buffersize>) is filled up with constant values

function saveOldTraj(oldTraj::OldTrajectory,zCurr::Array{Float64},uCurr::Array{Float64},lapStatus::LapStatus,posInfo::PosInfo,buffersize::Int64)
                println("Starting function")
                i               = lapStatus.currentIt-1         # i = number of points for 0 <= s < s_target (= cost of this lap)
                prebuf          = oldTraj.prebuf                # so many points of the end of the previous old traj will be attached to the beginning
                zCurr_export    = zeros(buffersize,6)
                uCurr_export    = zeros(buffersize,2)

                println("Creating exports")
                zCurr_export    = cat(1,oldTraj.oldTraj[oldTraj.oldCost[1]+1:oldTraj.oldCost[1]+prebuf,:,1],
                                        zCurr[1:i,:], NaN*ones(buffersize-i-prebuf,6))#[ones(buffersize-i-prebuf,1)*zCurr[i,1:5] zCurr[i,6]+collect(1:buffersize-i-prebuf)*dt*zCurr[i,1]])
                uCurr_export    = cat(1,oldTraj.oldInput[oldTraj.oldCost[1]+1:oldTraj.oldCost[1]+prebuf,:,1],
                                        uCurr[1:i,:], NaN*ones(buffersize-i-prebuf,2))#zeros(buffersize-i-prebuf,2))

                zCurr_export[1:prebuf,6] -= posInfo.s_target       # make the prebuf-values below zero
                costLap                   = i                      # the cost of the current lap is the time it took to reach the finish line
                println("Saving")
                # Save all data in oldTrajectory:
                if lapStatus.currentLap <= 2                        # if it's the first or second lap
                    oldTraj.oldTraj[:,:,1]  = zCurr_export          # ... just save everything
                    oldTraj.oldInput[:,:,1] = uCurr_export
                    oldTraj.oldTraj[:,:,2]  = zCurr_export
                    oldTraj.oldInput[:,:,2] = uCurr_export
                    oldTraj.oldCost = [costLap,costLap]
                else                                                # idea: always copy the new trajectory in the first array!
                    if oldTraj.oldCost[1] < oldTraj.oldCost[2]      # if the first old traj is better than the second
                        oldTraj.oldTraj[:,:,2]  = oldTraj.oldTraj[:,:,1]    # ... copy the first in the second
                        oldTraj.oldInput[:,:,2] = oldTraj.oldInput[:,:,1]   # ... same for the input
                        oldTraj.oldCost[2] = oldTraj.oldCost[1]
                    end
                    oldTraj.oldTraj[:,:,1]  = zCurr_export                 # ... and write the new traj in the first
                    oldTraj.oldInput[:,:,1] = uCurr_export
                    oldTraj.oldCost[1] = costLap
                end
end

function InitializeParameters(mpcParams::MpcParams,mpcParams_pF::MpcParams,trackCoeff::TrackCoeff,modelParams::ModelParams,
                              posInfo::PosInfo,oldTraj::OldTrajectory,mpcCoeff::MpcCoeff,lapStatus::LapStatus,buffersize::Int64,
                              obstacle::Obstacle,selectedStates::SelectedStates,oldSS::SafeSetData)
    selectedStates.simulator = false     # set this to TRUE if SIMULATOR is in use, set this to FALSE if BARC is in use

    if selectedStates.simulator == false   # if the BARC is in use

        selectedStates.Np           = 15                           # Number of points to take from each previous trajectory to build the convex hull
        selectedStates.Nl           = 2                             # Number of previous laps to include in the convex hull
        selectedStates.shift        = 8
        Nl                          = selectedStates.Nl
        selectedStates.selStates    = zeros(Nl*selectedStates.Np,6)
        selectedStates.statesCost   = zeros(Nl*selectedStates.Np)
        selectedStates.version      = false

        mpcParams.N                 = 14
        mpcParams.Q                 = [5.0,0.0,0.0,0.1,50.0,0.0]   # Q (only for path following mode)
        mpcParams.vPathFollowing    = 1                           # reference speed for first lap of path following
        mpcParams.Q_term            = 1.0*[20.0,1.0,10.0,20.0,50.0]   # weights for terminal constraints (LMPC, for xDot,yDot,psiDot,ePsi,eY).Not used if using convex hull
        mpcParams.R                 = 0*[10.0,10.0]                 # put weights on a and d_f
        mpcParams.QderivZ           = 10.0*[1,1,1,1,1,1]             # cost matrix for derivative cost of states
        mpcParams.QderivU           = 1.0*[5.0,1.0] #NOTE Set this to [5.0, 0/40.0]              # cost matrix for derivative cost of inputs
        mpcParams.Q_term_cost       = 3                        # scaling of Q-function
        mpcParams.delay_df          = 3                             # steering delay
        mpcParams.delay_a           = 1                             # acceleration delay
        mpcParams.Q_lane            = 1                      # weight on the soft constraint for the lane
        mpcParams.Q_vel             = 1                    # weight on the soft constraint for the maximum velocity
        mpcParams.Q_slack           = 1*[5*20.0,20.0,10.0,30.0,80.0,50.0]#[20.0,10.0,10.0,30.0,80.0,50.0]  #vx,vy,psiDot,ePsi,eY,s
        mpcParams.Q_obs             = ones(Nl*selectedStates.Np)# weight to esclude some of the old trajectories

    elseif selectedStates.simulator == true  # if the simulator is in use

        selectedStates.Np           = 15                           # Number of points to take from each previous trajectory to build the convex hull
        selectedStates.Nl           = 2                             # Number of previous laps to include in the convex hull
        selectedStates.shift        = 5
        Nl                          = selectedStates.Nl
        selectedStates.selStates    = zeros(Nl*selectedStates.Np,6)
        selectedStates.statesCost   = zeros(Nl*selectedStates.Np)
        selectedStates.version      = false

        mpcParams.N                 = 13
        mpcParams.Q                 = [5.0,0.0,0.0,0.1,50.0,0.0]   # Q (only for path following mode)
        mpcParams.vPathFollowing    = 1                           # reference speed for first lap of path following
        mpcParams.Q_term            = 1.0*[20.0,1.0,10.0,20.0,50.0]   # weights for terminal constraints (LMPC, for xDot,yDot,psiDot,ePsi,eY).Not used if using convex hull
        mpcParams.R                 = 0*[10.0,10.0]                 # put weights on a and d_f
        mpcParams.QderivZ           = 10.0*[1,1,1,1,1,1]             # cost matrix for derivative cost of states
        mpcParams.QderivU           = 1.0*[1,1.0] #NOTE Set this to [5.0, 0/40.0]              # cost matrix for derivative cost of inputs
        mpcParams.Q_term_cost       = 1                        # scaling of Q-function
        mpcParams.delay_df          = 3                             # steering delay
        mpcParams.delay_a           = 1                             # acceleration delay
        mpcParams.Q_lane            = 3                      # weight on the soft constraint for the lane
        mpcParams.Q_vel             = 1                    # weight on the soft constraint for the maximum velocity
        mpcParams.Q_slack           = 1*[20.0,1.0,10.0,30.0,80.0,50.0]#[20.0,10.0,10.0,30.0,80.0,50.0]  #vx,vy,psiDot,ePsi,eY,s
        mpcParams.Q_obs             = ones(Nl*selectedStates.Np)# weight to esclude some of the old trajectories
    end

    mpcParams_pF.N              = 16
    mpcParams_pF.Q              = [0.0,50.0,1.0,10.0]
    mpcParams_pF.R              = 0*[1.0,1.0]               # put weights on a and d_f
    mpcParams_pF.QderivZ        = 0.0*[0.0,0,1.0,0]           # cost matrix for derivative cost of states
    mpcParams_pF.QderivU        = 10*[1,1]                # cost matrix for derivative cost of inputs
    mpcParams_pF.vPathFollowing = 1                       # reference speed for first lap of path following
    mpcParams_pF.delay_df       = 3                         # steering delay (number of steps)
    mpcParams_pF.delay_a        = 1                         # acceleration delay

    trackCoeff.nPolyCurvature   = 8                         # 4th order polynomial for curvature approximation
    trackCoeff.coeffCurvature   = zeros(trackCoeff.nPolyCurvature+1)         # polynomial coefficients for curvature approximation (zeros for straight line)
    trackCoeff.width            = 0.6                       # width of the track (60cm)

    modelParams.l_A             = 0.125
    modelParams.l_B             = 0.125
    modelParams.dt              = 0.1                   # sampling time, also controls the control loop, affects delay_df and Qderiv
    modelParams.m               = 1.98
    modelParams.I_z             = 0.03
    modelParams.c_f             = 0.5                   # friction coefficient: xDot = - c_f*xDot (aerodynamic+tire)

    posInfo.s_target            = 19.11

    oldTraj.oldTraj             = NaN*ones(buffersize,7,30)
    oldTraj.oldInput            = zeros(buffersize,2,30)
    oldTraj.oldTimes            = NaN*ones(buffersize,30)
    oldTraj.count               = ones(30)*2
    oldTraj.oldCost             = ones(Int64,30)                   # dummies for initialization
    oldTraj.prebuf              = 30
    oldTraj.postbuf             = 30
    oldTraj.idx_start           = zeros(30)
    oldTraj.idx_end             = zeros(30)

    oldSS.oldSS                 = NaN*ones(buffersize,7,30)          # contains data from previous laps usefull to build the safe set
    oldSS.oldSS_xy              = NaN*ones(buffersize,4,30)
    oldSS.cost2target           = zeros(buffersize,30)     # cost to arrive at the target, i.e. how many iterations from the start to the end of the lap
    oldSS.oldCost               = ones(Int64,30)              # contains costs of laps
    oldSS.count                 = ones(30)*2               # contains the counter for each lap
    oldSS.prebuff                = 30
    oldSS.postbuff               = 50
    oldSS.idx_start             = ones(30)            # index of the first measurement with s > 0
    oldSS.idx_end               = zeros(30)              # index of the last measurement with s < s_target

    mpcCoeff.order              = 5
    mpcCoeff.coeffCost          = zeros(mpcCoeff.order+1,2)
    mpcCoeff.coeffConst         = zeros(mpcCoeff.order+1,2,5)
    mpcCoeff.pLength            = 5*2*mpcParams.N        # small values here may lead to numerical problems since the functions are only approximated in a short horizon
    mpcCoeff.c_Vx               = zeros(3)
    mpcCoeff.c_Vy               = zeros(4)
    mpcCoeff.c_Psi              = zeros(3)

    lapStatus.currentLap        = 1         # initialize lap number
    lapStatus.currentIt         = 0         # current iteration in lap

    obstacle.obstacle_active    = false     # true if we have to consider the obstacles in the optimization problem
    obstacle.lap_active         = 1000         # number of the first lap in which the obstacles are used
    obstacle.obs_detect         = 10         # maximum distance at which we can detect obstacles (in terms of s!!)
    obstacle.n_obs              = 1         # number of obstacles
    obstacle.s_obs_init         = [14]    # initial s coordinate of each obstacle
    obstacle.ey_obs_init        = [-0.3]       # initial ey coordinate of each obstacle
    obstacle.v_obs_init         = [0]       # initial velocity of each obstacles
    obstacle.r_s                = 0.2#0.5
    obstacle.r_ey               = 0.1#0.2
    obstacle.inv_step           = 0         # number of step of invariance required for the safe set
end

# Use this function to smooth data (moving average)
function smooth(x,n)
    y = zeros(size(x))
    for i=1:size(x,1)
        start = max(1,i-n)
        fin = min(size(x,1),start + 2*n)
        y[i,:] = mean(x[start:fin,:],1)
    end
    return y
end

max_stars_repo_name: RoyaFiroozi/barc

max_issues_repo_name: RoyaFiroozi/barc

max_forks_repo_name: RoyaFiroozi/barc


text:
module ScalarProducts

using Plots
using LinearAlgebra
using BenchmarkTools

"""
    scalar_product_v1(x, y)

Compute the scalar product of two vectors using the for loop.
"""
function scalar_product_v1(x::Array{Float64,1}, y::Array{Float64,1})
    @assert length(x) == length(y)

    tmp::Float64 = 0.

    for i in 1:length(x)
        tmp += x[i] * y[i]
    end

    return sum(tmp)
end

"""
    scalar_product_v2(x, y)

Compute the scalar product of two vectors using the list comprehension.
"""
function scalar_product_v2(x::Array{Float64,1}, y::Array{Float64,1})
    return sum([i * j for (i,j) in [x, y]])
end

"""
    scalar_product_v3(x, y)

Compute the scalar product of two vectors using element-wise `.*` operation.
"""
function scalar_product_v3(x::Array{Float64,1}, y::Array{Float64,1})
    return sum(x .* y)
end

"""
    scalar_product_v4(x, y)

Compute the scalar product of two vectors using `LinearAlgebra.dot()` function.
"""
function scalar_product_v4(x::Array{Float64,1}, y::Array{Float64,1})
    return dot(x, y)
end

function scalar_product_v5(x::Array{Float64,1}, y::Array{Float64,1})
    s = 0.
    foreach((v,w) -> s += v*w, x, y)
    return s
end

"""
    test(F, N)

Generates a couple of random vectors of length `N:100:M` and plot the time necessary
for a functions from your module to calculate the result.
"""
function test(N::Int64, M::Int64)
    # Don't call @time not in global scope if you
    # like to get proper allocation info from time macro, see
    # https://stackoverflow.com/a/46263352 and https://stackoverflow.com/a/46660228.

    ns = N:(10^4):M

    fs = Dict(
        "v1" => scalar_product_v1,
        "v2" => scalar_product_v2,
        "v3" => scalar_product_v3,
        "v4" => scalar_product_v4,
        # scalar_product_v5,
    )

    gr()    # Empty plot with GR backend.
    plot()

    for (name, f) in fs
        # The function `f` is compiled here because
        # you like to exclude compile time from test.

        f(rand(Float64, N), rand(Float64, N))
        println( f([2.0, 1.0], [1.0, 2.0])) # simple test
        ts = Vector{Float64}()

        for n in ns
            x = rand(Float64, n)
            y = rand(Float64, n)
            t = @elapsed f(x, y)
            append!(ts, t)
        end
        plot!(ns, ts, title = "Benchmarks", label=name)

    end

    savefig("Benchmarks.png")
end

end # module

max_stars_repo_name: groundf/course-julia

max_issues_repo_name: groundf/course-julia

max_forks_repo_name: groundf/course-julia


text:
#using Tracker
using Random, Flux
using MLDataUtils: kfolds

struct MLPClassifier
    net::Flux.Chain
    labels::Vector{<:Union{Int, Float64}}
    nfeatures::Int
    function MLPClassifier(number_features::Int, labels::Vector{<:Union{Int, Float64}},
                           hidden_layers::Vector{Int}=Int[])
        model = build_classifier_model(number_features, length(labels), hidden_layers, tanh)
        new(model, labels, number_features)
    end
    function MLPClassifier(number_features::Int, labels::Vector{<:Union{Int, Float64}},
                           weights, activations, hidden_layers::Vector{Int}=Int[])
        model = build_classifier_model(number_features, length(labels),  weights,
                                       activations, hidden_layers)
        new(model, labels, number_features)
    end
end

function build_classifier_model(n_inputs::Int, n_outputs::Int, hidden_layers::Vector{Int}=Int[],
                                activation::Function=identity)::Flux.Chain
    layer_sizes = [n_inputs; hidden_layers; n_outputs]
    nlayers = length(layer_sizes)
    layers = map(2:nlayers) do i
        size_below = layer_sizes[i-1]
        size_above = layer_sizes[i]
        if i != nlayers
            Dense(size_below, size_above, activation)
        else
            Dense(size_below, size_above)
        end
    end
    Chain(layers..., softmax)
end

function build_classifier_model(n_inputs::Int, n_outputs::Int,  weights, activations,
                                hidden_layers::Vector{Int}=Int[])::Flux.Chain
    layer_sizes = [n_inputs; hidden_layers; n_outputs]
    nlayers = length(layer_sizes)
    layers = map(2:nlayers) do i
        (w, b) = weights[i-1]
        act = activations[i-1]
        size_below = layer_sizes[i-1]
        size_above = layer_sizes[i]
        Dense(size_below, size_above, get_activation(act); initW = (out, in) -> w, initb = out -> b)
    end
    Chain(layers..., softmax)
end

function get_activation(f::Symbol)
    act_functions = [:tanh, :sigmoid, :σ, :relu, :identity]
    if any(map((s) -> f == s, act_functions))
        return eval(f)
    else
        return identity
    end
end

(m::MLPClassifier)(input::Vector{<:Real}) = m.net(input)

function calculate_topology(model::MLPClassifier)
    net = model.net
    number_layers = length(net)
    topology = Vector{Int}(undef, number_layers)
    output_layer = 0
    for l = 1:(number_layers-1)
        output, input = size(net[l].W)
        topology[l] = input
        if l == (number_layers-1) output_layer = output end
    end
    topology[number_layers] = output_layer
    return topology
end

number_hidden_layers(m::MLPClassifier) = length(m.net) - 2

function fit!(model::MLPClassifier, inputs, outputs ; max_epochs=1)
    function lossMSE(x, y)
        ŷ = model.net(x)
        #Flux.mse(ŷ, Flux.onehot(y[1], model.labels))
        Flux.mse(ŷ, y)
    end

    function loss(x, y)
        ŷ = model.net(x)
        Flux.crossentropy(ŷ, y)
    end

    function loss_all(data, model)
        l = 0f0
        for (x,y) in data
            l += loss(x, y)
        end
        l/length(data)
    end

    #cb = function ()
    #    accuracy() > 0.9 && Flux.stop()
    #end

    η = 0.2
    opt = Descent(η)
    params = Flux.params(model.net)
    new_outputs = convert_to_onehot(outputs, model.labels)
    #println("Outputs:")
    #println(new_outputs)
    dataset = zip(inputs, new_outputs) # (shuffle ∘ collect ∘ zip)

    for ep = 1:max_epochs
        Flux.train!(loss, params, (shuffle ∘ collect)(dataset), opt,  cb = Flux.throttle(() -> print(""), 1800))
    end
    #Flux.@epochs max_epochs Flux.train!(loss, params, (shuffle ∘ collect)(dataset), opt,  cb = Flux.throttle(() -> println("training"), 10))

    model
end

function early_stop_fit!(model::MLPClassifier, inputs, outputs ; max_epochs::Int=1, k_folds::Int=10)
    function lossMSE(x, y)
        ŷ = model.net(x)
        Flux.mse(ŷ, y)
    end

    function loss(x, y)
        ŷ = model.net(x)
        Flux.crossentropy(ŷ, y)
    end

    η = 0.2
    opt = Descent(η)
    #params = Flux.params(model.net)
    new_outputs = convert_to_onehot(outputs, model.labels)
    dataset = zip(inputs, new_outputs)

    # Early-stopping training using K-fold cross validation
    folds = kfolds(dataset, k = k_folds)
    record_loss_n_train = []
    record_loss_n_valid = []
    fold_select = 1
    early_stop = 0
    evaluate = (f, d) -> sum(map((s) -> f(s...), d))
    for epoch_idx in 1:max_epochs

        train, valid = folds[fold_select] # selection of the datasets

        evalcb = () -> (push!(record_loss_n_train, evaluate(loss, train)),
        push!(record_loss_n_valid, evaluate(loss, valid)))

        Flux.train!(loss, params(model), train, opt, cb = throttle(evalcb, 1))

        fold_select += 1 # for selecting the K-fold between 1 and the total number of folds
        if fold_select >= (k_folds+1)
            fold_select = 1
        end

        # for early stop
        if record_loss_n_valid[epoch_idx] > record_loss_n_valid[epoch_idx-1]
            early_stop += 1
        end
        if early_stop > 100
            break
        end
    end

    model
end

function k_fold_crossvalidation(model::MLPClassifier, inputs, outputs ; k_folds::Int=10)
    η = 0.2
    opt = Descent(η)

    function loss(m, x, y)
        ŷ = m.net(x)
        Flux.crossentropy(ŷ, y)
    end

    new_outputs = convert_to_onehot(outputs, model.labels)
    dataset = zip(inputs, new_outputs)
    folds = kfolds(dataset, k = k_folds)
    record_loss_n_valid = Float64[]
    for fold = 1:k_folds
        train, valid = folds[fold]
        m = deepcopy(model)
        evaluate = (f, d) -> sum(map((s) -> f(m, s...), d))
        Flux.train!(loss, params(m), train, opt, cb = throttle(evalcb, 1))
        push!(record_loss_n_valid, evaluate(loss, valid))
    end
    sum(record_loss_n_valid)/k_folds
end

function convert_to_onehot(output, labels)
    return [1*(Flux.onehot(elem[1], labels)) for elem in output]
    #=for i = 1:length(output)
        push!(new_output, Flux.onehot(output[i][1], labels))
    end
    new_output=#
end

function accuracy(data, model::MLPClassifier)
    acc = 0
    for (x,y) in data
        acc += sum(Flux.onecold(model.net(x)) .== Flux.onecold(Flux.onehot(y[1], model.labels)))*1 / size(x,2)
    end
    acc/length(data)
end

function accuracy(data, model::Chain, labels)
    acc = 0
    for (x,y) in data
        acc += sum(Flux.onecold(model(x)) .== Flux.onecold(Flux.onehot(y[1], labels)))*1 / size(x,2)
    end
    acc/length(data)
end

function error_rate(data, model::MLPClassifier)
    acc = 0
    for (x,y) in data
        acc += (predict(model, x) != y[1])
    end
    acc/length(data)
end

function predict(model::MLPClassifier, features)
    r = model.net(features)
    Flux.onecold(r, model.labels)
end

function my_custom_train!(loss, ps, data, opt)
  # training_loss is declared local so it will be available for logging outside the gradient calculation.
  local training_loss
  ps = Flux.Params(ps)
  for d in data
    gs = gradient(ps) do
      training_loss = loss(d...)
      # Code inserted here will be differentiated, unless you need that gradient information
      # it is better to do the work outside this block.
      return training_loss
    end
    # Insert whatever code you want here that nier.jl:88eeds training_loss, e.g. logging.
    # logging_callback(training_loss)
    # Insert what ever code you want here that needs gradient.
    # E.g. logging with TensorBoardLogger.jl as histogram so you can see if it is becoming huge.
    update!(opt, ps, gs)
    # Here you might like to check validation set accuracy, and break out to do early stopping.
  end
end

#=
model = Chain(Dense(5,10, σ), Dense(10,10, σ), Dense(10, 1))
function loss(x, y)
    x1 = deepcopy(x)
    y1 = deepcopy(y)
    Flux.Losses.mse(model(x1), y1)
end
params = Flux.params(model)
x_train, y_train, x_test, y_test = get_dataset1("../datasets/Drug-Classification")
#println("Acuracy before training: ", accuracy(zip(x_test, y_test), model, dataset1_labels))

η = 0.2
opt = Descent(η)

#new_outputs = convert_to_onehot(y_train, dataset1_labels)
#const x = x_train
#const y = y_train
dataset = zip(x_train, y_train)
#my_custom_train!(loss, params, dataset, opt)
#Flux.train!(loss, params, dataset, opt)
=#

max_stars_repo_name: Lucas-Muniz/MNEAT

max_issues_repo_name: Lucas-Muniz/MNEAT

max_forks_repo_name: Lucas-Muniz/MNEAT


text:
export get_bound_constraints

"""
sets up bound constraints, based on global min and max values,
possibly specific bounds for the water layer in marine seismic imaging
Only returns a vector of the vectorized model size if necessary,
because projections onto bounds using a scalar is faster
"""
function get_bound_constraints(comp_grid,constraint)

if typeof(constraint["m_min"])<:Real && typeof(constraint["m_max"])<:Real && haskey(constraint,"water_depth")==false
  LB = constraint["m_min"]
  UB = constraint["m_max"]
else
  TF=typeof(constraint["m_min"][1])
  if length(comp_grid.n)==3 && comp_grid.n[3]>1 #use 3D version
    LB = ones(TF,comp_grid.n[1],comp_grid.n[2],comp_grid.n[3])*constraint["m_min"];
    UB = ones(TF,comp_grid.n[1],comp_grid.n[2],comp_grid.n[3])*constraint["m_max"];

    if haskey(constraint,"water_depth")==1
      water_bottom_index = convert(Int,floor(constraint["water_depth"]/comp_grid.d[3]));
      if water_bottom_index==0
        water_bottom_index=1
      end
      LB_water = ones(TF,comp_grid.n[1],comp_grid.n[2],comp_grid.n[3]) .* constraint["water_min"];
      UB_water = ones(TF,comp_grid.n[1],comp_grid.n[2],comp_grid.n[3]) .* constraint["water_max"];
      LB_water[:,:,water_bottom_index+1:end]=LB[:,:,water_bottom_index+1:end];
      UB_water[:,:,water_bottom_index+1:end]=UB[:,:,water_bottom_index+1:end];

      #now find tightest combination
      LB=max.(vec(LB),vec(LB_water))
      UB=min.(vec(UB),vec(UB_water))
    end #end setting water layer bounds 3D
  else # use 2D version
    LB = ones(TF,comp_grid.n[1],comp_grid.n[2])*constraint["m_min"]
    UB = ones(TF,comp_grid.n[1],comp_grid.n[2])*constraint["m_max"]

    if haskey(constraint,"water_depth")==true
      water_bottom_index = convert(Int,floor(constraint["water_depth"]/comp_grid.d[2]));
      if water_bottom_index==0
        water_bottom_index=1
      end
      LB_water = ones(TF,comp_grid.n[1],comp_grid.n[2]) .* constraint["water_min"];
      UB_water = ones(TF,comp_grid.n[1],comp_grid.n[2]) .* constraint["water_max"];
      LB_water[:,water_bottom_index+1:end]=LB[:,water_bottom_index+1:end];
      UB_water[:,water_bottom_index+1:end]=UB[:,water_bottom_index+1:end];

      #now find tightest combination
      LB=max.(vec(LB),vec(LB_water))
      UB=min.(vec(UB),vec(UB_water))
    end #end setting water layer bounds 2D
  end #end 2D part of the function
  LB = vec(LB);
  UB = vec(UB);
end

return LB,UB
end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
"""
    adversarialProblem(C, c, d, Γ, X, α, pd)

Computes ``Adv`` with accuracy ϵ.

Check section 5.1 _Adversarial lower bound_ of [publication](https://arxiv.org/abs/1811.06719)
for more information about this algorithm.

# Arguments
- `C`: is a vector of nonnegative first stage costs.
- `c`: is a vector of a nonnegative nominal second stage costs.
- `d`: is a vector of maximal deviations of the costs from their nominal values.
- `Γ`: is a budget, or the amount of uncertainty, which can be allocated to the second stage costs.
- `X`: is a set of feasible solutions represented as a list functions, each of which accepts a list of JuMP variables as an argument and returns a JuMP linear constraint.
- `α`: fixed number belonging to ``[0, 1]``
- `pd`: an instance of [`ProblemDescriptor`](@ref)
"""
function adversarialProblem(C, c, d, Γ, X, α, pd)
    ϵ = getProperty("adversarialProblem.epsilon", parameterType = Float64)
    timeLimit = getProperty("adversarialProblem.timeLimit")
    numConstraints = 1

    Δt = @elapsed begin
        c₀ = initialScenario(c, d, Γ)
        (x, y, lb) = recoverableProblem(C, c₀, X, α, pd)

        (model, c̃ᵥ, c̃, t̃ᵥ, t̃) = relaxedAdversarialProblem(C, c, d, Γ, x, y)
        ub = t̃
    end
    while abs(ub - lb)/lb > ϵ && Δt <= timeLimit
        Δt += @elapsed begin

            (x, y, nlb) = recoverableProblem(C, c̃, X, α, pd)
            if lb < nlb
                lb = nlb
            end

            (model, c̃, t̃) = addConstraintAndSolve(model, C, c̃ᵥ, t̃ᵥ, x, y)
            ub = t̃

            numConstraints += 1
        end
    end

    @debug "$(numConstraints) constraints was added to this adversarial problem"

    lb
end

function relaxedAdversarialProblem(C, c, d, Γ, x, y)
    n = size(c, 1)

    model = Model(solver=CplexSolver(CPX_PARAM_TILIM = getProperty("adversarialProblem.timeLimit"),
        CPXPARAM_ScreenOutput = getProperty("adversarialProblem.cplexLogging")))

    @variable(model, t̃)
    if ndims(c) == 1
        @variable(model, c̃[1:n])
    elseif ndims(c) == 2
        m = size(c, 2)
        @variable(model, c̃[1:n, 1:m])
    else
        throw(ArgumentError("Unsupported number of dimensions ($(ndims(c)))"))
    end

    @objective(model, Max, t̃)

    @constraint(model, t̃ <= vecdot(C, x) + vecdot(c̃, y))

    @constraint(model, c̃ .>= c)
    @constraint(model, c̃ .<= c + d)
    @constraint(model, sum(c̃ - c) <= Γ)

    status = solve(model)
    return (model, c̃, getvalue(c̃), t̃, getvalue(t̃))
end

function addConstraintAndSolve(model, C, c̃, t̃, x, y)
    @constraint(model, t̃ <= vecdot(C, x) + vecdot(c̃, y))

    status = solve(model)

    return (model, getvalue(c̃), getvalue(t̃))
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: iot_jobs_data_plane
using AWS.Compat
using AWS.UUIDs

"""
    describe_job_execution(job_id, thing_name)
    describe_job_execution(job_id, thing_name, params::Dict{String,<:Any})

Gets details of a job execution.

# Arguments
- `job_id`: The unique identifier assigned to this job when it was created.
- `thing_name`: The thing name associated with the device the job execution is running on.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"executionNumber"`: Optional. A number that identifies a particular job execution on a
  particular device. If not specified, the latest job execution is returned.
- `"includeJobDocument"`: Optional. When set to true, the response contains the job
  document. The default is false.
"""
describe_job_execution(jobId, thingName; aws_config::AbstractAWSConfig=global_aws_config()) = iot_jobs_data_plane("GET", "/things/$(thingName)/jobs/$(jobId)"; aws_config=aws_config)
describe_job_execution(jobId, thingName, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_jobs_data_plane("GET", "/things/$(thingName)/jobs/$(jobId)", params; aws_config=aws_config)

"""
    get_pending_job_executions(thing_name)
    get_pending_job_executions(thing_name, params::Dict{String,<:Any})

Gets the list of all jobs for a thing that are not in a terminal status.

# Arguments
- `thing_name`: The name of the thing that is executing the job.

"""
get_pending_job_executions(thingName; aws_config::AbstractAWSConfig=global_aws_config()) = iot_jobs_data_plane("GET", "/things/$(thingName)/jobs"; aws_config=aws_config)
get_pending_job_executions(thingName, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_jobs_data_plane("GET", "/things/$(thingName)/jobs", params; aws_config=aws_config)

"""
    start_next_pending_job_execution(thing_name)
    start_next_pending_job_execution(thing_name, params::Dict{String,<:Any})

Gets and starts the next pending (status IN_PROGRESS or QUEUED) job execution for a thing.

# Arguments
- `thing_name`: The name of the thing associated with the device.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"statusDetails"`: A collection of name/value pairs that describe the status of the job
  execution. If not specified, the statusDetails are unchanged.
- `"stepTimeoutInMinutes"`: Specifies the amount of time this device has to finish
  execution of this job. If the job execution status is not set to a terminal state before
  this timer expires, or before the timer is reset (by calling UpdateJobExecution, setting
  the status to IN_PROGRESS and specifying a new timeout value in field stepTimeoutInMinutes)
  the job execution status will be automatically set to TIMED_OUT. Note that setting this
  timeout has no effect on that job execution timeout which may have been specified when the
  job was created (CreateJob using field timeoutConfig).
"""
start_next_pending_job_execution(thingName; aws_config::AbstractAWSConfig=global_aws_config()) = iot_jobs_data_plane("PUT", "/things/$(thingName)/jobs/$next"; aws_config=aws_config)
start_next_pending_job_execution(thingName, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_jobs_data_plane("PUT", "/things/$(thingName)/jobs/$next", params; aws_config=aws_config)

"""
    update_job_execution(job_id, status, thing_name)
    update_job_execution(job_id, status, thing_name, params::Dict{String,<:Any})

Updates the status of a job execution.

# Arguments
- `job_id`: The unique identifier assigned to this job when it was created.
- `status`: The new status for the job execution (IN_PROGRESS, FAILED, SUCCESS, or
  REJECTED). This must be specified on every update.
- `thing_name`: The name of the thing associated with the device.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"executionNumber"`: Optional. A number that identifies a particular job execution on a
  particular device.
- `"expectedVersion"`: Optional. The expected current version of the job execution. Each
  time you update the job execution, its version is incremented. If the version of the job
  execution stored in Jobs does not match, the update is rejected with a VersionMismatch
  error, and an ErrorResponse that contains the current job execution status data is
  returned. (This makes it unnecessary to perform a separate DescribeJobExecution request in
  order to obtain the job execution status data.)
- `"includeJobDocument"`: Optional. When set to true, the response contains the job
  document. The default is false.
- `"includeJobExecutionState"`: Optional. When included and set to true, the response
  contains the JobExecutionState data. The default is false.
- `"statusDetails"`:  Optional. A collection of name/value pairs that describe the status
  of the job execution. If not specified, the statusDetails are unchanged.
- `"stepTimeoutInMinutes"`: Specifies the amount of time this device has to finish
  execution of this job. If the job execution status is not set to a terminal state before
  this timer expires, or before the timer is reset (by again calling UpdateJobExecution,
  setting the status to IN_PROGRESS and specifying a new timeout value in this field) the job
  execution status will be automatically set to TIMED_OUT. Note that setting or resetting
  this timeout has no effect on that job execution timeout which may have been specified when
  the job was created (CreateJob using field timeoutConfig).
"""
update_job_execution(jobId, status, thingName; aws_config::AbstractAWSConfig=global_aws_config()) = iot_jobs_data_plane("POST", "/things/$(thingName)/jobs/$(jobId)", Dict{String, Any}("status"=>status); aws_config=aws_config)
update_job_execution(jobId, status, thingName, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_jobs_data_plane("POST", "/things/$(thingName)/jobs/$(jobId)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("status"=>status), params)); aws_config=aws_config)


text:
#=
The fraction 49/98 is a curious fraction, as an inexperienced mathematician in attempting to simplify it may incorrectly believe that 49/98 = 4/8, which is correct, is obtained by cancelling the 9s.

We shall consider fractions like, 30/50 = 3/5, to be trivial examples.

There are exactly four non-trivial examples of this type of fraction, less than one in value, and containing two digits in the numerator and denominator.

If the product of these four fractions is given in its lowest common terms, find the value of the denominator.
=#
function isdigitcancelling(a,b)
  da = digits(a)
  db = digits(b)
  if a / b > 1 || da[1] == db[1]
    return false
  end

  ra,rb = Int64[],Int64[]
  for i in intersect(da,db)
    if da[1] == db[2]
      push!(ra,da[2])
      push!(rb,db[1])
    end
  end
  return (length(ra) == 1 && length(rb) == 1 && (a/b) == (ra[1]/rb[1])) ? true : false
end

function calc()
  numerators = Int64[]
  denominators = Int64[]
  for i in 10:99
    for j in 10:99
      if isdigitcancelling(i,j)
        push!(numerators, i)
        push!(denominators, j)
      end
    end
  end
  f1,f2 = prod(numerators),prod(denominators)
  d = gcd(f1,f2)
  f1 /= d
  f2 /= d
  Int64(f2)
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
using TimedLoop

nl = 20
c = rand(nl)

@tml for j = 1:nl
    c[j] = j
    sleep(.1)
    println("looping $j")
end

@tml true for j = 1:nl
    c[j] = j
    sleep(.1)
    println("looping $j")
end

@tmlv for j = 1:nl
    c[j] = j
    sleep(.1)
    println("looping $j")
end

max_stars_repo_name: daviehh/TimedLoop

max_issues_repo_name: daviehh/TimedLoop

max_forks_repo_name: daviehh/TimedLoop


text:
module Utils

export comma_separated_to_histogram, parse_array

# Creates a histogram of the counts, one-indexed.
comma_separated_to_histogram(str::String) = begin
    ints = [parse(Int, a) for a in split(str, ",")]
    max = maximum(ints)
    foldl(ints, init=zeros(Int, max+1)) do acc, age
        acc[age+1] += 1
        acc
    end
end

# Parse a matrix
parse_array(input) = begin
    array = nothing
    for row in filter(!isempty, split(input, "\n"))
        numbers = [parse(Int8, digit) for digit in split(row, "")]
        if isnothing(array)
            array = zeros(Int8, length(numbers), 0)
        end
        array = hcat(array, numbers)
    end
    array
end

end


text:
# Run ttv computations:
#include("stepsize_precision.jl")

# Now run "high-resolution" simulations:
#dq = ttv_elements!(n,t0,h/16.,tmax,elements,tt1,count,0.0,0,0,rstar;fout="test_output_highres.txt",iout=1)
#dqbig = ttv_elements!(n,t0big,hbig/16,tmaxbig,elements_big,tt1big,count,big(0.0),0,0,rstarbig;fout="test_output_highres_big.txt",iout=1)

# Fits a 2D Keplerian to a set of data
using PyPlot
using Optim
using DelimitedFiles
using Statistics

# Kepler solver:
include("../../src/kepler.jl")

function return_xy(t,elements)
  P,ecc,tp,semi,omega = elements
  if ecc >= 1.0
    return 0.0,0.0
  end
  # First, compute M:
  M = (t-tp)*2pi/P
  # Now, compute eccentric anomaly:
  ekep = ekepler(M,ecc)
  # Now compute xp & yp (the planar values aligned with pericenter):
  xp = semi*(cos(ekep)-ecc)
  yp = semi*sqrt(1-ecc^2)*sin(ekep)
  # Now rotate:
  x = xp*cos(omega)-yp*sin(omega)
  y = xp*sin(omega)+yp*cos(omega)
return x,y
end

function fit_keplerian(t,x,y,elements0)

  # elements are given by: [P,ecc,tp,semi,omega]
  # Chi-square function:
  function chi_square(elements)
    # Call model:
    chisq = 0.0
    for i=1:length(t)
      xmod,ymod = return_xy(t[i],elements)
      chisq += (x[i]-xmod)^2+(y[i]-ymod)^2
    end
    #println("elements: ",elements," chi^2: ",chisq)
  return chisq
  end

  # Optimize the fit:
  result =  optimize(chi_square,elements0,LBFGS(),autodiff=:forward)
#  result =  optimize(chi_square,elements0,LBFGS())
#  result =  optimize(chi_square,elements0)
  elements_best = result.minimizer
  nt = length(t)
  xmod = zeros(nt)
  ymod = zeros(nt)
  for i=1:nt
     xmod[i],ymod[i] = return_xy(t[i],elements_best)
  end
return xmod,ymod, result
end

data_hr  = readdlm("test_output_highres.txt")
elements_all = readdlm("elements.txt",',')
semi = median(sqrt.(data_hr[:,5].^2 .+data_hr[:,7].^2))
period = elements_all[2,2]
elements_best = zeros(5)
chi_best = Inf
chisq = zeros(20,20)
tobs = data_hr[:,1]; xobs  = data_hr[:,5]-data_hr[:,2]; yobs = -(data_hr[:,7]-data_hr[:,4])
nobs = length(tobs)
xmod = zeros(nobs)
ymod = zeros(nobs)
k=1
for j=1:20
  tp = elements_all[2,3]+(j-0.5)/20*elements_all[2,2]
  M = (tobs[1]-tp)*2*pi/period
#  omega = atan2(-(yobs[1]*cos(M)-xobs[1]*sin(M))/semi,(yobs[1]*sin(M)+xobs[1]*cos(M))/semi)
  omega = atan(-(yobs[1]*cos(M)-xobs[1]*sin(M))/semi,(yobs[1]*sin(M)+xobs[1]*cos(M))/semi)
  for k=1:20
    omega = (k-0.5)/20*2*pi
    elements = [period,sqrt(elements_all[2,4]^2+elements_all[2,5]^2),tp,semi,omega]
  # P,ecc,tp,semi,omega = elements
    for i=1:length(tobs)
      xmod[i],ymod[i] = return_xy(tobs[i],elements)
    end
    chisq[j,k] =sum((xobs .- xmod).^2+(yobs .- ymod).^2)
    if chisq[j,k] < chi_best
      global chi_best = chisq[j,k]
      ibest = j
      global elements_best = elements
      println(j," chi_best: ",chi_best," elements_best: ",elements_best)
      clf()
      plot(xobs,xmod)
      plot(yobs,ymod)
    end
  end
end

xmod, ymod , res = fit_keplerian(tobs,xobs,yobs,elements_best)

# Now fit big-float data:
data_hr_big  = readdlm("test_output_highres_big.txt")
tobs_big = data_hr_big[:,1]; xobs_big  = data_hr_big[:,5]-data_hr_big[:,2]; yobs_big = -(data_hr_big[:,7]-data_hr_big[:,4])
elements_best_big = copy(res.minimizer)
xmod_big, ymod_big , res_big = fit_keplerian(tobs_big,xobs_big,yobs_big,elements_best_big)

clf()
#plot(tobs,xobs-xmod)
#plot(tobs,yobs-ymod)

# Plot the difference in angle between the n-body and Keplerian angles:
#dtheta = atan2(yobs.*xmod .- ymod.*xobs,xobs.*xmod+yobs.*ymod)
dtheta = atan.(yobs.*xmod .- ymod.*xobs,xobs.*xmod+yobs.*ymod)
xdot = data_hr[:,14]-data_hr[:,11]; ydot = -(data_hr[:,16]-data_hr[:,13])
# Compute the change in theta as a function of time:
dthetadt = (ydot.*xobs .-yobs.*xdot)./(xobs.^2 .+yobs.^2)
plot(tobs,-dtheta./dthetadt*24*3600)
#dtheta_big = atan2(yobs_big.*xmod_big .- ymod_big.*xobs_big,xobs_big.*xmod_big .+yobs_big.*ymod_big)
dtheta_big = atan.(yobs_big.*xmod_big .- ymod_big.*xobs_big,xobs_big.*xmod_big .+yobs_big.*ymod_big)
plot(tobs,-dtheta_big./dthetadt*24*3600)

max_forks_repo_name: langfzac/nbg-papers


text:
import Base: show

mutable struct SampleModel <: CmdStanModels
    name::AbstractString;              # Name of the Stan program
    model::AbstractString;             # Stan language model program
    n_chains::Vector{Int64};           # Number of chains
    seed::StanBase.RandomSeed;         # Seed section of cmd to run cmdstan
    init::StanBase.Init;               # Init section of cmd to run cmdstan
    output::StanBase.Output;           # Output section of cmd to run cmdstan
    tmpdir::AbstractString;            # Holds all below copied/created files
    output_base::AbstractString;       # Used for naming failes to be created
    exec_path::AbstractString;         # Path to the cmdstan excutable
    data_file::Vector{String};         # Array of data files input to cmdstan
    init_file::Vector{String};         # Array of init files input to cmdstan
    cmds::Vector{Cmd};                 # Array of cmds to be spawned/pipelined
    sample_file::Vector{String};       # Sample file array created by cmdstan
    log_file::Vector{String};          # Log file array created by cmdstan
    diagnostic_file::Vector{String};   # Diagnostic file array created by cmdstan
    summary::Bool;                     # Store cmdstan's summary as a .csv file
    printsummary::Bool;                # Print the summary
    cmdstan_home::AbstractString;      # Directory where cmdstan can be found
    method::Sample
    end

"""

Create a SampleModel

$(SIGNATURES)

# Extended help

### Required arguments
```julia
* `name::AbstractString`               : Name for the model
* `model::AbstractString`              : Stan model source
```

### Optional arguments
```julia
* `n_chains::Vector{Int64}=[4]`        : Optionally updated in stan_sample()
* `seed::RandomSeed`                   : Random seed settings
* `init::Init`                         : Interval bound for parameters
* `output::Output`                     : File output options
* `tmpdir::AbstractString`             : Directory where output files are stored
* `summary::Bool=true`                 : Create computed stan summary
* `printsummary::Bool=false`           : Print the summary
* `method::Sample                      : Will be Sample()
```
"""
function SampleModel(
    name::AbstractString,
    model::AbstractString,
    n_chains = [4];
    seed = -1,
    init = StanBase.Init(),
    output = StanBase.Output(),
    tmpdir = mktempdir(),
    summary::Bool = true,
    printsummary::Bool = false,
    method = Sample(),
)

    !isdir(tmpdir) && mkdir(tmpdir)

    StanBase.update_model_file(joinpath(tmpdir, "$(name).stan"), strip(model))

    output_base = joinpath(tmpdir, name)
    exec_path = StanBase.executable_path(output_base)
    cmdstan_home = get_cmdstan_home()

    error_output = IOBuffer()
    is_ok = cd(cmdstan_home) do
        success(pipeline(`make -f $(cmdstan_home)/makefile -C $(cmdstan_home) $(exec_path)`;
            stderr = error_output))
    end
    if !is_ok
        throw(StanModelError(name, String(take!(error_output))))
    end

    if typeof(seed) == RandomSeed
        nseed = seed
    elseif typeof(seed) == Int64
        if seed !== -1
            nseed = RandomSeed(seed)
        else
            nseed = RandomSeed()
        end
    end

    SampleModel(name, model, n_chains, nseed, init, output,
        tmpdir, output_base, exec_path, String[], String[],
        Cmd[], String[], String[], String[], summary, printsummary,
        cmdstan_home, method)
end

function Base.show(io::IO, ::MIME"text/plain", m::SampleModel)
    println(io, "  name =                    \"$(m.name)\"")
    println(io, "  n_chains =                $(StanBase.get_n_chains(m))")
    println(io, "  output =                  Output()")
    println(io, "    refresh =                 $(m.output.refresh)")
    println(io, "  tmpdir =                  \"$(m.tmpdir)\"")
    println(io, "  method =                  Sample()")
    println(io, "    num_samples =             ", m.method.num_samples)
    println(io, "    num_warmup =              ", m.method.num_warmup)
    println(io, "    save_warmup =             ", m.method.save_warmup)
    println(io, "    thin =                    ", m.method.thin)
    if isa(m.method.algorithm, Hmc)
        println(io, "    algorithm =               HMC()")
    if isa(m.method.algorithm.engine, Nuts)
        println(io, "      engine =                  NUTS()")
        println(io, "        max_depth =               ", m.method.algorithm.engine.max_depth)
    elseif isa(m.method.algorithm.engine, Static)
        println(io, "      engine =                  Static()")
        println(io, "        int_time =                ", m.method.algorithm.engine.int_time)
    end
    println(io, "      metric =                  ", typeof(m.method.algorithm.metric))
    println(io, "      stepsize =                ", m.method.algorithm.stepsize)
    println(io, "      stepsize_jitter =         ", m.method.algorithm.stepsize_jitter)
    else
        if isa(m.method.algorithm, Fixed_param)
            println(io, "    algorithm =               Fixed_param()")
        else
            println(io, "    algorithm =               Unknown")
        end
    end
    println(io, "    adapt =                   Adapt()")
    println(io, "      gamma =                   ", m.method.adapt.gamma)
    println(io, "      delta =                   ", m.method.adapt.delta)
    println(io, "      kappa =                   ", m.method.adapt.kappa)
    println(io, "      t0 =                      ", m.method.adapt.t0)
    println(io, "      init_buffer =             ", m.method.adapt.init_buffer)
    println(io, "      term_buffer =             ", m.method.adapt.term_buffer)
    println(io, "      window =                  ", m.method.adapt.window)
end


text:
# Tests for Cost Reverse DDF DEA Models
@testset "CostReverseDDFDEAModel" begin

    # Test using Book data
    X = [2 2; 1 4; 4 1; 4 3; 5 5; 6 1; 2 5; 1.6 8];
    Y = [1; 1; 1; 1; 1; 1; 1; 1];
    W = [1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1];

    # Cost Reverse DDF
    costrddf = deacostrddf(X, Y, W, :ERG)

    @test typeof(costrddf) == CostReverseDDFDEAModel

    @test nobs(costrddf) == 8
    @test ninputs(costrddf) == 2
    @test noutputs(costrddf) == 1
    @test ismonetary(costrddf) == false

    @test efficiency(costrddf) == efficiency(costrddf, :Economic)
    @test efficiency(costrddf, :Economic)   ≈ [0.0; 0.5; 0.5; 5/12; 0.6; 0.25; 0.525; 0.532609] atol = 1e-5
    @test efficiency(costrddf, :Technical)  ≈ [0.0; 0.0; 0.0; 5/12; 0.6; 1/6; 0.35; 7/16]
    @test efficiency(costrddf, :Allocative) ≈ [0.0; 0.5; 0.5; 0.0; 0.0; 1/12; 0.175; 0.095109] atol = 1e-5
    @test normfactor(costrddf) ≈ [4.0; 2.0; 2.0; 7.2; 10.0; 12.0; 40/7; 368/35]

    # With monetary = true
    costrddfmon = deacostrddf(X, Y, W, :ERG, monetary = true)

    @test ismonetary(costrddfmon) == true

    @test efficiency(costrddfmon) ≈ efficiency(costrddf) .* normfactor(costrddf)
    @test efficiency(costrddfmon, :Technical) ≈ efficiency(costrddf, :Technical) .* normfactor(costrddf)
    @test efficiency(costrddfmon, :Allocative) ≈ efficiency(costrddf, :Allocative) .* normfactor(costrddf)

    # With CRS (Same result in this example)
    costrddfcrs = deacostrddf(X, Y, W, :ERG, rts = :CRS)

    @test efficiency(costrddfcrs, :Economic)   ≈ [0.0; 0.5; 0.5; 5/12; 0.6; 0.25; 0.525; 0.532609] atol = 1e-5
    @test efficiency(costrddfcrs, :Technical)  ≈ [0.0; 0.0; 0.0; 5/12; 0.6; 1/6; 0.35; 7/16]
    @test efficiency(costrddfcrs, :Allocative) ≈ [0.0; 0.5; 0.5; 0.0; 0.0; 1/12; 0.175; 0.095109] atol = 1e-5

    # Print
    show(IOBuffer(), costrddf)

    # Test errors
    @test_throws DimensionMismatch deacostrddf([1; 2 ; 3], [4 ; 5], [1; 1; 1], :ERG) #  Different number of observations
    @test_throws DimensionMismatch deacostrddf([1; 2; 3], [4; 5; 6], [1; 2; 3; 4], :ERG) # Different number of observation in prices
    @test_throws DimensionMismatch deacostrddf([1 1; 2 2; 3 3 ], [4; 5; 6], [1 1 1; 2 2 2; 3 3 3], :ERG) # Different number of input prices and inputs
    @test_throws ArgumentError deacostrddf([1; 2; 3], [4; 5; 6], [1; 2; 3], :Error) # Invalid associated efficiency measure
    @test_throws ArgumentError deacostrddf([1; 2; 3], [4; 5; 6], [1; 2; 3], :ERG, rts = :Error) # Invalid returns to scale

end


text:
function iterate!(job::BasicMCJob, ::Type{AM}, ::Type{Univariate})
  job.sstate.count += 1

  if job.tuner.verbose
    job.sstate.tune.proposed += 1
  end

  if job.sstate.count <= job.sampler.t0
    set_normal!(job.sstate, job.sampler, job.pstate)
  else
    job.sstate.C = covariance(
      job.sstate.C, job.sstate.count-2, job.pstate.value, job.sstate.lastmean, job.sstate.secondlastmean
    )

    set_gmm!(job.sstate, job.sampler, job.pstate)
  end

  job.sstate.pstate.value = rand(job.sstate.proposal)

  job.parameter.logtarget!(job.sstate.pstate)

  job.sstate.ratio = job.sstate.pstate.logtarget-job.pstate.logtarget

  if job.sstate.count > job.sampler.t0
    job.sstate.ratio -= logpdf(job.sstate.proposal, job.sstate.pstate.value)
    set_gmm!(job.sstate, job.sampler, job.sstate.pstate)
    job.sstate.ratio += logpdf(job.sstate.proposal, job.pstate.value)
  end

  if (job.sstate.ratio > 0 || (job.sstate.ratio > log(rand())))
    job.pstate.value = job.sstate.pstate.value

    job.pstate.logtarget = job.sstate.pstate.logtarget
    if in(:loglikelihood, job.outopts[:monitor]) && job.parameter.loglikelihood! != nothing
      job.pstate.loglikelihood = job.sstate.pstate.loglikelihood
    end
    if in(:logprior, job.outopts[:monitor]) && job.parameter.logprior! != nothing
      job.pstate.logprior = job.sstate.pstate.logprior
    end

    if !isempty(job.sstate.diagnosticindices)
      if haskey(job.sstate.diagnosticindices, :accept)
        job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:accept]] = true
      end
    end
    if job.tuner.verbose
      job.sstate.tune.accepted += 1
    end
  else
    if !isempty(job.sstate.diagnosticindices)
      if haskey(job.sstate.diagnosticindices, :accept)
        job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:accept]] = false
      end
    end
  end

  job.sstate.secondlastmean = job.sstate.lastmean
  job.sstate.lastmean = recursive_mean(job.sstate.lastmean, job.sstate.count, job.pstate.value)

  if job.tuner.verbose
    if job.sstate.tune.totproposed <= job.range.burnin && mod(job.sstate.tune.proposed, job.tuner.period) == 0
      rate!(job.sstate.tune)
      println(
        "Burnin iteration ",
        job.fmt_iter(job.sstate.tune.totproposed),
        " of ",
        job.range.burnin,
        ": ",
        job.fmt_perc(100*job.sstate.tune.rate),
        " % acceptance rate"
      )
      reset_burnin!(job.sstate.tune)
    end
  end
end

function iterate!(job::BasicMCJob, ::Type{AM}, ::Type{Multivariate})
  job.sstate.count += 1

  if job.tuner.verbose
    job.sstate.tune.proposed += 1
  end

  if job.sstate.count <= job.sampler.t0
    set_normal!(job.sstate, job.sampler, job.pstate)
  else
    covariance!(
      job.sstate.C, job.sstate.C, job.sstate.count-2, job.pstate.value, job.sstate.lastmean,job.sstate.secondlastmean
    )

    job.sstate.C[:, :] = Hermitian(job.sstate.C)

    set_gmm!(job.sstate, job.sampler, job.pstate)
  end

  job.sstate.pstate.value[:] = rand(job.sstate.proposal)

  job.parameter.logtarget!(job.sstate.pstate)

  job.sstate.ratio = job.sstate.pstate.logtarget-job.pstate.logtarget

  if job.sstate.count > job.sampler.t0
    job.sstate.ratio -= logpdf(job.sstate.proposal, job.sstate.pstate.value)
    set_gmm!(job.sstate, job.sampler, job.sstate.pstate)
    job.sstate.ratio += logpdf(job.sstate.proposal, job.pstate.value)
  end

  if (job.sstate.ratio > 0 || (job.sstate.ratio > log(rand())))
    job.pstate.value = copy(job.sstate.pstate.value)
    job.pstate.logtarget = job.sstate.pstate.logtarget
    if in(:loglikelihood, job.outopts[:monitor]) && job.parameter.loglikelihood! != nothing
      job.pstate.loglikelihood = job.sstate.pstate.loglikelihood
    end
    if in(:logprior, job.outopts[:monitor]) && job.parameter.logprior! != nothing
      job.pstate.logprior = job.sstate.pstate.logprior
    end

    if !isempty(job.sstate.diagnosticindices)
      if haskey(job.sstate.diagnosticindices, :accept)
        job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:accept]] = true
      end
    end
    if job.tuner.verbose
      job.sstate.tune.accepted += 1
    end
  else
    if !isempty(job.sstate.diagnosticindices)
      if haskey(job.sstate.diagnosticindices, :accept)
        job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:accept]] = false
      end
    end
  end

  job.sstate.secondlastmean = copy(job.sstate.lastmean)
  recursive_mean!(job.sstate.lastmean, job.sstate.lastmean, job.sstate.count, job.pstate.value)

  if job.tuner.verbose
    if job.sstate.tune.totproposed <= job.range.burnin && mod(job.sstate.tune.proposed, job.tuner.period) == 0
      rate!(job.sstate.tune)
      println(
        "Burnin iteration ",
        job.fmt_iter(job.sstate.tune.totproposed),
        " of ",
        job.range.burnin,
        ": ",
        job.fmt_perc(100*job.sstate.tune.rate),
        " % acceptance rate"
      )
      reset_burnin!(job.sstate.tune)
    end
  end
end


text:
"""
Utilities for using AlphaZero.jl on RL environments that implement CommonRLInterface.jl.
"""
module CommonRLInterfaceWrapper

using ..AlphaZero
import CommonRLInterface
using Setfield
using Random: AbstractRNG
const RL = CommonRLInterface

#####
##### Wrappers
#####

mutable struct Env{E, H} <: AbstractGameEnv
  rlenv :: E
  last_reward :: Float64
  two_players :: Bool
  heuristic_value :: H
  vectorize_state :: Function
  symmetries :: Function
  render :: Function
  action_string :: Function
  parse_action :: Function
  read_state :: Function
end

mutable struct Spec{E, H} <: AbstractGameSpec
  env :: Env{E, H}
end

#####
##### Default functions and constructors
#####

function GI.vectorize_state(::RL.AbstractEnv, state)
  @assert isa(state, Array{<:Number}) "
  Your state is not a vector and so you should define `vectorize_state`"
  return convert(Array{Float32}, state)
end

default_vectorize_state(rl, state) = GI.vectorize_state(rl, state)

GI.symmetries(::RL.AbstractEnv, state) = Tuple{typeof(state), Vector{Int}}[]

default_symmetries(rl, state) = GI.symmetries(rl, state)

default_heurisic_value(rl)   = GI.heuristic_value(rl)
default_render(rl)           = GI.render(rl)
default_action_string(rl, a) = GI.action_string(rl, a)
default_parse_action(rl, s)  = GI.parse_action(rl, s)
default_read_state(rl)       = GI.read_state(rl)

"""
    Env(rlenv::CommonRLInterface.AbstractEnv; <kwargs>) <: AbstractGameEnv

Wrap an environment implementing the interface defined in CommonRLInterface.jl into
an `AbstractGameEnv`.

# Requirements

The following optional methods must be implemented for `rlenv`:

  - `clone`
  - `state`
  - `setstate!`
  - `valid_action_mask`
  - `player`
  - `players`

# Keyword arguments

The following optional functions from `GameInterface` are not present in
CommonRLInterface.jl and can be provided as keyword arguments:

  - `vectorize_state`: must be provided unless states already have type `Array{<:Number}`
  - `heuristic_value`
  - `symmetries`
  - `render`
  - `action_string`
  - `parse_action`
  - `read_state`

If `f` is not provided, the default implementation calls
`GI.f(::CommonRLInterface.AbstractEnv, ...)`.
"""
function Env(
    rlenv :: RL.AbstractEnv;
    heuristic_value = default_heurisic_value,
    vectorize_state = default_vectorize_state,
    symmetries      = default_symmetries,
    render          = default_render,
    action_string   = default_action_string,
    parse_action    = default_parse_action,
    read_state      = default_read_state)

  nplayers = length(RL.players(rlenv))
  @assert nplayers == 1 || nplayers == 2 "
  AlphaZero.jl only supports games with one or two players."
  two_players = nplayers == 2

  return Env(
    rlenv, 0., two_players,
    heuristic_value, vectorize_state, symmetries,
    render, action_string, parse_action, read_state)
end

"""
    Spec(rlenv::RL.AbstractEnv; kwargs...) = spec(Env(rlenv; kwargs...))
"""
Spec(rlenv::RL.AbstractEnv; funs...) = Spec(Env(rlenv; funs...))

#####
##### GameInterface API
#####

GI.clone(env::Env) = @set env.rlenv = RL.clone(env.rlenv)

GI.set_state!(env::Env, state) = RL.setstate!(env.rlenv, state)

function GI.init(spec::Spec, rng::AbstractRNG)
  env = GI.clone(spec.env)
  RL.reset!(env.rlenv, rng)
  return env
end

function GI.init(spec::Spec, state)
  env = GI.clone(spec.env)
  RL.setstate!(env.rlenv, state)
  return env
end

GI.spec(env::Env) = Spec(env)

# Queries on specifications

GI.two_players(spec::Spec) = spec.env.two_players

GI.actions(spec::Spec) = RL.actions(spec.env.rlenv)

GI.vectorize_state(spec::Spec, state) = spec.env.vectorize_state(spec.env.rlenv, state)

# Operations on environments

GI.current_state(env::Env) = RL.state(env.rlenv)

GI.game_terminated(env::Env) = RL.terminated(env.rlenv)

GI.white_playing(env::Env) = RL.player(env.rlenv) == 1

GI.actions_mask(env::Env) = RL.valid_action_mask(env.rlenv)

function GI.play!(env::Env, action)
  r = RL.act!(env.rlenv, action)
  env.last_reward = r
  return
end

GI.white_reward(env::Env) = env.last_reward

GI.heuristic_value(env::Env) = env.heuristic_value(env.rlenv)

# Symmetries

GI.symmetries(spec::Spec, state) = spec.env.symmetries(spec.env.rlenv, state)

# Interactive utilities

GI.render(env::Env) = env.render(env.rlenv)

GI.action_string(spec::Spec, a) = spec.env.action_string(spec.env.rlenv, a)

GI.parse_action(spec::Spec, s) = spec.env.parse_action(spec.env.rlenv, s)

GI.read_state(spec::Spec) = spec.env.read_state(spec.env.rlenv)

end


text:
using HDF5Logger
using Test
using HDF5 # For reading in wholesale time histories

# Set up some data to work with.
num_samples = 5
vector_data = [1., 2., 3.]
matrix_data = transpose([1 2 3; 4 5 6])
scalar_data = 1.
c           = 2

# Create the log.
file_name = "test_log.h5"
log = Log(file_name)
try

    # Test a vector
    add!(log, "/vector", vector_data, num_samples)
    log!(log, "/vector", vector_data)
    log!(log, "/vector", c * vector_data)

    # Test a matrix
    add!(log, "/group/matrix", matrix_data, num_samples, true) # Keep data.
    # log!(log, "/group/matrix", matrix_data) # First sample is kept this time!
    log!(log, "/group/matrix", c * matrix_data)

    # Test a scalar
    add!(log, "/group/scalar", scalar_data, num_samples)
    log!(log, "/group/scalar", scalar_data)
    log!(log, "/group/scalar", c * scalar_data)

catch err
    close!(log)
    rethrow(err)
end

close!(log)

# Now test that the right stuff happened. Use the regular HDF5 library to open
# the files we created. With the `do` block, the file will be automatically
# closed if anything goes wrong.
h5open(file_name, "r") do file

    vector_result = read(file, "/vector")
    matrix_result = read(file, "/group/matrix")
    scalar_result = read(file, "/group/scalar")

    @test(vector_result[:,1]   == vector_data)
    @test(matrix_result[:,:,1] == matrix_data)
    @test(scalar_result[1]     == scalar_data)

    @test(vector_result[:,2]   == c * vector_data)
    @test(matrix_result[:,:,2] == c * matrix_data)
    @test(scalar_result[2]     == c * scalar_data)

    @test(size(vector_result) == (3, num_samples))
    @test(size(matrix_result) == (3, 2, num_samples))
    @test(size(scalar_result) == (1, num_samples))

    @test(eltype(vector_result) == Float64)
    @test(eltype(matrix_result) == Int64)
    @test(eltype(scalar_result) == Float64)

end

# Clean up after ourselves.
rm(file_name)


text:
using Mimi

# Define a simple component
# =========================

@defcomp component1 begin

    # First define the state this component will hold
    savingsrate = Parameter()

    # Second, define the (optional) init function for the component
    function init(p, v, d)
    end

    # Third, define the run_timestep function for the component
    function run_timestep(p, v, d, t)
    end

end

# Create a model uses the component
# =================================

@defmodel begin

    m = Model()

    component(component1)
end

# Run model
# =========

run(m)

# Explore the model variables and parameters with the explorer UI
# =================================

explore(m)

# Access the variables in the model
# =================================