text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
function add_electric_tariff_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")
    r = Dict{String, Any}()
    m[Symbol("Year1UtilityEnergy"*_n)] = p.hours_per_timestep *
        sum(m[Symbol("dvGridPurchase"*_n)][ts, tier] for ts in p.time_steps, tier in 1:p.s.electric_tariff.n_energy_tiers)

    r["lifecycle_energy_cost"] = round(value(m[Symbol("TotalEnergyChargesUtil"*_n)]) * (1 - p.s.financial.offtaker_tax_pct), digits=2)
    r["year_one_energy_cost"] = round(value(m[Symbol("TotalEnergyChargesUtil"*_n)]) / p.pwf_e, digits=2)

    r["lifecycle_demand_cost"] = round(value(m[Symbol("TotalDemandCharges"*_n)]) * (1 - p.s.financial.offtaker_tax_pct), digits=2)
    r["year_one_demand_cost"] = round(value(m[Symbol("TotalDemandCharges"*_n)]) / p.pwf_e, digits=2)

    r["lifecycle_fixed_cost"] = round(m[Symbol("TotalFixedCharges"*_n)] * (1 - p.s.financial.offtaker_tax_pct), digits=2)
    r["year_one_fixed_cost"] = round(m[Symbol("TotalFixedCharges"*_n)] / p.pwf_e, digits=0)

    r["lifecycle_min_charge_adder"] = round(value(m[Symbol("MinChargeAdder"*_n)]) * (1 - p.s.financial.offtaker_tax_pct), digits=2)
    r["year_one_min_charge_adder"] = round(value(m[Symbol("MinChargeAdder"*_n)]) / p.pwf_e, digits=2)

    r["year_one_bill"] = r["year_one_energy_cost"] + r["year_one_demand_cost"] +
                                    r["year_one_fixed_cost"]  + r["year_one_min_charge_adder"]

    r["lifecycle_export_benefit"] = -1 * round(value(m[Symbol("TotalExportBenefit"*_n)]) * (1 - p.s.financial.offtaker_tax_pct), digits=2)
    r["year_one_export_benefit"] = -1 * round(value(m[Symbol("TotalExportBenefit"*_n)]) / p.pwf_e, digits=0)

    r["lifecycle_coincident_peak_cost"] = round(value(m[Symbol("TotalCPCharges"*_n)]), digits=2)
    r["year_one_coincident_peak_cost"] = round(r["lifecycle_coincident_peak_cost"] / p.pwf_e, digits=2)

    d["ElectricTariff"] = r
    nothing
end

function add_electric_tariff_results(m::JuMP.AbstractModel, p::MPCInputs, d::Dict; _n="")
    r = Dict{String, Any}()
    m[Symbol("energy_purchased"*_n)] = p.hours_per_timestep *
        sum(m[Symbol("dvGridPurchase"*_n)][ts] for ts in p.time_steps)

    r["energy_cost"] = round(value(m[Symbol("TotalEnergyChargesUtil"*_n)]), digits=2)

    r["demand_cost"] = round(value(m[Symbol("TotalDemandCharges"*_n)]), digits=2)

    r["export_benefit"] = -1 * round(value(m[Symbol("TotalExportBenefit"*_n)]), digits=0)

    d["ElectricTariff"] = r
    nothing
end

max_stars_repo_name: NREL/REoptLite

max_issues_repo_name: NREL/REoptLite

max_forks_repo_name: NREL/REoptLite


text:
using SafeTestsets
using SymPy
@safetestset "Basic function test" begin include("basic_test.jl") end
@safetestset "SEq for one variable" begin include("SEq1_test.jl") end
@safetestset "Polynomial" begin include("Polynomial_test.jl") end

ext: jl


text:
"""
A Column is a vector (or similar store) of data annotated by a column name
"""
immutable Column{Name, StorageType}
    data::StorageType
    function Column{T}(x::T)
        check_Column(Val{Name}, StorageType)
        new(convert(StorageType, x))
    end
end

@compat @inline (::Type{Column{Name}}){Name, StorageType}(x::StorageType) = Column{Name, StorageType}(x)
@compat @inline (::Type{Column{Name,StorageType}}){Name, StorageType}() = Column{Name, StorageType}(StorageType())

Base.convert{Name, T1, T2}(::Type{Column{Name, T2}}, x::Column{Name,T1}) = Column{Name,T2}(x.data)
Base.convert{Name, T1, T2}(::Type{Column{Name, T2}}, x::Cell{Name,T1}) = Column{Name,T2}([x.data]) # ??

@generated function check_Column{Name, StorageType}(::Type{Val{Name}}, ::Type{StorageType})
    if !isa(Name, Symbol)
        return :( error("Field name $Name should be a Symbol") )
    elseif Name == :Row
        return :( error("Field name cannot be :Row") )
    elseif eltype(StorageType) == StorageType
        warn("Column :$Name storage type $StorageType doesn't appear to be a storage container")
    end
    return nothing
end

@compat Base.:(==){Name}(col1::Column{Name}, col2::Column{Name}) = (col1.data == col2.data)

@inline rename{Name1, Name2, ElType}(x::Column{Name1, ElType}, ::Type{Val{Name2}}) = Column{Name2, ElType}(x.data)

@inline name{Name, StorageType}(::Column{Name, StorageType}) = Name
@inline name{Name}(::Type{Column{Name}}) = Name
@inline name{Name, StorageType}(::Type{Column{Name, StorageType}}) = Name
@inline Base.eltype{Name, StorageType}(::Column{Name, StorageType}) = eltype(StorageType)
@inline Base.eltype{Name, StorageType}(::Type{Column{Name, StorageType}}) = eltype(StorageType)
@inline storagetype{Name, StorageType}(::Column{Name, StorageType}) = StorageType
@inline storagetype{Name, StorageType}(::Type{Column{Name, StorageType}}) = StorageType

@inline nrow{Name, StorageType}(col::Column{Name, StorageType}) = length(col.data)
@inline ncol{Name, StorageType}(col::Column{Name, StorageType}) = 1
@inline ncol{C <: Column}(::Type{C}) = 1

Base.getindex{Name}(c::Column{Name}, ::Type{Val{Name}}) = c.data
Base.getindex{Name1, Name2}(c::Column{Name1}, ::Type{Val{Name2}}) = error("Tried to index column of name :$Name1 with name :$Name2")
Base.getindex{Name}(c::Column{Name}, ::Type{Val{:Row}}) = 1:length(c.data)

@inline Base.length{Name, StorageType}(col::Column{Name, StorageType}) = length(col.data)
Base.ndims{Name, StorageType}(col::Column{Name, StorageType}) = 1
Base.size{Name, StorageType}(col::Column{Name, StorageType}) = (length(col.data),)
Base.size{Name, StorageType}(col::Column{Name, StorageType}, i::Int) = i == 1 ? length(col.data) : error("Columns are one-dimensional")
Base.isempty{Name, StorageType}(col::Column{Name, StorageType}) = isempty(col.data)
Base.endof{Name, StorageType}(col::Column{Name, StorageType}) = endof(col.data)

# Iterators
Base.start{Name, StorageType}(col::Column{Name, StorageType}) = start(col.data)
Base.next{Name, StorageType}(col::Column{Name, StorageType}, i) = next(col.data, i)
Base.done{Name, StorageType}(col::Column{Name, StorageType}, i) = done(col.data, i)

# get/set index
Base.getindex{Name, StorageType}(col::Column{Name, StorageType}, idx::Int) = getindex(col.data, idx)
Base.getindex{Name, StorageType}(col::Column{Name, StorageType}, idx) = Column{Name}(getindex(col.data, idx))

Base.setindex!{Name, StorageType}(col::Column{Name, StorageType}, val::Cell{Name}, idx::Integer) = setindex!(col.data, val.data, idx)
Base.setindex!{Name, StorageType}(col::Column{Name, StorageType}, val::Column{Name}, idx) = setindex!(col.data, val.data, idx)
Base.setindex!{Name, StorageType}(col::Column{Name, StorageType}, val, idx) = setindex!(col.data, val, idx)

Base.unsafe_getindex{Name, StorageType}(col::Column{Name, StorageType}, idx::Int) = unsafe_getindex(col.data, idx)
Base.unsafe_getindex{Name, StorageType}(col::Column{Name, StorageType}, idx) = Column{Name}(unsafe_getindex(col.data, idx))

Base.unsafe_setindex!{Name, StorageType}(col::Column{Name, StorageType}, val::Cell{Name}, idx::Integer) = unsafe_setindex!(col.data, val.data, idx)
Base.unsafe_setindex!{Name, StorageType}(col::Column{Name, StorageType}, val::Column{Name}, idx) = unsafe_setindex!(col.data, val.data, idx)
Base.unsafe_setindex!{Name, StorageType}(col::Column{Name, StorageType}, val, idx) = unsafe_setindex!(col.data, val, idx)

# Mutators: push!, append!, pop!, etc
Base.pop!(col::Column) = pop!(col.data)
Base.shift!(col::Column) = shift!(col.data)

Base.push!(col::Column, data_in) = (push!(col.data, data_in); col)
Base.push!{Name}(col::Column{Name}, cell::Cell{Name}) = (push!(col.data, cell.data); col)
Base.push!{Name}(col::Column{Name}, i::Integer, cell::Cell) = error("Column with name :$Name don't match cell with name :$(name(cell))")
Base.unshift!(col::Column, data_in) = (unshift!(col.data, data_in); col)
Base.unshift!{Name}(col::Column{Name}, cell::Cell{Name}) = (unshift!(col.data, cell.data); col)
Base.unshift!{Name}(col::Column{Name}, i::Integer, cell::Cell) = error("Column with name :$Name don't match cell with name :$(name(cell))")

Base.append!(col::Column, data_in) = (append!(col.data, data_in); col)
Base.append!{Name}(col::Column{Name}, cell::Column{Name}) = (append!(col.data, cell.data); col)
Base.append!{Name}(col::Column{Name}, i::Integer, cell::Cell) = error("Column with name :$Name don't match column with name :$(name(cell))")
Base.prepend!(col::Column, data_in) = (prepend!(col.data, data_in); col)
Base.prepend!{Name}(col::Column{Name}, cell::Column{Name}) = (prepend!(col.data, cell.data); col)
Base.prepend!{Name}(col::Column{Name}, i::Integer, cell::Cell) = error("Column with name :$Name don't match column with name :$(name(cell))")

Base.insert!(col::Column, i::Integer, v) = (insert!(col.data, i, v); col)
Base.insert!{Name}(col::Column{Name}, i::Integer, cell::Cell{Name}) = (insert!(col.data, i, cell.data); col)
Base.insert!{Name}(col::Column{Name}, i::Integer, cell::Cell) = error("Column with name :$Name don't match cell with name :$(name(cell))")
Base.deleteat!(col::Column, i) = (deleteat!(col.data, i); col)
Base.splice!{Name}(col::Column{Name}, i::Integer) = splice!(col.data, i)
Base.splice!{Name}(col::Column{Name}, i::Integer, r) = splice!(col.data, i, r)
Base.splice!{Name}(col::Column{Name}, i) = Column{Name}(splice!(col.data, i))
Base.splice!{Name}(col::Column{Name}, i, r) = Column{Name}(splice!(col.data, i, r))

Base.empty!(col::Column) = (empty!(col.data); col)

# unique/unique! (union, etc??)

Base.sort{Name, StorageType}(col::Column{Name, StorageType}; kwargs...) = Column{Name, StorageType}(sort(col.data; kwargs...))
Base.sort!(col::Column; kwargs...) = (sort!(col.data; kwargs...); col)

