text:
# JuMP:         for DenseAxisArray
# Random:       for random strings
# JLD2:          to save data to Harddrive
module DataCreator

using JuMP, Random, JLD, Dictionaries, NamedDims

export createAllData

# create array data used to generate all other data
function createBasicData(n::Int64)

    # create random data used for all tests
    arrayCoeff = rand(n,n)
    arrayBound = (ones(n)+rand(n))*sqrt(n)/2

    return arrayCoeff, arrayBound
end

# create array of strings to use for all models with string indexing
function createArrayString(n::Int64)
    str_obj = Array{String,1}(undef,n)
    str_con = similar(str_obj)

    #string have a length of 20 and are numbered
    for p in 1:n
		str_obj[p] = string("obj_", randstring(15-length(string(p))), "_", string(p))
		str_con[p] = string("con_", randstring(15-length(string(p))), "_", string(p))
    end

    return str_obj, str_con
end

# denseAxisArray data
function createDenseAxisArrayData(coeff::Array{Float64,2}, bound::Array{Float64,1})

	# random n by n matrix
	denseAxisArrayCoeff = JuMP.Containers.DenseAxisArray(coeff,1:length(coeff[:,1]),1:length(coeff[1,:]))

    # vector used defining contraints
    denseAxisArrayBound = JuMP.Containers.DenseAxisArray(bound,1:length(bound))

    return denseAxisArrayCoeff, denseAxisArrayBound
end

# denseAxisArray data
function createDenseAxisArrayDataStr(coeff::Array{Float64,2}, bound::Array{Float64,1}, str_obj::Array{String,1}, str_con::Array{String,1})

	# random n by n matrix
	denseAxisArrayCoeffStr = JuMP.Containers.DenseAxisArray(coeff, str_con, str_obj)

    # vector used defining contraints
    denseAxisArrayBoundStr = JuMP.Containers.DenseAxisArray(bound, str_con)

    return denseAxisArrayCoeffStr, denseAxisArrayBoundStr
end

# denseAxisArray data
function createDictData(coeff::Array{Float64,2}, bound::Array{Float64,1})

    DictCoeff = Dict()
    DictBound = Dict()

    for i in 1:length(coeff[:,1])
        for j in 1:length(coeff[1,:])
            DictCoeff[i,j] = coeff[i,j]
        end
        DictBound[i] = bound[i]
    end

    return DictCoeff, DictBound
end

# denseAxisArray data
function createDictDataStr(coeff::Array{Float64,2}, bound::Array{Float64,1}, str_obj::Array{String,1}, str_con::Array{String,1})

    DictCoeff = Dict{Tuple{String,String},Float64}()
    DictBound = Dict{String,Float64}()

    for i in 1:length(coeff[:,1])
        for j in 1:length(coeff[1,:])
            DictCoeff[str_con[i],str_obj[j]] = coeff[i,j]
        end
        DictBound[str_con[i]] = bound[i]
    end

    return DictCoeff, DictBound
end

function createDictionaryData(coeff::Array{Float64,2}, bound::Array{Float64,1})

    dictCoeff = Dictionary(coeff)
    dictBound = Dictionary(bound)

    return dictCoeff, dictBound
end

function createDictionaryData(coeff::Array{Float64,2}, bound::Array{Float64,1})

    dictCoeff = Dictionary(coeff)
    dictBound = Dictionary(bound)

    return dictCoeff, dictBound
end

# denseAxisArray data
function createNamedDimsData(coeff::Array{Float64,2}, bound::Array{Float64,1})

    NamedDimsCoeff = NamedDimsArray{(:x1,:x2)}(coeff)
    NamedDimsBound = NamedDimsArray{(:x1,)}(bound)

    return NamedDimsCoeff, NamedDimsBound
end

# create all data and save to harddrive for specified problem size n
function createAllData(n::Int)
    arrayCoeff, arrayBound = createBasicData(n)
    str_obj, str_con = createArrayString(n)
    denseAxisArrayCoeff, denseAxisArrayBound = createDenseAxisArrayData(arrayCoeff, arrayBound)
    denseAxisArrayCoeffStr, denseAxisArrayBoundStr = createDenseAxisArrayDataStr(arrayCoeff, arrayBound,str_obj, str_con)
    dictCoeff, dictBound = createDictData(arrayCoeff, arrayBound)
    dictCoeffStr, dictBoundStr = createDictDataStr(arrayCoeff, arrayBound, str_obj, str_con)
    dictionaryCoeff, dictionaryBound = createDictionaryData(arrayCoeff, arrayBound)
    NamedDimsCoeff, NamedDimsBound = createNamedDimsData(arrayCoeff, arrayBound)

    #save data to .jld file
    save("Data/arrayData.jld", "n", n, "arrayCoeff", arrayCoeff, "arrayBound", arrayBound)

    save("Data/denseAxisArrayData.jld", "n", n, "arrayCoeff", denseAxisArrayCoeff, "arrayBound", denseAxisArrayBound)
    save("Data/denseAxisArrayDataStr.jld", "n", n, "str_obj", str_obj, "str_con", str_con, "arrayCoeff", denseAxisArrayCoeffStr, "arrayBound", denseAxisArrayBoundStr)

    save("Data/dictData.jld", "n", n, "arrayCoeff", dictCoeff, "arrayBound", dictBound)
    #save("Data/dictDataStr.jld", "n", n, "str_obj", str_obj, "str_con", str_con, "arrayCoeff", dictCoeffStr, "arrayBound", dictBoundStr)

    save("Data/dictionaryData.jld", "n", n, "arrayCoeff", dictionaryCoeff, "arrayBound", dictionaryBound)

    save("Data/namedDimsData.jld", "n", n, "arrayCoeff", NamedDimsCoeff, "arrayBound", NamedDimsBound)

    println("Data Saved succesfully")
end

end

max_stars_repo_name: SonntagK/JuMPContainerBenchmarks

max_issues_repo_name: SonntagK/JuMPContainerBenchmarks

max_forks_repo_name: SonntagK/JuMPContainerBenchmarks


text:
#! /usr/bin/env julia

# raise exception when receiving SIGINT
ccall(:jl_exit_on_sigint, Void, (Cint,), 0)

using PipeLayout
using PipeLayout.GndStr

function solve_with(instance::AbstractString, solver::GroundStructureSolver)
    inst = PipeLayout.read_instance(instance)
    topo = PipeLayout.read_topology(instance)

    optimize(inst, topo, solver)
end

if length(ARGS) != 2
    println("usage: run.jl <config> <instname>")
    exit(1)
end

config = ARGS[1]
instname = ARGS[2]

include(config) # creates solver

println("solving $instname with $(basename(config))")

try
    tic()
    println("vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv")
    result = solve_with(instname, solver)
    println("^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^")
    toc()
    println("Status: $(result.status)")
    if result.solution != nothing
        println("Solution: $(findn(result.solution.zsol))")
        println("Primal bound: $(result.primalbound)")
        println("Dual bound: $(result.dualbound)")
    else
        println("No solution")
    end
catch ex
    if isa(ex, InterruptException)
        println("-- was interrupted --")
    else
        rethrow()
    end
end

ext: jl


text:
export GlobalToLocal, getGlobalToLocal
export interpLocalToGlobal, interpGlobalToLocal
export prepareGlobalToLocal

"""
mutable struct jInv.InverseSolve.GlobalToLocal

Maps global model to local model

sigLocal = PForInv*sigGlobal + sigmaBackground

Fields:
	PForInv::SparseMatrixCSC - linear operator between inverse mesh and
	                           forward mesh, e.g., linear interpolation matrix
                               and projection on active set.
	sigmaBackground::Vector  - background model

Constructors:
    getGlobalToLocal(P::SparseMatrixCSC)
	getGlobalToLocal(P::SparseMatrixCSC,sigmaBack::Vector)

Example:
	Mesh2Mesh = getInterpolationMatrix(Minv,Mfwd)   # Minv and Mfwd have different resolutions
	sigmaBack = 1.2*ones(Minv.nc)                   # put background conductivity
	gloc      = getGlobalToLocal(Mesh2Mesh,sigmaBack)

	sigLocal     = gloc.PForInv\' * sigGlobal + gloc.sigmaBack
	# this call is equivalent, but needed in case the Mesh2Mesh matrix is compressed
	sigLocalFast = interpGlobalToLocal(sigGlobal,gloc.PForInv,gloc.sigmaBack)
"""
struct GlobalToLocal
	PForInv::Union{SparseMatrixCSC,AbstractFloat,AbstractModelTransform} # interpolation matrix from fwd mesh to inv mesh
	sigmaBackground::Union{Vector{Float64},AbstractFloat,AbstractModel} #  (# of cells fwd mesh)
end # mutable struct GlobalToLocal

# Constructors
getGlobalToLocal(P) = GlobalToLocal(P,1e-8)
getGlobalToLocal(P,sigBack) = GlobalToLocal(P,sigBack)
getGlobalToLocal(P,sigBack,fname) = GlobalToLocal(P,sigBack)