# Concatenate cells and columns into colums
Base.vcat{Name}(c1::Union{Cell{Name}, Column{Name}}) = Column{Name}(vcat(c1.data))
Base.vcat{Name}(c1::Union{Cell{Name}, Column{Name}}, c2::Union{Cell{Name}, Column{Name}}) = Column{Name}(vcat(c1.data, c2.data))
@generated function Base.vcat{Name}(c1::Union{Cell{Name}, Column{Name}}, c2::Union{Cell{Name}, Column{Name}}, cs::Union{Cell{Name}, Column{Name}}...)
    # Do our best to help inference here
    exprs = [:(cs[$j].data) for j = 1:length(cs)]
    vcat_expr = Expr(:call, :vcat, :(c1.data), :(c2.data), exprs...)
    return Expr(:call, Column{Name}, vcat_expr)
end

# copy
Base.copy{Name, StorageType}(col::Column{Name, StorageType}) = Column{Name, StorageType}(copy(col.data))

# @Column and @Cell are very similar
macro Column(expr)
    if expr.head != :(=) && expr.head != :(kw) # strange Julia bug, see issue 7669
        error("A Expecting expression like @Column(name::Type = value) or @Column(name = value)")
    end
    local field
    value = expr.args[2]
    if isa(expr.args[1], Symbol)
        name = expr.args[1]
        return :( TypedTables.Column{$(QuoteNode(name))}($(esc(value))) )
    elseif isa(expr.args[1],Expr)
        if expr.args[1].head != :(::) || length(expr.args[1].args) != 2 || !isa(expr.args[1].args[1], Symbol)
            error("B Expecting expression like @Column(name::Type = value) or @Column(name = value)")
        end
        name = expr.args[1].args[1]
        eltype = expr.args[1].args[2]
        field = :( TypedTables.Column{$(QuoteNode(name)), $(esc(eltype))}($(esc(value))) )
    else
        error("C Expecting expression like @Column(name::Type = value) or @Column(name = value)")
    end
end


text:
#!/usr/bin/env julia
include("sedor20.jl")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepackage{fontspec}")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\setmainfont{Asana-Math}")
cas = [Dict(:label => "tumör", :caexpr => "C00-D48"),
    Dict(:label => "COVID-19", :caexpr => "U07-U07"),
    Dict(:label => "cirkulation", :caexpr=>"I00-I99|F01")]
allca = Dict(:caexpr => "A00-Y98", :label => "alla orsaker")
popca = Dict(:caexpr => "pop", :label => "folkmängd")
pgfsave("../../images/se_cirk_cov_tum_prop_2020h1.svg",
    cas_sexes_prop_plot(cas, allca))
pgfsave("../../images/se_cirk_cov_tum_rate_2020h1.svg",
    cas_sexes_prop_plot(cas, popca, "log"))

max_stars_repo_name: klpn/static-dust

max_issues_repo_name: klpn/static-dust

max_forks_repo_name: klpn/static-dust


text:
using TrajOptPlots
using RobotZoo: Cartpole
using MeshCat
using StaticArrays

# Start visualizer
vis = Visualizer()
open(vis)

# Display the model
model = Cartpole()
TrajOptPlots.set_mesh!(vis, model)

# Visualize a single state
x = SA[0.5, pi/3, 0, 0]
visualize!(vis, model, x)

# Visualize a trajectory
tf = 2.0
X = [SA[sin(t), cos(t), 0, 0] for t in range(0, 2.0, length=21)]
visualize!(vis, model, tf, X)


text:
# CUFFT API function return values
const cufftStatus_t = UInt32
const CUFFT_STATUS_SUCCESS        = 0  #  The cuFFT operation was successful
const CUFFT_STATUS_INVALID_PLAN   = 1  #  cuFFT was passed an invalid plan handle
const CUFFT_STATUS_ALLOC_FAILED   = 2  #  cuFFT failed to allocate GPU or CPU memory
const CUFFT_STATUS_INVALID_TYPE   = 3  #  No longer used
const CUFFT_STATUS_INVALID_VALUE  = 4  #  User specified an invalid pointer or parameter
const CUFFT_STATUS_INTERNAL_ERROR = 5  #  Driver or internal cuFFT library error
const CUFFT_STATUS_EXEC_FAILED    = 6  #  Failed to execute an FFT on the GPU
const CUFFT_STATUS_SETUP_FAILED   = 7  #  The cuFFT library failed to initialize
const CUFFT_STATUS_INVALID_SIZE   = 8  #  User specified an invalid transform size
const CUFFT_STATUS_UNALIGNED_DATA = 9  #  No longer used
const CUFFT_STATUS_INCOMPLETE_PARAMETER_LIST = 10 #  Missing parameters in call
const CUFFT_STATUS_INVALID_DEVICE = 11 #  Execution of a plan was on different GPU than plan creation
const CUFFT_STATUS_PARSE_ERROR    = 12 #  Internal plan database error
const CUFFT_STATUS_NO_WORKSPACE   = 13  #  No workspace has been provided prior to plan execution
const CUFFT_STATUS_NOT_IMPLEMENTED = 14 # Function does not implement functionality for parameters given.
const CUFFT_STATUS_LICENSE_ERROR  = 15 # Used in previous versions.
const CUFFT_STATUS_NOT_SUPPORTED  = 16  # Operation is not supported for parameters given.

const cufftReal = Float32
const cufftDoubleReal = Float64

const cufftComplex = ComplexF32
const cufftDoubleComplex = ComplexF64

# CUFFT transform directions
const CUFFT_FORWARD = -1 # Forward FFT
const CUFFT_INVERSE =  1 # Inverse FFT

# CUFFT supports the following transform types
const cufftType = Cint
const CUFFT_R2C = 0x2a     # Real to Complex
const CUFFT_C2R = 0x2c     # Complex to Real
const CUFFT_C2C = 0x29     # Complex to Complex
const CUFFT_D2Z = 0x6a     # Double to Double-Complex
const CUFFT_Z2D = 0x6c     # Double-Complex to Double
const CUFFT_Z2Z = 0x69     # Double-Complex to Double-Complex

const cufftCompatibility = Cint
const   CUFFT_COMPATIBILITY_NATIVE          = 0x00
const   CUFFT_COMPATIBILITY_FFTW_PADDING    = 0x01
const   CUFFT_COMPATIBILITY_FFTW_ASYMMETRIC = 0x02
const   CUFFT_COMPATIBILITY_FFTW_ALL        = 0x03

const cufftHandle_t = Cint


text: f (1)

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %% [markdown]
# Straightforward translation of https://github.com/rmeinl/apricot-julia/blob/5f130f846f8b7f93bb4429e2b182f0765a61035c/notebooks/python_reimpl.ipynb
#
# See also https://github.com/genkuroki/public/blob/main/0016/apricot/python_reimpl.ipynb

# %%
using Seaborn
using ScikitLearn: @sk_import
@sk_import datasets: fetch_covtype
using Random
using StatsBase: sample

# %%
digits_data = fetch_covtype()

# %%
X_digits = permutedims(abs.(digits_data["data"]))
summary(X_digits)

# %%
"""`calculate_gains!(X, gains, current_values, idxs, current_concave_values_sum)` mutates `gains` only"""
function calculate_gains!(X, gains, current_values, idxs, current_concave_values_sum)
    Threads.@threads for i in eachindex(idxs)
        @inbounds idx = idxs[i]
        @inbounds gains[i] = sum(j -> sqrt(current_values[j] + X[j, idx]), axes(X, 1))
    end
    gains .-= current_concave_values_sum
end

@doc calculate_gains!

# %%
function fit(X, k; calculate_gains! = calculate_gains!)
    d, n = size(X)

    cost = 0.0

    ranking = Int[]
    total_gains = Float64[]

    current_values = zeros(d)
    current_concave_values_sum = sum(sqrt, current_values)

    idxs = collect(1:n)

    gains = zeros(n)
    while cost < k
        calculate_gains!(X, gains, current_values, idxs, current_concave_values_sum)

        idx = argmax(gains)
        best_idx = idxs[idx]
        curr_cost = 1.0

        cost + curr_cost > k && break

        cost += curr_cost
        # Calculate gains
        gain = gains[idx] * curr_cost

        # Select next
        current_values .+= @view X[:, best_idx]
        current_concave_values_sum = sum(sqrt, current_values)

        push!(ranking, best_idx)
        push!(total_gains, gain)

        popat!(idxs, idx)
    end
    return ranking, total_gains
end

# %%
k = 1000

# %%
@time ranking0, gains0 = fit(X_digits, k; calculate_gains! = calculate_gains!);

# %%
@time ranking0, gains0 = fit(X_digits, k; calculate_gains! = calculate_gains!);

# %%
tic = time()
ranking0, gains0 = fit(X_digits, k; calculate_gains! = calculate_gains!)
toc0 = time() - tic
toc0

# %%
@time begin
idxs = sample(axes(X_digits, 2), k; replace=false)
X_subset = X_digits[:, idxs]
gains1 = cumsum(X_subset; dims=2)
gains1 = vec(sum(sqrt, gains1; dims=1))
end;

# %%
@time begin
idxs = sample(axes(X_digits, 2), k; replace=false)
X_subset = X_digits[:, idxs]
gains1 = cumsum(X_subset; dims=2)
gains1 = vec(sum(sqrt, gains1; dims=1))
end;

# %%
tic = time()
idxs = sample(axes(X_digits, 2), k; replace=false)
X_subset = X_digits[:, idxs]
gains1 = cumsum(X_subset; dims=2)
gains1 = vec(sum(sqrt, gains1; dims=1))
toc1 = time() - tic
toc1

# %%
plt.figure(figsize=(9, 4.5))

plt.subplot(121)
plt.plot(cumsum(gains0), label="Naive")
plt.plot(gains1, label="Random")
plt.ylabel("F(S)")
plt.xlabel("Subset Size")
plt.legend()
plt.grid(lw=0.3)

plt.subplot(122)
plt.bar(1:2, [toc0,  toc1])
plt.ylabel("Time (s)")
plt.xticks(1:2, ["Naive", "Random"], rotation=90)
plt.grid(lw=0.3)
plt.title("Julia 1.6.2")

plt.tight_layout()
#plt.show()

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
try
    using FLoopsTests
    true
catch
    false
end || begin
    push!(LOAD_PATH, @__DIR__)
    using FLoopsTests
end


text:
function rarity_code(x)
    x == "COMMON" && return 1
    x == "RARE" && return 2
    x == "EPIC" && return 3
    x == "LEGENDARY" && return 4
    x == "FREE" && return 0
    error("Unkown rarity $x")
end

const card_db = eval(Meta.parse(read(joinpath(@__DIR__,"..", "data", "card_db.txt"), String)))

max_stars_repo_name: asieev/HearthstoneSim

max_issues_repo_name: asieev/HearthstoneSim

max_forks_repo_name: asieev/HearthstoneSim


text:
#=
run tests using Pardiso sparse linear system solver caches
requires that Pardiso.jl be installed and built successfully; Requires.jl handles this optional dependency
=#

using Test
import Pardiso
import Hypatia.Solvers

include(joinpath(@__DIR__, "nativeinstances.jl"))
include(joinpath(@__DIR__, "nativesets.jl"))

options = (verbose = true,)

@testset "Pardiso cache tests" begin

@testset "cache setup: $cache_type" for cache_type in [Hypatia.PardisoSymCache,
    Hypatia.PardisoNonSymCache]
    cache = cache_type()
    @test !cache.analyzed
    @test_throws Exception cache_type{Float32}()
end

T = Float64
@testset "NaiveSparse tests: $inst_name" for inst_name in inst_cones_many
    inst_function = eval(Symbol(inst_name))
    inst_function(T, solver = Solvers.Solver{T}(syssolver =
        Solvers.NaiveSparseSystemSolver{T}(fact_cache =
        Hypatia.PardisoNonSymCache{T}()); options...))
end
@testset "NaiveElimSparse tests: $inst_name" for inst_name in inst_cones_many
    inst_function = eval(Symbol(inst_name))
    inst_function(T, solver = Solvers.Solver{T}(syssolver =
        Solvers.NaiveElimSparseSystemSolver{T}(fact_cache =
        Hypatia.PardisoNonSymCache{T}()); options...))
end
@testset "SymIndefSparse tests: $inst_name" for inst_name in inst_cones_many
    inst_function = eval(Symbol(inst_name))
    inst_function(T, solver = Solvers.Solver{T}(syssolver =
        Solvers.SymIndefSparseSystemSolver{T}(fact_cache =
        Hypatia.PardisoSymCache{T}()); options...))
end

end
;


text:
using Documenter, Mag2Dpoly #, PyPlot

makedocs(sitename="Mag2Dpoly.jl",
         modules = [Mag2Dpoly],
         authors = "Andrea Zunino, Alessandro Ghirotto",
         format = Documenter.HTML(prettyurls=get(ENV,"CI",nothing)=="true"),
         pages = [
             "Home" => "index.md",
         ]
         )

deploydocs(
    repo="github.com/inverseproblem/Mag2Dpoly.jl.git",
    devbranch = "main",
    branch = "gh-pages"
)


text:
module Catlab

include("core/Core.jl")
include("theories/Theories.jl")
include("categorical_algebra/CategoricalAlgebra.jl")
include("wiring_diagrams/WiringDiagrams.jl")
include("graphics/Graphics.jl")
include("programs/Programs.jl")
include("linear_algebra/LinearAlgebra.jl")

end


text:
include("../test_fdm.jl")
nparts = (2,2,2)
prun(test_fdm,mpi,nparts)


text:
using Talkon
using Documenter

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

makedocs(;
    modules=[Talkon],
    authors="TBD",
    repo="https://github.com/tsurkys/Talkon.jl/blob/{commit}{path}#{line}",
    sitename="Talkon.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://tsurkys.github.io/Talkon.jl",
        siteurl="https://github.com/tsurkys/Talkon.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/tsurkys/Talkon.jl",
)


text:
##############################################################################
#
# REFERENCES: "Statistical Distributions"
#
##############################################################################

immutable Erlang <: ContinuousUnivariateDistribution
    shape::Int
    scale::Float64
    nested_gamma::Gamma
end

Erlang(scale::Real) = Erlang(1, 1.0)
Erlang() = Erlang(1, 1.0)

cdf(d::Erlang, x::Real) = cdf(d.nested_gamma, x)

entropy(d::Erlang) = entropy(d.nested_gamma)

insupport(d::Erlang, x::Number) = isreal(x) && isfinite(x) && 0.0 <= x

kurtosis(d::Erlang) = kurtosis(d.nested_gamma)

mean(d::Erlang) = d.shape * d.scale

median(d::Erlang) = median(d.nested_gamma)

mgf(d::Erlang, t::Real) = mgf(d.nested_gamma, t)
cf(d::Erlang, t::Real) = cf(d.nested_gamma, t)

modes(d::Erlang) = modes(d.nested_gamma)

function pdf(d::Erlang, x::Real)
    b, c = d.scale, d.shape
    return ((x / b)^(c - 1.0) * exp(-x / b)) / (b * gamma(c))
end

quantile(d::Erlang, p::Real) = quantile(d.nested_gamma, p)

function rand(d::Erlang)
    b, c = d.scale, d.shape
    z = 0.0
    for i in 1:c
        z *= rand()
    end
    return -b * log(z)
end

skewness(d::Erlang) = skewness(d.nested_gamma)

var(d::Erlang) = d.scale^2 * d.shape


text:
@deprecate einsumexp(args...) loop_einsum(args...)
@deprecate einsumexp!(args...) loop_einsum!(args...)

@deprecate einsumopt einsum

# deprecated manually
function Iterators.flatten(c::Union{NestedEinsum,EinCode})
    @warn "`Iterators.flatten(eincode)` has been deprecated, use `OMEinsum.flatten` instead."
    flatten(c)
end


text:
using RAFF

#
# This example is the "Basic Usage" example in the documentation
#

# Matrix with data
A=[-2.0  5.0;
   -1.5  3.25;
   -1.0  2.0 ;
   -0.5  1.25;
    0.0  1.0 ;
    0.5  2.55;
    1.0  2.0 ;
    1.5  3.25;
    2.0  5.0 ;]

# Define the model to fit data
model(x, θ) = θ[1] * x[1]^2 + θ[2]

# Number of parameters in the model
n = 2

# Run RAFF
output = raff(model, A, n)

# Print output
println(output)


text:
type Min <: ContinuousUnivariateStreamStat
    m::Float64
    n::Int
end

Min() = Min(Inf, 0)

function update!(stat::Min, x::Real)
    stat.n += 1
    stat.m = min(stat.m, x)
    return
end

Base.minimum(stat::Min) = stat.m

state(stat::Min) = Base.minimum(stat)

nobs(stat::Min) = stat.n

Base.copy(stat::Min) = Min(stat.m, stat.n)

function Base.merge(a::Min, b::Min)
    m = min(a.m, b.m)
    n = a.n + b.n
    return Min(m, n)
end

function Base.empty!(stat::Min)
    stat.m = Inf
    stat.n = 0
    return
end

function Base.show(io::IO, stat::Min)
    m = minimum(stat)
    n = nobs(stat)
    @printf(io, "Online Min\n")
    @printf(io, " * Min:  %f\n", m)
    @printf(io, " * N:    %d\n", n)
    return
end


text:
#=
  Demos for gradient descent glm algos
=#

using glmSolverjl
using Random: seed!

seed!(0);
gammaLogX, gammaLogY = simulateData(Float64, GammaDistribution(), LogLink(), 30, 100_000);
gammaLogBlockX = matrixToBlock(gammaLogX, 100);
gammaLogBlockY = matrixToBlock(gammaLogY, 100);