function prepareGlobalToLocal(Mesh2Mesh,Iact,sigmaBackground,fname)
	return getGlobalToLocal(Iact'*Mesh2Mesh,interpGlobalToLocal(sigmaBackground,Mesh2Mesh),fname)
end


text:
# ------------------------------------------------------------------------------------------
# # Julia 绘图简介
#
# Julia 有包括调用 `PyPlot` 在内的、多种不同的绘图方式。
#
# 这里将介绍如何使用 `Plots.jl`。
# ------------------------------------------------------------------------------------------

using Pkg; Pkg.add("Plots")
using Plots

# ------------------------------------------------------------------------------------------
# `Plots.jl` 的优势之一是可以无缝地切换后端（backends）。
# 在这个 notebook 中，我们将尝试 `gr()` 和 `plotlyjs()` 后端。
#
# 以科学调查之名，我们用这个 notebook 调查一下大概1860到2000年的全球温度和海盗数量的关系。
# ------------------------------------------------------------------------------------------

globaltemperatures = [14.4, 14.5, 14.8, 15.2, 15.5, 15.8]
numpirates = [45000, 20000, 15000, 5000, 400, 17]

# ------------------------------------------------------------------------------------------
# **要绘制出图线，先载入GR后端**
# ------------------------------------------------------------------------------------------

gr()

# ------------------------------------------------------------------------------------------
# 现在我们可以调用如 `plot` 和 `scatter` 的函数来绘制图像。
# ------------------------------------------------------------------------------------------

plot(numpirates, globaltemperatures, label="line")
scatter!(numpirates, globaltemperatures, label="points")

# ------------------------------------------------------------------------------------------
# `scatter!` 函数名后面的 `!` 意味着它是一个原地修改传入变量的函数，表示散点图会被添加到已存在的图像上。
#
# 与此对应的，你可以试试换成 `scatter` 看看会发生什么。
#
# 接着，我们通过 `xlabel!`，`ylabel!` 和 `title!` 函数来给图像加上更多信息。
# ------------------------------------------------------------------------------------------

xlabel!("Number of Pirates [Approximate]")
ylabel!("Global Temperature (C)")
title!("Influence of pirate population on global warming")

# ------------------------------------------------------------------------------------------
# 这看上去仍不对劲。自1860以来海盗数量是减少的，而从左往右看其实时间上是倒序的。
# 我们来把 X 轴反过来，可以更清楚地看到时间顺序下海盗人口是如何导致全球温度的变化！
# ------------------------------------------------------------------------------------------

xflip!()

# ------------------------------------------------------------------------------------------
# 绘图就是这样简单！
#
# 注意：这是一个关于人们是如何经常结合相关性和因果性的笑话。
#
# **不需要修改语法，我们可以在 `unicodeplots()` 后端中绘制同样的图像**
# ------------------------------------------------------------------------------------------

using Pkg; Pkg.add("UnicodePlots")
unicodeplots()

plot(numpirates, globaltemperatures, label="line")
scatter!(numpirates, globaltemperatures, label="points")
xlabel!("Number of Pirates [Approximate]")
ylabel!("Global Temperature (C)")
title!("Influence of pirate population on global warming")

# ------------------------------------------------------------------------------------------
# 注意到第二幅图和第一幅的区别！
# ------------------------------------------------------------------------------------------

max_stars_repo_name: NahsiN/JuliaTutorials

max_issues_repo_name: NahsiN/JuliaTutorials

max_forks_repo_name: NahsiN/JuliaTutorials


text:
module TestUnivariateFiniteMethods

using Test
using MLJBase
using CategoricalArrays
import Distributions:pdf, logpdf, support
import Distributions
using StableRNGs
import Random
rng = StableRNG(123)

v = categorical(collect("asqfasqffqsaaaa"), ordered=true)
V = categorical(collect("asqfasqffqsaaaa"))
a, s, q, f = v[1], v[2], v[3], v[4]
A, S, Q, F = V[1], V[2], V[3], V[4]

@testset "set 1" begin

    # ordered (OrderedFactor)
    dict = Dict(s=>0.1, q=> 0.2, f=> 0.7)
    d    = UnivariateFinite(dict)
    @test classes(d) == [a, f, q, s]
    @test classes(d) == classes(s)
    @test levels(d) == levels(s)
    @test support(d) == [f, q, s]
    @test MLJBase.sample_scitype(d) == OrderedFactor{4}
    # levels!(v, reverse(levels(v)))
    # @test classes(d) == [s, q, f, a]
    # @test support(d) == [s, q, f]
    @test ismissing(pdf(d, missing))
    @test pdf(d, s) ≈ 0.1
    @test pdf(d, 's') ≈ 0.1
    @test pdf(d, q) ≈ 0.2
    @test pdf(d, 'q') ≈ 0.2
    @test pdf(d, f) ≈ 0.7
    @test pdf(d, 'f') ≈ 0.7
    @test pdf(d, a) ≈ 0.0
    @test logpdf(d, s) ≈ log(0.1)
    @test logpdf(d, 's') ≈ log(0.1)
    @test logpdf(d, q) ≈ log(0.2)
    @test logpdf(d, 'q') ≈ log(0.2)
    @test logpdf(d, f) ≈ log(0.7)
    @test logpdf(d, 'f') ≈ log(0.7)
    @test isinf(logpdf(d, a))
    @test mode(d) == f

    @test UnivariateFinite(support(d), [0.7, 0.2, 0.1]) ≈ d

    N = 50
    rng = StableRNG(125)
    samples = [rand(rng,d) for i in 1:50];
    rng = StableRNG(125)
    @test samples == [rand(rng, d) for i in 1:N]

    N = 10000
    samples = rand(StableRNG(123), d, N);
    @test Set(samples) == Set(support(d))
    freq = Distributions.countmap(samples)
    @test isapprox(freq[f]/N, 0.7, atol=0.05)
    @test isapprox(freq[s]/N, 0.1, atol=0.05)
    @test isapprox(freq[q]/N, 0.2, atol=0.05)

    # test unnormalized case gives same answer:
    dd = UnivariateFinite(support(d), [70, 20, 10])
    samples = rand(StableRNG(123), dd, N);
    @test Set(samples) == Set(support(d))
    ffreq = Distributions.countmap(samples)
    @test isapprox(freq[f]/N, ffreq[f]/N)
    @test isapprox(freq[s]/N, ffreq[s]/N)
    @test isapprox(freq[q]/N, ffreq[q]/N)

    #
    # unordered (Multiclass):
    dict = Dict(S=>0.1, Q=> 0.2, F=> 0.7)
    d    = UnivariateFinite(dict)
    @test classes(d) == [a, f, q, s]
    @test classes(d) == classes(s)
    @test levels(d) == levels(s)
    @test support(d) == [f, q, s]
    @test MLJBase.sample_scitype(d) == Multiclass{4}
    # levels!(v, reverse(levels(v)))
    # @test classes(d) == [s, q, f, a]
    # @test support(d) == [s, q, f]

    @test pdf(d, S) ≈ 0.1
    @test pdf(d, 's') ≈ 0.1
    @test pdf(d, Q) ≈ 0.2
    @test pdf(d, 'q') ≈ 0.2
    @test pdf(d, F) ≈ 0.7
    @test pdf(d, 'f') ≈ 0.7
    @test pdf(d, A) ≈ 0.0
    @test logpdf(d, S) ≈ log(0.1)
    @test logpdf(d, 's') ≈ log(0.1)
    @test logpdf(d, Q) ≈ log(0.2)
    @test logpdf(d, 'q') ≈ log(0.2)
    @test logpdf(d, F) ≈ log(0.7)
    @test logpdf(d, 'f') ≈ log(0.7)
    @test isinf(logpdf(d, A))
    @test mode(d) == F

    @test UnivariateFinite(support(d), [0.7, 0.2, 0.1]) ≈ d

    N = 50
    rng = StableRNG(661)
    samples = [rand(rng,d) for i in 1:50];
    rng = StableRNG(661)
    @test samples == [rand(rng, d) for i in 1:N]

    N = 10000
    samples = rand(rng, d, N);
    @test Set(samples) == Set(support(d))
    freq = Distributions.countmap(samples)
    @test isapprox(freq[F]/N, 0.7, atol=0.05)
    @test isapprox(freq[S]/N, 0.1, atol=0.05)
    @test isapprox(freq[Q]/N, 0.2, atol=0.05)

    # corner case:
    d = UnivariateFinite(["ying", "yang"], 0.3, augment=true,
                         ordered=true, pool=missing)
    @test pdf(d, "yang") == 0.3
    @test classes(d)[1] == "ying"
    d = UnivariateFinite(classes(d), 0.3, augment=true)
    @test pdf(d, "yang") == 0.3

    # no support specified:
    @test_logs (:warn, r"No ") UnivariateFinite([0.7, 0.3])
    d = UnivariateFinite([0.7, 0.3], pool=missing)
    @test pdf(d, "class_1") == 0.7
end

@testset "constructor arguments not categorical values" begin
    @test_throws ArgumentError UnivariateFinite(Dict('f'=>0.7, 'q'=>0.2))
    @test_throws ArgumentError UnivariateFinite(Dict('f'=>0.7, 'q'=>0.2),
                                                pool=missing)
    @test_logs((:warn, r"Ignoring"),
               UnivariateFinite(Dict('f'=>0.7, 'q'=>0.3),
                                pool=f, ordered=true))

    d = UnivariateFinite(Dict('f'=>0.7, 'q'=>0.2, 's'=>0.1), pool=v)
    @test pdf(d, 'f') ≈ 0.7
    @test pdf(d, 's') ≈ 0.1
    @test pdf(d, 'q') ≈ 0.2
    @test_logs((:warn, r"No "),
               UnivariateFinite(['f', 'q', 's'],  [0.7, 0.2, 0.1]))
    # @test_throws(MethodError,
    #              UnivariateFinite(['f', 'q', 's'],  ["junk", 0.2, 0.1],
    #                               pool=missing))
    d = UnivariateFinite(['f', 'q', 's'],  [0.7, 0.2, 0.1], pool=missing)
    @test pdf(d, 'f') ≈ 0.7
    @test pdf(d, 's') ≈ 0.1
    @test pdf(d, 'q') ≈ 0.2

    d = UnivariateFinite(Dict('f'=>0.7, 'q'=>0.2, 's'=>0.1), pool=f)
    @test pdf(d, 'f') ≈ 0.7
    @test pdf(d, 's') ≈ 0.1
    @test pdf(d, 'q') ≈ 0.2

    d = UnivariateFinite(['f', 'q', 's'],  [0.7, 0.2, 0.1], pool=a)
    @test pdf(d, 'f') ≈ 0.7
    @test pdf(d, 's') ≈ 0.1
    @test pdf(d, 'q') ≈ 0.2

    d = UnivariateFinite(Dict('f'=>0.7, 'q'=>0.2, 's'=>0.1), pool=v)
    @test pdf(d, 'f') ≈ 0.7
    @test pdf(d, 's') ≈ 0.1
    @test pdf(d, 'q') ≈ 0.2

    d = UnivariateFinite(['f', 'q', 's'],  [0.7, 0.2, 0.1], pool=v)
    @test pdf(d, 'f') ≈ 0.7
    @test pdf(d, 's') ≈ 0.1
    @test pdf(d, 'q') ≈ 0.2
end

@testset "Univariate mode" begin
    v = categorical(1:101)
    p = rand(rng,101)
    s = sum(p) - p[42]
    p[42] = 0.5
    p = p/2/s
    p[42] = 0.5
    d = UnivariateFinite(v, p)
    @test mode(d) == 42
end

@testset "UnivariateFinite methods" begin
    y = categorical(["yes", "no", "yes", "yes", "maybe"])
    yes = y[1]
    no = y[2]
    maybe = y[end]
    prob_given_class = Dict(yes=>0.7, no=>0.3)
    d = UnivariateFinite(prob_given_class)
    @test pdf(d, yes) ≈ 0.7
    @test pdf(d, no) ≈ 0.3
    @test pdf(d, maybe) ≈ 0

    v = categorical(collect("abcd"))
    d = UnivariateFinite(v, [0.2, 0.3, 0.1, 0.4])
    sample = rand(rng,d, 10^4)
    freq_given_class = Distributions.countmap(sample)
    pairs = collect(freq_given_class)
    sort!(pairs, by=pair->pair[2], alg=QuickSort)
    sorted_classes = first.(pairs)
    @test sorted_classes == ['c', 'a', 'b', 'd']

    junk = categorical(['j',])
    j = junk[1]
    v = categorical(['a', 'b', 'a', 'b', 'c', 'b', 'a', 'a', 'f'])
    a = v[1]
    f = v[end]
    # remove f from sample:
    v = v[1 : end - 1]
    d = Distributions.fit(UnivariateFinite, v)
    @test pdf(d, a) ≈ 0.5
    @test pdf(d, 'a') ≈ 0.5
    @test pdf(d, 'b') ≈ 0.375
    @test pdf(d, 'c') ≈ 0.125
    @test pdf(d, 'f') == 0
    @test pdf(d, f) == 0
    @test_throws DomainError pdf(d, 'j')

    d2 = Distributions.fit(UnivariateFinite, v, nothing)
    @test d2 ≈ d

    # with weights:
    w = [2, 3, 2, 3, 5, 3, 2, 2]
    d = Distributions.fit(UnivariateFinite, v, w)
    @test pdf(d, a) ≈ 8/22
    @test pdf(d, 'a') ≈ 8/22
    @test pdf(d, 'b') ≈ 9/22
    @test pdf(d, 'c') ≈ 5/22
    @test pdf(d, 'f') == 0
    @test pdf(d, f) == 0
    @test_throws DomainError pdf(d, 'j')

    # regression test
    # https://github.com/alan-turing-institute/MLJBase.jl/pull/432/files#r502299301
    d = UnivariateFinite(classes(categorical(UInt32[0, 1])), [0.4, 0.6])
    @test pdf(d, UInt32(1)) == 0.6
end

@testset "approx for UnivariateFinite" begin
    y = categorical(["yes", "no", "maybe"])
    yes   = y[1]
    no    = y[2]
    maybe = y[3]
    @test(UnivariateFinite([yes, no, maybe], [0.1, 0.2, 0.7]) ≈
          UnivariateFinite([maybe, yes, no], [0.7, 0.1, 0.2]))
    @test(!(UnivariateFinite([yes, no, maybe], [0.1, 0.2, 0.7]) ≈
          UnivariateFinite([maybe, yes, no], [0.7, 0.2, 0.1])))
end

@testset "UnivariateFinite arithmetic" begin
    v = categorical(collect("abc"))
    a , b, c = v[1], v[2], v[3]

    d1  = UnivariateFinite([a, b], [0.2, 0.8])
    d2  = UnivariateFinite([b, c], [0.3, 0.7])
    dvec = [d1, d2]
    d = average(dvec)
    @test pdf(d, 'a') ≈ 0.1
    @test pdf(d, 'b') ≈ 0.55
    @test pdf(d, 'c') ≈ 0.35
    w = [4, 6]
    d = average(dvec, weights=w)
    @test pdf(d, 'a') ≈ 0.4*0.2
    @test pdf(d, 'b') ≈ 0.4*0.8 + 0.6*0.3
    @test pdf(d, 'c') ≈ 0.6*0.7
end

end # module

true


text:
################################################################################
#  Copyright 2020, Marta Vanin, Tom Van Acker                                  #
################################################################################
# PowerModelsDistributionStateEstimation.jl                                    #
# An extention package of PowerModels(Distribution).jl for Static Power System #
# State Estimation.                                                            #
################################################################################

"solves the AC state estimation in polar coordinates (ACP formulation)"
function solve_acp_mc_se(data::Union{Dict{String,<:Any},String}, solver; kwargs...)
    return solve_mc_se(data, _PMD.ACPUPowerModel, solver; kwargs...)
end

"solves the AC state estimation in rectangular coordinates (ACR formulation)"
function solve_acr_mc_se(data::Union{Dict{String,<:Any},String}, solver; kwargs...)
    return solve_mc_se(data, _PMD.ACRUPowerModel, solver; kwargs...)
end

"solves state estimation in current and voltage rectangular coordinates (IVR formulation)"
function solve_ivr_mc_se(data::Union{Dict{String,<:Any},String}, solver; kwargs...)
    return solve_mc_se(data, _PMD.IVRUPowerModel, solver; kwargs...)
end

"solves state estimation with a positive semi-definite fomrulation of the power flow equations (SDP formulation)"
function solve_sdp_mc_se(data::Union{Dict{String,<:Any},String}, solver; kwargs...)
    return solve_mc_se(data, _PMD.SDPUBFPowerModel, solver; kwargs...)
end

"solves the reduced AC state estimation in polar coordinates (ReducedACP formulation)"
function solve_acp_red_mc_se(data::Union{Dict{String,<:Any},String}, solver; kwargs...)
    return solve_mc_se(data, ReducedACPUPowerModel, solver; kwargs...)
end

"solves the reduced AC state estimation in rectangular coordinates (ReducedACR formulation)"
function solve_acr_red_mc_se(data::Union{Dict{String,<:Any},String}, solver; kwargs...)
    return solve_mc_se(data, ReducedACRUPowerModel, solver; kwargs...)
end

"solves the reduced state estimation in current and voltage rectangular coordinates (ReducedIVR formulation)"
function solve_ivr_red_mc_se(data::Union{Dict{String,<:Any},String}, solver; kwargs...)
    return solve_mc_se(data, ReducedIVRUPowerModel, solver; kwargs...)
end

"solves state estimation with a linear approximation of the power flow equations (LinDist3Flow formulation)"
function solve_linear_mc_se(data::Union{Dict{String,<:Any},String}, solver; kwargs...)
    return solve_mc_se(data, _PMD.LinDist3FlowPowerModel, solver; kwargs...)
end

function solve_mc_se(data::Union{Dict{String,<:Any},String}, model_type::Type, solver; kwargs...)
    if haskey(data["se_settings"], "criterion")
        _PMDSE.assign_unique_individual_criterion!(data)
    end
    if !haskey(data["se_settings"], "rescaler")
        data["se_settings"]["rescaler"] = 1
        @warn "Rescaler set to default value, edit data dictionary if you wish to change it."
    end
    if !haskey(data["se_settings"], "number_of_gaussian")
        data["se_settings"]["number_of_gaussian"] = 10
        @warn "Estimation criterion set to default value, edit data dictionary if you wish to change it."
    end
    return _PMD.solve_mc_model(data, model_type, solver, build_mc_se; kwargs...)
end

"specification of the state estimation problem for a bus injection model - ACP and ACR formulations"
function build_mc_se(pm::_PMD.AbstractUnbalancedPowerModel)

    # Variables
    _PMDSE.variable_mc_bus_voltage(pm; bounded = true)
    _PMD.variable_mc_branch_power(pm; bounded = true)
    _PMD.variable_mc_transformer_power(pm; bounded = true)
    _PMD.variable_mc_generator_power(pm; bounded = true)
    variable_mc_load(pm; report = true)
    variable_mc_residual(pm; bounded = true)
    variable_mc_measurement(pm; bounded = false)

    # Constraints
    for (i,gen) in _PMD.ref(pm, :gen)
        _PMD.constraint_mc_generator_power(pm, i)
    end
    for (i,bus) in _PMD.ref(pm, :ref_buses)
        @assert bus["bus_type"] == 3
        _PMD.constraint_mc_theta_ref(pm, i)
    end
    for (i,bus) in _PMD.ref(pm, :bus)
        _PMDSE.constraint_mc_power_balance_se(pm, i)
    end
    for (i,branch) in _PMD.ref(pm, :branch)
        _PMD.constraint_mc_ohms_yt_from(pm, i)
        _PMD.constraint_mc_ohms_yt_to(pm,i)
    end
    for (i,meas) in _PMD.ref(pm, :meas)
        constraint_mc_residual(pm, i)
    end

    for i in _PMD.ids(pm, :transformer)
        _PMD.constraint_mc_transformer_power(pm, i)
    end

    # Objective
    objective_mc_se(pm)

end

"specification of the state estimation problem for the IVR Flow formulation"
function build_mc_se(pm::_PMD.AbstractUnbalancedIVRModel)
    # Variables

    _PMD.variable_mc_bus_voltage(pm, bounded = true)
    _PMDSE.variable_mc_branch_current(pm, bounded = true)
    variable_mc_generator_current_se(pm, bounded = true)
    _PMD.variable_mc_transformer_current(pm, bounded = false)
    variable_mc_load_current(pm, bounded = false)#TODO bug in the bounds assignment
    variable_mc_residual(pm, bounded = true)
    variable_mc_measurement(pm, bounded = false)

    # Constraints
    for (i,bus) in _PMD.ref(pm, :ref_buses)
        @assert bus["bus_type"] == 3
        _PMD.constraint_mc_theta_ref(pm, i)
    end

    # gens should be constrained before KCL, or Pd/Qd undefined
    for id in _PMD.ids(pm, :gen)
        constraint_mc_generator_power_se(pm, id)
    end

    for (i,bus) in _PMD.ref(pm, :bus)
        constraint_mc_current_balance_se(pm, i)
    end

    if typeof(pm) <: ReducedIVRUPowerModel
        for i in _PMD.ids(pm, :branch)
            constraint_current_to_from(pm, i)
            constraint_mc_bus_voltage_drop(pm, i)
        end
    else
        for i in _PMD.ids(pm, :branch)
            _PMD.constraint_mc_current_from(pm, i)
            _PMD.constraint_mc_current_to(pm, i)
            _PMD.constraint_mc_bus_voltage_drop(pm, i)
        end
    end

    for i in _PMD.ids(pm, :transformer)
        _PMD.constraint_mc_transformer_power(pm, i)
    end

    for (i,meas) in _PMD.ref(pm, :meas)
        constraint_mc_residual(pm, i)
    end

    objective_mc_se(pm)

end

"specification of the state estimation problem for a branch flow model - SDP and LinDist3Flow formulations"
function build_mc_se(pm::_PMD.AbstractUBFModels)

    # Variables
    _PMD.variable_mc_bus_voltage(pm)
    _PMD.variable_mc_branch_current(pm)
    _PMD.variable_mc_branch_power(pm)
    _PMD.variable_mc_transformer_power(pm; bounded=true)
    _PMD.variable_mc_generator_power(pm; bounded=true)
    variable_mc_load(pm; report = true)
    variable_mc_residual(pm, bounded = true)
    variable_mc_measurement(pm, bounded = false)

    #Constraints
    _PMD.constraint_mc_model_current(pm)

    for (i,bus) in _PMD.ref(pm, :ref_buses)
        @assert bus["bus_type"] == 3
        _PMD.constraint_mc_theta_ref(pm, i)
    end

    for id in _PMD.ids(pm, :gen)
        _PMD.constraint_mc_generator_power(pm, id)
    end

    for (i,bus) in _PMD.ref(pm, :bus)
        _PMDSE.constraint_mc_power_balance_se(pm,i)
    end

    for i in _PMD.ids(pm, :branch)
        _PMD.constraint_mc_power_losses(pm, i)
        _PMD.constraint_mc_model_voltage_magnitude_difference(pm, i)
        _PMD.constraint_mc_voltage_angle_difference(pm, i)
    end

    for i in _PMD.ids(pm, :transformer)
        _PMD.constraint_mc_transformer_power(pm, i)
    end

    for (i,meas) in _PMD.ref(pm, :meas)
        constraint_mc_residual(pm, i)
    end

    objective_mc_se(pm)

end


text:
using Documenter, CardioModels

makedocs(
    modules = [CardioModels],
    format = Documenter.HTML(; prettyurls = get(ENV, "CI", nothing) == "true"),
    authors = "Marius Pille",
    sitename = "CardioModels.jl",
    pages = [
        "General" => "index.md",
        "DeBoer" => "DeBoer.md",
        "Karemaker" => "Karemaker.md",
        "Baselli" => "Baselli.md",
    ]
    # strict = true,
    # clean = true,
    # checkdocs = :exports,
)

deploydocs(
    repo = "github.com/mapi1/CardioModels.jl.git",
    push_preview = true
)


text:
mutable struct Chord
    notes::Vector{Interval}
end

const MAJOR_CHORD = Chord([ROOT, MAJOR_THIRD, PERFECT_FIFTH])
const MINOR_CHORD = Chord([ROOT, MINOR_THIRD, PERFECT_FIFTH])
const POWER_CHORD = Chord([ROOT, PERFECT_FIFTH, OCTAVE])


text:
#= Atom
Implements static and mutable types for representing atoms and associated atom features.

Depreciated: StaticAtom implements properties that do not change depending on atomic context, but OxygenAtom contains properties that do.
Current: OxygenAtom is a static type, and other variable properties are always calculated with respect to their environment.

=#

abstract type Atomic end

struct OxygenAtom <: Atomic
    number::Int
    symbol::String
    mass::Float64
    orbitals::Electrons
    electronegativity::Float64
end

"""
    OxygenAtom

Immutable atom type.

# Examples
```
julia> OxygenAtom("C")
"C" atom (number 6)

julia> OxygenAtom(6)
"C" atom (number 6)
```
"""
OxygenAtom

OxygenAtom(symbol::String) = begin
    number = periodic_table[symbol]
    OxygenAtom(number, symbol, periodic_table_mass[number], Electrons(number), periodic_table_electronegativity[number])
end

OxygenAtom(number::Int) = begin
    symbol = periodic_table_reverse[number]
    OxygenAtom(number, symbol, periodic_table_mass[number], Electrons(number), periodic_table_electronegativity[number])
end

Base.:(==)(a::Atomic, b::Atomic) = isequal(a.number, b.number) && isequal(a.symbol, b.symbol) && isequal(a.mass, b.mass) && isequal(a.orbitals, b.orbitals)

Base.show(io::IO, a::Atomic) = print(io, "\"", a.symbol, "\" atom (number ", a.number, ")")

# Base.show(io::IO, ::MIME"text/plain", a::Atomic) =
#            print(io, "$(typeof(a)):\n   ", a)


text:
function emax_outneighbors(g::AbstractGraph{T}, v::S, w::AbstractMatrix{U}) where {S,T <: Integer} where {U <: Real}
    nb = outneighbors(g, v)
    wnb = w[v, nb]
    length(nb) == 0 && return nb
    nb[wnb .== maximum(wnb)]
end

emax_outneighbors(g::AbstractMetaGraph, v::Integer) = emax_outneighbors(g, v, weights(g))

function emax_inneighbors(g::AbstractGraph{T}, v::S, w::AbstractMatrix{U}) where {S,T <: Integer} where {U <: Real}
    nb = inneighbors(g, v)
    wnb = w[nb, v]
    length(nb) == 0 && return nb
    nb[wnb .== maximum(wnb)]
end

emax_inneighbors(g::AbstractMetaGraph, v::Integer) = emax_inneighbors(g, v, weights(g))

function max_outweight(g::AbstractGraph{T}, v::S, w::AbstractMatrix{U}) where {S,T <: Integer} where {U <: Real}
    nb = outneighbors(g, v)
    length(nb) > 0 || return 0
    maximum(w[v, nb])
end

function max_inweight(g::AbstractGraph{T}, v::S, w::AbstractMatrix{U}) where {S,T <: Integer} where {U <: Real}
    nb = inneighbors(g, v)
    length(nb) > 0 || return 0
    maximum(w[nb, v])
end

abstract type MainPathAlgorithm <: AbstractGraphAlgorithm end
abstract type MainPathResult <: AbstractGraphResult end

struct StandardGlobal <: MainPathAlgorithm end

struct StandardGlobalResult{T<:Integer, U<:Real} <: MainPathResult
    edges::Vector{SimpleEdge{T}}
    weights::Vector{U}
    totalweight::U
end

function mainpath(
    g::AbstractGraph{T},
    weights::AbstractMatrix{U},
    ::StandardGlobal) where T <: Integer where U <: Real

    add_source_target!(g)
    vt = topological_sort(g)
    source = vt[1]; sink = vt[end]
    n = nv(g)
    lpath = zeros(U, n)
    lpath[vt[1]] = 1.0

    @inbounds for v in vt
        (v == n-1 || v == n-2) && continue
        nb = inneighbors(g, v)
        if length(nb) == 0
            lpath[v] = 0
        else
            val = lpath[nb]
            lpath[v] = maximum(val) + weights[nb[findmax(val)[2]], v]
        end
    end

    # This is adapted from LightGraphs breadthfirst

    start = sink
    visited = falses(n)
    cur_level = Vector{T}()
    sizehint!(cur_level, n)
    next_level = Vector{T}()
    sizehint!(next_level, n)
    edgesmp = Vector{SimpleEdge{T}}()
    push!(cur_level, start)

    while !isempty(cur_level)
        @inbounds for v in cur_level
            nb = inneighbors(g, v)
            length(nb) == 0 && continue
            val = lpath[nb]
            maxnb = nb[val .== maximum(val)]
            for u in maxnb
                if !visited[u]
                    (v != sink && u != source) && push!(edgesmp, SimpleEdge(u, v))
                    push!(next_level, u)
                    visited[u] = true
                else
                    (v != sink && u != source) && push!(edgesmp, SimpleEdge(u, v))
                end
            end
        end
        empty!(cur_level)
        cur_level, next_level = next_level, cur_level
        sort!(cur_level, alg=QuickSort)
    end

    for s in [sink, source] rem_vertex!(g, s) end
    w = [weights[src(e), dst(e)] for e in edgesmp]
    return StandardGlobalResult(edgesmp, w, sum(w))
end

mainpath(g::AbstractGraph{T}, ::StandardGlobal) where {T <: Integer} = mainpath(g, weights(g), StandardGlobal())


text:
include("./options.jl")
using ArgParse

function parse_commandline()

	settings = ArgParseSettings()

	@add_arg_table! settings begin
		"file"
			help = "audio file to be classified or dataset config file"
			required = false

		"--shufflesamples", "-s"
			help = "shuffles the samples in the dataset and save the result in " *
						"the csv file (do not use with previously started training " *
						"or the test set will be tainted with training and " *
						"validation samples)"
			action = :store_true

		"--trainnetwork", "-t"
			help = "train a network with the dataset described by the csv file"
			action = :store_true

		"--preprocess", "-p"
			help = "extract features from mp3 audio files to be used during the " *
					"network training, overwriting existing extracted feature " *
					"files (note: can be used with the -t flag to start training " *
					"after pre-processing)"
			action = :store_true

		"--selectedfeature", "-f"
			help = "define the features that will be extracted from mp3 audio " *
					"files and/or that will be used during the network training. " *
					"Can be set to 'chroma', 'melspectrogram' or 'mfcc'"
			arg_type = String
			default = "melspectrogram"

		"--audiofeaturesdir", "-F"
			help = "the directory to save extracted features after preprocessing " *
					"and to read features before training, will be appended to the " *
					"paths found in the CSV file to generate a filepath for each " *
					"item in the dataset (default: \"./preprocessed_audio/SELECTEDFEATURE\")"
			arg_type = String

		"--model", "-m"
			help = "the trained model to be used to classify a music, or to be " *
					"trained"
			arg_type = String

		"--modeloutputdir", "-M"
			help = "the directory to save model parameters learned during the " *
					"training process"
			arg_type = String
			default = "./models"

		"--dumpoutputdir", "-D"
			help = "the directory to save dump files of loss/accuracy data"
			arg_type = String
			default = "./data_dump"

		"--epochs", "-e"
			help = "the number of training epochs"
			arg_type = Int
			default = 20

		"--batchsize", "-b"
			help = "the number of samples per batch"
			arg_type = Int
			default = 5

		"--eta", "-n"
			help = "the initial learning rate to be used when training"
			arg_type = Float64
			default = 1E-3

	end

	args = parse_args(settings)

	if args["audiofeaturesdir"] === nothing
		args["audiofeaturesdir"] = "extracted_features/$(args["selectedfeature"])"
	end

	if args["selectedfeature"] ∉ ["chroma", "melspectrogram", "mfcc"]
		println("The SELECTEDFEATURE must be set to 'chroma', 'melspectrogram' or 'mfcc'.")
		println(usage_string(settings))
		exit(1)
	end

	return args

end

merge!(options, parse_commandline())

max_stars_repo_name: gabriel-ss/music-autotagger

max_issues_repo_name: gabriel-ss/music-autotagger

max_forks_repo_name: gabriel-ss/music-autotagger


text:
pprint(it...;sep=" ",ed="\n")=length(it)>0 ? (@inbounds for i in 1:length(it)-1;print(it[i],sep);end;print(last(it),ed)) : print(ed)

k,n=parse.(Int64, split(readline()))
a=parse.(Int64, split(readline()))
mx=0
for i in 1:n-1
    global mx
    mx=max(a[i+1]-a[i],mx)
end
mx=max(mx,(k-a[n])+a[1])
pprint(k-mx)

max_stars_repo_name: sugitanishi/competitive-programming

max_issues_repo_name: sugitanishi/competitive-programming

max_forks_repo_name: sugitanishi/competitive-programming


text:
#=
Created on 09/01/2021 12:54:49
Last update: -

@author: Michiel Stock
michielfmstock@gmail.com

Estimate pi through Monte Carlo sampling.
Do this by simuting to throw `n` pebbles in the [-1, 1] x [-1, 1] square
and track the fraction that land in the unit circle.
=#

function estimatepi(n)
    hits = 0
    for i in 1:n
        x = rand()
        y = rand()
        if x^2 + y^2 ≤ 1.0
            hits += 1
        end
    end
    return 4hits/n
end

estimatepi2(n) = 4count(x-> x ≤ 1.0, sum(rand(n, 2).^2, dims=2)) / n


text:
module MinkReduction

using LinearAlgebra, Printf
export GaussReduce, RandUnimodMat, RandLowerTri, minkReduce, DeviousMat, isMinkReduced

"""
    minkReduce(U, V, W, debug=false)

Find the shortest equivalent basis of that lattice formed by {`U`, `V`, `W`}

```jldoctest
julia> U = [1, 2, 3]; V = [-1, 2, 3]; W = [3, 0, 4]; minkReduce(U,V,W)
[1.0 2.0 -1.0; 2.0 -2.0 2.0; 3.0 1.0 3.0]
([1.0, 2.0, 3.0], [2.0, -2.0, 1.0], [-2.0, 0.0, 0.0])
```
"""
function minkReduce(U, V, W; debug=false)
    i = 0
    origMat = hcat(U,V,W)
    println("Input basis vectors:")
    println("U: ",U, "  norm: ",round(norm(U),digits=2))
    println("V: ",V, "  norm: ",round(norm(V),digits=2))
    println("W: ",W, "  norm: ",round(norm(W),digits=2))
    while true
        i +=1
        println(" --< iteration $i >--")
        norms = [norm(U), norm(V), norm(W)]
        println("sort vectors into ascending order by ℓ-2 norms")
        p = sortperm(norms)
        U,V,W = (U,V,W)[p] # sort into ascending order
        U,V,W = shortenW_in_UVW(U, V, W)
        #println("integer transform? ", inv(origMat)*hcat(U,V,W))
        i > 10 && error("minkReduce: Too many iterations")
        norm(W) ≥ norm(V) && break
    end
    println("final vectors:")
    if !debug return U, V, W end
    return U, V, W, i
end

"""
    shortenW_in_UVW

Reduce vector W so that it is as close to the origin as possible.

Gauss reduce U and V. Subtract multiples of new U, V from W to bring its
U-V projection as close as possible to the origin. The projection of the new
W will be contained in the parallelogram formed by U, V.
(See Lecture notes in computer science, ISSN 0302-974, ANTS - VI : algorithmic
number theory, 2004, vol. 3076, pp. 338-357 ISBN 3-540-22156-5)
"""
function shortenW_in_UVW(U,V,W)
    U, V = GaussReduce(U,V)
    println("Next: Shorten W in U-V-W")
    println("U:",U)
    println("V:",V)
    println("W:",W)
    @printf("norms: U, V, W:\n %7.3f\n %7.3f\n %7.3f\n ", norm(U), norm(V), norm(W))
    # n̂ is a unit vector ⟂ to U-V plane. Subtract W⋅n̂ multiples of n̂
    # from W to get W's projection, T, into the U-V plane.
    n̂ = (U×V)/norm(U×V)
    T = W - W⋅n̂ * n̂ # Get the projection of W into U-V plane
    𝕄 = hcat(U,V,W) # Matrix with U, V, W as columns
    # Get the number of multiples of U, V needed to move T inside the U-V cell
    println("unrounded/unfloored multiples of U,V,W: ", inv(𝕄)*T)
    latCoords = floor.(Int,inv(𝕄)*T)
     println("U,V multiples needed to shift W into first cell: ",latCoords)
     println("Shift in Cartesian coords: ", 𝕄*latCoords)
    # Find the corner of the U-V cell closest to shifted T and shift T to be closest to origin
    Wnew = W - 𝕄*latCoords # Try the corner at the origin first
        # Now try the other three corners, keep the shortest that is found
    norm(Wnew) > norm(W - U) && (Wnew = W - U)
    norm(Wnew) > norm(W - V) && (Wnew = W - V)
    norm(Wnew) > norm(W - U - V) && (Wnew = W - U -V)
    W = Wnew
    println("shortened W:",W)
    @printf("norms: U, V, Wnew:\n %7.3f\n %7.3f\n %7.3f\n ", norm(U), norm(V), norm(Wnew))
    return U, V, Wnew
end

"""
    GaussReduce(U, V)

Reduce the basis vectors {`U`, `V`} to the shortest possible basis.

# Examples
```julia-repl
julia> GaussReduce([5 8], [8 13])
([-1.0 0.0], [0.0 -1.0])
```
"""
function GaussReduce(U, V)
    println("GaussReduce begins")
    #maxval = max(abs.(U)...,abs.(V)...)
    if norm(U) > norm(V) U, V = V, U end
    i = 0
    while true
        println("U: ",U,"  V: ",V)
        V, U = U, V - round(Int,(U⋅V)/(U⋅U))*U
        i += 1
        if norm(U) > norm(V) || norm(U)≈norm(V) break; end
        i > 50 && error("GaussReduce: Too many iterations") # failsafe to break out if not converging
    end
    println("U: ",U,"  V: ",V)
    println("GaussReduce complete")
    return V, U
end

"""
    RandUnimodMat(n)

Generate a random unimodular 2x2 matrix. `n` is a small integer (number of row and column operations).

See also: `RandLowerTri(n)`, `FibonacciMat(n)`, `DeviousMat(n)`

# Examples
```jldoctest
julia> RandUnimodMat(5)
2×2 Array{Int64,2}:
 3   8
 7  19
```
"""
function RandUnimodMat(n)
    mat = RandLowerTri(1)
    for i ∈ 1:n
        mat = mat*RandLowerTri(1)
        mat = mat*transpose(RandLowerTri(1))
    end
    return mat
end

"""
    RandLowerTri(n)
Generate a random 2x2 matrix of the form [1 0; 0 ±n].

See also: `RandUnimodMat(n)`, `FibonacciMat(n)`, `DeviousMat(n)`

# Examples
```jldoctest
julia> RandLowerTri(4)
2×2 Array{Int64,2}:
  1  0
 -3  1
```
"""
function RandLowerTri(n)
    return [1 0; rand(-n:n) 1]
end

"""
    FibonacciMat(k)
Generate a 2x2 matrix of the form [f2 f3; f1 f2] where f1, f2, f3 consecutive Fibonacci-like

See also: `RandUnimodMat(n)`, `FibonacciMat(n)`, `DeviousMat(n)`

"""
function FibonacciMat(k)
    f1 = round(Int32,1.61803398875^k/sqrt(5))
    f2 = round(Int32,1.61803398875^(k+1)/sqrt(5))
    f3 = f1 + f2
    any(i -> i < 1, [f1 f2 f3]) && error("Overflow in FibonacciMat function")
    return [f2 f3; f1 f2]
end

"""
    DeviousMat(n)

Make a unimodular 3x3 matrix that requires a large number of steps to reduce
(See email from Rod Feb 1 2020)

`n` dictates the size of the entries
"""
function DeviousMat(n)
    n < 3 && error("for DeviousMat, n ≥ 3")
    u,v = round(Int64,(2+√3)^n/(2*√3)), round(Int64,(2+√3)^n/2)
    a,b = convert(Int64,(u+v+1)/2), -u
    c,d = a-1, v-u
    return [a;b;c], [b;d;b], [c;b;a]
end

"""
    isMinkReduced(U,V,W)

Check if the basis {`U`,`V`,`W`} is Minkoswki reduced.

"""
function isMinkReduced(U,V,W)
    if norm(U) > norm(V)+eps()     println("Condition 1  failed"); return false end
    if norm(V) > norm(W)+eps()     println("Condition 2  failed"); return false end
    if norm(V) > norm(U+V)+eps()   println("Condition 3  failed"); return false end
    if norm(V) > norm(U-V)+eps()   println("Condition 4  failed"); return false end
    if norm(W) > norm(U+W)+eps()   println("Condition 5  failed"); return false end
    if norm(W) > norm(U-W)+eps()   println("Condition 6  failed"); return false end
    if norm(W) > norm(V+W)+eps()   println("Condition 7  failed"); return false end
    if norm(W) > norm(V-W)+eps()   println("Condition 8  failed"); return false end
    if norm(W) > norm(U+V+W)+eps() println("Condition 9  failed"); return false end
    if norm(W) > norm(U-V+W)+eps() println("Condition 10 failed"); return false end
    if norm(W) > norm(U+V-W)+eps() println("Condition 11 failed"); return false end
    if norm(W) > norm(U-V-W)+eps() println("Condition 12 failed"); return false end
    return true
end

end

max_stars_repo_name: glwhart/JuSym

max_issues_repo_name: glwhart/JuSym

max_forks_repo_name: glwhart/JuSym


text:
######################################################################################################
# Copy from https://bitbucket.org/wmkoolen/tidnabbil/src/master/purex_games_paper/ by Wouter Koolen  #
######################################################################################################
using Plots;
using Distributed;
using Printf;
@everywhere include("runit.jl");
include("experiment_helpers.jl");

# exercise all code paths
function test_all()
    δ = 10^-2; # confidence
    problems = (((μs, pep)
        for dist in (Gaussian(), Bernoulli(), Exponential())
        for (μs, pep) in (([.5, .6, .7, .8, .9],    MinimumThreshold(dist, .08)),
                          ([.5, .6, .7, .8, .9],    BestArm(dist)),
                          ([.5, .6, .5, .6, .55],   SameSign(dist, .45)),
                          ([.1 .2 .3; .6 .5 .4][:], LargestProfit(dist))
                          ))
                ...,
                ([.4, .4], HeteroSkedastic((Gaussian(1), Gaussian(2)), .5))
                );
    configs = ((μs, sr, pep)
               for (μs, pep) in problems
               for sr in everybody(pep, μs)
               );
    β = (t) -> -log(δ); # gross cheat here, no log(log(1/\δ)) or log(log(t)) terms
    for (μs, sr, pep) in configs
        Tstar, _ = oracle(pep, μs);
        ((irec, N, time),) = runit(1, sr, μs, pep, (β,));
        ⋆ = istar(pep, μs);
        println("problem $pep\nsampling $sr\ninstance $μs\nTstar $(Tstar == nothing ? "nope" : Tstar*d(Bernoulli(), δ, 1-δ))  τ $(sum(N))\n⋆ $(⋆)  rec $irec\ntime $(time/1000000.0)\n");
    end
end

function launch()
    δ = 10^-4; # confidence
    dist = Gaussian(1);
    case = 1;
    if case == 1
        μs = [.3, .6]; # parameters
        if false
            γ = .08;
        else
            γ = .6;
        end
        pep = MinimumThreshold(dist, γ);
    elseif case == 2
        μs = [.5, .6, .7, .8, .9]; # parameters
        pep = BestArm(dist);
    elseif case == 3
        μs = [.5, .6, .5, .6, .55]; # parameters
        pep = SameSign(dist, .45);
    elseif case == 4
        μs = [.1 .2 .3; .6 .5 .4][:]
        pep = LargestProfit(dist);
    elseif case == 5
        μs = [.5, .5]
        pep = HeteroSkedastic((Gaussian(1), Gaussian(5)), .3);
    else
        @assert false "unknown case $case"
    end

    K = length(μs)
    Tstar, wstar = oracle(pep, μs);
    β = (t) -> -log(δ); # gross cheat here, no log(K), log(log(1/\δ)) or log(log(t)) terms
    # correct answer
    ⋆ = istar(pep, μs);

    nreps = 1000;
    @time runs = pmap(seed -> begin
                      ((irec, N, H, A),) = runit(seed, DaBomb(DTracking, nanswers(pep, μs)), μs, pep, (β,));
                      if irec != ⋆
                      println("error $irec != $(⋆) at τ $(sum(N))  N $N");
                      end
                      irec, N
                      end,
                      1:nreps);

    if Tstar != nothing
        lbd = Tstar*d(Bernoulli(), δ, 1-δ);
        println("w⋆ = $(wstar)");
        println("lbd = $lbd");
    end

    println("E[N] = ", [sum(x->x[2][k], runs)/nreps for k in 1:K])

    Eτ = sum(sum(x[2] for x in runs))/nreps;
    perr = sum(first.(runs) .!= ⋆)/length(runs);

    if perr > δ
        @warn "too many errors: $perr > $δ";
    end

    println("E[τ] ≈ $Eτ");
    println("P(err) ≈ $(perr)");

    # histogram of run times
    stephist([sum(x[2]) for x in runs], normalize=:pdf, bins = 50, label="tau");
    if Tstar != nothing
        vline!([lbd], label="lower bd");
    end
    vline!([Eτ], label="emp. mean");

    gui();
    runs
end


text:
#=

The four adjacent digits in the 1000-digit number that have the greatest product are 9 × 9 × 8 × 9 = 5832.

73167176531330624919225119674426574742355349194934
96983520312774506326239578318016984801869478851843
85861560789112949495459501737958331952853208805511
12540698747158523863050715693290963295227443043557
66896648950445244523161731856403098711121722383113
62229893423380308135336276614282806444486645238749
30358907296290491560440772390713810515859307960866
70172427121883998797908792274921901699720888093776
65727333001053367881220235421809751254540594752243
52584907711670556013604839586446706324415722155397
53697817977846174064955149290862569321978468622482
83972241375657056057490261407972968652414535100474
82166370484403199890008895243450658541227588666881
16427171479924442928230863465674813919123162824586
17866458359124566529476545682848912883142607690042
24219022671055626321111109370544217506941658960408
07198403850962455444362981230987879927244284909188
84580156166097919133875499200524063689912560717606
05886116467109405077541002256983155200055935729725
71636269561882670428252483600823257530420752963450

Find the thirteen adjacent digits in the 1000-digit number that have the greatest product. What is the value of this product?

=#

str ="""73167176531330624919225119674426574742355349194934
        96983520312774506326239578318016984801869478851843
        85861560789112949495459501737958331952853208805511
        12540698747158523863050715693290963295227443043557
        66896648950445244523161731856403098711121722383113
        62229893423380308135336276614282806444486645238749
        30358907296290491560440772390713810515859307960866
        70172427121883998797908792274921901699720888093776
        65727333001053367881220235421809751254540594752243
        52584907711670556013604839586446706324415722155397
        53697817977846174064955149290862569321978468622482
        83972241375657056057490261407972968652414535100474
        82166370484403199890008895243450658541227588666881
        16427171479924442928230863465674813919123162824586
        17866458359124566529476545682848912883142607690042
        24219022671055626321111109370544217506941658960408
        07198403850962455444362981230987879927244284909188
        84580156166097919133875499200524063689912560717606
        05886116467109405077541002256983155200055935729725
        71636269561882670428252483600823257530420752963450"""
long_string = replace(str, "\n" => "")
len_adjacent = 13

function largest_adjacent_product(string, len_adjacent)
    num_digits = length(string)
    largest_product = 0
    largest_piece = 0
    for offset in 0:num_digits-len_adjacent
        number_piece = string[1+offset:len_adjacent+offset]
        #println(number_piece)
        product = prod(digits(parse(Int, number_piece)))
        if product > largest_product
            largest_product = product
            largest_piece = number_piece
        end

    end
    return largest_product, largest_piece
end

println()
print(largest_adjacent_product(long_string, len_adjacent))

max_stars_repo_name: pedvide/project-euler

max_issues_repo_name: pedvide/project-euler

max_forks_repo_name: pedvide/project-euler


text:
struct Falsi <: AbstractBracketingAlgorithm
end

function alg_cache(alg::Falsi, left, right, p, ::Val{true})
  nothing
end

function alg_cache(alg::Falsi, left, right, p, ::Val{false})
  nothing
end

function perform_step(solver, alg::Falsi, cache)
  @unpack f, p, left, right, fl, fr = solver

  fzero = zero(fl)
  fl * fr > fzero && error("Bracket became non-containing in between iterations. This could mean that "
  + "input function crosses the x axis multiple times. Bisection is not the right method to solve this.")

  mid = (fr * left - fl * right) / (fr - fl)

  if right == mid || right == mid
    @set! solver.force_stop = true
    @set! solver.retcode = FLOATING_POINT_LIMIT
    return solver
  end

  fm = f(mid, p)

  if iszero(fm)
    # todo: phase 2 bisection similar to the raw method
    @set! solver.force_stop = true
    @set! solver.left = mid
    @set! solver.fl = fm
    @set! solver.retcode = EXACT_SOLUTION_LEFT
  else
    if sign(fm) == sign(fl)
      @set! solver.left = mid
      @set! solver.fl = fm
    else
      @set! solver.right = mid
      @set! solver.fr = fm
    end
  end
  return solver
end


text:
 """
```
pseudo_measurement(m::Model1010_1spread{T}, TTT::Matrix{T}, RRR::Matrix{T},
                   CCC::Vector{T}) where {T<:AbstractFloat}
```

Assign pseudo-measurement equation (a linear combination of states):

```
x_t = ZZ_pseudo*s_t + DD_pseudo
```
"""
function pseudo_measurement(m::Model1010_1spread{T},
                            TTT::Matrix{T},
                            RRR::Matrix{T},
                            CCC::Vector{T}) where {T<:AbstractFloat}
    endo      = m.endogenous_states
    endo_addl = m.endogenous_states_augmented
    pseudo    = m.pseudo_observables

    _n_states = n_states_augmented(m)
    _n_pseudo = n_pseudo_observables(m)

    # Compute TTT^10, used for ExpectedAvg10YearRateGap, ExpectedAvg10YearRate, and ExpectedAvg10YearNaturalRate
    TTT5  = (1/20)*((UniformScaling(1.) - TTT)\(UniformScaling(1.) - TTT^20))
    TTT10 = (1/40)*((UniformScaling(1.) - TTT)\(UniformScaling(1.) - TTT^40))
    TTT20 = (1/80)*((UniformScaling(1.) - TTT)\(UniformScaling(1.) - TTT^80))

	# Initialize pseudo ZZ and DD matrices
    ZZ_pseudo = zeros(_n_pseudo, _n_states)
    DD_pseudo = zeros(_n_pseudo)

	##########################################################
	## PSEUDO-OBSERVABLE EQUATIONS
	##########################################################

	## Output
	ZZ_pseudo[pseudo[:y_t],endo[:y_t]] = 1.

	## Flexible Output
	ZZ_pseudo[pseudo[:y_f_t],endo[:y_f_t]] = 1.

	## Output Gap
	ZZ_pseudo[pseudo[:OutputGap],endo[:y_t]] = 1
	ZZ_pseudo[pseudo[:OutputGap],endo[:y_f_t]] = -1

    if subspec(m) in ["ss2", "ss4", "ss22"]
	    ## π_t
	    ZZ_pseudo[pseudo[:π_t],endo[:π_t]] = 1.
	    DD_pseudo[pseudo[:π_t]]            = 100*(m[:π_star]-1)

	    ## Long Run Inflation
	    ZZ_pseudo[pseudo[:LongRunInflation],endo[:π_star_t]] = 1.
	    DD_pseudo[pseudo[:LongRunInflation]]                 = 100.0*(m[:π_star]-1.)

        ## Marginal Cost
        ZZ_pseudo[pseudo[:MarginalCost],endo[:mc_t]] = 1.

	    ## Wages
	    ZZ_pseudo[pseudo[:Wages],endo[:w_t]] = 1.

	    ## Flexible Wages
	    ZZ_pseudo[pseudo[:FlexibleWages],endo[:w_f_t]] = 1.

	    ## z_t
	    ZZ_pseudo[pseudo[:z_t], endo[:z_t]] = 1.
    end

	## Hours
	ZZ_pseudo[pseudo[:Hours],endo[:L_t]] = 1.

    if subspec(m) in ["ss2", "ss4", "ss22"]
	    ## Flexible Hours
	    ZZ_pseudo[pseudo[:FlexibleHours],endo[:L_f_t]] = 1.
    end

	## Natural Rate
	ZZ_pseudo[pseudo[:RealNaturalRate],endo[:r_f_t]] = 1.
	DD_pseudo[pseudo[:RealNaturalRate]]              = 100. * (m[:rstar]-1.)

	## Ex Ante Real Rate
	ZZ_pseudo[pseudo[:ExAnteRealRate],endo[:R_t]]  = 1
	ZZ_pseudo[pseudo[:ExAnteRealRate],endo[:Eπ_t]] = -1
	DD_pseudo[pseudo[:ExAnteRealRate]]             = m[:Rstarn] - 100. * (m[:π_star]-1.)

	## Nominal FFR
	ZZ_pseudo[pseudo[:NominalFFR], endo[:R_t]] = 1.
	DD_pseudo[pseudo[:NominalFFR]] = m[:Rstarn]

	## Ex Ante Real Deposit Rate
	ZZ_pseudo[pseudo[:RealDepositRate],endo[:Rd_t]] = 1
	ZZ_pseudo[pseudo[:RealDepositRate],endo[:Eπ_t]] = -1
	DD_pseudo[pseudo[:RealDepositRate]]             = m[:Rstard] - 100. * (m[:π_star]-1.)

	## Nominal Deposit Rate
	ZZ_pseudo[pseudo[:NominalDepRate],endo[:Rd_t]]  = 1
	DD_pseudo[pseudo[:NominalDepRate]]             = m[:Rstard]

	## Output: GDP
	ZZ_pseudo[pseudo[:gdp], endo[:y_t]]          = 1.0
    ZZ_pseudo[pseudo[:gdp], endo_addl[:y_t1]]    = -1.0
    ZZ_pseudo[pseudo[:gdp], endo[:z_t]]          = 1.0
    DD_pseudo[pseudo[:gdp]]                      = 100*(exp(m[:z_star])-1)

	## Output: GDI
	ZZ_pseudo[pseudo[:gdi], endo[:y_t]]          = m[:γ_gdi]
    ZZ_pseudo[pseudo[:gdi], endo_addl[:y_t1]]    = -m[:γ_gdi]
    ZZ_pseudo[pseudo[:gdi], endo[:z_t]]          = m[:γ_gdi]
    DD_pseudo[pseudo[:gdi]]                      = 100*(exp(m[:z_star])-1) + m[:δ_gdi]

	## PCE Inflation
	ZZ_pseudo[pseudo[:InflationPCE], endo[:π_t]] = 1.0
    DD_pseudo[pseudo[:InflationPCE]]             = 100*(m[:π_star]-1)

	## Inflation (GDP Deflator)
	ZZ_pseudo[pseudo[:GDPdeflator], endo[:π_t]] = m[:Γ_gdpdef]
	DD_pseudo[pseudo[:GDPdeflator]]             = 100*(m[:π_star]-1) + m[:δ_gdpdef]

	## TFP
	ZZ_pseudo[pseudo[:tfp], endo[:z_t]]       = (1-m[:α])*m[:Iendoα] + 1*(1-m[:Iendoα])
	ZZ_pseudo[pseudo[:tfp], endo[:u_t]]       = m[:α]/( (1-m[:α])*(1-m[:Iendoα]) + 1*m[:Iendoα] )
	ZZ_pseudo[pseudo[:tfp], endo_addl[:u_t1]]  = -(m[:α]/( (1-m[:α])*(1-m[:Iendoα]) + 1*m[:Iendoα]) )

	## 10-year nominal bond yield
	ZZ_pseudo[pseudo[:LongRate], :] = ZZ_pseudo[pseudo[:NominalFFR], :]' * TTT10
	DD_pseudo[pseudo[:LongRate]]    = m[:Rstarn]

	## 20 yrs forward transition matrix
    TTT20 = if subspec(m) == "ss16"
        eye(size(TTT,1))
    else
        (1/80)*((UniformScaling(1.) - TTT) \ (UniformScaling(1.) - TTT^80))
    end

    ## BAA or BBB Spread
    ZZ_pseudo[pseudo[:BBBspread], endo[:ERtil_k_t]]    = 1.0
    ZZ_pseudo[pseudo[:BBBspread], endo[:R_t]]          = -1.0
    ZZ_pseudo[pseudo[:BBBspread], :]                   = ZZ_pseudo[pseudo[:BBBspread], :]' * TTT20
    DD_pseudo[pseudo[:BBBspread]]                      = 100*log(m[:spr]*m[:lnb_safe]*m[:lnb_liq])

    if subspec(m) in ["ss2", "ss4", "ss22"]
        ## Expected Average Nominal Natural Rate
        ZZ_pseudo[pseudo[:ExpectedAvgNominalNaturalRate], endo[:r_f_t]] = 1.
        ZZ_pseudo[pseudo[:ExpectedAvgNominalNaturalRate], endo[:Eπ_t]]  = 1.
        DD_pseudo[pseudo[:ExpectedAvgNominalNaturalRate]]               = m[:Rstarn]

        ## Nominal Rate Gap
        ZZ_pseudo[pseudo[:NominalRateGap], endo[:R_t]]   = 1.
        ZZ_pseudo[pseudo[:NominalRateGap], endo[:r_f_t]] = -1.
        ZZ_pseudo[pseudo[:NominalRateGap], endo[:Eπ_t]]  = -1.
    end

    ## Expected Average 10-Year Real Interest Rate
    ZZ_pseudo[pseudo[:ExpectedAvg10YearRealRate], :] = TTT10[endo[:R_t], :] - TTT10[endo[:Eπ_t], :]
    DD_pseudo[pseudo[:ExpectedAvg10YearRealRate]]    = m[:Rstarn] - 100*(m[:π_star]-1)

    ## Expected Average 10-Year Real Natural Rate
    ZZ_pseudo[pseudo[:ExpectedAvg10YearRealNaturalRate], :] = TTT10[endo[:r_f_t], :]
    DD_pseudo[pseudo[:ExpectedAvg10YearRealNaturalRate]]    = m[:Rstarn] - 100*(m[:π_star]-1)

    if subspec(m) in ["ss2", "ss4", "ss22"]
        ## Expected Average 10-Year Nominal Interest Rate
        ZZ_pseudo[pseudo[:ExpectedAvg10YearNominalRate], :] = TTT10[endo[:R_t], :]
        DD_pseudo[pseudo[:ExpectedAvg10YearNominalRate]]    = m[:Rstarn]

        ## Expected Average 10-Year Nominal Natural Rate
        ZZ_pseudo[pseudo[:ExpectedAvg10YearNominalNaturalRate], :] = TTT10[endo[:r_f_t], :] + TTT10[endo[:Eπ_t], :]
        DD_pseudo[pseudo[:ExpectedAvg10YearNominalNaturalRate]]    = m[:Rstarn]

        ## Expected Average 10-Year Rate Gap
        ZZ_pseudo[pseudo[:ExpectedAvg10YearRateGap], :] = TTT10[endo[:R_t], :] - TTT10[endo[:r_f_t], :] - TTT10[endo[:Eπ_t], :]
    end

    ## Expected Average 5-Year Real Interest Rate
    ZZ_pseudo[pseudo[:ExpectedAvg5YearRealRate], :] = TTT5[endo[:R_t], :] - TTT5[endo[:Eπ_t], :]
    DD_pseudo[pseudo[:ExpectedAvg5YearRealRate]]    = m[:Rstarn] - 100*(m[:π_star]-1)

    ## Expected Average 5-Year Real Natural Rate
    ZZ_pseudo[pseudo[:ExpectedAvg5YearRealNaturalRate], :] = TTT5[endo[:r_f_t], :]
    DD_pseudo[pseudo[:ExpectedAvg5YearRealNaturalRate]]    = m[:Rstarn] - 100*(m[:π_star]-1)

    if subspec(m) in ["ss2", "ss4", "ss22"]
        ## Expected Average 5-Year Nominal Interest Rate
        ZZ_pseudo[pseudo[:ExpectedAvg5YearNominalRate], :] = TTT5[endo[:R_t], :]
        DD_pseudo[pseudo[:ExpectedAvg5YearNominalRate]]    = m[:Rstarn]

        ## Expected Average 5-Year Nominal Natural Rate
        ZZ_pseudo[pseudo[:ExpectedAvg5YearNominalNaturalRate], :] = TTT5[endo[:r_f_t], :] + TTT5[endo[:Eπ_t], :]
        DD_pseudo[pseudo[:ExpectedAvg5YearNominalNaturalRate]]    = m[:Rstarn]

        ## Expected Average 5-Year Rate Gap
        ZZ_pseudo[pseudo[:ExpectedAvg5YearRateGap], :] = TTT5[endo[:R_t], :] -
        	TTT5[endo[:r_f_t], :] - TTT5[endo[:Eπ_t], :]
    end

    ## Expected Average 20-Year Real Interest Rate
    ZZ_pseudo[pseudo[:ExpectedAvg20YearRealRate], :] = TTT20[endo[:R_t], :] - TTT20[endo[:Eπ_t], :]
    DD_pseudo[pseudo[:ExpectedAvg20YearRealRate]]    = m[:Rstarn] - 100*(m[:π_star]-1)

    ## Expected Average 20-Year Real Natural Rate
    ZZ_pseudo[pseudo[:ExpectedAvg20YearRealNaturalRate], :] = TTT20[endo[:r_f_t], :]
    DD_pseudo[pseudo[:ExpectedAvg20YearRealNaturalRate]]    = m[:Rstarn] - 100*(m[:π_star]-1)

    if subspec(m) in ["ss2", "ss4", "ss22"]
        ## Expected Average 20-Year Nominal Interest Rate
        ZZ_pseudo[pseudo[:ExpectedAvg20YearNominalRate], :] = TTT20[endo[:R_t], :]
        DD_pseudo[pseudo[:ExpectedAvg20YearNominalRate]]    = m[:Rstarn]

        ## Expected Average 20-Year Nominal Natural Rate
        ZZ_pseudo[pseudo[:ExpectedAvg20YearNominalNaturalRate], :] = TTT20[endo[:r_f_t], :] + TTT20[endo[:Eπ_t], :]
        DD_pseudo[pseudo[:ExpectedAvg20YearNominalNaturalRate]]    = m[:Rstarn]

        ## Expected Average 20-Year Rate Gap
        ZZ_pseudo[pseudo[:ExpectedAvg20YearRateGap], :] = TTT20[endo[:R_t], :] -
        	TTT20[endo[:r_f_t], :] - TTT20[endo[:Eπ_t], :]

        ## 5-year forward rate gap
        TTT5_fwd = TTT^20
        ZZ_pseudo[pseudo[:Forward5YearRateGap], :] = TTT5_fwd[endo[:R_t], :] -
        	TTT5_fwd[endo[:r_f_t], :] - TTT5_fwd[endo[:Eπ_t], :]

        ## 10-year forward rate gap
        TTT10_fwd = TTT^40
        ZZ_pseudo[pseudo[:Forward10YearRateGap], :] = TTT10_fwd[endo[:R_t], :] - TTT10_fwd[endo[:r_f_t], :] - TTT10_fwd[endo[:Eπ_t], :]
    end

    ## 20-Year Forward Real Rate
    TTT20_fwd = TTT^80
    ZZ_pseudo[pseudo[:Forward20YearRealRate], :] = ZZ_pseudo[pseudo[:ExAnteRealRate], :]' * TTT20_fwd
    DD_pseudo[pseudo[:Forward20YearRealRate]]    = m[:Rstarn] - 100*(m[:π_star]-1)

    ## 20-Year Real Natural Forward Rate
    ZZ_pseudo[pseudo[:Forward20YearRealNaturalRate], :] = ZZ_pseudo[pseudo[:RealNaturalRate], :]' * TTT20_fwd
    DD_pseudo[pseudo[:Forward20YearRealNaturalRate]]    = m[:Rstarn] - 100*(m[:π_star]-1)

    ## 30-Year Real Forward Rate
    TTT30_fwd = TTT^120
    ZZ_pseudo[pseudo[:Forward30YearRealRate], :] = ZZ_pseudo[pseudo[:ExAnteRealRate], :]' * TTT30_fwd
    DD_pseudo[pseudo[:Forward30YearRealRate]]    = m[:Rstarn] - 100*(m[:π_star]-1)

    ## 30-Year Real Natural Forward Rate
    ZZ_pseudo[pseudo[:Forward30YearRealNaturalRate], :] = ZZ_pseudo[pseudo[:RealNaturalRate], :]' * TTT30_fwd
    DD_pseudo[pseudo[:Forward30YearRealNaturalRate]]    = m[:Rstarn] - 100*(m[:π_star]-1)

    if !(subspec(m) in ["ss2", "ss4", "ss22"])
        ## Nominal Natural Rate
        ZZ_pseudo[pseudo[:NominalNaturalRate], endo[:r_f_t]] = 1.
        ZZ_pseudo[pseudo[:NominalNaturalRate], endo[:Eπ_t]]  = 1.
        DD_pseudo[pseudo[:NominalNaturalRate]]               = m[:Rstarn]

        ## 5-Year Forward Real Rate
        TTT5_fwd = TTT^20
        ZZ_pseudo[pseudo[:Forward5YearRealRate], :] = ZZ_pseudo[pseudo[:ExAnteRealRate], :]' * TTT5_fwd
        DD_pseudo[pseudo[:Forward5YearRealRate]]    = m[:Rstarn] - 100*(m[:π_star]-1)

        ## 5-Year Real Natural Forward Rate
        ZZ_pseudo[pseudo[:Forward5YearRealNaturalRate], :] = ZZ_pseudo[pseudo[:RealNaturalRate], :]' * TTT5_fwd
        DD_pseudo[pseudo[:Forward5YearRealNaturalRate]]    = m[:Rstarn] - 100*(m[:π_star]-1)

        ## 10-Year Forward Real Rate
        TTT10_fwd = TTT^40
        ZZ_pseudo[pseudo[:Forward10YearRealRate], :] = ZZ_pseudo[pseudo[:ExAnteRealRate], :]' * TTT10_fwd
        DD_pseudo[pseudo[:Forward10YearRealRate]]    = m[:Rstarn] - 100*(m[:π_star]-1)

        ## 10-Year Real Natural Forward Rate
        ZZ_pseudo[pseudo[:Forward10YearRealNaturalRate], :] = ZZ_pseudo[pseudo[:RealNaturalRate], :]' * TTT10_fwd
        DD_pseudo[pseudo[:Forward10YearRealNaturalRate]]    = m[:Rstarn] - 100*(m[:π_star]-1)
    end

    ## Exogenous processes
    if subspec(m) in ["ss22"]
        to_add = [:g_t, :b_t, :μ_t, :z_t, :λ_f_t, :λ_w_t, :rm_t, :σ_ω_t, :μ_e_t,
                  :γ_t, :π_star_t]
        to_add_addl = [:lr_t, :tfp_t, :e_gdpdef_t, :e_corepce_t, :e_gdp_t, :e_gdi_t]
        for i in to_add
            ZZ_pseudo[pseudo[i], endo[i]] = 1.
        end
        for i in to_add_addl
            ZZ_pseudo[pseudo[i], endo_addl[i]] = 1.
        end
    end

    # Collect indices and transforms
    return PseudoMeasurement(ZZ_pseudo, DD_pseudo)
end


text:
module FlippingModel

using DataDeps, CSV
using DataFrames, ShiftedArrays
using Optim
using Distributions, StructArrays, StatsBase

export pokes, streaks, process_bouts, fit
export AbstractModel, Log_LikeliHood, params
export PoissonLapseUniform, PoissonLapseExponential

include("datadeps.jl")
include("loaddata.jl")
include("bouts.jl")
include("Models.jl")
include("LogLikeliHood.jl")
include("survival.jl")

greet() = print("Hello World!")

end # module


text:
# Autogenerated wrapper script for Cbc_jll for aarch64-apple-darwin-libgfortran5
export cbc, libCbc, libOsiCbc, libcbcsolver

using ASL_jll
using Cgl_jll
using Clp_jll
using Osi_jll
using CoinUtils_jll
using OpenBLAS32_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("Cbc")
JLLWrappers.@declare_library_product(libCbc, "@rpath/libCbc.3.dylib")
JLLWrappers.@declare_library_product(libOsiCbc, "@rpath/libOsiCbc.3.dylib")
JLLWrappers.@declare_library_product(libcbcsolver, "@rpath/libCbcSolver.3.dylib")
JLLWrappers.@declare_executable_product(cbc)
function __init__()
    JLLWrappers.@generate_init_header(ASL_jll, Cgl_jll, Clp_jll, Osi_jll, CoinUtils_jll, OpenBLAS32_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libCbc,
        "lib/libCbc.3.10.5.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libOsiCbc,
        "lib/libOsiCbc.3.10.5.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcbcsolver,
        "lib/libCbcSolver.3.10.5.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        cbc,
        "bin/cbc",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
pages = [
        "GlobalSensitivity.jl: Global Sensitivity Analysis (GSA)" => "index.md",
        "Tutorials" => Any[
          "tutorials/parallelized_gsa.md",
          "tutorials/juliacon21.md"
        ],
        "Methods" => Any[
          "methods/morris.md",
          "methods/sobol.md",
          "methods/regression.md",
          "methods/efast.md",
          "methods/delta.md",
          "methods/dgsm.md",
          "methods/easi.md",
          "methods/fractional.md",
          "methods/rbdfast.md"
        ],
    ]


text:
function getimagepreprocessing(data, computestats::Bool; kwargs...)
    if isnothing(data) && computestats
        error("If `computestats` is `true`, you have to pass in a data container `data`.")
    end
    return if computestats
        setup(ImagePreprocessing, Image{2}(), data; kwargs...)
    else
        ImagePreprocessing(; kwargs...)
    end
end


text:
# This file is a part of SimilaritySearch.jl

using SimilaritySearch, Random
using Test

#
# This file contains a set of tests for SearchGraph over databases of vectors (of Float32)
#

@testset "vector indexing with SearchGraph" begin
    # NOTE: The following algorithms are complex enough to say we are testing it doesn't have syntax errors, a more grained test functions are requiered
    Random.seed!(0)
    ksearch = 10
    n, m, dim = 1000, 30, 3

    db = [rand(Float32, dim) for i in 1:n]
    queries = [rand(Float32, dim) for i in 1:m]

    dist = SqL2Distance()
    seq = ExhaustiveSearch(dist, db)
    perf = Performance(seq, queries, ksearch)

    for search_algo_fun in [() -> IHCSearch(restarts=4), () -> BeamSearch(bsize=2)]
        search_algo = search_algo_fun()
        @info "=================== $search_algo"
        graph = SearchGraph(; dist, search_algo=search_algo)
        graph.neighborhood.reduce = SatNeighborhood()
        append!(graph, db)
        @time p = probe(perf, graph)
        @info "testing search_algo: $(string(graph.search_algo)), p: $(p)"
        @test p.macrorecall >= 0.6
        @info "queries per second: $(1/p.searchtime)"
        @info "===="
    end
end


text:
"""
    StateIterator{H,S}

An iterator for enumerating all the states in a basis
"""
struct StateIterator{T,H}
    basis::H
end

"""
    states(hilb) -> iterator

Returns an iterator to iterate all states in the hilbert space
"""
states(T::Type{<:Number}, h::AbstractHilbert) = StateIterator{T, typeof(h)}(h)
states(h::AbstractHilbert) = states(STD_REAL_PREC, h)

Base.length(iter::StateIterator) = spacedimension(iter.basis)
Base.eltype(iter::StateIterator) = typeof(state(iter.basis))
Base.getindex(iter::StateIterator{T}, i::Int) where T =
    state_i(T, iter.basis, i)

function Base.iterate(iter::StateIterator{T}, idx = 1) where T
    idx > spacedimension(iter.basis) && return nothing

    @inbounds σ = state_i(T, iter.basis, idx)
    return (σ, idx + 1)
end


text:
using GeoPhyInv
using Test
using Random
@init_parallel_stencil(3, false, Float32, 4)

mgrid = [range(0.0, stop=10.,step=0.1), range(0.0, stop=30.,step=0.2)];
ageom=AGeom(mgrid, SSrcs(10), Srcs(10), Recs(10));

tgrid=range(0, stop=1.0, step=0.1);

data=Records(tgrid, ageom, [:p,:vz]);

Random.randn!(data[3][:p]);


text:
export mapelites_step!,mapelites_run!
function mapelites_step!(e::Evolution,
                             map_el::MapElites,
                             map_ind_to_b::Function;
                             mutation::Function=Cambrian.uniform_mutation,
                             crossover::Function=Cambrian.uniform_crossover,
                             evaluate::Function=Cambrian.random_evaluate)
    e.gen += 1
    if (e.gen == 1)
        Cambrian.fitness_evaluate!(e;fitness=evaluate)
        for i in eachindex(e.population)
            add_to_map(map_el,map_ind_to_b(e.population[i]),e.population[i],e.population[i].fitness)
        end

    else

        new_pop = Array{Individual}(undef,0)

        if e.cfg["n_elite"] > 0
            sort!(e.population)
            append!(new_pop,
                    e.population[(length(e.population)-e.cfg["n_elite"]+1):end])
        end

        for i in (e.cfg["n_elite"]+1):e.cfg["n_population"]
            p1 = select_random(map_el)
            child = deepcopy(p1)

            if e.cfg["p_crossover"] > 0 && rand() < e.cfg["p_crossover"]
                parents = vcat(p1, select_random(map_el))
                child = crossover(parents...)
            end

            if e.cfg["m_rate"] > 0 && rand() < e.cfg["m_rate"]
                child = mutation(child)
            end

            push!(new_pop, child)
        end
        e.population = new_pop
        Cambrian.fitness_evaluate!(e;fitness=evaluate)
        for i in eachindex(e.population)
            add_to_map(map_el,map_ind_to_b(e.population[i]),e.population[i],e.population[i].fitness)
        end
    end
    if ((e.cfg["log_gen"] > 0) && mod(e.gen, e.cfg["log_gen"]) == 0)
        Cambrian.log_gen(e)
    end
    if ((e.cfg["save_gen"] > 0) && mod(e.gen, e.cfg["save_gen"]) == 0)
        Cambrian.save_gen(e)
    end
end

function mapelites_run!(e::Evolution,
                        mapel :: MapElites,
                        map_ind_to_b::Function;
                        mutation::Function=Cambrian.uniform_mutation,
                        crossover::Function=Cambrian.uniform_crossover,
                        evaluate::Function=Cambrian.random_evaluate)

    for i in (e.gen+1):e.cfg["n_gen"]
        mapelites_step!(e,mapel,map_ind_to_b;
                        mutation = mutation,
                        crossover = crossover,
                        evaluate = evaluate)
    end
end


text:
# https://adventofcode.com/2020/day/17

# Each cube only ever considers its neighbors: any of the 26 other cubes where any of their
# coordinates differ by at most 1. For example, given the cube at x=1,y=2,z=3, its neighbors
# include the cube at x=2,y=2,z=2, the cube at x=0,y=2,z=3, and so on.
#
# During a cycle, all cubes simultaneously change their state according to the following rules:
#
#  * If a cube is active and exactly 2 or 3 of its neighbors are also active, the cube remains
#    active. Otherwise, the cube becomes inactive.
#  * If a cube is inactive but exactly 3 of its neighbors are active, the cube becomes active.
#    Otherwise, the cube remains inactive.

import Base: ==, iterate, length

example = "
.#.
..#
###
"

Coords = Tuple{Vararg{Int}}

mutable struct Grid
    coords::Vector{Coords}
end

(==)(x::Grid, y::Grid) = all(x.coords .== y.coords)
iterate(grid::Grid, i...) = iterate(grid.coords, i...)
length(grid::Grid) = length(grid.coords)

function read_init(string::AbstractString, other_dims::Tuple{Vararg{Int}}=(1,))
    rows = split(string, '\n', keepempty=false)
    coords = Coords[]
    for (i, row) in enumerate(rows)
        for (j, ch) in enumerate(collect(row))
            if ch == '#'
                # using rotated coordinates, so they are in same order
                # as when looking at the ASCII examples
                push!(coords, (other_dims..., i, j))
            end
        end
    end
    return Grid(coords)
end

@assert read_init(example) == Grid([(1, 1, 2), (1, 2, 3), (1, 3, 1), (1, 3, 2), (1, 3, 3)])
@assert length(read_init(example)) == 5

function isneighbour(x::Coords, y::Coords)
    return !all(x .== y) && all(abs.(x .- y) .<= 1)
end

@assert isneighbour((1, 2, 3), (2, 2, 2))
@assert isneighbour((1, 2, 3), (0, 2, 3))
@assert !isneighbour((1, 2, 3), (3, 4, 1))

function find_neighbours(point::Coords, grid::Grid)
    return filter(x -> isneighbour(point, x), grid.coords)
end

@assert find_neighbours((1, 1, 1), read_init(example)) == Coords[(1, 1, 2)]

function cycle!(grid::Grid)
    new_coords = Coords[]
    for point in grid.coords
        # explore neighbourhood of every activated point
        for neighbourhood in Iterators.product([-1:1 for _ = 1:length(point)]...)
            new_point = point .+ neighbourhood
            if new_point in new_coords
                continue
            end
            neighbours_count = length(find_neighbours(new_point, grid))
            if new_point in grid.coords
                if 2 <= neighbours_count <= 3
                    push!(new_coords, new_point)
                end
            else
                if neighbours_count == 3
                    push!(new_coords, new_point)
                end
            end
        end
    end
    grid.coords = new_coords
    return grid
end

function boot!(init::Grid, num_cycles=6)
    for _ = 1:num_cycles
        cycle!(init)
    end
    return init
end

@assert length(cycle!(read_init(example))) == 11
@assert length(cycle!(cycle!(read_init(example)))) == 21

function part1(input)
    init = read_init(input)
    return length(boot!(init))
end

@assert part1(example) == 112

function part2(input)
    init = read_init(input, (1, 1))
    return length(boot!(init))
end

@assert part2(example) == 848

test = read("data/day-17.txt", String)
println("Part 1: $(result1 = part1(test))")
println("Part 2: $(result2 = part2(test))")

@assert result1 == 304
@assert result2 == 1868

max_stars_repo_name: twolodzko/advent-of-code

max_issues_repo_name: twolodzko/advent-of-code

max_forks_repo_name: twolodzko/advent-of-code


text:
using firstTask, Test

@testset "goldenSection" begin
  f(x) = x^4 + 3x^3 + x^2 + sin(x)
  a = -10
  b = 10
  tol = 1e-6
  @test goldenSection(f, a, b, tol) == (-1.9712755845857126, -4.915216432826815)

  f(x) = 3(x - 4)^2
  a = 0
  b = 10
  tol = 1e-6
  @test goldenSection(f, a, b, tol) == (4.000000004370129, 5.729408629607859e-17)
end

max_stars_repo_name: KK-studia/firstTask

max_issues_repo_name: KK-studia/firstTask

max_forks_repo_name: KK-studia/firstTask


text:
using CUDA
ff = CUDA.zeros(3)
function gggg(ff)

  shmemblockDataa = @cuDynamicSharedMem(Float32,(1,2))
    return
end
@cuda threads=(32,2) blocks=2 shmem = 3 gggg(ff)
ff[1]


text:
function pulaymix!(vin::Array{Float64,1}, vout::Array{Float64,1},
                   beta::Float64, dvmat, vmat, iter::Int64, mixdim::Int64)
    # output [vnew,df,dv]
    Npoints = size(vin)[1]
    #
    # construct and modify vmat and dvmat
    #
    if iter <= mixdim
        vmat[:,iter] = vin
        dvmat[:,iter] = vout - vin
    else
        # delete the first leading column  ???
        #vmat[:,1] = []; dvmat[:,1] = []
        # append vin and vout-vin at the end of vmat and dvmat;
        vmat[:,mixdim] = vin
        dvmat[:,mixdim] = vout - vin
    end

    thresh = 1e-6
    if iter > 1
        ibeg = 2
        iend = min(iter,mixdim)
        A = dvmat[:,ibeg:iend]
        b = dvmat[:,ibeg-1]
        ncols = size(A)[2]
        A = b*ones(1,ncols) - A
        U,S,V = svd(A,thin=true)
        while S[ncols] <  thresh*S[1]
            ibeg = ibeg + 1
            if ibeg > iend
                @printf("warning: no mixing...\n")
                vnew = vout
                return
            end
            A = dvmat[:,ibeg:iend]
            b = dvmat[:,ibeg-1]
            ncols = size(A)[2]
            A = b*ones(1,ncols) - A
            U,S,V = svd(A,thin=true)
        end
        g2 = V*(Matrix(Diagonal(S))\(U'*b))
        g1 = 1 - sum(g2)
        g = [g1;g2]
        vnew = ( vmat[:,ibeg-1:iend] + beta*dvmat[:,ibeg-1:iend] )*g
    else
        vnew = vout
    end
    return vnew
end


text:
using RDatasets, DataArrays, Gadfly

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


text:
""" [Grabowski2006](@cite) """
function TRMM_LBA_z(::Type{FT}) where {FT}
    z_in = FT[0.130,  0.464,  0.573,  1.100,  1.653,  2.216,  2.760,
              3.297,  3.824,  4.327,  4.787,  5.242,  5.686,  6.131,
              6.578,  6.996,  7.431,  7.881,  8.300,  8.718,  9.149,
              9.611, 10.084, 10.573, 11.008, 11.460, 11.966, 12.472,
             12.971, 13.478, 13.971, 14.443, 14.956, 15.458, 16.019,
             16.491, 16.961, 17.442, 17.934, 18.397, 18.851, 19.331,
             19.809, 20.321, 20.813, 21.329, 30.000] .* 1000 .- 130 #LES z is in meters
    return z_in
end
""" [Grabowski2006](@cite) """
function TRMM_LBA_p(::Type{FT}) where {FT}
    z_in = TRMM_LBA_z(FT)
    p_in = FT[991.3, 954.2, 942.0, 886.9, 831.5, 778.9, 729.8,
              684.0, 641.7, 603.2, 570.1, 538.6, 509.1, 480.4,
              454.0, 429.6, 405.7, 382.5, 361.1, 340.9, 321.2,
              301.2, 281.8, 263.1, 246.1, 230.1, 213.2, 197.0,
              182.3, 167.9, 154.9, 143.0, 131.1, 119.7, 108.9,
              100.1,  92.1,  84.6,  77.5,  71.4,  65.9,  60.7,
               55.9,  51.3,  47.2,  43.3,  10.3] .* 100 # LES pres is in pasc
    return Dierckx.Spline1D(z_in, p_in; k = 1)
end
""" [Grabowski2006](@cite) """
function TRMM_LBA_T(::Type{FT}) where {FT}
    z_in = TRMM_LBA_z(FT)
    T_in = FT[23.70,  23.30,  22.57,  19.90,  16.91,  14.09,  11.13,
               8.29,   5.38,   2.29,  -0.66,  -3.02,  -5.28,  -7.42,
             -10.34, -12.69, -15.70, -19.21, -21.81, -24.73, -27.76,
             -30.93, -34.62, -38.58, -42.30, -46.07, -50.03, -54.67,
             -59.16, -63.60, -67.68, -70.77, -74.41, -77.51, -80.64,
             -80.69, -80.00, -81.38, -81.17, -78.32, -74.77, -74.52,
             -72.62, -70.87, -69.19, -66.90, -66.90] .+ FT(273.15) # LES T is in deg K
    return Dierckx.Spline1D(z_in, T_in; k = 1)
end
""" [Grabowski2006](@cite) """
function TRMM_LBA_RH(::Type{FT}) where {FT}
  z_in = TRMM_LBA_z(FT)
    RH_in = FT[98.00,  86.00,  88.56,  87.44,  86.67,  83.67,  79.56,
               84.78,  84.78,  89.33,  94.33,  92.00,  85.22,  77.33,
               80.11,  66.11,  72.11,  72.67,  52.22,  54.67,  51.00,
               43.78,  40.56,  43.11,  54.78,  46.11,  42.33,  43.22,
               45.33,  39.78,  33.78,  28.78,  24.67,  20.67,  17.67,
               17.11,  16.22,  14.22,  13.00,  13.00,  12.22,   9.56,
                7.78,   5.89,   4.33,   3.00,   3.00]
    return Dierckx.Spline1D(z_in, RH_in; k = 1)
end
""" [Grabowski2006](@cite) """
function TRMM_LBA_u(::Type{FT}) where {FT}
    z_in = TRMM_LBA_z(FT)
    u_in = FT[0.00,   0.81,   1.17,   3.44,   3.53,   3.88,   4.09,
              3.97,   1.22,   0.16,  -1.22,  -1.72,  -2.77,  -2.65,
             -0.64,  -0.07,  -1.90,  -2.70,  -2.99,  -3.66,  -5.05,
             -6.64,  -4.74,  -5.30,  -6.07,  -4.26,  -7.52,  -8.88,
             -9.00,  -7.77,  -5.37,  -3.88,  -1.15,  -2.36,  -9.20,
             -8.01,  -5.68,  -8.83, -14.51, -15.55, -15.36, -17.67,
            -17.82, -18.94, -15.92, -15.32, -15.32]
    return Dierckx.Spline1D(z_in, u_in; k = 1)
end
""" [Grabowski2006](@cite) """
function TRMM_LBA_v(::Type{FT}) where {FT}
    z_in = TRMM_LBA_z(FT)
    v_in = FT[-0.40,  -3.51,  -3.88,  -4.77,  -5.28,  -5.85,  -5.60,
              -2.67,  -1.47,   0.57,   0.89,  -0.08,   1.11,   2.15,
               3.12,   3.22,   3.34,   1.91,   1.15,   1.01,  -0.57,
              -0.67,   0.31,   2.97,   2.32,   2.66,   4.79,   3.40,
               3.14,   3.93,   7.57,   2.58,   2.50,   6.44,   6.84,
               0.19,  -2.20,  -3.60,   0.56,   6.68,   9.41,   7.03,
               5.32,   1.14,  -0.65,   5.27,   5.27]
    return Dierckx.Spline1D(z_in, v_in; k = 1)
end

""" [Grabowski2006](@cite) """
function TRMM_LBA_tke(::Type{FT}) where {FT}
    z -> if z <= 2500.0
        FT(1) - z / 3000
    else
        FT(0)
    end
end

""" [Grabowski2006](@cite) """
function TRMM_LBA_radiation(::Type{FT}) where {FT}
    rad_time = range(FT(10), FT(360); length = 36) .* 60
    z_in = FT[42.5, 200.92, 456.28, 743, 1061.08, 1410.52, 1791.32, 2203.48, 2647,3121.88, 3628.12,
              4165.72, 4734.68, 5335, 5966.68, 6629.72, 7324.12,
              8049.88, 8807, 9595.48, 10415.32, 11266.52, 12149.08, 13063, 14008.28,
              14984.92, 15992.92, 17032.28, 18103, 19205.08, 20338.52, 21503.32, 22699.48]
    rad_in = [[-1.386, -1.927, -2.089, -1.969, -1.805, -1.585, -1.406, -1.317, -1.188, -1.106, -1.103, -1.025,
              -0.955, -1.045, -1.144, -1.119, -1.068, -1.092, -1.196, -1.253, -1.266, -1.306,  -0.95,  0.122,
               0.255,  0.258,  0.322,  0.135,      0,      0,      0,      0,      0],
             [ -1.23, -1.824, -2.011, -1.895, -1.729, -1.508, -1.331, -1.241, -1.109, -1.024, -1.018,  -0.94,
              -0.867, -0.953, -1.046, -1.018, -0.972, -1.006, -1.119, -1.187, -1.209, -1.259, -0.919,  0.122,
               0.264,  0.262,  0.326,  0.137,      0,      0,      0,      0,     0],
             [-1.043, -1.692, -1.906, -1.796,  -1.63,  -1.41, -1.233, -1.142,  -1.01,  -0.92, -0.911, -0.829,
              -0.754, -0.837, -0.923,  -0.89, -0.847, -0.895, -1.021, -1.101, -1.138, -1.201,  -0.88,  0.131,
               0.286,  0.259,  0.332,   0.14,      0,      0,      0,      0,      0],
             [-0.944, -1.613, -1.832,  -1.72, -1.555, -1.339, -1.163, -1.068, -0.935, -0.846, -0.835,  -0.75,
              -0.673, -0.751, -0.833, -0.798,  -0.76, -0.817, -0.952, -1.042, -1.088, -1.159, -0.853,  0.138,
               0.291,  0.265,  0.348,  0.136,      0,      0,      0,      0,      0],
             [-0.833, -1.526, -1.757, -1.648, -1.485,  -1.27, -1.093, -0.998, -0.867, -0.778, -0.761, -0.672,
              -0.594, -0.671, -0.748, -0.709, -0.676, -0.742, -0.887, -0.986, -1.041, -1.119, -0.825,  0.143,
               0.296,  0.271,  0.351,  0.138,      0,      0,      0,      0,      0],
             [-0.719, -1.425, -1.657,  -1.55, -1.392, -1.179, -1.003, -0.909, -0.778, -0.688, -0.667, -0.573,
              -0.492, -0.566, -0.639, -0.596, -0.568, -0.647, -0.804, -0.914, -0.981,  -1.07, -0.793,  0.151,
               0.303,  0.279,  0.355,  0.141,      0,      0,      0,      0,      0],
             [-0.724, -1.374, -1.585, -1.482, -1.328, -1.116, -0.936, -0.842, -0.715, -0.624, -0.598, -0.503,
              -0.421, -0.494, -0.561, -0.514,  -0.49,  -0.58, -0.745, -0.863, -0.938, -1.035, -0.764,  0.171,
               0.291,  0.284,  0.358,  0.144,      0,      0,      0,      0,      0],
             [-0.587,  -1.28, -1.513, -1.416, -1.264, -1.052, -0.874, -0.781, -0.655, -0.561, -0.532, -0.436,
              -0.354, -0.424, -0.485, -0.435, -0.417, -0.517, -0.691, -0.817, -0.898,     -1,  -0.74,  0.176,
               0.297,  0.289,   0.36,  0.146,      0,      0,      0,      0,      0],
             [-0.506, -1.194, -1.426, -1.332, -1.182, -0.972, -0.795, -0.704, -0.578,  -0.48, -0.445, -0.347,
              -0.267, -0.336, -0.391, -0.337, -0.325, -0.436,  -0.62, -0.756, -0.847,  -0.96, -0.714,   0.18,
               0.305,  0.317,  0.348,  0.158,      0,      0,      0,      0,      0],
             [-0.472,  -1.14, -1.364, -1.271, -1.123, -0.914, -0.738, -0.649, -0.522, -0.422, -0.386, -0.287,
              -0.207, -0.273, -0.322, -0.267,  -0.26, -0.379, -0.569, -0.712, -0.811, -0.931, -0.696,  0.183,
               0.311,   0.32,  0.351,   0.16,      0,      0,      0,      0,     0],
             [-0.448, -1.091, -1.305, -1.214, -1.068, -0.858, -0.682, -0.594, -0.469, -0.368, -0.329, -0.229,
              -0.149, -0.213, -0.257,   -0.2, -0.199, -0.327, -0.523, -0.668, -0.774, -0.903, -0.678,  0.186,
               0.315,  0.323,  0.355,  0.162,      0,      0,      0,      0,      0],
             [-0.405, -1.025, -1.228, -1.139, -0.996, -0.789, -0.615, -0.527, -0.402,   -0.3, -0.256, -0.156,
              -0.077, -0.136, -0.173, -0.115, -0.121, -0.259, -0.463, -0.617, -0.732, -0.869, -0.656,   0.19,
               0.322,  0.326,  0.359,  0.164,      0,      0,      0,      0,      0],
             [-0.391, -0.983, -1.174, -1.085, -0.945, -0.739, -0.566, -0.478, -0.354, -0.251, -0.205, -0.105,
              -0.027, -0.082, -0.114, -0.056, -0.069, -0.213,  -0.42, -0.579, -0.699,  -0.84, -0.642,  0.173,
               0.327,  0.329,  0.362,  0.165,      0,      0,      0,      0,      0],
             [-0.385, -0.946, -1.121, -1.032, -0.898, -0.695, -0.523, -0.434, -0.307, -0.203, -0.157, -0.057,
               0.021, -0.031, -0.059, -0.001, -0.018, -0.168, -0.381, -0.546, -0.672, -0.819, -0.629,  0.176,
               0.332,  0.332,  0.364,  0.166,      0,      0,      0,      0,      0],
             [-0.383, -0.904, -1.063, -0.972, -0.834, -0.632, -0.464, -0.378, -0.251, -0.144, -0.096,  0.001,
               0.079,  0.032,  0.011,  0.069,  0.044, -0.113, -0.332, -0.504, -0.637, -0.791, -0.611,  0.181,
               0.338,  0.335,  0.367,  0.167,      0,      0,      0,      0,      0],
             [-0.391, -0.873, -1.016, -0.929, -0.794, -0.591, -0.423, -0.337, -0.212, -0.104, -0.056,  0.043,
               0.121,  0.077,  0.058,  0.117,  0.088, -0.075, -0.298, -0.475, -0.613, -0.772, -0.599,  0.183,
               0.342,  0.337,   0.37,  0.168,      0,      0,      0,      0,      0],
             [-0.359, -0.836, -0.976, -0.888, -0.755, -0.554, -0.386,   -0.3, -0.175, -0.067, -0.018,  0.081,
                0.16,  0.119,  0.103,  0.161,  0.129, -0.039, -0.266, -0.448, -0.591, -0.755, -0.587,  0.187,
               0.345,  0.339,  0.372,  0.169,      0,      0,      0,      0,     0],
             [-0.328, -0.792, -0.928, -0.842, -0.709, -0.508, -0.341, -0.256, -0.131, -0.022,  0.029,  0.128,
               0.208,   0.17,  0.158,  0.216,  0.179,  0.005, -0.228, -0.415, -0.564, -0.733, -0.573,   0.19,
               0.384,  0.313,  0.375,   0.17,      0,      0,      0,      0,      0],
             [-0.324, -0.767, -0.893, -0.807, -0.676, -0.476,  -0.31, -0.225, -0.101,  0.008,   0.06,  0.159,
               0.239,  0.204,  0.195,  0.252,  0.212,  0.034, -0.203, -0.394, -0.546, -0.719, -0.564,  0.192,
               0.386,  0.315,  0.377,  0.171,      0,      0,      0,      0,      0],
             [ -0.31,  -0.74,  -0.86, -0.775, -0.647, -0.449, -0.283, -0.197, -0.073,  0.036,  0.089,  0.188,
               0.269,  0.235,  0.229,  0.285,  0.242,  0.061, -0.179, -0.374,  -0.53, -0.706, -0.556,  0.194,
               0.388,  0.317,  0.402,  0.158,      0,      0,      0,      0,      0],
             [-0.244, -0.694, -0.818,  -0.73, -0.605, -0.415, -0.252, -0.163, -0.037,  0.072,  0.122,   0.22,
               0.303,  0.273,  0.269,  0.324,  0.277,  0.093, -0.152,  -0.35,  -0.51, -0.691, -0.546,  0.196,
               0.39,   0.32,  0.403,  0.159,      0,      0,      0,      0,      0],
             [-0.284, -0.701, -0.803, -0.701, -0.568, -0.381, -0.225, -0.142, -0.017,  0.092,  0.143,  0.242,
               0.325,  0.298,  0.295,   0.35,    0.3,  0.112, -0.134, -0.334, -0.497,  -0.68,  -0.54,  0.198,
               0.392,  0.321,  0.404,   0.16,      0,      0,      0,      0,      0],
             [-0.281, -0.686, -0.783,  -0.68, -0.547, -0.359, -0.202, -0.119,  0.005,  0.112,  0.163,  0.261,
               0.345,  0.321,  0.319,  0.371,  0.319,   0.13, -0.118, -0.321, -0.486, -0.671, -0.534,  0.199,
               0.393,  0.323,  0.405,  0.161,      0,      0,      0,      0,      0],
             [-0.269, -0.667,  -0.76, -0.655, -0.522, -0.336, -0.181, -0.096,  0.029,  0.136,  0.188,  0.286,
                0.37,  0.346,  0.345,  0.396,  0.342,   0.15, -0.102, -0.307, -0.473, -0.661, -0.528,    0.2,
               0.393,  0.324,  0.405,  0.162,      0,      0,      0,      0,      0],
             [-0.255, -0.653, -0.747, -0.643, -0.511, -0.325, -0.169, -0.082,  0.042,  0.149,  0.204,  0.304,
               0.388,  0.363,  0.36 ,  0.409,  0.354,  0.164, -0.085, -0.289, -0.457, -0.649, -0.523,  0.193,
               0.394,  0.326,  0.406,  0.162,      0,      0,      0,      0,      0],
             [-0.265,  -0.65, -0.739, -0.634,   -0.5, -0.314, -0.159, -0.072,  0.052,  0.159,  0.215,  0.316,
               0.398,  0.374,  0.374,  0.424,   0.37,  0.181, -0.065, -0.265, -0.429, -0.627, -0.519,   0.18,
               0.394,  0.326,  0.406,  0.162,      0,      0,      0,      0,      0],
             [-0.276, -0.647, -0.731, -0.626, -0.492, -0.307, -0.152, -0.064,  0.058,  0.166,  0.227,  0.329,
               0.411,  0.389,   0.39,  0.441,  0.389,  0.207, -0.032, -0.228, -0.394, -0.596, -0.494,  0.194,
               0.376,  0.326,  0.406,  0.162,      0,      0,      0,      0,      0],
             [-0.271, -0.646,  -0.73, -0.625, -0.489, -0.303, -0.149, -0.061,  0.062,  0.169,  0.229,  0.332,
               0.412,  0.388,  0.389,  0.439,  0.387,  0.206, -0.028, -0.209, -0.347, -0.524, -0.435,  0.195,
               0.381,  0.313,  0.405,  0.162,      0,      0,      0,      0,      0],
             [-0.267, -0.647, -0.734, -0.628,  -0.49, -0.304, -0.151, -0.062,  0.061,  0.168,  0.229,  0.329,
               0.408,  0.385,  0.388,  0.438,  0.386,  0.206, -0.024, -0.194, -0.319,  -0.48,  -0.36,  0.318,
               0.405,  0.335,  0.394,  0.162,      0,      0,      0,      0,      0],
             [-0.274, -0.656, -0.745,  -0.64,   -0.5, -0.313, -0.158, -0.068,  0.054,  0.161,  0.223,  0.325,
               0.402,  0.379,  0.384,  0.438,  0.392,  0.221,  0.001, -0.164, -0.278, -0.415, -0.264,  0.445,
               0.402,  0.304,  0.389,  0.157,      0,      0,      0,      0,      0],
             [-0.289, -0.666, -0.753, -0.648, -0.508,  -0.32, -0.164, -0.073,  0.049,  0.156,   0.22,  0.321,
               0.397,  0.374,  0.377,   0.43,  0.387,  0.224,  0.014, -0.139, -0.236, -0.359, -0.211,  0.475,
                 0.4,  0.308,  0.375,  0.155,      0,      0,      0,      0,      0],
             [-0.302, -0.678, -0.765, -0.659, -0.517, -0.329, -0.176, -0.085,  0.038,  0.145,  0.208,   0.31,
               0.386,  0.362,  0.366,  0.421,  0.381,  0.224,  0.022, -0.119, -0.201,   -0.3, -0.129,  0.572,
               0.419,  0.265,  0.364,  0.154,      0,      0,      0,      0,      0],
             [-0.314, -0.696, -0.786, -0.681, -0.539, -0.349, -0.196, -0.105,  0.019,  0.127,  0.189,  0.289,
               0.364,   0.34,  0.346,  0.403,   0.37,  0.222,  0.036, -0.081, -0.133, -0.205, -0.021,  0.674,
               0.383,  0.237,  0.359,  0.151,      0,      0,      0,      0,      0],
             [-0.341, -0.719, -0.807, -0.702, -0.558, -0.367, -0.211,  -0.12,  0.003,  0.111,  0.175,  0.277,
               0.351,  0.325,  0.331,   0.39,   0.36,  0.221,  0.048, -0.046, -0.074, -0.139,  0.038,  0.726,
               0.429,  0.215,  0.347,  0.151,      0,      0,      0,      0,      0],
             [ -0.35, -0.737, -0.829, -0.724, -0.577, -0.385, -0.229, -0.136, -0.011,  0.098,  0.163,  0.266,
               0.338,   0.31,  0.316,  0.378,  0.354,  0.221,  0.062, -0.009, -0.012, -0.063,  0.119,  0.811,
               0.319,  0.201,  0.343,  0.148,      0,      0,      0,      0,      0],
             [-0.344,  -0.75, -0.856, -0.757, -0.607, -0.409,  -0.25, -0.156, -0.033,  0.076,  0.143,  0.246,
               0.316,  0.287,  0.293,  0.361,  0.345,  0.225,  0.082,  0.035,  0.071,  0.046,  0.172,  0.708,
               0.255,   0.21,  0.325,  0.146,      0,      0,      0,      0,      0]] ./ 86400

    rad_in = hcat(rad_in...)'
    profile = Dierckx.Spline2D(rad_time, z_in, rad_in; kx = 1, ky = 1)
    return profile
end


text:
module TestLoading

# using Revise
using Test
using MLJModels

@loadcode RandomForestClassifier pkg=DecisionTree verbosity=0

@load RidgeRegressor pkg=MultivariateStats verbosity=0

@test isdefined(TestLoading, :RidgeRegressor)
@test MLJModels.info("RidgeRegressor", pkg="MultivariateStats") in
    localmodels(modl=TestLoading)

# if we load the same model again:
program, _ = @test_logs((:info, r"For silent"),
           (:info, r"Model code"),
           MLJModels._load(TestLoading,
                           :(RidgeRegressor),
                           :(pkg=MultivariateStats)))
eval(program)

@test !isdefined(TestLoading, :RidgeRegressor2)

# load the same model again, with a different binding:
@load RidgeRegressor pkg=MultivariateStats name=Foo

@test Foo() == RidgeRegressor()

# try to use the name of an existing object for new type name
program, _ = MLJModels._load(TestLoading,
                :(DecisionTreeClassifier),
                :(pkg=DecisionTree),
                             :(name=RidgeRegressor))

@test_throws Exception eval(program)

@test_throws Exception load("model", pkg = "pkg")
@test_throws Exception load(models()[1])

@testset "scope=:local inside a @testset" begin
    @load RidgeRegressor pkg=MultivariateStats verbosity=0 scope=:local
end

@testset "install_pkgs=true" begin
    @load KMeans pkg=Clustering verbosity=0 scope=:local install_pkgs=true
    @load KMeans pkg=Clustering verbosity=0 scope=:local install_pkgs=true
end

end # module

true


text:
import Registrator.CommentBot: accept_regex, parse_comment

using Test

@testset "Server" begin

@testset "accept_regex" begin
    r = Regex("@JuliaRegistrator $accept_regex")
    @test match(r, "@JuliaRegistrator register()")[1] == "register()"
    @test match(r, "@JuliaRegistrator register()\r\n\r\n")[1] == "register()"
    @test match(r, """@JuliaRegistrator register()
                Patch notes:
                Release v0.1.0
                Fixes bugs
                """)[1] == "register()"
end

@testset "parse_comment" begin
    @test parse_comment("register()") == ("register", [], Dict{Symbol,String}())
    @test parse_comment("register(qux)") == ("register", ["qux"], Dict{Symbol,String}())
    @test parse_comment("register(qux, baz)") == ("register", ["qux", "baz"], Dict{Symbol,String}())

    @test parse_comment("approved()") == ("approved", [], Dict{Symbol,String}())

    @test parse_comment("register(target=qux)") == ("register", String[], Dict(:target=>"qux"))
    @test parse_comment("register(target=qux, branch=baz)") == ("register", String[], Dict(:branch=>"baz",:target=>"qux"))
    @test parse_comment("register(qux, baz, target=foo, branch=bar)") == ("register", ["qux", "baz"], Dict(:branch=>"bar",:target=>"foo"))

    @test parse_comment("register(branch=\"release-1.0\")") == ("register", String[], Dict(:branch=>"release-1.0"))
end

end


text:
################################################################################
# Function that writes tabular data into table
# Latex dependencies:
# - \usepackage[flushleft]{threeparttable}
# - \usepacakge{booktabs}
# Needs threeparttable package (\usepackage[flushleft]{threeparttable}) in Latex
function writeLatexTable(d,fname,caption,colnames,rownames,comment;
	label = "", midrules=[],format="%.2f", pos = "htb",percentcols = [])
	f = open(fname,"w+")

	n = size(d,1)
	k = size(d,2)

	tt = ["\\begin{table}[$pos] \\centering ";
			"\\caption{$caption} \\label{$label}" ;
			"\\begin{threeparttable}" ;
			"\\begin{tabular}{$(repeat("c",k+2))} ";
			"\\midrule"]
	for ii = 1: size(tt,1)
		println(f,tt[ii])
	end
	print(f,"&")
	for ii = 1:k
		print(f,"& $(colnames[ii])")
	end
	println(f,"\\\\")
	println(f,"\\midrule")
	for ii = 1:n
		print(f, rownames[ii])
		print(f,"&")
		for jj = 1:k
			print(f,"&")
			if !isnan(d[ii,jj])
				fmttxt = @eval @sprintf $format $d[$ii,$jj]
				if jj in percentcols
					fmttxt *= "\\%"
				end
				print(f,fmttxt)
			else
				print(f," ")
			end
		end
		println(f,"\\\\")
		if !isempty(midrules) && midrules[ii]==1
			println(f,"\\midrule")
		end

	end
	println(f,"\\midrule ")
	println(f,"\\end{tabular}")
	if !(isempty(comment))
		println(f,"\\begin{tablenotes}")
		println(f,"\\item \\footnotesize{\\emph{Notes:} $comment }")
		println(f,"\\end{tablenotes} ")
	end
	println(f,"\\end{threeparttable}")
	println(f,"\\end{table}")

	close(f)
end

#####
# Overload functions in Distributions
import Distributions: cov
# Required since cov(::UnivariateDistribution) not defined in Distributions.jl
function cov(d::UnivariateDistribution)
	return var(d)
end

#####################################################################
# Printing

function printD(d,i)
	println("Purch = $(d[i].purch) ")
	println("Stop = $(d[i].stop) ")
	println("Zd = $(round.(d[i].zd,digits=2))")
	println("Pid | Pos  | Path | Zs   | u ")
	display(round.([d[i].pid d[i].pos d[i].path d[i].zs d[i].u ],digits=2))
end

#####################################################################
# Max/min functions used in code (for performance)
@inline function findmax_subsel(A,S)
	val, ind = typemin(eltype(A)), -1
	for ii in eachindex(A)
		val_ = A[ii]
		if val_ > val && S[ii] == true
			val, ind = val_, ii
		end
	end
	return val, ind
end

@inline function findmax_range(A,R)
	val, ind = typemin(eltype(A)), -1
	for ii in R
		val_ = A[ii]
		if val_ > val
			val, ind = val_, ii
		end
	end
	return val, ind
end

@inline function findmin_range(A,R)
	val, ind = typemax(eltype(A)), -1
	for ii in R
		val_ = A[ii]
		if val_ < val
			val, ind = val_, ii
		end
	end
	return val, ind
end

@inline function findmax_range(A,R,ee)
	val, ind = typemin(eltype(A)), -1
	for ii in R
		val_ = A[ii,ee]
		if val_ > val
			val, ind = val_, ii
		end
	end
	return val, ind
end

#####################################################################
# Get range within thread
@inline function getrange(n)
	tid = Threads.threadid()
	nt = Threads.nthreads()
	d , r = divrem(n, nt)
	from = (tid - 1) * d + min(r, tid - 1) + 1
	to = from + d - 1 + (tid ≤ r ? 1 : 0)
	from:to
end

#####################################################################
# Translate string cost spec into a function

function getCfun(s::String)
	if s == "linear"
		(cs::Float64) -> cs
	elseif s == "exp"
		(cs::Float64) -> exp(cs)
	elseif s == "quadratic"
		(cs::Float64) -> cs^2
	else
		error("cfun not correctly specified")
	end
end

function getZdfun(s::String)
	if s == "linear"
		(z,dze,pos) -> z + dze * pos
	elseif s.zdfun == "exp"
		(z,dze,pos) -> z - exp(dze*pos) + 1.0
	elseif s.zdfun == "quadratic"
		(z,dze,pos) -> z + dze * pos^2
	else
		error("zdfun not correctly specified")
	end
end

#####################################################################
# Draw products
function genProducts(m::Model,seli,nProdTot,nProdCons;
						weights::AbstractWeights= fweights(fill(1,nProdTot-1)))
		# Derived values
		nChar = length(m.dChars) # number of characteristics (excl. outside option dummy)
		nCons = length(seli)
		# Product ids, 0 will be outside option
		pid = 1:nProdTot-1 ;

		# Product characteristics
		if nChar == 1
			chars = [0.0 ;rand(m.dChars,nProdTot-1,1)]
		else
			chars = [fill(0.0,1,nChar) ; permutedims(rand(m.dChars,nProdTot-1),[2,1])]
		end

		# Empty Array for consumer-products
		charsC = zeros(seli[end][end],nChar + m.outDummy)

		# Draw products per consumer randomly from set of products for each consumer
		pidC = [ [0 ; sample(pid,weights,nProdCons-1,replace=false)] for i = 1:nCons]

		# Fill into array for consumer-products
		for ii = 1:nCons
			charsC[seli[ii],1:end-m.outDummy] = chars[pidC[ii] .+ 1,:]
			if m.outDummy
				charsC[seli[ii][1],end] = 1.0
			end
		end
	return pidC,charsC
end


text:
module OpenHPLJulia

using DataStructures
include("parameters.jl")
include("fluid.jl")
include("unitdesign.jl")
export penstock, headraceNsurgetank,f2x1y
end # module