function gdDataDemo(niter::Int64 = 20, learningRate::Float64 = 2E-8)

  #= Number of parameters =#
  p = size(gammaLogX)[2]

  gammaModel = glm(Block1DParallel(), gammaLogBlockX,
        gammaLogBlockY, GammaDistribution(), LogLink(),
        #= solver =# GESVSolver(), inverse = GETRIInverse());
  println("Full GLM Solve\n", gammaModel.coefficients)

  println("The outputs for all these models should be the same.");

  gammaModel = glm(RegularData(), gammaLogX,
        gammaLogY, GammaDistribution(), LogLink(),
        #= solver =# GradientDescentSolver(learningRate), inverse = GETRIInverse(),
        control = Control{Float64}(maxit = niter),
        calculateCovariance = true,
        doStepControl = false);
  println("Gradient Descent With Regular Data\n", gammaModel.coefficients)
  #===================================================================#
  #= Gradient Descent Block Model =#
  gammaModel = glm(Block1D(), gammaLogBlockX,
        gammaLogBlockY, GammaDistribution(), LogLink(),
        #= solver =# GradientDescentSolver(learningRate),
        inverse = GETRIInverse(), control = Control{Float64}(maxit = niter),
        calculateCovariance = true,
        doStepControl = false);
  println("Gradient Descent With Block Data \n", gammaModel.coefficients)

  #= Gradient Descent Block1DParallel Model =#
  gammaModel = glm(Block1DParallel(), gammaLogBlockX,
        gammaLogBlockY, GammaDistribution(), LogLink(),
        #= solver =# GradientDescentSolver(learningRate),
        inverse = GETRIInverse(), control = Control{Float64}(maxit = niter),
        calculateCovariance = true,
        doStepControl = false);
  println("Gradient Descent With Parallel Block Data \n", gammaModel.coefficients)

  return
end

function momentumDataDemo(niter::Int64 = 20, learningRate::Float64 = 2E-8, momentum::Float64 = 0.7)

  #= Number of parameters =#
  p = size(gammaLogX)[2]

  gammaModel = glm(Block1DParallel(), gammaLogBlockX,
      gammaLogBlockY, GammaDistribution(), LogLink(),
      #= solver =# GESVSolver(), inverse = GETRIInverse());
  println("Full GLM Solve\n", gammaModel.coefficients)

  println("The outputs for all these models should be the same.");

  gammaModel = glm(RegularData(), gammaLogX,
      gammaLogY, GammaDistribution(), LogLink(),
      #= solver =# MomentumSolver(learningRate, momentum, p),
      inverse = GETRIInverse(),
      control = Control{Float64}(maxit = niter),
      calculateCovariance = true,
      doStepControl = false);
  println("Gradient Descent With Regular Data\n", gammaModel.coefficients)
  #===================================================================#
  #= Gradient Descent Block Model =#
  gammaModel = glm(Block1D(), gammaLogBlockX,
      gammaLogBlockY, GammaDistribution(), LogLink(),
      #= solver =# MomentumSolver(learningRate, momentum, p),
      inverse = GETRIInverse(),
      control = Control{Float64}(maxit = niter),
      calculateCovariance = true,
      doStepControl = false);
  println("Gradient Descent With Block Data \n", gammaModel.coefficients)

  #= Gradient Descent Block1DParallel Model =#
  gammaModel = glm(Block1DParallel(), gammaLogBlockX,
      gammaLogBlockY, GammaDistribution(), LogLink(),
      #= solver =# MomentumSolver(learningRate, momentum, p),
      inverse = GETRIInverse(), control = Control{Float64}(maxit = niter),
      calculateCovariance = true,
      doStepControl = false);
  println("Gradient Descent With Parallel Block Data \n", gammaModel.coefficients)

  return
end

function nesterovDataDemo(niter::Int64 = 30, learningRate::Float64 = 2E-8, momentum::Float64 = 0.70)

  #= Number of parameters =#
  p = size(gammaLogX)[2]

  gammaModel = glm(Block1DParallel(), gammaLogBlockX,
      gammaLogBlockY, GammaDistribution(), LogLink(),
      #= solver =# GESVSolver(), inverse = GETRIInverse());
  println("Full GLM Solve\n", gammaModel.coefficients)

  println("The outputs for all these models should be the same.");

  gammaModel = glm(RegularData(), gammaLogX,
      gammaLogY, GammaDistribution(), LogLink(),
      #= solver =# NesterovSolver(learningRate, momentum, p),
      inverse = GETRIInverse(),
      control = Control{Float64}(maxit = niter),
      calculateCovariance = true,
      doStepControl = false);
  println("Gradient Descent With Regular Data\n", gammaModel.coefficients)
  #===================================================================#
  #= Gradient Descent Block Model =#
  gammaModel = glm(Block1D(), gammaLogBlockX,
      gammaLogBlockY, GammaDistribution(), LogLink(),
      #= solver =# NesterovSolver(learningRate, momentum, p),
      inverse = GETRIInverse(),
      control = Control{Float64}(maxit = niter),
      calculateCovariance = true,
      doStepControl = false);
  println("Gradient Descent With Block Data \n", gammaModel.coefficients)

  #= Gradient Descent Block1DParallel Model =#
  gammaModel = glm(Block1DParallel(), gammaLogBlockX,
      gammaLogBlockY, GammaDistribution(), LogLink(),
      #= solver =# NesterovSolver(learningRate, momentum, p),
      inverse = GETRIInverse(), control = Control{Float64}(maxit = niter),
      calculateCovariance = true,
      doStepControl = false);
  println("Gradient Descent With Parallel Block Data \n", gammaModel.coefficients)

  return
end

function adagradDataDemo(niter::Int64 = 10, learningRate::Float64 = 1E-10, epsilon::Float64 = 1E-8)

  #= Number of parameters =#
  p = size(gammaLogX)[2]

  gammaModel = glm(Block1DParallel(), gammaLogBlockX,
      gammaLogBlockY, GammaDistribution(), LogLink(),
      #= solver =# GESVSolver(), inverse = GETRIInverse());
  println("Full GLM Solve\n", gammaModel.coefficients)

  println("The outputs for all these models should be the same.");

  gammaModel = glm(RegularData(), gammaLogX,
      gammaLogY, GammaDistribution(), LogLink(),
      #= solver =# AdagradSolver(learningRate, p, epsilon),
      inverse = GETRIInverse(),
      control = Control{Float64}(maxit = niter),
      calculateCovariance = true,
      doStepControl = false);
  println("Gradient Descent With Regular Data\n", gammaModel.coefficients)
  #===================================================================#
  #= Gradient Descent Block Model =#
  gammaModel = glm(Block1D(), gammaLogBlockX,
      gammaLogBlockY, GammaDistribution(), LogLink(),
      #= solver =# AdagradSolver(learningRate, p, epsilon),
      inverse = GETRIInverse(),
      control = Control{Float64}(maxit = niter),
      calculateCovariance = true,
      doStepControl = false);
  println("Gradient Descent With Block Data \n", gammaModel.coefficients)

  #= Gradient Descent Block1DParallel Model =#
  gammaModel = glm(Block1DParallel(), gammaLogBlockX,
      gammaLogBlockY, GammaDistribution(), LogLink(),
      #= solver =# AdagradSolver(learningRate, p, epsilon),
      inverse = GETRIInverse(), control = Control{Float64}(maxit = niter),
      calculateCovariance = true,
      doStepControl = false);
  println("Gradient Descent With Parallel Block Data \n", gammaModel.coefficients)

  return
end

#= Poor convergence performance for this one =#
function adadeltaDataDemo(niter::Int64 = 10, momentum::Float64 = 0.9, epsilon::Float64 = 1E-8)

  #= Number of parameters =#
  p = size(gammaLogX)[2]

  gammaModel = glm(Block1DParallel(), gammaLogBlockX,
      gammaLogBlockY, GammaDistribution(), LogLink(),
      #= solver =# GESVSolver(), inverse = GETRIInverse());
  println("Full GLM Solve\n", gammaModel.coefficients)

  println("The outputs for all these models should be the same.");

  gammaModel = glm(RegularData(), gammaLogX,
      gammaLogY, GammaDistribution(), LogLink(),
      #= solver =# AdadeltaSolver(momentum, p, epsilon),
      inverse = GETRIInverse(),
      control = Control{Float64}(maxit = niter),
      calculateCovariance = true,
      doStepControl = false);
  println("Gradient Descent With Regular Data\n", gammaModel.coefficients)
  #===================================================================#
  #= Gradient Descent Block Model =#
  gammaModel = glm(Block1D(), gammaLogBlockX,
      gammaLogBlockY, GammaDistribution(), LogLink(),
      #= solver =# AdadeltaSolver(momentum, p, epsilon),
      inverse = GETRIInverse(),
      control = Control{Float64}(maxit = niter),
      calculateCovariance = true,
      doStepControl = false);
  println("Gradient Descent With Block Data \n", gammaModel.coefficients)

  #= Gradient Descent Block1DParallel Model =#
  gammaModel = glm(Block1DParallel(), gammaLogBlockX,
      gammaLogBlockY, GammaDistribution(), LogLink(),
      #= solver =# AdadeltaSolver(momentum, p, epsilon),
      inverse = GETRIInverse(), control = Control{Float64}(maxit = niter),
      calculateCovariance = true,
      doStepControl = false);
  println("Gradient Descent With Parallel Block Data \n", gammaModel.coefficients)

  return
end

function RMSPropDataDemo(niter::Int64 = 10, learningRate::Float64 = 1E-2, epsilon::Float64 = 1E-8)

  #= Number of parameters =#
  p = size(gammaLogX)[2]

  gammaModel = glm(Block1DParallel(), gammaLogBlockX,
        gammaLogBlockY, GammaDistribution(), LogLink(),
        #= solver =# GESVSolver(), inverse = GETRIInverse());
  println("Full GLM Solve\n", gammaModel.coefficients)

  println("The outputs for all these models should be the same.");

  gammaModel = glm(RegularData(), gammaLogX,
        gammaLogY, GammaDistribution(), LogLink(),
        #= solver =# RMSpropSolver(learningRate, p, epsilon),
        inverse = GETRIInverse(),
        control = Control{Float64}(maxit = niter),
        calculateCovariance = true,
        doStepControl = false);
  println("Gradient Descent With Regular Data\n", gammaModel.coefficients)
  #===================================================================#
  #= Gradient Descent Block Model =#
  gammaModel = glm(Block1D(), gammaLogBlockX,
        gammaLogBlockY, GammaDistribution(), LogLink(),
        #= solver =# RMSpropSolver(learningRate, p, epsilon),
        inverse = GETRIInverse(), control = Control{Float64}(maxit = niter),
        calculateCovariance = true,
        doStepControl = false);
  println("Gradient Descent With Block Data \n", gammaModel.coefficients)

  #= Gradient Descent Block1DParallel Model =#
  gammaModel = glm(Block1DParallel(), gammaLogBlockX,
        gammaLogBlockY, GammaDistribution(), LogLink(),
        #= solver =# RMSpropSolver(learningRate, p, epsilon),
        inverse = GETRIInverse(),
        control = Control{Float64}(maxit = niter),
        calculateCovariance = true,
        doStepControl = false);
  println("Gradient Descent With Parallel Block Data \n", gammaModel.coefficients)

  return
end

function adamDataDemo(niter::Int64 = 50, learningRate::Float64 = 1E-3, b1::Float64 = 0.9, b2::Float64 = 0.999, epsilon::Float64 = 1E-8)

  #= Number of parameters =#
  p = size(gammaLogX)[2]

  gammaModel = glm(Block1DParallel(), gammaLogBlockX,
        gammaLogBlockY, GammaDistribution(), LogLink(),
        #= solver =# GESVSolver(), inverse = GETRIInverse());
  println("Full GLM Solve\n", gammaModel.coefficients)

  println("The outputs for all these models should be the same.");

  gammaModel = glm(RegularData(), gammaLogX,
        gammaLogY, GammaDistribution(), LogLink(),
        #= solver =# AdamSolver(learningRate, b1, b2, p, epsilon),
        inverse = GETRIInverse(),
        control = Control{Float64}(maxit = niter),
        calculateCovariance = true,
        doStepControl = false);
  println("Gradient Descent With Regular Data\n", gammaModel.coefficients)
  #===================================================================#
  #= Gradient Descent Block Model =#
  gammaModel = glm(Block1D(), gammaLogBlockX,
        gammaLogBlockY, GammaDistribution(), LogLink(),
        #= solver =# AdamSolver(learningRate, b1, b2, p, epsilon),
        inverse = GETRIInverse(), control = Control{Float64}(maxit = niter),
        calculateCovariance = true,
        doStepControl = false);
  println("Gradient Descent With Block Data \n", gammaModel.coefficients)

  #= Gradient Descent Block1DParallel Model =#
  gammaModel = glm(Block1DParallel(), gammaLogBlockX,
        gammaLogBlockY, GammaDistribution(), LogLink(),
        #= solver =# AdamSolver(learningRate, b1, b2, p, epsilon),
        inverse = GETRIInverse(),
        control = Control{Float64}(maxit = niter),
        calculateCovariance = true,
        doStepControl = false);
  println("Gradient Descent With Parallel Block Data \n", gammaModel.coefficients)

  return
end

function adaMaxDataDemo(niter::Int64 = 10, learningRate::Float64 = 1E-1, b1::Float64 = 0.9, b2::Float64 = 0.999)

  #= Number of parameters =#
  p = size(gammaLogX)[2]

  gammaModel = glm(Block1DParallel(), gammaLogBlockX,
        gammaLogBlockY, GammaDistribution(), LogLink(),
        #= solver =# GESVSolver(), inverse = GETRIInverse());
  println("Full GLM Solve\n", gammaModel.coefficients)

  println("The outputs for all these models should be the same.");

  gammaModel = glm(RegularData(), gammaLogX,
        gammaLogY, GammaDistribution(), LogLink(),
        #= solver =# AdaMaxSolver(learningRate, b1, b2, p),
        inverse = GETRIInverse(),
        control = Control{Float64}(maxit = niter),
        calculateCovariance = true,
        doStepControl = false);
  println("Gradient Descent With Regular Data\n", gammaModel.coefficients)
  #===================================================================#
  #= Gradient Descent Block Model =#
  gammaModel = glm(Block1D(), gammaLogBlockX,
        gammaLogBlockY, GammaDistribution(), LogLink(),
        #= solver =# AdaMaxSolver(learningRate, b1, b2, p),
        inverse = GETRIInverse(),
        control = Control{Float64}(maxit = niter),
        calculateCovariance = true,
        doStepControl = false);
  println("Gradient Descent With Block Data \n", gammaModel.coefficients)

  #= Gradient Descent Block1DParallel Model =#
  gammaModel = glm(Block1DParallel(), gammaLogBlockX,
        gammaLogBlockY, GammaDistribution(), LogLink(),
        #= solver =# AdaMaxSolver(learningRate, b1, b2, p),
        inverse = GETRIInverse(),
        control = Control{Float64}(maxit = niter),
        calculateCovariance = true,
        doStepControl = false);
  println("Gradient Descent With Parallel Block Data \n", gammaModel.coefficients)

  return
end

function nadamDataDemo(niter::Int64 = 30, learningRate::Float64 = 3E-3, b1::Float64 = 0.9, b2::Float64 = 0.999, epsilon::Float64 = 1E-8)

  #= Number of parameters =#
  p = size(gammaLogX)[2]

  gammaModel = glm(Block1DParallel(), gammaLogBlockX,
        gammaLogBlockY, GammaDistribution(), LogLink(),
        #= solver =# GESVSolver(), inverse = GETRIInverse());
  println("Full GLM Solve\n", gammaModel.coefficients)

  println("The outputs for all these models should be the same.");

  gammaModel = glm(RegularData(), gammaLogX,
        gammaLogY, GammaDistribution(), LogLink(),
        #= solver =# NAdamSolver(learningRate, b1, b2, p, epsilon),
        inverse = GETRIInverse(),
        control = Control{Float64}(maxit = niter),
        calculateCovariance = true,
        doStepControl = false);
  println("Gradient Descent With Regular Data\n", gammaModel.coefficients)
  #===================================================================#
  #= Gradient Descent Block Model =#
  gammaModel = glm(Block1D(), gammaLogBlockX,
        gammaLogBlockY, GammaDistribution(), LogLink(),
        #= solver =# NAdamSolver(learningRate, b1, b2, p, epsilon),
        inverse = GETRIInverse(), control = Control{Float64}(maxit = niter),
        calculateCovariance = true,
        doStepControl = false);
  println("Gradient Descent With Block Data \n", gammaModel.coefficients)

  #= Gradient Descent Block1DParallel Model =#
  gammaModel = glm(Block1DParallel(), gammaLogBlockX,
        gammaLogBlockY, GammaDistribution(), LogLink(),
        #= solver =# NAdamSolver(learningRate, b1, b2, p, epsilon),
        inverse = GETRIInverse(),
        control = Control{Float64}(maxit = niter),
        calculateCovariance = true,
        doStepControl = false);
  println("Gradient Descent With Parallel Block Data \n", gammaModel.coefficients)

  return
end

function amsGradDataDemo(niter::Int64 = 50, learningRate::Float64 = 9E-3, b1::Float64 = 0.9, b2::Float64 = 0.999, epsilon::Float64 = 1E-8)

  #= Number of parameters =#
  p = size(gammaLogX)[2]

  gammaModel = glm(Block1DParallel(), gammaLogBlockX,
        gammaLogBlockY, GammaDistribution(), LogLink(),
        #= solver =# GESVSolver(), inverse = GETRIInverse());
  println("Full GLM Solve\n", gammaModel.coefficients)

  println("The outputs for all these models should be the same.");

  gammaModel = glm(RegularData(), gammaLogX,
        gammaLogY, GammaDistribution(), LogLink(),
        #= solver =# AMSGradSolver(learningRate, b1, b2, p, epsilon),
        inverse = GETRIInverse(),
        control = Control{Float64}(maxit = niter),
        calculateCovariance = true,
        doStepControl = false);
  println("Gradient Descent With Regular Data\n", gammaModel.coefficients)
  #===================================================================#
  #= Gradient Descent Block Model =#
  gammaModel = glm(Block1D(), gammaLogBlockX,
        gammaLogBlockY, GammaDistribution(), LogLink(),
        #= solver =# AMSGradSolver(learningRate, b1, b2, p, epsilon),
        inverse = GETRIInverse(), control = Control{Float64}(maxit = niter),
        calculateCovariance = true,
        doStepControl = false);
  println("Gradient Descent With Block Data \n", gammaModel.coefficients)

  #= Gradient Descent Block1DParallel Model =#
  gammaModel = glm(Block1DParallel(), gammaLogBlockX,
        gammaLogBlockY, GammaDistribution(), LogLink(),
        #= solver =# AMSGradSolver(learningRate, b1, b2, p, epsilon),
        inverse = GETRIInverse(),
        control = Control{Float64}(maxit = niter),
        calculateCovariance = true,
        doStepControl = false);
  println("Gradient Descent With Parallel Block Data \n", gammaModel.coefficients)

  return
end

#===================================================================#

# Simple Gradient Descent
gdDataDemo(30, 2E-8)
# Gradient Descent with Momentum
momentumDataDemo(20, 2E-8, 0.7)
# Nesterov
nesterovDataDemo(30, 2E-8, 0.7)
# Adagrad
adagradDataDemo(10, 1E-10, 1E-8)
# Adadelta (Poor performance)
adadeltaDataDemo(10, 0.9, 1E-8)
# RMSProp
RMSPropDataDemo(10, 1E-2, 1E-8)
# Adam
adamDataDemo(50, 1E-3, 0.9, 0.999, 1E-8)
# AdaMax
adaMaxDataDemo(10, 1E-1, 0.9, 0.999)
# NAdam
nadamDataDemo(30, 3E-3, 0.9, 0.999, 1E-8)
# AMSGrad
amsGradDataDemo(40, 5E-4, 0.9, 0.999, 1E-8)

max_stars_repo_name: dataPulverizer/glmSolver

max_issues_repo_name: dataPulverizer/glmSolver

max_forks_repo_name: dataPulverizer/glmSolver


text:
using StochasticDiffEq, Test, Random
using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems; importsdeproblems()
using DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_linear, prob_sde_2Dlinear

Random.seed!(100)

for prob in (prob_sde_linear, prob_sde_2Dlinear)
  prob2 = remake(prob; tspan = (1.0, 0.0))

  solEM = solve(prob2, EM(), dt = -1//2^(4), tstops = [0.33])
  @test solEM.t[15] > solEM.t[16]

  solSRIW1 = solve(prob2, SRIW1(), tstops = [0.33])
  @test solSRIW1.t[15] > solSRIW1.t[16]
end


text:
const r = -3:6
const v = collect(r)
const v2 = 1.:10.
const m = reshape(1:100, 10, 10)
const m2 = reshape(1:35, 5, 7)

const args = [
    :(),                                # nil
    :(r,),                              # vec (range)
    :(v,),                              # vec (array)
    :(v',),                             # vec (transposed array)
    :(m,),                              # mat
    :(complex.(v,v),),                   # complex vec
    :(complex.(m,m),),                   # complex mat
    :(v,v2),                            # vec/vec
    :(m,m),                             # mat/mat (square)
    :(m2,m2),                           # mat/mat (non-square)
    :(v,m),                             # vec/mat
    :(2:6,m2),                          # vec/mat (cols)
    :(2:8,m2),                          # vec/mat (rows)
    :(m2,2:6),                          # mat/vec (cols)
    :(m2,2:8),                          # mat/vec (rows)
]

for i = 1:length(args)
    f = Symbol(@sprintf("plot%03d",i))
    body = copy(args[i])
    body.head = :call
    pushfirst!(body.args, :plot)
    push!(body.args, Expr(:kw, :title, string(body)))
    eval(Expr(:toplevel, Expr(:export, f), :(($f)() = $body)))
end


text:
struct QuadFloat64 <: Real
    a1::Float64
    a2::Float64
    a3::Float64
    a4::Float64
end


text:
module DynamicIterators
using Trajectories

export collectfrom, DynamicIterator, dyniterate

export Key, NextKey # keywords

export Evolution, evolve, timelift_evolve, from # evolution
export attach, Bind, bind, bindafter, Evolve, TimeLift, mix, synchronize, mixture # combinators

export trace, endtime, lastiterate # trajectories

export control, timed # control

# random
export WhiteNoise, Randn, InhomogeneousPoisson,
    MetropolisHastings #,Sampled

using Random, Base.Iterators

using Base.Iterators
using Base: SizeUnknown, HasEltype
import Base: iterate, IteratorSize, @propagate_inbounds, IsInfinite, eltype, IteratorEltype,
    rand, EltypeUnknown, HasEltype, bind

# keyword arguments:
# jump, key, t

"""
    DynamicIterator

`DynamicIterator`s which extend the iterator protocol by
keywords for the `iterate` function.

"""
abstract type DynamicIterator
end

abstract type Message
end

macro returnnothing(exp)
    quote let ϕ = $(esc(exp)); if ϕ === nothing; return nothing; end; ϕ end end
end
function evolve
end
function dyniterate
end

include("messages.jl")

Sample(x) = Sample(x, Random.GLOBAL_RNG)

Base.getindex(M::Message1) = getfield(M, 1)

iterate(M::Message2) = getfield(M, 1), 1
iterate(M::Message2, Any) = getfield(M, 2), nothing
iterate(M::Message2, ::Nothing) = nothing

iterate(M::Message1) = getfield(M, 1), nothing
iterate(M::Message1, Any) = nothing

export Message, Message1, Message2, Start, Control, Value, Steps, Sample,
    NewKey, Key, NextKey, NextKeys, BindOnce

dub(x) = x === nothing ? nothing : (x, x)
dedub(x) = x === nothing ? nothing : x[1]
dubwith(x, c) = x === nothing ? nothing : (x, c(x))

# keyword functions shouldn't shadow non-keyword functions
# when keywords are absent
iteratefallback(iter, state) = iterate(iter, state)
iteratefallback(iter, ::Nothing) = iterate(iter)
iteratefallback(iter::DynamicIterator, state) = throw(ArgumentError("Trying to use `iteration` fallback for DynamicIterator"))
iteratefallback(iter::DynamicIterator, ::Nothing) = throw(ArgumentError("Trying to use `iteration` fallback for DynamicIterator"))
dyniterate(iter, state) = iteratefallback(iter, state)

iterate(iter::DynamicIterator) = dyniterate(iter, nothing)
iterate(iter::DynamicIterator, state) = dyniterate(iter, state)

IteratorSize(::DynamicIterator) = SizeUnknown()

include("evolution.jl")
include("time.jl")
include("combinators.jl")
include("trajectories.jl")
include("random.jl")
include("control.jl")
# Examples

end


text:
@testset "lasing" begin

# Create a system.
Ngrid = [3,3,3]
N = 3prod(Ngrid)
ind_in = 5:7
εc = ones(ComplexF64, N)
εc[ind_in] .= 12+0.1im
m = 1
ωₘ = 100.0
ωₐ = 1.0
γ˔ = 1.0
D₀ = fill(0.01, N)
D₀[ind_in] .= 1.0
gp = GainProfile(gen_gain_2lv(ωₐ,γ˔), D₀)

# Create a solution.
M = 3
ω = rand(M)
ω[m] = ωₘ
a² = rand(M)
abs2gain = SALTBase.gen_abs2gain_2lv(ωₐ,γ˔)
abs2gain′ = SALTBase.gen_abs2gain′_2lv(ωₐ,γ˔)
Ψ = randn(ComplexF64,N,M)
ψ = [Ψ[:,j] for j = 1:M]
iarray = rand(1:N, M)
vtemp = similar(ψ[1])
sol = LasingSol(ω, a², ψ, [abs2.(ψ[m]) for m = 1:M], iarray, vtemp)  # ω = rand(M); this is fine because init_reduced_var does not use ω
sol.active = fill(true, M)
sol.m_active = collect(1:M)

# Check reduced variables.
hb = Vector{Float64}(undef, N)
hole_burning!(hb, abs2gain, ω, a², [abs2.(ψ[m]) for m = 1:M])
D = D₀ ./ hb
D′ = -D₀ ./ abs2.(hb)
rvar = SALTBase.LasingReducedVar(N, 1)  # K = 1
SALTBase.init_reduced_var!(rvar, sol, gp)
@testset "reduced variables" begin
    @test rvar.D[1] ≈ D
    @test rvar.D′[1] ≈ D′
end

mvar = SALTBase.LasingModalVar(DefaultMXS(), N)
SALTBase.init_modal_var!(mvar, m, sol, rvar.D, gp, εc)

γ = gen_gain_2lv(ωₐ, γ˔)(ωₘ)
γ′ = SALTBase.gen_gain′_2lv(ωₐ, γ˔)(ωₘ)
ε = εc + γ*D
@testset "modal variables" begin
    @test mvar.∂f∂ω ≈ (2ωₘ*ε + ωₘ^2*γ′*D) .* ψ[m]
end

end  # @testset "lasing"

# @test mvar.rowA⁻¹ᵢₐ ≈ transpose(transpose(eᵢₐ)/A)


text:
#=  A cell, single conceptual entity on the terminal screen. The terminal screen
 * is basically a 2d array of cells. It has the following fields:
 *  - 'ch' is a unicode character
 *  - 'fg' foreground color and attributes
 *  - 'bg' background color and attributes
  =#
#=struct tb_cell {
  uint32_t ch;
  uint16_t fg;
  uint16_t bg;
};=#
export tb_cell
mutable struct tb_cell
  ch::Cuint
  fg::Cushort
  bg::Cushort
  tb_cell() = new()
  tb_cell(ch, fg=0x0001, bg = 0x0000) = new(Cuint(ch), fg, bg)
end

#=  This struct represents a termbox event. The 'mod', 'key' and 'ch' fields are
 * valid if 'type' is TB_EVENT_KEY. The 'w' and 'h' fields are valid if 'type'
 * is TB_EVENT_RESIZE. The 'x' and 'y' fields are valid if 'type' is
 * TB_EVENT_MOUSE.
  =#
#=struct tb_event {
  uint8_t type;
  uint8_t mod;
  uint16_t key;
  uint32_t ch;
  int32_t w;
  int32_t h;
  int32_t x;
  int32_t y;
  };=#
export tb_event
mutable struct tb_event
  etype::Cuchar
  mod::Cuchar
  key::Cushort
  ch::Cuint
  w::Cint
  h::Cint
  x::Cint
  y::Cint
  tb_event() = new()
end


text:
abstract type ReturnLevelModel{T<:fittedEVA} end

"""
    ReturnLevel

ReturnLevel type constructed by the function [`returnlevel`](@ref).
"""
struct ReturnLevel{T<:fittedEVA}
      model::ReturnLevelModel{T}
      returnperiod::Real
      value::Array{<:Real}
end

struct BlockMaximaModel{T<:fittedEVA{BlockMaxima}} <: ReturnLevelModel{T}
  fm::T
end

struct PeakOverThreshold{T<:fittedEVA{ThresholdExceedance}} <: ReturnLevelModel{T}
  fm::T
  threshold::Real
  nobservation::Int
  nobsperblock::Int
end

"""
    Base.show(io::IO, obj::ReturnLevel)

Override of the show function for the objects of type ReturnLevel.
"""
function Base.show(io::IO, obj::ReturnLevel)

    println(io, "ReturnLevel")
    println(io, "returnperiod :\t", obj.returnperiod)
    println(io, "value :\t\t", typeof(obj.value), "[", length(obj.value), "]")

end


text:
# SPDX-License-Identifier: MIT

using ArgParse
using HDF5

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "features"
            required = true
            help = "HDF5 archive containing the features"
        "out"
            required = true
            help = "output HDF5 file containing the statistics"
    end
    s.description = """
    Extract statistics from the a set of utterances. The list
    of utterances to use is read from the standard input.
    """
    parse_args(s)
end

function main(args)
    h5open(args["features"], "r") do f

        # Initialize the accumulators with the 1st utterance
        data = read(f[readline(stdin)])
        D, N = size(data)
        x = dropdims(sum(data, dims = 2), dims = 2)
        xxᵀ = sum([c*c' for c in eachcol(data)])

        uttcount = 0
        for line in eachline(stdin)
            uttcount += 1
            @debug "reading $(uttcount)th utterance: $line"
            data = read(f[line])
            N += size(data, 2)
            x += dropdims(sum(data, dims = 2), dims = 2)
            xxᵀ += sum([c*c' for c in eachcol(data)])
        end
        μ = x ./ N
        Σ = (xxᵀ .- μ*μ') ./ N

        h5open(args["out"], "w") do fout
            fout["N"] = N
            fout["dim"] = D
            fout["mean"] = μ
            fout["cov"] = Σ
        end
    end
end

args = parse_commandline()
main(args)

max_stars_repo_name: MartinKocour/SpeechLab

max_issues_repo_name: MartinKocour/SpeechLab

max_forks_repo_name: MartinKocour/SpeechLab


text:
using FileWatching

dir = "."
fm = FolderMonitor(dir)

function reload_chrome()
   run(`osascript chrome.scpt`)
end

watching_files = ["test.css"]

loop = Task() do
    while fm.open
        (fname, events) = wait(fm)::Pair
        println("fname ", fname)
        fname in watching_files && reload_chrome()
    end
end
schedule(loop)

max_stars_repo_name: wookay/TestBukdu

max_issues_repo_name: wookay/TestBukdu

max_forks_repo_name: wookay/TestBukdu


text:
"""
    jgraph_walk(g::JGraph, walk_func)

This function will iteratively change the colors of nodes and edges in a graph
to draw a path along them of the graph.

# Arguments
- `g::JGraph` the `JGraph` on which the walk should be animated
- `walk_func` a function that takes as argument and `AbstractGraph` and returns
a set of nodes linked each to the following one by and edge present in the graph.
One can use `randomwalk` and `non_backtracking_randomwalk` from `LightGraphs`.
"""
function jgraph_walk(g::JGraph, walk_func)

    # TODO add color palette for sequential coloring of repeatedly transversed edges
    Jnodes = jnodes(g)
    Jedges = jedges(g)
    random_walk = walk_func(g.data.graph)
    step_length = g.data.frames[end] ÷ length(random_walk)
    acted_nodes = Dict()
    acted_edges = Dict()
    act!(Jnodes[random_walk[1]], Action(1:step_length, change(:color, colorant"black"=>colorant"red")))
    for (idx, v) in enumerate(random_walk[2:end])
        from = (idx - 1) * step_length + 1
        to = idx * step_length

        if !haskey(acted_nodes, v)
            act!(Jnodes[v], Action(from:to, change(:color, colorant"black"=>colorant"red")))
            acted_nodes[v] = 1
        else
            act!(Jnodes[v], Action(from:to, change(:color, colorant"red"=>colorant"red")))
            acted_nodes[v] += 1
        end

        ce = if haskey(Jedges, v => random_walk[idx])
            v => random_walk[idx]
        else
            random_walk[idx] => v
        end

        if  !haskey(acted_edges, ce)
            act!(Jedges[ce], Action(from:to, change(:color, colorant"black"=>colorant"red")))
            acted_edges[ce] = 1
        else
            act!(Jedges[ce], Action(from:to, change(:color, colorant"red"=>colorant"red")))
            acted_edges[ce] += 1
        end
    end
end


text:
@testset "Frame" begin
    @testset "Basic usage" begin
        frame = Frame()
        @test size(frame) == 0

        resize!(frame, 4)
        @test size(frame) == 4

        pos = positions(frame)
        expected = Array{Float64}(undef, 3, 4)
        for i=1:3, j=1:4
            pos[i, j] = i*j
            expected[i, j] = i*j
        end
        @test positions(frame) == expected

        @test has_velocities(frame) == false
        add_velocities!(frame)
        @test has_velocities(frame) == true

        vel = velocities(frame)
        for i=1:3, j=1:4
            vel[i, j] = i*j
        end
        @test velocities(frame) == expected

        set_cell!(frame, UnitCell(3, 4, 5))
        cell = UnitCell(frame)
        @test lengths(cell) == [3.0, 4.0, 5.0]

        topology = Topology()
        add_atom!(topology, Atom("H"))
        add_atom!(topology, Atom("O"))
        add_atom!(topology, Atom("Zn"))
        add_atom!(topology, Atom("H"))
        set_topology!(frame, topology)

        @test name(Atom(frame, 0)) == "H"
        @test name(Atom(frame, 2)) == "Zn"

        @test step(frame) == 0
        set_step!(frame, 42)
        @test step(frame) == 42

        @test size(frame) == 4

        position = Float64[0.2, 0.8, 0]
        velocity = Float64[1, 2, 1]
        add_atom!(frame, Atom("H"), position, velocity)
        @test size(frame) == 5

        remove_atom!(frame, 4)
        @test size(frame) == 4
    end

    @testset "Frame copy" begin
        frame = Frame()
        @test size(frame) == 0

        resize!(frame, 4)
        @test size(frame) == 4

        copy = deepcopy(frame)
        @test size(copy) == 4

        resize!(copy, 10)
        @test size(copy) == 10
        @test size(frame) == 4
    end

    @testset "Geometry tests" begin
        frame = Frame()
        add_atom!(frame, Atom(""), [0.0, 0.0, 0.0])
        add_atom!(frame, Atom(""), [1.0, 2.0, 3.0])
        @test distance(frame, 0, 1) ≈ sqrt(14.0) atol=1e-10

        frame = Frame()
        add_atom!(frame, Atom(""), [1.0, 0.0, 0.0])
        add_atom!(frame, Atom(""), [0.0, 0.0, 0.0])
        add_atom!(frame, Atom(""), [0.0, 1.0, 0.0])
        @test angle(frame, 0, 1, 2) ≈ pi/2 atol=1e-10

        frame = Frame()
        add_atom!(frame, Atom(""), [1.0, 0.0, 0.0])
        add_atom!(frame, Atom(""), [0.0, 0.0, 0.0])
        add_atom!(frame, Atom(""), [0.0, 1.0, 0.0])
        add_atom!(frame, Atom(""), [0.0, 1.0, 1.0])
        @test dihedral(frame, 0, 1, 2, 3) ≈ pi/2 atol=1e-10

        frame = Frame()
        add_atom!(frame, Atom(""), [0.0, 0.0, 0.0])
        add_atom!(frame, Atom(""), [0.0, 0.0, 2.0])
        add_atom!(frame, Atom(""), [1.0, 0.0, 0.0])
        add_atom!(frame, Atom(""), [0.0, 1.0, 0.0])
        @test out_of_plane(frame, 0, 1, 2, 3) ≈ 2 atol=1e-10
    end

    @testset "Frame add/remove bonds" begin
        frame = Frame();
        add_atom!(frame, Atom("H"), [1.0, 0.0, 0.0])
        add_atom!(frame, Atom("O"), [0.0, 0.0, 0.0])
        add_atom!(frame, Atom("H"), [0.0, 1.0, 0.0])

        add_bond!(frame, 0, 1)
        add_bond!(frame, 1, 2, Chemfiles.TripleBond)

        # the bonds are actually stored inside the topology
        topology = Topology(frame)

        @test bonds(topology) == reshape(UInt64[0,1, 1,2], (2,2))
        @test angles(topology) == reshape(UInt64[0, 1, 2], (3,1))

        @test bond_order(topology, 0, 1) == Chemfiles.UnknownBond
        @test bond_order(topology, 1, 2) == Chemfiles.TripleBond

        remove_bond!(frame, 1, 0)
        @test bonds(Topology(frame)) == reshape(UInt64[1,2], (2,1))
    end

    @testset "Frame add residues" begin
        frame = Frame()
        add_atom!(frame, Atom("Zn"), [0.0, 0.0, 0.0])
        add_atom!(frame, Atom("Fe"), [1.0, 2.0, 3.0])

        residue = Residue("first")
        add_atom!(residue,0)
        add_residue!(frame,residue)

        # residues are actually stored in the topology
        @test count_residues(Topology(frame)) == 1
    end

    @testset "Frame iteration" begin
        frame = Frame()
        add_atom!(frame, Atom(""), [0.0, 0.0, 0.0])
        add_atom!(frame, Atom(""), [0.0, 0.0, 0.0])
        add_atom!(frame, Atom(""), [0.0, 0.0, 0.0])
        add_atom!(frame, Atom(""), [0.0, 0.0, 0.0])

        count = 0
        for atom in frame
            count += 1
        end

        @test count == 4
    end
end


text:
include("SingleIntegrator.jl")
include("LinearSystem.jl")
include("MultipleEnvs.jl")
include("ReferenceModel.jl")
include("TwoDimensionalNonlinearOscillator.jl")
include("TwoDimensionalNonlinearPolynomialSystem.jl")
include("TwoDimensionalNonlinearDTSystem.jl")


text:
using EllZeroTrendFiltering
using Test, Random, LinearAlgebra
import Printf: @printf

tests = [
    "quadratic_polynomials.jl",
    "quadratic_forms.jl",
    "piecewise_quadratics.jl",
    "continuous_transition_costs.jl",
    "discrete_transition_costs.jl",
    "brute_force_search.jl",
    "pwq_dp_constrained.jl",
    "exact_tests.jl",
    "fit_pwl_constrained.jl",
    "fit_pwl_regularized.jl",
    "subsample_test.jl",
    "examples.jl"]

@testset "All Tests" begin
    @testset "Testfile: $test" for test in tests
        include(test)
    end
end


text:
abstract type AbstractUniqueIndex <: AbstractIndex; end

struct UniqueIndex <: AbstractUniqueIndex; end

accelerate(a::AbstractArray, ::Type{UniqueIndex}) = AcceleratedArray(a, UniqueIndex())

Base.unique(a::AcceleratedArray{<:Any, <:Any, <:Any, <:AbstractUniqueIndex}) = a


text:
using Yao, Zygote, YaoFlux
import Zygote: Context

function mapply!(r::ArrayReg, rb::RotationGate)
    # v0 = copy(r.state)
    # apply!(r, rb.block)
    # NOTE: we should not change register's memory address,
    # or batch operations may fail
    # r.state .= -im*sin(rb.theta/2)*r.state + cos(rb.theta/2)*r.state
    t = -im*sin(rb.theta/2)*r.state
    copyto!(r.state, t)
    return r
end

t = ArrayReg(bit"0") + ArrayReg(bit"1")
normalize!(t)

function loss(x)
    r = mapply!(zero_state(1), Rx(x))
    return fidelity(t, r)[]
end

Zygote.refresh()
_, back = Zygote._forward(loss, 0.1)

back(1)
(loss(0.1+1e-8) - loss(0.1)) / 1e-8

r = ComplexF64[1, 0]
t = rand(ComplexF64, 2)

function f1!(r, x)
    tt = -im * sin(x) * r
    copyto!(r, tt)
    return abs(t' * r)
end

function f2!(r, x)
    tt = -im * sin(x) * r
    return abs(t' * r)
end

Zygote.refresh()
gradient(f1!, r, 0.1)[2]

(f1!(statevec(zero_state(1)), 0.1+sqrt(eps())/2) - f1!(statevec(zero_state(1)), 0.1-sqrt(eps())/2)) / sqrt(eps())

gradient(f2!, r, 0.1)[2]
(f2!(statevec(zero_state(1)), 0.1+1e-8) - f2!(statevec(zero_state(1)), 0.1)) / 1e-8


text:
using PkgSRA
using Plots, Random, Distributions, Interpolations
pyplot()
Random.seed!(11)
using BSON: @save
include("../scripts/paper_settings.jl");

#####
##### Define the controlled dynamical system
#####
# Load example problem
include(EXAMPLE_FOLDERNAME*"example_falling_ball.jl")
U_starts = [6.0]
    U_widths = 0.2;
    amplitude = 75.0;
U_func_kick(t, u) = U_func_time(t, u,
                                U_widths, U_starts,
                                F_dim=2,
                                amplitude=amplitude)

# Define the spatial forcing function

U_func_wall(X) = U_func_spring(X, k=1e3, r=1.0);

# Solve the system
prob = ODEProblem((du, u, p, t)->
        ball_system(du, u, p, t,
             U_func_time=U_func_kick,
             U_func_space=U_func_wall),
         u1, tspan, p1)
dt = ts[2] - ts[1]

dyn_with_ctr = solve(prob, Tsit5(), saveat=ts, adaptive=false, dt = dt);
# dyn_with_ctr = solve_ball_system(U_func_space=U_func_wall,
#                                  U_func_time=U_func_kick)
dat = Array(dyn_with_ctr)
dat_with_const = vcat(dat, ones(1,size(dat,2)))
numerical_grad = numerical_derivative(dat, ts)

# Get just the control signals
dyn_control_wall = zeros(size(dat))
for i in 1:size(dat, 2)
    dyn_control_wall[:,i] = U_func_wall(dat[:, i])
end

dyn_control_kick = zeros(size(dat))
for (i, t) in enumerate(ts)
    dyn_control_kick[:,i] = U_func_kick(t, dat[:, i])
end

U_true = dyn_control_kick .+ dyn_control_wall

#####
##### Fit a naive linear model (with constant)
#####
naive_model = dmdc(dat_with_const, numerical_grad)[1]

f(u,p,t) = vcat(naive_model * u, [1])
prob = ODEProblem(f, vcat(u1,[1]), tspan, p1)
dat_naive = Array(solve(prob, Tsit5(), saveat=ts))

# plot(dat_naive[1,:]); plot!(dat[1,:])

#####
##### Build SRA object and analyze
#####
# Initialize
this_model = sra_stateful_object(ts, tspan, dat, u1, numerical_grad)
prams = this_model.parameters
prams.sindyc_ensemble_parameters[:selection_criterion] =
    sindy_cross_validate;
# Several changes are required because there are only 2 variables here
# prams.variable_names = ["x", "v"];
prams.sindy_terms_list = Iterators.product(1:2, 1:2)
prams.sindy_library["cross_terms"] = [2, 3] # Also include cubic terms

fit_first_model(this_model, 1);

# U_true = vcat(dyn_control_kick, dyn_control_wall)
ind = 1:500
i = 2
    plot(residual_unctr[i,ind], ribbon=sample_noise)
    plot!(U_true[:,ind]', label="True")
    title!("Residual and true control")

## Subsample
accepted_ind = subsample_using_residual(
            residual_unctr, sample_noise,
            min_length=1, shorten_length=2)
num_pts = 400; start_ind = 101
subsample_ind = accepted_ind[start_ind:num_pts+start_ind-1]

U_sub = U_true[:,subsample_ind]
plot(U_sub')
    # plot!(residual[:,subsample_ind]', ribbon=sample_trajectory_noise)
    title!("Control signals in the subsampled dataset")

chain_ctr, best_sindy_ctr = calc_distribution_of_models(
    dat[:,subsample_ind],
    numerical_grad[:,subsample_ind],
    sindy_library,
    val_list = Iterators.product(1:3,1:3), # DEBUG
    # val_list = combinations(1:3, 2),
    # val_list = calc_permutations(3,2),
    chain_opt = (iterations=300, train_ind=1:num_pts)
)

# print_equations(best_sindy_ctr)
# plot(chain_ctr)

# Calculate updated control signal
(residual_ctr, _, noise_ctr, _) =
        calc_distribution_of_residuals(
        dat, numerical_grad, chain_ctr, 1:length(ts), best_sindy_ctr)
ctr_final = process_residual(residual_ctr, noise_ctr)

plot(ctr_final', label="Learned", lw=3, color=:black)
    plot!(U_true', label="True")

# Get the actually controlled model
# val_list = calc_permutations(3,2)
# val_list = combinations(1:3, 2)
val_list = Iterators.product(1:31:3) # DEBUG
(final_sindy_model,best_criterion,all_criteria,_, _) =
     sindyc_ensemble(dat, numerical_grad, sindy_library, val_list,
                     U=ctr_final, ts=ts,
                     selection_criterion=my_aicc,
                     sparsification_mode="num_terms",
                     selection_dist=Normal(0,1),
                     use_clustering_minimization=false)

scatter(sum.(val_list), all_criteria)

prob = ODEProblem(final_sindy_model, dat[:,1], tspan)
# ts2 = range(tspan[1], tspan[end], length=10000)
dt = ts[2] - ts[1]
sol = solve(prob, AB3(), dt=dt, saveat=ts)
dat_final = Array(sol)

plot(sol)
#####
##### Save
#####
this_dat_name = DAT_FOLDERNAME*"dat_ball_"

# Raw data
fname = this_dat_name*"raw.bson";
@save fname dat numerical_grad dyn_control_kick dyn_control_wall

# Naive model
fname = this_dat_name*"naive_model.bson";
@save fname dat_naive naive_model

# Controlled model
fname = this_dat_name*"controlled_model.bson";
@save fname chain_ctr ctr_final dat_final


text:
module GtkGUI

    if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@optlevel"))
        @eval Base.Experimental.@optlevel 1
    end

    using MappedArrays, IndirectArrays, ColorSchemes, Images, GaussianBeamLens, GtkReactive, Gtk, Gtk.ShortNames, ImageView

    export gui

    const cache = Dict{Tuple{Int64,Int64,Int64,Int64},IndirectArray{RGB{Float64}, 2}}()

    """
        color_me_scaleminmax(A, cmap)

    Apply colormap `cmap` to an array `A`.
    See https://discourse.julialang.org/t/how-to-convert-a-matrix-to-an-rgb-image-using-images-jl/7265/9 by Tim Holy
    """
    function color_me_scaleminmax(A, cmap)
        n = length(cmap)
        scale = takemap(scaleminmax, A)
        f = s->clamp(round(Int, (n-1)*scale(s))+1, 1, n)  # safely convert 0-1 to 1:n
        Ai = mappedarray(f, A)       # like f.(A) but does not allocate significant memory
        IndirectArray(Ai, cmap)      # colormap array
    end

    """
        GaussianBeamLensPropagateColored(args...; colormap = ColorSchemes.inferno.colors)

    Apply a colormap to the output of `GaussianBeamLens.GaussianBeamLensPropagate`. See that function for args.
    """
    GaussianBeamLensPropagateColored(args...; colormap = ColorSchemes.inferno.colors) =
        haskey(cache, args) ?
            cache[args] :
            cache[args] = color_me_scaleminmax(GaussianBeamLens.GaussianBeamLensPropagateUncached(args...)', colormap)

    function gui()
        win = Window("GaussianBeamLens.jl") |> (bx = Box(:v))
        I = GaussianBeamLensPropagateColored(0, 25, 1000, 1000)'
        c = canvas(UserUnit, size(I,2), size(I,1))
        imshow!(c,I)
        x0sl = slider(-240:10:240)
        w0sl = slider(2:50)
        fsl = slider(0:10:2000)
        Lens_zsl = slider(0:10:2000)
        push!(bx, c)
        push!(bx, label("x_0, Lateral offset of the input beam"))
        push!(bx, x0sl)
        push!(bx, label("w_0, Width of the Gaussian beam in object space"))
        push!(bx, w0sl)
        push!(bx, label("f, Focal length of the lens"))
        push!(bx, fsl)
        push!(bx, label("Lens_z, Axial position of the lens"))
        push!(bx, Lens_zsl)

        signals = map(signal, [x0sl, w0sl, fsl, Lens_zsl])
        redraw = draw(c, signals...) do c, x0, w0, f, Lens_z
            copy!(c, GaussianBeamLensPropagateColored(x0, w0, f, Lens_z))
        end

        Gtk.showall(win)
        (win = win, c = c, x0sl = x0sl, w0sl = w0sl, fsl = fsl, Lens_zsl = Lens_zsl, redraw = redraw)
    end

    function __init__()
        @time gui()
        readline()
    end

    precompile(color_me_scaleminmax, (Matrix{Float64}, Vector{RGB{Float64}}))
    precompile(GaussianBeamLensPropagateColored,(Int64,Int64,Int64,Int64))
    precompile(gui,())
end

max_stars_repo_name: mkitti/GaussianBeamLens

max_issues_repo_name: mkitti/GaussianBeamLens

max_forks_repo_name: mkitti/GaussianBeamLens


text:
function read_d_drp(drp)
    @info "  - Dutch data training set"
    dp = DataFrame(ID          = String[],  #  1. Interbull ID
                   Milk        = Float64[], #  2. milk yield DRP
                   Milk_ERC    = Float64[], #  3. milk yield ERC
                   Fat         = Float64[], #  4. fat yield DRP
                   Fat_ERC     = Float64[], #  5. fat yield ERC
                   Protein     = Float64[], #  6. protein yield DRP
                   Protein_ERC = Float64[], #  7. protein yield ERC
                   SCS         = Float64[], #  8. somatic cell score DRP
                   SCS_ERC     = Float64[], #  9. somatic cell score ERC
                   NR56        = Float64[], # 10. NR56 DRP
                   NR56_ERC    = Float64[], # 11. NR56 ERC
                   x12         = Float64[], # 12. calving to 1st ins DRP
                   x13         = Float64[], # 13. calving to 1st ins ERC
                   x14         = Float64[], # 14. 1st-Last ins. Cows DRP
                   x15         = Float64[], # 15. 1st-Last ins. Cows ERC
                   x16         = Float64[], # 16. calving interval DRP
                   x17         = Float64[], # 17. calving interval ERC
                   x18         = Int[]      # 18. Numerical ID (as in “all.ped”)
                   )
    for line in eachline(drp)
        v = split(line)
        id = v[1]
        y = parse.(Float64, v[2:end-1])
        n = parse(Int, v[end])
        push!(dp, [id; y; n])
    end
    dp
end

function read_g_drp(drp)
    @info "  - German data training set"
    gp = DataFrame(ID          = String[],  #  1. Interbull ID
                   Milk        = Float64[], #  2. milk yield DRP
                   Milk_ERC    = Float64[], #  3. milk yield ERC
                   Fat         = Float64[], #  4. fat yield DRP
                   Fat_ERC     = Float64[], #  5. fat yield ERC
                   Protein     = Float64[], #  6. protein yield DRP
                   Protein_ERC = Float64[], #  7. protein yield ERC
                   SCS         = Float64[], #  8. somatic cell score DRP
                   SCS_ERC     = Float64[], #  9. somatic cell score ERC
                   NR56        = Float64[], # 10. NR56 heifers DRP
                   NR56_ERC    = Float64[], # 11. NR56 heifers ERC
                   x12         = Float64[], # 12. calving to 1st ins DRP
                   x13         = Float64[], # 13. calving to 1st ins ERC
                   x14         = Float64[], # 14. NR56 cows DRP
                   x15         = Float64[], # 15. NR56 cows ERC
                   x16         = Float64[], # 16. 1st-Last ins. Cows DRP
                   x17         = Float64[], # 17. 1st-Last ins. Cows ERC
                   x18         = Float64[], # 18. Days open cows DRP
                   x19         = Float64[], # 19. Days open cows ERC
                   N           = Int[]      # 20. Numerical ID (as in “all.ped”)
                   )
    for line in eachline(drp)
        v = split(line)
        id = v[1]
        y = parse.(Float64, v[2:end-1])
        n = parse(Int, v[end])
        push!(gp, [id; y; n])
    end
    gp
end

function read_n_drp(drp)
    @info "  - Norwegian data training set"
    np = DataFrame(ID       = String[],
                   Milk     = Float64[],
                   Milk_ERC = Float64[],
                   SCS      = Float64[],
                   SCS_ERC  = Float64[],
                   N        = Int[])
    for line in eachline(drp)
        v = split(line)
        length(v) != 6 && continue # because two lines have only 2 columns in trnn set
        id = v[1]
        y = parse.(Float64, v[2:end-1])
        n = parse(Int, v[end])
        push!(np, [id; y; n])
    end
    np
end

"""
    function read_training_20210311()
---
## For bulls
I checked the DRP files: the weights in the files are ERC’s (Effective Record
Contributions) and not reliabilities, so that does seem to be correct.

The values do look somewhat funny, because I rescaled all heritabilities to be
0.5 (details are in the document attached, which I sent last year).

With `h² = 0.5, λ =(1 - h²)/h² = (1 - 0.5)/0.5 = 1`

With a reliability of 99.0 (on a 0-100% scale), the ERC becomes:
`ERC = λ⋅rel/(1-rel) = 1 × 0.99/(1-0.99) = 99.0` !
(note that `rel` is the reliability on a 0-1 scale).

With a reliability of 98.0 (on a 0-100% scale), the ERC becomes:
`ERC = λ⋅rel/(1-rel) = 1×0.98/(1-0.98) = 49.0`

## And for cows:

With a reliability of 63.0 (on a 0-100% scale), the ERC becomes:
`ERC = λ⋅rel/(1-rel) = 1×0.63/(1-0.63) = 1.7

So, I think we need to revert to simpler deregression (as we discussed last
time!) which simply is: `DRP = EBV/rel`, after first subtracting the mean EBV.

I have attached files with DRPs computed using this simpler deregression; the
files are in the same format as before.  The weights are still ERC’s.

## Two additional points:

I think there may have been a mistake in the file with Norwegian DRPs, where I
linked IDs to the wrong data.  That should be fixed now.

The difference between the Dutch bulls and cows (at least this is what I think
we saw from one of Xijiang’s plots), may be due to a different base pertaining
to the cow and bull EBVs.  So I wonder if I should do the deregression
separately for cows and bulls.  Note that with the deregression the mean of the
DRP is removed.  Doing it separately for cows and bulls would effectively give
both the same mean.

Anyway, I’m curious to see the results with these new DRPs.

-- Mario, 2021.3
"""
function read_training_20210311()
    trn_dir = joinpath(dat_dir, "mario/20210311")
    dp = read_n_drp(joinpath(trn_dir, "dutch_DRPs.txt"))
    gp = read_g_drp(joinpath(trn_dir, "german_DRPs.txt"))
    np = read_n_drp(joinpath(trn_dir, "norway_DRPs.txt"))
    return (dp, gp, np)
end

function refresh_drp_20210311()
    @info "Read phenotypes, EBV and CV sets"
    dts, gts, nts = read_training() # ts = training set
    @save joinpath(dat_dir, "jld/drp-training-20210311.jld") {compress=true} dts gts nts
end

function refresh_norsk_drp_210408()
    @info "Read phnotypes, EBV and CV sets"
    @load "$dat_dir/jld/drp-training-2021-03-11.jld" dts gts nts
    nts = read_n_drp(joinpath(dat_dir, "mario/20210408/norway_DRPs.txt"))
    @save joinpath(dat_dir, "jld/drp-training-2021-04-08.jld") {compress=true} dts gts nts
end

"""
The results for SCS look very plausible. These results also seem to confirm that for milk, the results where the German population is involved, do not match expectations.

The only reasonable explanation that I can come up with, is that the EBV (for milk) in the different countries refer to a different base. Before the deregression I make the scale of EBV the same in all countries (by dividing by the genetic standard deviation within population). During the deregression the mean is removed within each country. This may cause inconsistencies across populations, for instance, if the Dutch population has on average higher BV than the German population (on the “true” scale). Because the German data includes cows only, perhaps subtracting the mean EBV is not really appropriate, because the mean EBV in the data may not be representative for the mean in the entire population (see also histograms below; these include EBV that were used to compute deregressed EBV).

To test: could you re-do the analyses using the attached data for Germany? The only difference, is that in this case I did not remove the mean from the EBV for milk yield (during the deregression).

This function is together with
`RD2.cv_21_04_22("drp-training-2021-04-21.jld")`
"""
function refresh_german_drp_210421()
    @info "Read phenotypes, EBV and CV sets"
    @load "$dat_dir/jld/drp-training-2021-04-08.jld" dts gts nts
    gts = read_g_drp(joinpath(dat_dir, "mario/20210421/german_DRPs.txt"))
    @save joinpath(dat_dir, "jld/drp-training-2021-04-21.jld") {compress=true} dts gts nts
end


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

export
    av_hash_alloc,
    av_hash_names,
    av_hash_get_name,
    av_hash_get_size,
    av_hash_init,
    av_hash_update,
    av_hash_final,
    av_hash_final_bin,
    av_hash_final_hex,
    av_hash_final_b64,
    av_hash_freep

function av_hash_alloc(ctx, name)
    ccall((:av_hash_alloc, libavutil), Cint, (Ptr{Ptr{AVHashContext}}, Cstring), ctx, name)
end

function av_hash_names(i::Integer)
    ccall((:av_hash_names, libavutil), Cstring, (Cint,), i)
end

function av_hash_get_name(ctx)
    ccall((:av_hash_get_name, libavutil), Cstring, (Ptr{AVHashContext},), ctx)
end

function av_hash_get_size(ctx)
    ccall((:av_hash_get_size, libavutil), Cint, (Ptr{AVHashContext},), ctx)
end

function av_hash_init(ctx)
    ccall((:av_hash_init, libavutil), Cvoid, (Ptr{AVHashContext},), ctx)
end

function av_hash_update(ctx, src, len::Integer)
    ccall((:av_hash_update, libavutil), Cvoid, (Ptr{AVHashContext}, Ptr{UInt8}, Cint), ctx, src, len)
end

function av_hash_final(ctx, dst)
    ccall((:av_hash_final, libavutil), Cvoid, (Ptr{AVHashContext}, Ptr{UInt8}), ctx, dst)
end

function av_hash_final_bin(ctx, dst, size::Integer)
    ccall((:av_hash_final_bin, libavutil), Cvoid, (Ptr{AVHashContext}, Ptr{UInt8}, Cint), ctx, dst, size)
end

function av_hash_final_hex(ctx, dst, size::Integer)
    ccall((:av_hash_final_hex, libavutil), Cvoid, (Ptr{AVHashContext}, Ptr{UInt8}, Cint), ctx, dst, size)
end

function av_hash_final_b64(ctx, dst, size::Integer)
    ccall((:av_hash_final_b64, libavutil), Cvoid, (Ptr{AVHashContext}, Ptr{UInt8}, Cint), ctx, dst, size)
end

function av_hash_freep(ctx)
    ccall((:av_hash_freep, libavutil), Cvoid, (Ptr{Ptr{AVHashContext}},), ctx)
end


text:
module DevilSpawn

import Base.Threads

"""
    @spawn_background expr

Like `@spawn`, but ensures that the spawned task isn't scheduled on the main thread, except
if `nthreads() == 1`.

If all background threads are currently blocked, we retry until one becomes available using
a small exponential backoff. A `yield()` call or automatically placed yield point will free
up the thread and allow another task to be scheduled.
"""
macro spawn_background(exprs...)
    esc(quote
        $Threads.@spawn $spawn_background(()->$(exprs...))
    end)
end

function spawn_background(f, attempt_num=0)
    if Threads.threadid() == 1 && Threads.nthreads() > 1
        sleep(1E-9 * (1 << min(32, attempt_num)))
        # Maintain demand on the new Task, so that errors will propagate.
        return fetch(Threads.@spawn spawn_background(f, attempt_num+1))
    else
        f()
    end
end

export @spawn_background

end # module


text:
using GLVisualize, GeometryTypes, Colors

function visu(x, sym=:default)
    tic()
    @profile r = visualize(x, sym)
    t = toq()
    println("visualize $t")
    tic()
    @profile _view(r)
    t = toq()
    println("view $t")
    tic()
    r = visualize(x, sym, boundingbox=nothing, color=RGBA{Float32}(0,0,0,0))
    t = toq()
    println("no bb $t")
end

function test(w)
    x = rand(Float32, 10, 10)
    visu(x)
    visu(x, :surface)
    visu((x,x,x), :surface)

    y = rand(Point3f0, 10)
    visu(y)
    visu((GLVisualize.CIRCLE, y))
    visu(y, :lines)
    visu(y, :line_segments)
    cat = loadasset("cat.obj")
    visu(cat)
end

function profileit()
    w = GLVisualize.current_screen()
    test(w)
    test(w)
    empty!(w)
    Profile.init()
    Profile.clear()
    for i=1:100
        test(w)
        empty!(w)
        println("#################################################")
    end
end
w = glscreen()
profileit()
using ProfileView
ProfileView.view()


text:
using Test

@testset "RRTMGP gray radiation tests" begin

    include("gray_atm.jl")

end

@testset "RRTMGP clear-sky tests" begin
    include("clear_sky.jl")
end

@testset "RRTMGP clear-sky tests" begin
    include("all_sky.jl")
end

max_stars_repo_name: climate-machine/JRRTMGP

max_forks_repo_name: climate-machine/JRRTMGP