text:
#-------------------------------------------------------------------------------
# File: greedy.jl
# Description: This file contains all functions related to the solution.
# Date: October 31, 2019
# Author: Jonathan Fontaine, Killian Fretaud, Rémi Garcia,
#         Boualem Lamraoui, Benoît Le Badezet, Benoit Loger
#-------------------------------------------------------------------------------
"""
    Batch

Initialize the set of batches
"""
mutable struct Batch
    width::Int
        # Width of the batch
    start::Int
        # First position of the batch

    Batch(w::Int,s::Int) = new(w, s)
end

"""
    Solution

Create a empty sequence with matrices of zeros.
"""
mutable struct Solution
    sequence::Array{Int, 1}
        # vector of cars (sequence pi or s in the following algorithms)
    M1::Array{Int, 2}
        # M1_ij is the number of cars with oj in the subsequence of
        # q(oj) cars starting at position i of sequence pi
    M2::Array{Int, 2}
        # M2_ij is the number of subsequences of q(oj) cars starting
        # at positions 1 up to i in which the number of cars that
        # require oj is greater than p(oj)
    M3::Array{Int, 2}
        # M3_ij is the number of subsequences in which the number of
        # cars that require oj is greater than or equal to p(oj)

    length::Int
    colors::Union{Nothing, Array{Batch, 1}}

    Solution(nC::Int,nO::Int) = new(
        zeros(Int,nC),#collect(1:nC),
        zeros(Int,nO,nC),
        zeros(Int,nO,nC),
        zeros(Int,nO,nC),
        0,
        nothing
    )
end

# Build an initial
"""

"""
function init_solution(instance::Instance)
    n = instance.nb_cars
    m = instance.nb_HPRC + instance.nb_LPRC # number of ratio
    solution = Solution(n, m)

    for i in 1:instance.nb_late_prec_day
        solution.sequence[i] = i
    end
    solution.length = instance.nb_late_prec_day
    return solution
end

init_solution(nom_fichier::String, type_fichier::String) =
    init_solution(parser(nom_fichier, type_fichier))

"""

"""
function is_solution_valid(solution::Solution, instance::Instance)
    if solution.length != instance.nb_cars
        return false
    end

    batch = instance.nb_late_prec_day+1
    valid = true
    while batch <= solution.length
        if solution.colors[batch].width <= instance.nb_paint_limitation
            batch += solution.colors[batch].width
        else
            valid = false
            break
        end
    end
    return valid
end

"""
    update_matrices!(solution::Solution, instance::Instance)

Updates `solution.M1`, `solution.M2` and `solution.M3` for known cars at positions 1 to `solution.length`.
"""
function update_matrices!(solution::Solution, instance::Instance)
    nb_RC = instance.nb_HPRC + instance.nb_LPRC
    nb = solution.length

    # Last column has just one car
    for option in 1:nb_RC
        car = solution.sequence[nb]
        if instance.RC_flag[car, option]
            solution.M1[option, nb] = 1
        else
            solution.M1[option, nb] = 0
        end
    end

    # Dynamic (right to left)
    for counter in 1:(nb-1)
        index = nb - counter
        car = solution.sequence[index]
        for option in 1:nb_RC
            # Has option -> next sequence + 1
            if instance.RC_flag[car, option]
                solution.M1[option, index] = solution.M1[option, index+1] + 1
            else
                solution.M1[option, index] = solution.M1[option, index+1]
            end

            # Is there one car not reach anymore ?
            index_first_out = index + instance.RC_q[option]
            if index_first_out <= nb
                car_first_out = solution.sequence[index_first_out]
                # It had option -> sequence - 1
                if instance.RC_flag[car_first_out, option]
                    solution.M1[option, index] = solution.M1[option, index] - 1
                end
            end
        end
    end

    # Update M2 and M3 (left to right)
    for option in 1:nb_RC
        for index in 1:nb
            if index > 1
                solution.M2[option, index] = solution.M2[option, index-1]
                solution.M3[option, index] = solution.M3[option, index-1]
            else
               solution.M2[option, index] = 0
               solution.M3[option, index] = 0
            end
            # M3 is >=
            if solution.M1[option, index] >= instance.RC_p[option]
                solution.M3[option, index] += 1
                # M2 is >
                if solution.M1[option, index] > instance.RC_p[option]
                    solution.M2[option, index] += 1
                end
            end
        end
    end
end

"""
    update_matrices_new_car!(solution::Solution, position::Int, instance::Instance)

Updates `solution.M1`, `solution.M2` and `solution.M3` for a new car add at a given position.
DOES NOT update columns for sequences excluding the position. Assert made: the function is called
when we add a car in a tail-less sequence.
"""
function update_matrices_new_car!(solution::Solution, position::Int, instance::Instance)
    nb_RC = instance.nb_HPRC + instance.nb_LPRC
    car = solution.sequence[position]

    for option in 1:nb_RC
        for index in (position - instance.RC_q[option] + 1):position
            if index > 0
                # TODO: flag not raise -> index loop skipped
                # new car has option -> update M1
                if instance.RC_flag[car, option]
                    solution.M1[option, index] = solution.M1[option, index] + 1
                end

                # Columns of M2 and M3 can be update
                if index == 1
                    solution.M2[option, index] = 0
                    solution.M3[option, index] = 0
                else
                    solution.M2[option, index] = solution.M2[option, index-1]
                    solution.M3[option, index] = solution.M3[option, index-1]
                end

                # M3 is >=
                if solution.M1[option, index] >= instance.RC_p[option]
                    solution.M3[option, index] += 1
                    # M2 is >
                    if solution.M1[option, index] > instance.RC_p[option]
                        solution.M2[option, index] += 1
                    end
                end
            end
        end
    end
end

"""
    initialize_batches!(solution::Solution, instance::Instance)

"""
function initialize_batches!(solution::Solution, instance::Instance)
    solution.colors = Array{Batch, 1}(undef, 0)
    b0 = instance.nb_late_prec_day+1
    batch = Batch(b0-1,1)
    for i in 1:(b0-1)
        push!(solution.colors, batch)
    end
    batch = Batch(1, b0)
    current_color = instance.color_code[solution.sequence[b0]]
    push!(solution.colors, batch)
    for position in (b0+1):solution.length
        car = solution.sequence[position]
        if instance.color_code[car] == current_color
            batch.width += 1
            push!(solution.colors, batch)
        else
            current_color = instance.color_code[car]
            batch = Batch(1, position)
            push!(solution.colors, batch)
        end
    end
end


text:
function getPost(testResult)

    # Create a Dictionary containing the possible combinations:
    P = Dict()

    # Idea: solve by recursive calls of function,
    # constantly updating the prior --> not implemented

    deseaseFrequency = 0.001
    truePositives    = 0.99
    falsePositives   = 0.05

    ### Priors:
    P[":("] = deseaseFrequency
    P[":)"] = 1.0 - P[":("]

    ### Likelihoods
    # Positive test restuls, given status:
    P["+|:("] = truePositives
    P["+|:)"] = falsePositives
    # Negative test restuls, given status:
    P["-|:("] = 1.0 - truePositives
    P["-|:)"] = 1.0 - falsePositives

    # Update for each new measurement:
    for i in 1:length(testResult)
        # Marginal Likelihoods:
        P["+"] = P["+|:("] * P[":("] + P["+|:)"] * P[":)"]
        P["-"] = P["-|:("] * P[":("] + P["-|:)"] * P[":)"]

        ### Posteriors:
        P[":(|+"] = P["+|:("] * P[":("] / P["+"]
        P[":(|-"] = P["-|:("] * P[":("] / P["-"]
        P[":)|+"] = 1.0 - P[":(|+"]
        P[":)|-"] = 1.0 - P[":(|-"]

        # Update priors to the posteriors:
        P[":("] = P[string(":(|",testResult[i])]
        P[":)"] = 1.0 - P[":("]
    end

    # Desired Posteriors are the new Priors:
    # If no test result, then this also works, because the posterior = prior
    posteriors = [P[":("], P[":)"]]

    return posteriors
end

# Enter the received results:
myResult = "+-"

# Show the posteriors for each step when receiving test results,
# including "no result" at the start:
for i in 1:length(myResult)+1
    print("\nAccumulated Results: ")
    println(myResult[1:i-1])
    println("[P(:(|T), P(:)|T)] = ",getPost(myResult[1:i-1]))
    print("Sanity check: Sum(P) = ")
    println(sum(getPost(myResult[1:i-1])))
end

ext: jl


text:
# ---
# title: 500. Keyboard Row
# id: problem500
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Hash Table
# link: <https://leetcode.com/problems/keyboard-row/description/>
# hidden: true
# ---
#
# Given a List of words, return the words that can be typed using letters of
# **alphabet** on only one row's of American keyboard like the image below.
#
#
#
# ![](https://assets.leetcode.com/uploads/2018/10/12/keyboard.png)
#
#
#
# **Example:**
#
#
#
#     Input: ["Hello", "Alaska", "Dad", "Peace"]
#     Output: ["Alaska", "Dad"]
#
#
#
#
# **Note:**
#
#   1. You may use one character in the keyboard more than once.
#   2. You may assume the input string will only contain letters of alphabet.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*--------------------------------------------------------------------
#* Models to simulate a power plant.
#*--------------------------------------------------------------------
#* Author: Argimiro R. Secchi
#* $Id: power_plant.mso 195 2007-03-07 20:30:12Z arge $
#*-------------------------------------------------------------------
# Declaracao de tipos
type CorrenteZ
	CorrenteZ()=begin
		Ncomp=outers.Ncomp
		new(
			Corrente(),
			DanaInteger(),
			fill(fraction()),
			[:Ncomp,],
			[:z,]
		)
	end
	_base_1::Corrente
	Ncomp::DanaInteger
	z::Array{fraction}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
end
export CorrenteZ


text:
## gaussian_random_fields.jl : Gaussian random field types and generator types

"""
Abstract type `GaussianRandomFieldGenerator`

The following Gaussian random field generators are implemented:
- `Cholesky`: Cholesky factorization of the covariance matrix, exact but expensive for random fields in dimension `d` > 1
- `Spectral`: spectral (eigenvalue) decomposition of the covariance matrix, exact but expensive for random fields in dimension `d` > 1
- `KarhunenLoeve`: Karhunen-Lo\u00e8ve expansion, inexact but very efficient for "smooth" random fields when used with a low truncation dimension
- `CirculantEmbedding`: circulant embedding method, exact and efficient, but can only be used for random fields on structured grids

See also: [`Cholesky`](@ref), [`Spectral`](@ref), [`KarhunenLoeve`](@ref), [`CirculantEmbedding`](@ref)
"""
abstract type GaussianRandomFieldGenerator end

"""
    GaussianRandomField([mean,] cov, generator, pts...)
    GaussianRandomField([mean,] cov, generator, nodes, elements)

Compute a Gaussian random field with mean `mean` and covariance structure `cov` defined in the points `pts`, and computed using the Gaussian random field generator `generator`.

# Examples
```jldoctest label2
julia> cov = CovarianceFunction(2, Matern(.3, 1))
2d Matérn covariance function (λ=0.3, ν=1.0, σ=1.0, p=2.0)

julia> pts = pts = range(0, stop=1, length=51)
0.0:0.02:1.0

julia> mean = fill(π, (51, 51))
[...]

julia> grf = GaussianRandomField(mean, cov, Cholesky(), pts, pts)
Gaussian random field with 2d Matérn covariance function (λ=0.3, ν=1.0, σ=1.0, p=2.0) on a 51×51 structured grid, using a Cholesky decomposition

```
If no `mean` is specified, a zero-mean Gaussian random field is assumed.
```jldoctest label2
julia> grf = GaussianRandomField(cov, Cholesky(), pts, pts)
Gaussian random field with 2d Matérn covariance function (λ=0.3, ν=1.0, σ=1.0, p=2.0) on a 51×51 structured grid, using a Cholesky decomposition

```
The  Gaussian random field generator `generator` can be `Cholesky()`, `Spectral()`, `KarhunenLoeve(n)` (where `n` is the number of terms in the expansion), or `CirculantEmbedding()`. The points `pts` can be specified as arguments of type `AbstractVector`, in which case a tensor (Kronecker) product is assumed, or as a Finite Element mesh with node table `nodes` and element table `elements`.
```jldoctest label2
julia> grf = GaussianRandomField(cov, KarhunenLoeve(500), pts, pts)
Gaussian random field with 2d Matérn covariance function (λ=0.3, ν=1.0, σ=1.0, p=2.0) on a 51×51 structured grid, using a KL expansion with 500 terms

julia> exponential_cov = CovarianceFunction(2, Exponential(.1))
2d exponential covariance function (λ=0.1, σ=1.0, p=2.0)

julia> grf = GaussianRandomField(exponential_cov, CirculantEmbedding(), pts, pts)
Gaussian random field with 2d exponential covariance function (λ=0.1, σ=1.0, p=2.0) on a 51×51 structured grid, using a circulant embedding

```
Separable Gaussian random fields can be defined using `SeparableCovarianceFunction`.
```jldoctest label2
julia> separable_cov = SeparableCovarianceFunction(Exponential(.1), Exponential(.1))
2d separable covariance function [ exponential (λ=0.1, σ=1.0, p=2.0), exponential (λ=0.1, σ=1.0, p=2.0) ]

julia> grf = GaussianRandomField(separable_cov, KarhunenLoeve(500), pts, pts)
Gaussian random field with 2d separable covariance function [ exponential (λ=0.1, σ=1.0, p=2.0), exponential (λ=0.1, σ=1.0, p=2.0) ] on a 51×51 structured grid, using a KL expansion with 500 terms

julia> plot_eigenfunction(grf, 3)
[...]

```
We also offer support for anisotropic random fields.
```jldoctest label2
julia> anisotropic_cov = CovarianceFunction(2, AnisotropicExponential([500 400; 400 500]))
2d anisotropic exponential covariance function (A=[500 400; 400 500], σ=1.0)

julia> grf = GaussianRandomField(anisotropic_cov, CirculantEmbedding() , pts, pts)
Gaussian random field with 2d anisotropic exponential covariance function (A=[500 400; 400 500], σ=1.0) on a 51×51 structured grid, using a circulant embedding

julia> heatmap(grf)
[...]

```
For irregular domains, specify the points as matrices containing the `nodes` and `elements` of a finite element mesh. To compute the value of the random field at the element centers, use the optional keyword `mode="center"`.
```jldoctest label2
julia> nodes, elements = Lshape()
[...]

julia> grf = GaussianRandomField(cov, Spectral(), nodes, elements)
Gaussian random field with 2d Matérn covariance function (λ=0.3, ν=1.0, σ=1.0, p=2.0) on a mesh with 998 points and 1861 elements, using a spectral decomposition

```

Samples from the random field can be computed using the `sample` function.
```jldoctest label2
julia> sample(grf)
[...]

```
See also: [`Cholesky`](@ref), [`Spectral`](@ref), [`KarhunenLoeve`](@ref), [`CirculantEmbedding`](@ref), [`sample`](@ref)
"""
struct GaussianRandomField{G<:GaussianRandomFieldGenerator,C<:AbstractCovarianceFunction,
                           P,M,D}
    mean::M
    cov::C
    pts::P
    data::D
end

function GaussianRandomField(mean::Array{<:Real}, cov::CovarianceFunction{d}, method::GaussianRandomFieldGenerator, pts::Vararg{AbstractVector,d}; kwargs...) where d
    size(mean) == length.(pts) || throw(DimensionMismatch("size of the mean does not correspond to the dimension of the points"))

    method isa CirculantEmbedding && !(pts isa NTuple{d,AbstractRange}) &&
        throw(ArgumentError("can only use circulant embedding on a regular grid, supply ranges for pts"))

    method isa Union{CirculantEmbedding,KarhunenLoeve} && any(pt -> length(pt) < 2, pts) &&
        throw(ArgumentError("must have at least 2 points in each direction to use circulant embedding or KL expansion"))

    _GaussianRandomField(mean, cov, method, pts...; kwargs...)
end

# zero-mean GRF
GaussianRandomField(cov::CovarianceFunction{d}, method::GaussianRandomFieldGenerator, pts::Vararg{AbstractVector,d}; kwargs...) where d = GaussianRandomField(zeros(eltype(cov), length.(pts)), cov, method, pts...; kwargs...)

# constant mean GRF
GaussianRandomField(mean::Real, cov::CovarianceFunction{d}, method::GaussianRandomFieldGenerator, pts::Vararg{AbstractVector,d}; kwargs...) where d = GaussianRandomField(fill(convert(eltype(cov), mean), length.(pts)), cov, method, pts...; kwargs...)

# obtain generator used to compute random field
generator(::GaussianRandomField{G}) where G = G()

"""
	sample(grf)
        sample(grf[, xi])

Take a sample from the Gaussian random field `grf` using the (optional) normally distributed random numbers `xi`. The vector`xi` must have appropriate length.

# Examples
```jldoctest
julia> cov = CovarianceFunction(2, Whittle(.1))
2d Whittle covariance function (λ=0.1, σ=1.0, p=2.0)

julia> pts = pts = range(0, stop=1, length=51)
0.0:0.02:1.0

julia> grf = GaussianRandomField(cov, CirculantEmbedding(), pts, pts)
Gaussian random field with 2d Whittle covariance function (λ=0.1, σ=1.0, p=2.0) on a 51×51 structured grid, using a circulant embedding

julia> sample(grf)
[...]

julia> sample(grf, xi=randn(randdim(grf)))
[...]
```
See also: [`GaussianRandomField`](@ref), [`Matern`](@ref), [`CovarianceFunction`](@ref), [`CirculantEmbedding`]
"""
function sample(grf::GaussianRandomField; xi::AbstractArray{<:Real}=randn(randdim(grf)))
    length(xi) == prod(randdim(grf)) || throw(DimensionMismatch("length of random points vector must be equal to $(randdim(grf))"))
    _sample(grf, xi)
end

function Base.show(io::IO, grf::GaussianRandomField)
    print(io, "Gaussian random field with ", grf.cov, " on a ")
    showpoints(io, grf.pts)
    print(io, ", using a ", generator(grf))
end

function showpoints(io::IO, points)
    d = length(points)
    if d == 1
        print(io, length(points[1]), "-point")
    else
        print(io, length(points[1]))
        @inbounds for i in 2:d
            print(io, "×", length(points[i]))
        end
    end
    print(io, " structured grid")
end

"""
    randdim(grf)

Returns the number of random numbers used to sample from the Gaussian random field `grf`.

# Examples
```jldoctest
julia> cov = CovarianceFunction(2, Whittle(.1))
2d Whittle covariance function (λ=0.1, σ=1.0, p=2.0)

julia> pts = pts = range(0, stop=1, length=51)
0.0:0.02:1.0

julia> grf = GaussianRandomField(cov, KarhunenLoeve(200), pts, pts)
Gaussian random field with 2d Whittle covariance function (λ=0.1, σ=1.0, p=2.0) on a 51×51 structured grid, using a KL expansion with 200 terms

julia> randdim(grf)
200

```
"""
randdim(grf::GaussianRandomField) = 0


text:
include("./0_setup.jl")

N = 250_000_000
K = 100

# srand(1);
# id6 = rand(Int32(1):Int32(round(N/K)), N);
# v1 =  rand(Int32(1):Int32(5), N);
# @elapsed sumby_multi_rg(id6, v1)

function bench_sumby_multi_van()
    srand(1);
    id6 = rand(Int32(1):Int32(round(N/K)), N);
    v1 =  rand(Int32(1):Int32(5), N);
    @elapsed sumby_multi_van(id6, v1)
end

function bench_sumby()
    srand(1)
    id6 = rand(Int32(1):Int32(round(N/K)), N)
    v1 =  rand(Int32(1):Int32(5), N)
    @elapsed sumby(id6, v1)
end

function bench_sumby_radixgroup()
    srand(1)
    id6 = rand(Int32(1):Int32(round(N/K)), N)
    v1 =  rand(Int32(1):Int32(5), N)
    @elapsed sumby_radixgroup(id6, v1)
end

function bench_sumby_radixsort()
    srand(1)
    id6 = rand(Int32(1):Int32(round(N/K)), N)
    v1 =  rand(Int32(1):Int32(5), N)
    @elapsed sumby_radixsort(id6,v1)
end

a  = [bench_sumby_multi_rs() for i=1:5]
b = [bench_sumby_radixgroup() for i=1:5]

1-mean(a)/mean(b)
mean(a)
mean(b)

srand(1)
id6 = rand(Int32(1):Int32(round(N/K)), N)
v1 =  rand(Int32(1):Int32(5), N)
@time res = sumby_multi_rs(by, v1)

#[bench_sumby_multi_van() for i=1:5]
#[bench_sumby() for i=1:5]
#[bench_sumby_radixsort() for i=1:5]


text:
using Poltergeist, Plots, ApproxFun; pyplot()
d = Interval(0,1);
M1(x) = 2x+sin(2pi*x)/6; M2(x) = 2-2x
M = MarkovMap([M1,M2],[0..0.5,0.5..1]);
L = Transfer(M)

scatter(eigvals(L,80),legend=true,grid=false,legend=false,
title = "Eigenvalues of transfer operator")
xlims!(-1.1,1.1)
ylims!(-0.5,0.5)
savefig("eigvals.pdf")
run(`sips -s format png eigvals.pdf --out eigvals.png`)

K = SolutionInv(L);
ρ = acim(K);
dρ = linearresponse(K,sinpi)
ϵ = 0.05
Mϵ = perturb(M,sinpi,ϵ)
plot(ρ,legend=true,grid=false,label="\$\\rho\_\{0\}\$",
title = "Linear response");
plot!(ρ+ϵ*dρ,label="\$\\rho\_\{\\epsilon\}\$ estimate")
plot!(acim(Mϵ), label="\$\\rho\_\{\\epsilon\}\$ actual")
xlabel!("\$x\$"); ylabel!("\$\\rho(x)\$")
# xlabel!("x"); ylabel!("ρ(x)")
ylims!(0.,1.8)
savefig("acim.pdf")
run(`sips -s format png acim.pdf --out acim.png`)


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/NodeNumbering.jl/blob/master/LICENSE

"""
    node_degrees(adjacency::Dict{Int, Vector{Int}})

Calculate the degrees for the adjacency graph nodes. Node degree = number of nodes a certain node is adjacent to.

# Examples

```jldoctest
julia> adjacency = Dict{Int, Vector{Int}}(
                   1 => [3, 8, 9],
                   2 => [3, 8, 7],
                   3 => [1, 2],
                   4 => [8, 9],
                   5 => [7, 8],
                   6 => [2, 7],
                   7 => [5, 2, 6],
                   8 => [1, 2, 4, 5],
                   9 => [1, 4]);

julia> node_degrees(adjacency::Dict{Int, Vector{Int}})
Dict{Int64,Int64} with 9 entries:
  7 => 3
  4 => 2
  9 => 2
  2 => 3
  3 => 2
  5 => 2
  8 => 4
  6 => 2
  1 => 3

```

# References

* Wikipedia contributors. "Degree (graph theory)". Wikipedia, The Free Encyclopedia. 24 Nov. 2016. Web. 18 Jul. 2017. https://en.wikipedia.org/wiki/Degree_(graph_theory)
"""
function node_degrees(adjacency::Dict{Int, Vector{Int}})
    degrees = Dict{Int, Int}()
    for (k, v) in adjacency
        degree = length(v)
        degrees[k] = degree
    end
    return degrees
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ApplicationGatewayAvailableRequestHeadersResult <: SwaggerModel

    function ApplicationGatewayAvailableRequestHeadersResult(;)
        o = new()
        o
    end
end # type ApplicationGatewayAvailableRequestHeadersResult

const _property_map_ApplicationGatewayAvailableRequestHeadersResult = Dict{Symbol,Symbol}()
const _property_types_ApplicationGatewayAvailableRequestHeadersResult = Dict{Symbol,String}()
Base.propertynames(::Type{ ApplicationGatewayAvailableRequestHeadersResult }) = collect(keys(_property_map_ApplicationGatewayAvailableRequestHeadersResult))
Swagger.property_type(::Type{ ApplicationGatewayAvailableRequestHeadersResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ApplicationGatewayAvailableRequestHeadersResult[name]))}
Swagger.field_name(::Type{ ApplicationGatewayAvailableRequestHeadersResult }, property_name::Symbol) =  _property_map_ApplicationGatewayAvailableRequestHeadersResult[property_name]

function check_required(o::ApplicationGatewayAvailableRequestHeadersResult)
    true
end

function validate_property(::Type{ ApplicationGatewayAvailableRequestHeadersResult }, name::Symbol, val)
end


text:
using LinearAlgebra
using SparseArrays

BLAS.set_num_threads(1)

function solvematrix(A::SparseMatrixCSC)
    N = size(A, 1)
    e = ones(N, 1)
    b = A * e

    x, time, space = @timed A \ b
    error = norm(x - e) / norm(e)

    return (N, time, error, space)
end

max_stars_repo_name: saiteki-kai/direct-methods

max_issues_repo_name: saiteki-kai/direct-methods

max_forks_repo_name: saiteki-kai/direct-methods


text:
using Documenter
using LaserTypes

makedocs(
    sitename = "LaserTypes",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    pages = [
        "index.md",
        "setup.md",
        "Laser profiles" => [
            "Spatial profiles" => [
                "Gaussian profile" => "gauss.md",
                "Laguerre-Gauss profile" => "laguerre-gauss.md"
            ],
            "temporal-profiles.md"
        ],
        "covariant.md",
        "advanced.md"
    ],
    modules = [LaserTypes]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/SebastianM-C/LaserTypes.jl",
    push_preview = true
)


text:
"""
    CPIDataBase

Librería base para tipos y funcionalidad básica para manejo de los datos
desagregados del IPC a nivel de república.
"""
module CPIDataBase

    using Dates
    using DataFrames
    using PrettyTables

    # Exportar tipos
    export IndexCPIBase, VarCPIBase, FullCPIBase
    export CountryStructure, UniformCountryStructure, MixedCountryStructure

    # Exportar funciones
    export capitalize, varinterm, varinteran,
        capitalize!, varinterm!, varinteran!,
        periods, infl_periods, infl_dates,
        getunionalltype

    # Exportar tipos para implementar nuevas funciones de inflación
    export InflationFunction, EnsembleInflationFunction
    export EnsembleFunction, CombinationFunction
    export InflationEnsemble, InflationCombination # alias de los 2 anteriores
    export components # componentes de una InflationCombination
    export num_measures, weights, measure_name, measure_tag, params

    # Exportar tipos necesarios para especificar tipos de los resultados
    export CPIIndex, CPIVarInterm

    # Función básica de inflación
    export InflationTotalCPI

    # Definición de tipos para bases del IPC
    include("CPIBase.jl")
    include("CountryStructure.jl")

    # Operaciones básicas
    include("utils/capitalize.jl")
    include("utils/varinterm.jl")
    include("utils/varinteran.jl")

    # Estructura básica para medidas de inflación
    include("inflation/InflationFunction.jl")
    include("inflation/EnsembleFunction.jl")
    include("inflation/CombinationFunction.jl")

    # Medida de inflación básica
    include("inflation/InflationTotalCPI.jl")

    # Funciones de utilidad
    export getdates
    include("utils/utils.jl")

    ##  ------------------------------------------------------------------------
    #   Funcionalidades de árboles de cómputo del IPC
    #   ------------------------------------------------------------------------
    using AbstractTrees
    import AbstractTrees: children, printnode
    # Tipos básicos, un gasto básico y una estructura de grupo
    export Item, Group
    # Funciones para construir a partir de estructura de códigos
    export get_cpi_tree, cpi_tree_nodes
    # Operaciones básicas:
    # find_tree : Encontrar un nodo en el árbol con un código especificado
    # compute_index : Computar el índice de cualquier nodo en el árbol con la fórmula del IPC
    # compute_index! : Similar a compute_index pero utiliza una caché de diccionario
    export find_tree, compute_index, compute_index!
    export children, print_tree  # reexport from AbstractTrees

    # Estructura envolvente CPITree de FullCPIBase y estructura anidada Group
    export CPITree

    include("tree/CPItree.jl")

    ##  ------------------------------------------------------------------------
    #   Submódulo de pruebas
    #   ------------------------------------------------------------------------
    # Submódulo con funciones relacionadas con los tipos de este paquete para
    # realizar pruebas en paquetes que extiendan la funcionalidad. Este módulo
    # no se exporta por defecto, requiere carga explícita (e.g using
    # CPIDataBase.TestHelpers)
    module TestHelpers
        using Dates, ..CPIDataBase

        export getrandomweights, getbasedates
        export getzerobase, getrandombase
        export getzerocountryst, getrandomcountryst

        include("helpers/test_helpers.jl")
    end

end


text:
using Test
using AbstractTrees

using FlatZincParser
using FlatZincParser: TokenStream, AbstractNode, DataNode, Node, next!, reset!, context, peek
using FlatZincParser: match!, match_many!, match_token, match_many

@testset "parser testes" begin
    @testset "match matchmany" begin
        ts = TokenStream(tokenize("1.2 2.6 0.2 ;"))
        childs = AbstractNode[]
        match!(childs, ts, :float_literal)
        match!(childs, ts, :float_literal)
        match!(childs, ts, :float_literal)
        @test_throws ParsingError match!(childs, ts, :float_literal)
        peek(ts).type === :semicolon

        ts = TokenStream(tokenize("1.2 2.6 0.2 ;"))
        @test match_many!(childs, ts, :float_literal) == 3
        peek(ts).type === :semicolon

        ts = TokenStream(tokenize("1.2, 2.6, 0.2 ;"))
        match_many(ts, :float_literal; delimiter=:comma)
    end

    @testset "matchtoken" begin
        ts = TokenStream("predicate foo bar;")
        match_token(ts, :keyword, "predicate")
        match_token(ts, :identifier, "foo")
        @test_throws ParsingError match_token(ts, :identifier, "foo")
        match_token(ts, :identifier)
        match_token(ts, :semicolon)
    end

    @testset "predicate_item" begin
        ts = TokenStream("predicate cumulativeChoco (array [int] of var int: s, array [int] of var int: d,array [int] of var int: r,var int: b);")
        match(ts, :predicate_item)
    end

    @testset "range and set literal" begin
        ts = TokenStream("1..2;")
        print_tree(match(ts, :set_literal))
        ts = TokenStream("1.2..2.5;")
        print_tree(match(ts, :set_literal))
        ts = TokenStream("{1.2,2.5,};")
        print_tree(match(ts, :set_literal))
        ts = TokenStream("{1,2,3,4,5,6,7};")
        print_tree(match(ts, :set_literal))
    end

    @testset "par_decl_item" begin
        streams = [
            TokenStream("array [1..100] of int: duration = [88,740,752,503,536,537,300,668,332,693,30,249,673,391,51,386,328,313,103,190,741,800,425,552,749,31,632,690,530,706,131,377,505,656,654,720,757,90,331,450,103,276,571,782,568,772,106,198,183,800,705,140,107,542,597,97,580,59,325,336,76,482,776,428,237,433,701,220,478,102,32,617,454,360,541,107,609,730,107,274,589,305,249,96,365,651,163,202,560,571,104,740,720,437,230,157,145,318,531,481];"),
            TokenStream("array [1..2] of int: X_INTRODUCED_1952_ = [1,-1];"),
        ]
        match(reset!(streams[1]), :par_decl_item)
        match(reset!(streams[2]), :par_decl_item)
    end

    @testset "array literal" begin
        str = "[X_INTRODUCED_0_,X_INTRODUCED_1_,X_INTRODUCED_2_,X_INTRODUCED_3_]"
        match(TokenStream(str), :array_literal)
        match(TokenStream("[1..4]"), :expr)
        match(TokenStream("[1..4]"), :array_literal)
        match(TokenStream("1..4"), :basic_literal_expr)
    end

    @testset "var_decl_item" begin
        streams = [
            TokenStream("var 1..4: X_INTRODUCED_0_;"),
            TokenStream("array [1..4] of var int: queens:: output_array([1..4]) = [X_INTRODUCED_0_,X_INTRODUCED_1_,X_INTRODUCED_2_,X_INTRODUCED_3_];"),
        ]
        match(reset!(streams[1]), :var_decl_item)
        match(reset!(streams[2]), :var_decl_item)
    end

    @testset "constraint_item" begin
        stream = TokenStream("constraint int_lin_ne(X_INTRODUCED_4_,[X_INTRODUCED_0_,X_INTRODUCED_1_],-1);")
        match(reset!(stream), :constraint_item) |> print_tree

        str = "X_INTRODUCED_4_,[X_INTRODUCED_0_,X_INTRODUCED_1_],-1"
        stream = TokenStream(str)
        match(stream, :expr)
        match_token(stream, :comma)
        match(stream, :expr) |> print_tree
        match_token(stream, :comma)
        match(stream, :expr) |> print_tree
    end

    @testset "solver_item" begin
        streams = [TokenStream("solve :: seq_search([int_search(X_INTRODUCED_2947_,first_fail,indomain_min,complete),int_search(X_INTRODUCED_2946_,input_order,indomain_min,complete)]) minimize objective;"),
                   TokenStream("solve :: int_search(queens,first_fail,indomain_min,complete) satisfy;"),
                   TokenStream("solve :: int_search(X_INTRODUCED_51094_,input_order,indomain_min,complete) satisfy;")]

        match(reset!(streams[1]), :solve_item)
        match(reset!(streams[2]), :solve_item)
        match(reset!(streams[3]), :solve_item)
    end

    @testset "parse files" begin
        small_file = joinpath(@__DIR__,"files", "queens3_4.fzn")
        big_file = joinpath(@__DIR__,"files", "2018_test-scheduling_t100m10r3-2.fzn")
        # huge_file = joinpath(@__DIR__,"files", "oocsp_racks_030_f7.fzn")

        println("Tokenize small file")
        @time tokens = open(tokenize, small_file);
        println("Parse small file")
        @time match(TokenStream(tokens), :model);

        println("Tokenize big file")
        @time tokens = open(tokenize, big_file);
        println("Parse big file")
        @time match(TokenStream(tokens), :model);

        @info "small file asyn=true"
        @time m1 = parsefile(small_file; async=true);
        @info "small file asyn=false"
        @time m2 = parsefile(small_file; async=false);
        @test m1 == m2

        @info "big file asyn=true, spawn=true"
        @time m1 = parsefile(big_file; async=true, spawn=true);
        @info "big file asyn=true, spawn=false"
        @time m2 = parsefile(big_file; async=true, spawn=false);
        @info "big file asyn=false"
        @time m3 = parsefile(big_file; async=false);
        @test m1 == m2 == m3

        @info "small file parallel=true"
        @time m1 = parsefile(small_file, parallel=true);
        @info "small file parallel=false"
        @time m2 = parsefile(small_file, parallel=false);
        @test m1 == m2

        @info "big file parallel=true"
        @time m1 = parsefile(big_file, parallel=true);
        @info "big file parallel=false"
        @time m2 = parsefile(big_file, parallel=false);
        @test m1 == m2
    end
end


text:
function CMM_driver(modeloutputpathandfilename,dtoutput)

include("user_settings.jl")
#import Pkg.add("user_settings.jl")
# Define User_Settings

#using user_settings

# Define Constants

include("constants.jl")

# Define Grid

[ zu,zw,xu,xs ] = griddef(nz,nx,dz,dx,zu,zw,xu,xs);

# Initialize the Base State
[ tb,qb,qbs,tbv,pisfc, pib, rhou, rhow,rhb,ub,um,u ] = ...
    base( profile_method,nx,nz,dz,psurf,qsurf,q4km,ztr,temptr,ttr,tsurf,p0,cp,g,...
    rd,xk,c_v,zu,rhow,rhou,tb,tbv,qb,qbs,rhb,pib,ub,um,u);

# Set Initial Conditions

[ th,thm,pim,pic,pprt] = cmm_init(xs,g,nx,nz,zu,dx,dz,zcnt,xcnt,radz,radx,trigpi,...
    cp,rhou,tb,qb,thermamp,th,thm,pim,pic,pprt,bubble_switch);

# Integrate the Model

integ(tbv,pib,p0,lv,rd,ub, g, cp, c_sound,rhow,rhou,tb, zu, zw,xu,xs,x_dist_in, ...
    latdmpcoef,raydmpcoef,raydmpz,trigpi,qb,um,u,up,wm,w,wp,thm,th,thp,pim,pic,pip,pprt,...
    qvtot,qvm,qv,qvp,qcm,qc,qcp,qrainm,qrain,qrainp,dx,dz,nt,nz,nx,dt,asscoef,cmixh,cmixv,...
    qc0,k1,k2,thvm,thv,thvp,modeloutputpathandfilename,dtoutput);

# Save Model Data to Output File

# save (modeloutputpathandfilename)

end

max_stars_repo_name: jeffreycunn/toy-atmospheric-model

max_issues_repo_name: jeffreycunn/toy-atmospheric-model

max_forks_repo_name: jeffreycunn/toy-atmospheric-model


text:
CLUST_FOR_OPT=ENV["CLUST_FOR_OPT"]
push!(LOAD_PATH, normpath(joinpath(CLUST_FOR_OPT,"src"))) #adds the location of ClustForOpt to the LOAD_PATH
using ClustForOpt
using TimeWarp
 using PyPlot
 plt = PyPlot
plt.close()
# example not working, somehow input is not in correct format
# x = [1,2,2,3,3,4]
# y = [1,3,4]
# z = [1,2,2,4]
# avg,result = dba([x,y,z],ClassicDTW(),iterations=15)  # make sure to put in ClassicDTW()
# println(avg)

 ######## DATA INPUT ##########

 # region
region = "GER"

# read in original data
data_orig_daily = load_pricedata(region)

println("data loaded")

seq1 = data_orig_daily[:,1]
seq2 = data_orig_daily[:,2]
seq = data_orig_daily[:,1:365]  # do not load as sequence

# test dtw
D = dtw(seq1,seq2) # dtw can be fed with normal array/vector or Sequence
println("D ",D)
E = fastdtw(seq1,seq2,24) # fastdtw can only be fed with array/vector
println("E ",E)

#dtwplot(Sequence(seq1),Sequence(seq2))  # dtwplot has to be fed with sequence, somehow does not work currently, shows empty plot
#show()

#avg, result = dba(seq,ClassicDTW(),iterations=15) #always include DTWMethod. Put initial one in example or actual code
#println(avg)

##
n_seq=10
tic()
mean_dtw,result = dba(seq[:,1:n_seq],ClassicDTW(),iterations=15,rtol=1e-5,show_progress=true)
toc()

mean_euc = mean(seq[:,1:n_seq],2)

figure()
plt.plot(seq[:,1:n_seq],color="0.75")
plt.plot(mean_dtw,label="dtw",color="red")
plt.plot(mean_euc,label="euclidean",color="blue")
plt.legend()

##########################
# normalized clustering hourly

seq_norm, hourly_mean, hourly_sdv = z_normalize(data_orig_daily[:,1:n_seq],hourly=true)
tic()
mean_dtw_norm,result_norm = dba(seq_norm[:,1:n_seq],ClassicDTW(),iterations=15,rtol=1e-5,show_progress=true)
toc()

mean_euc_norm = mean(seq_norm[:,1:n_seq],2)

mean_dtw = undo_z_normalize(reshape(mean_dtw_norm,size(mean_dtw_norm)[1],1),hourly_mean,hourly_sdv)
mean_euc = undo_z_normalize(mean_euc_norm,hourly_mean,hourly_sdv)

figure()
plt.plot(seq[:,1:n_seq],color="0.75")
plt.plot(mean_dtw,label="dtw",color="red")
plt.plot(mean_euc,label="euclidean",color="blue")
plt.title("normalized")
plt.legend()

#######################
# normalized and include warping window
rad_sc = 2 # sakoe chiba band radius
rmin,rmax = sakoe_chiba_band(rad_sc,24)

seq_norm, hourly_mean, hourly_sdv = z_normalize(data_orig_daily[:,1:n_seq],hourly=true)
tic()
mean_dtw_norm,result_norm = dba(seq_norm[:,1:n_seq],ClassicDTW(),iterations=15,rtol=1e-5,show_progress=true,i2min=rmin,i2max=rmax)
toc()

mean_euc_norm = mean(seq_norm[:,1:n_seq],2)

mean_dtw = undo_z_normalize(reshape(mean_dtw_norm,size(mean_dtw_norm)[1],1),hourly_mean,hourly_sdv)
mean_euc = undo_z_normalize(mean_euc_norm,hourly_mean,hourly_sdv)

figure()
plt.plot(seq[:,1:n_seq],color="0.75")
plt.plot(mean_dtw,label="dtw",color="red")
plt.plot(mean_euc,label="euclidean",color="blue")
plt.title("normalized")
plt.legend()

show()


text:
# Autogenerated wrapper script for GoogleTest_jll for x86_64-w64-mingw32-cxx03
export libgtest, libgtest_main

JLLWrappers.@generate_wrapper_header("GoogleTest")
JLLWrappers.@declare_library_product(libgtest, "libgtest.dll")
JLLWrappers.@declare_library_product(libgtest_main, "libgtest_main.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libgtest,
        "bin\\libgtest.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgtest_main,
        "bin\\libgtest_main.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
abstract type DiffusionMethod end
struct ExplicitDiffusion <: DiffusionMethod end
struct CrankNickolsonDiffusion <: DiffusionMethod end
struct ADIDiffusion <: DiffusionMethod end

function perform_step!(model::DiffusionModel, method::ExplicitDiffusion, cache::ExplicitDiffusionCache)
    u  = cache.u
    un = cache.uprev
    nx = model.grid.nx
    ny = model.grid.ny
    Δx = model.grid.Δx
    Δy = model.grid.Δy
    Δt = model.Δt
    ν  = model.parameters.ν

    copy!(un,u)
    @inbounds for j in 2:ny-1, i in 2:nx-1
        u[j,i] = un[j,i] +
            ν*Δt/(Δx*Δx)*(un[j  ,i+1]-2*un[j,i]+un[j  ,i-1]) +
            ν*Δt/(Δy*Δy)*(un[j+1,i  ]-2*un[j,i]+un[j-1,i  ])
    end

    return nothing
end
function perform_step!(model::DiffusionModel, method::CrankNickolsonDiffusion, cache::CNDiffusionCache)
    u   = cache.u
    un  = cache.uprev
    rhs = cache.rhs
    a1  = cache.a1
    a2  = cache.a2

    copy!(un,u)

    unV = @view un[2:end-1,2:end-1]
    unVC = @view unV[:]
    uV = @view u[2:end-1,2:end-1]
    uVC = @view uV[:]
    rhsV = @view rhs[2:end-1,2:end-1]
    rhsVC = @view rhsV[:]

    mul!(rhsVC, a2, unVC)
    uVC .= a1\rhsV[:]

    return nothing
end
function perform_step!(model::DiffusionModel, method::ADIDiffusion, cache::ADIDiffusionCache)
    u      = cache.u
    uhalf  = cache.uhalf
    un     = cache.uprev
    rhs    = cache.rhs
    a1x    = cache.a1x
    a1y    = cache.a1y
    a2x    = cache.a2x
    a2y    = cache.a2y
    nx, ny = model.grid.nx, model.grid.ny

    copy!(un,u)

    un_i(i) = @view un[2:end-1,i]
    u_i(i) = @view u[2:end-1,i]
    rhs_j(j) = @view rhs[j,2:end-1]
    rhs_i(i) = @view rhs[2:end-1,i]
    uhalf_j(j) = @view uhalf[j,2:end-1]

    # Step one
    for i = 2:nx-1
        mul!(rhs_i(i), a2y, un_i(i))
    end
    for j = 2:ny-1
        ldiv!(uhalf_j(j), a1x, rhs_j(j))
    end

    # Step two
    for j = 2:ny-1
        mul!(rhs_j(j), a2x, uhalf_j(j))
    end
    for i = 2:nx-1
        ldiv!(u_i(i), a1y, rhs_i(i))
    end
    return nothing
end


text:
using HaploSoup
using Documenter

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

makedocs(;
    modules=[HaploSoup],
    authors="Roland Laboulaye <rlaboulaye@gmail.com> and contributors",
    repo="https://github.com/rlaboulaye/HaploSoup.jl/blob/{commit}{path}#{line}",
    sitename="HaploSoup.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://rlaboulaye.github.io/HaploSoup.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/rlaboulaye/HaploSoup.jl",
)


text:
using Modia
using Base.Test

@testset "RunExamples" begin

Modia.ModiaLogging.resetTestStatus()

global figure = 1
include("CurrentController.jl")
include("Rectifier.jl")
include("CauerLowPassFilter.jl")
include("LinearSystems.jl")
include("SynchronousExamples.jl")
include("ElectricalVehicleAndCharger.jl")
include("CollidingBalls.jl")

Modia.ModiaLogging.printTestStatus()

end # testset
nothing


text:
function reconstruction(u0, f, p; alg, problem_kwargs, solver_kwargs, isshow=true)
    prob = GalacticOptim.OptimizationProblem(f, u0, p; problem_kwargs...)
    trace = Trace(prob; isshow)
    sol = solve(prob, alg; cb=trace, solver_kwargs...)
    trace(sol.u, f(sol.u, p))
    trace
end

function pp_reconstruction(fx, x0; alg=LBFGS(), problem_kwargs, solver_kwargs, isshow=true, Ns, σs)
    e = one(eltype(x0))

    ux = rand(eltype(x0), length(x0))
    traces_x = []
    for (N, σ) in zip(Ns, σs)
        J, L, K = 3, 8, 2

        wp_x = WaveletParams(e / 2, N, J, L, K, σ)

        p = get_params(x0, wp_x)

        t_x = reconstruction(ux, fx, p; alg, problem_kwargs, solver_kwargs, isshow)

        ux .= t_x[end]
        push!(traces_x, t_x)

        @info "x reconstruction" N σ improvement=t_x.way_l[end] / t_x.way_l[1]
    end

    traces_x
end

function scaled_colored_reconstruction(fy, x0, y0, ux, s::Int = 2; alg=LBFGS(), problem_kwargs, solver_kwargs, isshow=true, Ns = [128, 128, 128, 128], σs = [one(eltype(y0)) / 2^n for n in 4:7])
    e = one(eltype(y0))

    uy = s * rand(eltype(y0), s * s * length(y0))
    traces_y = []
    for (N, σ) in zip(Ns, σs)
        J, L, K = 3, 8, 2

        wp_x = WaveletParams(e / 2, N, J, L, K, σ)
        wp_y = WaveletParams(e / 2, N, J, L, K, σ)
        wp_y_cross = WaveletParams(e / 2, N, J, L, K, 4σ)

        p = get_colored_params(x0, y0, wp_x, wp_y, wp_y_cross; s)

        py = (; p, ux)
        t_y = reconstruction(uy, fy, py; alg, problem_kwargs, solver_kwargs, isshow)

        uy .= t_y[end]
        push!(traces_y, t_y)

        @info "y reconstruction" N σ improvement=t_y.way_l[end] / t_y.way_l[1]
    end

    traces_y
end


text:
# File: histogram.jl
# Author:
# Date: 29-June-2020

function histogram()
    # import XL file
    [dates,returns] = getXLData('hfma_matlab_data.xlsx','CTA Index');

    mn = abs(min(returns)); # get min percent absolute
    mx = abs(max(returns)); # get max percent absolute

    figure; # create figure
    y=returns; # set returns to y

    x=-4.5:0.25:4.5;  # bins -4.5 to 4.5 intervals 0.25
    freq = histc(y,x);  # count frequency per bin
    freq = y.count(x);
    rfreq = (freq/sum(freq))*100; # relative frequency

    bar(x,rfreq);  # plot bar chart
    xlabel('Bin'); # add x label
    ylabel('Rel. Freq. (# )'); # add y label
end


text:
module Access_Test

using Compat
using Compat.Test
using IntervalArithmetic
using EAGOBranchBound

S = BnBSolver()
B = BnBModel([Interval(1.0,2.0)])
#B1 = BnBModel()
tL = 3.0
tU = 5.0
X1 = [Interval(1.0,2.0)]
X2 = [Interval(3.0,5.0)]
X = [Interval(3.0,5.0)]
pos = 4
EAGOBranchBound.BM_breadth!(S,B,tL,tU,X1,X2,pos)
EAGOBranchBound.BM_Single!(S,B,tL,tU,X,pos)
B.soln = [7.0]
B.UBDg = 14.7

#@test B1.Init_Box == [Interval(0,1)]
@test getsolution(B) == [7.0]
@test getobjval(B) == 14.7
@test getobjbound(B) == 14.7
@test getfeasibility(B) == false
@test LBDtime(B) == 0.0
@test UBDtime(B) == 0.0
end


text:
# This file is part of Kpax3. License is MIT.

function computelocalmode!(v::Vector{Int},
                           n1s::Matrix{Float64},
                           C::Matrix{UInt8},
                           cl::Vector{Int},
                           k::Int,
                           logpC::Vector{Float64},
                           priorC::PriorColPartition)
  fill!(logpC, 0.0)

  m = size(n1s, 2)

  tmp1 = zeros(Float64, 4)
  tmp2 = zeros(Float64, 2, k)
  tmp3 = 0x00

  lγ = zeros(Float64, 3)

  g = 0

  for b in 1:m
    lγ[1] = priorC.logγ[1]
    lγ[2] = priorC.logγ[2]
    lγ[3] = priorC.logγ[3]

    tmp1[1] = SpecialFunctions.logbeta(priorC.A[1, b], priorC.B[1, b])
    tmp1[2] = SpecialFunctions.logbeta(priorC.A[2, b], priorC.B[2, b])
    tmp1[3] = SpecialFunctions.logbeta(priorC.A[3, b], priorC.B[3, b]) -
              priorC.logω[k][1]
    tmp1[4] = SpecialFunctions.logbeta(priorC.A[4, b], priorC.B[4, b]) -
              priorC.logω[k][2]

    for l in 1:k
      g = cl[l]

      # noise
      lγ[1] += SpecialFunctions.logbeta(priorC.A[1, b] + n1s[g, b],
                                        priorC.B[1, b] + v[g] - n1s[g, b]) -
               tmp1[1]

      # weak signal
      lγ[2] += SpecialFunctions.logbeta(priorC.A[2, b] + n1s[g, b],
                                        priorC.B[2, b] + v[g] - n1s[g, b]) -
               tmp1[2]

      # strong signal but not characteristic
      tmp2[1, l] = SpecialFunctions.logbeta(priorC.A[3, b] + n1s[g, b],
                                            priorC.B[3, b] + v[g] - n1s[g, b]) -
                   tmp1[3]

      # strong signal and characteristic
      tmp2[2, l] = SpecialFunctions.logbeta(priorC.A[4, b] + n1s[g, b],
                                            priorC.B[4, b] + v[g] - n1s[g, b]) -
                   tmp1[4]

      if tmp2[1, l] >= tmp2[2, l]
        tmp3 = log1p(exp(tmp2[2, l] - tmp2[1, l]))
        lγ[3] += tmp2[1, l] + tmp3
        tmp2[1, l] = tmp3
      else
        tmp3 = log1p(exp(tmp2[1, l] - tmp2[2, l]))
        lγ[3] += tmp2[2, l] + tmp3
        tmp2[2, l] = tmp3
      end
    end

    if lγ[3] > lγ[2]
      if lγ[3] > lγ[1]
        logpC[1] += priorC.logγ[3]
        logpC[2] -= log1p(exp(lγ[1] - lγ[3]) + exp(lγ[2] - lγ[3]))
        for l in 1:k
          if tmp2[1, l] >= tmp2[2, l]
            C[cl[l], b] = 0x03
            logpC[1] += priorC.logω[k][1]
            logpC[2] -= tmp2[1, l]
          else
            C[cl[l], b] = 0x04
            logpC[1] += priorC.logω[k][2]
            logpC[2] -= tmp2[2, l]
          end
        end
      else
        logpC[1] += priorC.logγ[1]
        logpC[2] -= log1p(exp(lγ[2] - lγ[1]) + exp(lγ[3] - lγ[1]))
        for l in 1:k
          C[cl[l], b] = 0x01
        end
      end
    elseif lγ[2] > lγ[1]
      logpC[1] += priorC.logγ[2]
      logpC[2] -= log1p(exp(lγ[1] - lγ[2]) + exp(lγ[3] - lγ[2]))
      for l in 1:k
        C[cl[l], b] = 0x02
      end
    else
      logpC[1] += priorC.logγ[1]
      logpC[2] -= log1p(exp(lγ[2] - lγ[1]) + exp(lγ[3] - lγ[1]))
      for l in 1:k
        C[cl[l], b] = 0x01
      end
    end
  end

  nothing
end

function computelocalmode!(state::State,
                           priorC::PriorColPartition)
  fill!(state.logpC, 0.0)

  m = size(state.n1s, 2)
  k = state.k

  tmp1 = zeros(Float64, 4)
  tmp2 = zeros(Float64, 2, k)
  tmp3 = 0x00

  lγ = zeros(Float64, 3)

  g = 0

  for b in 1:m
    lγ[1] = priorC.logγ[1]
    lγ[2] = priorC.logγ[2]
    lγ[3] = priorC.logγ[3]

    tmp1[1] = SpecialFunctions.logbeta(priorC.A[1, b], priorC.B[1, b])
    tmp1[2] = SpecialFunctions.logbeta(priorC.A[2, b], priorC.B[2, b])
    tmp1[3] = SpecialFunctions.logbeta(priorC.A[3, b], priorC.B[3, b]) -
              priorC.logω[k][1]
    tmp1[4] = SpecialFunctions.logbeta(priorC.A[4, b], priorC.B[4, b]) -
              priorC.logω[k][2]

    for l in 1:k
      g = state.cl[l]

      # noise
      lγ[1] += SpecialFunctions.logbeta(priorC.A[1, b] + state.n1s[g, b],
                                        priorC.B[1, b] + state.v[g] -
                                        state.n1s[g, b]) - tmp1[1]

      # weak signal
      lγ[2] += SpecialFunctions.logbeta(priorC.A[2, b] + state.n1s[g, b],
                                        priorC.B[2, b] + state.v[g] -
                                        state.n1s[g, b]) - tmp1[2]

      # strong signal but not characteristic
      tmp2[1, l] = SpecialFunctions.logbeta(priorC.A[3, b] + state.n1s[g, b],
                                            priorC.B[3, b] + state.v[g] -
                                            state.n1s[g, b]) - tmp1[3]

      # strong signal and characteristic
      tmp2[2, l] = SpecialFunctions.logbeta(priorC.A[4, b] + state.n1s[g, b],
                                            priorC.B[4, b] + state.v[g] -
                                            state.n1s[g, b]) - tmp1[4]

      if tmp2[1, l] >= tmp2[2, l]
        tmp3 = log1p(exp(tmp2[2, l] - tmp2[1, l]))
        lγ[3] += tmp2[1, l] + tmp3
        tmp2[1, l] = tmp3
      else
        tmp3 = log1p(exp(tmp2[1, l] - tmp2[2, l]))
        lγ[3] += tmp2[2, l] + tmp3
        tmp2[2, l] = tmp3
      end
    end

    if lγ[3] > lγ[2]
      if lγ[3] > lγ[1]
        state.logpC[1] += priorC.logγ[3]
        state.logpC[2] -= log1p(exp(lγ[1] - lγ[3]) + exp(lγ[2] - lγ[3]))
        for l in 1:k
          if tmp2[1, l] >= tmp2[2, l]
            state.C[state.cl[l], b] = 0x03
            state.logpC[1] += priorC.logω[k][1]
            state.logpC[2] -= tmp2[1, l]
          else
            state.C[state.cl[l], b] = 0x04
            state.logpC[1] += priorC.logω[k][2]
            state.logpC[2] -= tmp2[2, l]
          end
        end
      else
        state.logpC[1] += priorC.logγ[1]
        state.logpC[2] -= log1p(exp(lγ[2] - lγ[1]) + exp(lγ[3] - lγ[1]))
        for l in 1:k
          state.C[state.cl[l], b] = 0x01
        end
      end
    elseif lγ[2] > lγ[1]
      state.logpC[1] += priorC.logγ[2]
      state.logpC[2] -= log1p(exp(lγ[1] - lγ[2]) + exp(lγ[3] - lγ[2]))
      for l in 1:k
        state.C[state.cl[l], b] = 0x02
      end
    else
      state.logpC[1] += priorC.logγ[1]
      state.logpC[2] -= log1p(exp(lγ[2] - lγ[1]) + exp(lγ[3] - lγ[1]))
      for l in 1:k
        state.C[state.cl[l], b] = 0x01
      end
    end
  end

  nothing
end


text:
export @set

struct PolynomialInequality{PT<:AbstractPolynomialLike}
    p::PT
end
function Base.intersect(ineq::PolynomialInequality; kws...)
    return BasicSemialgebraicSet(FullSpace(), [ineq.p])
end

struct PolynomialEquality{PT<:AbstractPolynomialLike}
    p::PT
end
function equality(lhs, rhs)
    return PolynomialEquality(lhs - rhs)
end
function Base.intersect(eq::PolynomialEquality; lib_or_solver=nothing)
    if lib_or_solver === nothing
        return algebraicset([eq.p])
    else
        return algebraicset([eq.p], lib_or_solver)
    end
end

const Element = Union{PolynomialInequality, PolynomialEquality, FixedVariable}

function Base.intersect(el::Element,
                        args...; kws...)
    return intersect(intersect(el; kws...), args...; kws...)
end

function Base.intersect(set::AbstractBasicSemialgebraicSet,
                        el::Element, args...; kws...)
    return intersect(set, intersect(el; kws...), args...; kws...)
end

# Taken from JuMP/macros.jl
function _canonicalize_sense(sns::Symbol, _error)
    if sns == :(==)
        return (:(==),false)
    elseif sns == :(>=) || sns == :(≥)
        return (:(>=),false)
    elseif sns == :(<=) || sns == :(≤)
        return (:(<=),false)
    elseif sns == :(.==)
        return (:(==),true)
    elseif sns == :(.>=) || sns == :(.≥)
        return (:(>=),true)
    elseif sns == :(.<=) || sns == :(.≤)
        return (:(<=),true)
    else
        _error("Unrecognized sense $sns")
    end
end

function appendconstraints!(elements, expr, _error)
    if Base.Meta.isexpr(expr, :call)
        try
            sense, vectorized = _canonicalize_sense(expr.args[1], _error)
            @assert !vectorized
            if sense == :(>=)
                push!(elements, esc(:(SemialgebraicSets.PolynomialInequality($(expr.args[2]) - $(expr.args[3])))))
            elseif sense == :(<=)
                push!(elements, esc(:(SemialgebraicSets.PolynomialInequality($(expr.args[3]) - $(expr.args[2])))))
            elseif sense == :(==)
                push!(elements, esc(:(SemialgebraicSets.equality($(expr.args[2]), $(expr.args[3])))))
            else
                _error("Unrecognized sense $(string(sense)) in domain specification")
            end
        catch
            push!(elements, esc(expr))
        end
    elseif Base.Meta.isexpr(expr, :&&)
        map(t -> appendconstraints!(elements, t, _error), expr.args)
    else
        push!(elements, esc(expr))
    end
    return nothing
end

macro set(expr, library=nothing)
    elements = []
    appendconstraints!(elements, expr, msg -> error("In @set($expr: ", msg))
    if library === nothing
        return :( intersect($(elements...)) )
    else
        return :( intersect($(elements...); lib_or_solver=$(esc(library))) )
    end
end


text:
module UVH5Splicer

export uvh5_splice

using Pkg
using HDF5

#=const=# EQUALITY_EXEMPT_HEADERS = (
    "Nblts",
    "Nfreqs",
    "Ntimes",
    "ant_1_array",
    "ant_2_array",
    "channel_width",
    "flex_spw_id_array",
    "extra_keywords",
    "freq_array",
    "history",
    "integration_time",
    "lst_array",
    "time_array",
    "uvw_array"
)

function all_same(hdrs, name)
    values = getindex.(getindex.(hdrs, name))
    if all(Ref(values[1]) .== values[2:end])
        return true
    end
    @warn "mismatch detected in $name" _module=nothing _file=nothing
    return false
end

"""
    check_commensurate(hdrs, check_names)

Check that the datasets listed in `check_names` in the `hdrs` are equal.
Returns an Array of `check_names` datasets that are not equal across all `hdrs`,
which will be empty if all the datasets are equal across all headers.
"""
function check_commensurate(hdrs, check_names)
    oks = all_same.(Ref(hdrs), check_names)
    check_names[.!oks]
end

function splice_dataset_by_frequency(h5out, h5ins, name, outdims)
    # Create data type for virtual dataset
    vtype = datatype(h5ins[1][name])
    # Create dataspace for virtual dataset
    vspace = dataspace(outdims) # (Npols, Nfreqsout, Nbltsout)
    # Create dataset creation property list for virtual dataset
    dcpl = create_property(HDF5.H5P_DATASET_CREATE)
    # Create source dataspace based on Nfreqs of first dataset
    nfreqin = h5ins[1]["Header/Nfreqs"][]
    sspace = dataspace(h5ins[1][name])
    sspace = HDF5.hyperslab(sspace, 1:outdims[1], 1:nfreqin, 1:outdims[3])

    # Setup hyperslab parameters
    freq_offset = 0

    # Map regions of vspace to regions in source files/datasets
    for h5 in h5ins
        # Ensure source dataspace has proper Nfreqs dimension
        if nfreqin != h5["Header/Nfreqs"][]
            nfreqin = h5["Header/Nfreqs"][]
            close(sspace)
            sspace = dataspace(h5ins[1][name])
            sspace = HDF5.hyperslab(sspace, 1:outdims[1], 1:nfreqin, 1:outdims[3])
        end

        # Select hyperslab in virtual dataspace
        mapspace = HDF5.hyperslab(copy(vspace), 1:outdims[1], (1:nfreqin).+freq_offset, 1:outdims[3])

        # Map it!
        HDF5.h5p_set_virtual(dcpl, mapspace, HDF5.filename(h5), name, sspace)
        close(mapspace)

        # Adjust freq_offset
        freq_offset += nfreqin
    end
    # Close source dataspace
    close(sspace)
    # Since h5s_select_all is not currently availble in HDF5.jl, close and
    # recreate vspace.
    close(vspace)
    vspace = dataspace(outdims) # (Npols, Nfreqsout, Nbltsout)

    # Create virtual dataset
    # TODO Ensure links are created with UTF8 cset encoding?
    HDF5.h5d_create(h5out, name, vtype, vspace,
                    HDF5.DEFAULT_PROPERTIES, dcpl, HDF5.DEFAULT_PROPERTIES)

    # Close up
    close(vtype)
    close(vspace)
end

function uvh5_splice(outfile, infiles)
    # Open input files
    h5ins = h5open.(infiles)

    # Get Header groups
    hdrins = getindex.(h5ins, "Header")
    # Get all names from all input headers
    all_names = keys.(hdrins)
    # Get names common to all input headers
    common_names = reduce(intersect, all_names)
    # Reject names in the EQUALITY_EXEMPT_HEADERS list
    copy_names = filter(∉(EQUALITY_EXEMPT_HEADERS), common_names)

    # Check that the must-be-equal datasets are equal
    mismatches = check_commensurate(hdrins, copy_names)
    if !isempty(mismatches)
        error("input datasets $mismatches are not commensurate")
    end

    # Determine the minimum Nblts

    Nbltsout = minimum(getindex.(getindex.(hdrins, "Nblts")))
    # Check 1D arrays dimensioned by Nblts
    for name in ("ant_1_array", "ant_2_array", "integration_time",
                 "lst_array", "time_array")
        vals = getindex.(getindex.(hdrins, name))
        for i = 1:length(vals)
            if vals[1][1:Nbltsout] != vals[i][1:Nbltsout]
                error("input datasets $name are not commensurate")
            end
        end
    end
    # Check 1D arrays dimensioned by Nblts
    vals = getindex.(getindex.(hdrins, "uvw_array"))
    for i = 1:length(vals)
        if vals[1][:, 1:Nbltsout] != vals[i][:, 1:Nbltsout]
            error("input datasets uvw_array are not commensurate")
        end
    end

    # Everything is commensurate!  OK to create output file.
    h5out = h5open(outfile, "w")
    hdrout = create_group(h5out, "Header")

    # Copy the "easy" datasets
    for name in copy_names
        @info "copying $name"
        HDF5.h5o_copy(hdrins[1], name, hdrout, name, HDF5.H5P_DEFAULT, HDF5.H5P_DEFAULT)
    end

    # Copy the Nblts datasets, then adjust to Nbltsout
    # Scalar
    @info "storing Nblts"
    hdrout["Nblts", layout=HDF5.H5D_COMPACT] = Nbltsout
    # 1D
    for name in ("ant_1_array", "ant_2_array", "integration_time",
                 "lst_array", "time_array")
        @info "copying $name and setting dims"
        HDF5.h5o_copy(hdrins[1], name, hdrout, name, HDF5.H5P_DEFAULT, HDF5.H5P_DEFAULT)
        HDF5.set_extent_dims(hdrout[name], (Nbltsout,))
    end
    # 2D
    @info "copying uvw_array and setting dims"
    HDF5.h5o_copy(hdrins[1], "uvw_array", hdrout, "uvw_array", HDF5.H5P_DEFAULT, HDF5.H5P_DEFAULT)
    HDF5.set_extent_dims(hdrout["uvw_array"], (3, Nbltsout))
    # Recompute Ntimes
    time_array_out = hdrout["time_array"][]
    Ntimesout = length(unique(time_array_out))
    @info "storing Ntimes"
    hdrout["Ntimes", layout=HDF5.H5D_COMPACT] = Ntimesout

    # Concatenate frequency datasets
    freq_array_out = reduce(vcat, getindex.(getindex.(hdrins, "freq_array")))
    chan_width_out = reduce(vcat, getindex.(getindex.(hdrins, "channel_width")))
    Nfreqsout = length(freq_array_out)
    @info "storing concatenated freq_array"
    hdrout["freq_array"] = freq_array_out
    @info "storing concatenated channel_width"
    hdrout["channel_width"] = chan_width_out
    @info "storing Nfreqs"
    hdrout["Nfreqs", layout=HDF5.H5D_COMPACT] = Nfreqsout

    # Output history dataset
    toml = Pkg.TOML.parsefile(joinpath(pkgdir(@__MODULE__), "Project.toml"))
    history = reduce((a,b)->"$a\n- $b", HDF5.filename.(h5ins),
                     init="# Created with $(toml["name"]) version $(toml["version"])")

    # This is probably more complicated than it needs to be
    @info "storing history"
    dt = HDF5.create_datatype(HDF5.H5T_STRING, sizeof(history))
    HDF5.h5t_set_strpad(dt, HDF5.H5T_STR_NULLTERM)
    HDF5.h5t_set_cset(dt, HDF5.H5T_CSET_UTF8)
    history_dataset = create_dataset(hdrout, "history", dt, ())
    write_dataset(history_dataset, dt, history)
    close(dt)

    # Create virtual datasets in "Data" group
    dataout = create_group(h5out, "Data")

    # Create dataspace for spliced datasets
    Npols = hdrins[1]["Npols"][]
    for name in ("visdata", "flags", "nsamples")
        @info "splicing $name"
        splice_dataset_by_frequency(h5out, h5ins, "/Data/$name", (Npols, Nfreqsout, Nbltsout))
    end

    # Close input files and output file
    close.(h5ins)
    close(h5out)
end

end # module UVH5Splicer


text:
module HwlocViz

using Hwloc, Plots, AbstractTrees, GraphRecipes

const topo = topology_load()

AbstractTrees.children(x::Hwloc.Object) = x.children

function AbstractTrees.printnode(io::IO, obj::Hwloc.Object)
    idxstr = obj.type_ in (:Package, :Core, :PU) ? "L#$(obj.logical_index) P#$(obj.os_index) " : ""
    attrstr = string(obj.attr)

    if obj.type_ in (:L1Cache, :L2Cache, :L3Cache)
        tstr = first(string(obj.type_), 2)
        csize = obj.attr.size ÷ 1024
        if csize > 1000
            attrstr = "($(csize ÷ 1024) MB)"
        else
            attrstr = "($csize KB)"
        end
    elseif obj.type_ == :Package
        tstr = "Package"
        # TODO: add total memory info
    else
        tstr = string(obj.type_)
    end

    print(io, tstr, " ",
        idxstr,
        attrstr)
end

function print_topology()
    # print topology
    print_tree(topo, maxdepth=12)
    # print summary
    h = histmap(topo)
    caches = cachesize() .÷ 1024
    l1 = caches[1] > 1024 ? "$(caches[1] ÷ 1024) MB" : "$(caches[1]) KB"
    l2 = caches[2] > 1024 ? "$(caches[2] ÷ 1024) MB" : "$(caches[2]) KB"
    l3 = caches[3] > 1024 ? "$(caches[3] ÷ 1024) MB" : "$(caches[3]) KB"
    println("\nPackages: $(h[:Package]) \t L3: $l3")
    println("Cores: $(h[:Core]) \t L2: $l2")
    println("PUs: $(h[:PU]) \t L1: $l1")
    return nothing
end

function cachesize()
    l1 = l2 = l3 = zero(Int64)
    fl1 = fl2 = fl3 = false
    for obj in topo
        if obj.type_ == :L1Cache
            l1 = obj.attr.size::Int64
            fl1 = true
        elseif obj.type_ == :L2Cache
            l2 = obj.attr.size::Int64
            fl2 = true
        elseif obj.type_ == :L3Cache
            l3 = obj.attr.size::Int64
            fl3 = true
        end
        (fl1 && fl2 && fl3) && break
    end
    return (l1, l2, l3)
end

plot_topology(; kwargs...) = begin
    return plot(TreePlot(topo), fontsize=10, nodeshape=:circle, nodecolor=:white, size=(800,800), kwargs...)
end

plot_topology(filename::String; kwargs...) = begin
    p = plot_topology(; kwargs...)
    savefig(p, filename)
    return nothing
end

export plot_topology, print_topology

end

max_stars_repo_name: crstnbr/HwlocViz

max_issues_repo_name: crstnbr/HwlocViz

max_forks_repo_name: crstnbr/HwlocViz


text:
module PreProcess

include("../Types.jl")

export run

INTERPOLATIONTIME = 10 # in millis

scaletime = (time, starttime) -> (time - starttime)

function linearInterpolate(dp1:: Types.Types.Datapoint, dp2:: Types.Types.Datapoint, interpolatedtime:: Int)::Types.Types.Datapoint
    dt = dp2.time - dp1.time
    dv = dp2.magnitude - dp1.magnitude
    mag = (dv/dt) * (interpolatedtime - dp1.time) + dp1.magnitude
    Types.Types.Datapoint(interpolatedtime, mag)
end

function run(input:: Channel)
    window = []
    starttime = -1
    interpolationcount = 0
    return function (output:: Channel)
        for dp in input
            if starttime == -1
                starttime = copy(dp.time)
            end
            local_dp = Types.Types.Datapoint(scaletime(dp.time, starttime), dp.magnitude)
            push!(window, local_dp)
            if size(window)[1] >= 2
                dp1 = deepcopy(window[1])
                dp2 = deepcopy(window[2])
                numberofpoints = ceil(Int, (dp2.time - dp1.time) / INTERPOLATIONTIME)
                for i in range(1, numberofpoints; step=1)
                    interpolatedtime = interpolationcount * INTERPOLATIONTIME
                    if (dp1.time <= interpolatedtime < dp2.time)
                        interpolated = linearInterpolate(dp1, dp2, interpolatedtime)
                        put!(output, interpolated)
                        interpolationcount += 1
                    end
                    if !isempty(window)
                        popfirst!(window)
                    else
                        println("EMPTY!")
                    end
                end
            end
        end
    end
end

end

max_stars_repo_name: eliabiraschi/pedometer

max_issues_repo_name: eliabiraschi/pedometer

max_forks_repo_name: eliabiraschi/pedometer


text:
module MOModelingExtensions

using JuMP
import Base: max, min, maximum, minimum, &, |, any, all, ⊻, !, ~, ⇒

function _check_var_has_bounds(var::AbstractVariableRef)
    if ! has_lower_bound(var)
        error("Variable $(var) has no lower bound")
    end
    if ! has_upper_bound(var)
        error("Variable $(var) has no upper bound")
    end
end

function _check_var_is_binary(var::AbstractVariableRef)
    if ! is_binary(var)
        error("Variable $(var) is not binary")
    end
end

# Minimum and maximum.

function max(a::AbstractVariableRef, b::AbstractVariableRef)::AbstractVariableRef
    # Return a new variable, x. This function uses an internal binary variable, z.
    # Model: x <= a + z * M, x <= b + (1 - z) * M,   x >= a, x >= b
    #        \_________ x <= max{a, b} _________/ \_ x >= max{a, b} _/
    # Interpretation: z = 0 if a is the maximum, 1 if it is b.
    model = owner_model(a)
    check_belongs_to_model(b, model)

    _check_var_has_bounds(a)
    _check_var_has_bounds(b)

    x_lower_bound = max(lower_bound(a), lower_bound(b))
    x_upper_bound = max(upper_bound(a), upper_bound(b))
    # Maximum difference between a and b.
    M = max(abs(lower_bound(a) - upper_bound(b)), abs(lower_bound(b) - upper_bound(a)))

    x = @variable(model, lower_bound=x_lower_bound, upper_bound=x_upper_bound)
    z = @variable(model, binary=true)

    @constraint(model, x <= a + M * z)
    @constraint(model, x <= b + M * (1 - z))
    @constraint(model, x >= a)
    @constraint(model, x >= b)

    return x
end

function maximum(vars::AbstractVector{<:AbstractVariableRef})::AbstractVariableRef
    # Same model as max, but with one binary variable per input variable.
    model = owner_model(vars[1])
    for var in vars # TODO: Only check for first variable.
        check_belongs_to_model(var, model)
    end

    for var in vars
        _check_var_has_bounds(var)
    end

    x_lower_bound = maximum(lower_bound.(vars))
    x_upper_bound = maximum(upper_bound.(vars))
    # Maximum difference between two variables.
    M = maximum(abs.(lower_bound(v1) - upper_bound(v2) for v1 in vars, v2 in vars))

    x = @variable(model, lower_bound=x_lower_bound, upper_bound=x_upper_bound)
    z = @variable(model, [eachindex(vars)], binary=true)

    @constraint(model, [var in eachindex(vars)], x <= vars[var] + M * (1 - z[var]))
    @constraint(model, [var in eachindex(vars)], x >= vars[var])
    @constraint(model, sum(z) == 1)

    return x
end

function min(a::AbstractVariableRef, b::AbstractVariableRef)::AbstractVariableRef
    # Return a new variable, x. This function uses an internal binary variable, z.
    # Model: x >= a + z * M, x >= b + (1 - z) * M,   x <= a, x <= b
    #        \_________ x >= min{a, b} _________/ \_ x <= min{a, b} _/
    # Interpretation: z = 0 if a is the minimum, 1 if it is b.
    model = owner_model(a)
    check_belongs_to_model(b, model)

    _check_var_has_bounds(a)
    _check_var_has_bounds(b)

    x_lower_bound = min(lower_bound(a), lower_bound(b))
    x_upper_bound = min(upper_bound(a), upper_bound(b))
    # Maximum difference between a and b.
    M = max(abs(lower_bound(a) - upper_bound(b)), abs(lower_bound(b) - upper_bound(a)))

    x = @variable(model, lower_bound=x_lower_bound, upper_bound=x_upper_bound)
    z = @variable(model, binary=true)

    @constraint(model, x >= a - M * z)
    @constraint(model, x >= b - M * (1 - z))
    @constraint(model, x <= a)
    @constraint(model, x <= b)

    return x
end

function minimum(vars::AbstractVector{<:AbstractVariableRef})::AbstractVariableRef
    # Same model as min, but with one binary variable per input variable.
    model = owner_model(vars[1])
    for var in vars # TODO: Only check for first variable.
        check_belongs_to_model(var, model)
    end

    for var in vars
        _check_var_has_bounds(var)
    end

    x_lower_bound = maximum(lower_bound.(vars))
    x_upper_bound = maximum(upper_bound.(vars))
    # Maximum difference between two variables.
    M = maximum(abs.(lower_bound(v1) - upper_bound(v2) for v1 in vars, v2 in vars))

    x = @variable(model, lower_bound=x_lower_bound, upper_bound=x_upper_bound)
    z = @variable(model, [eachindex(vars)], binary=true)

    @constraint(model, [var in eachindex(vars)], x >= vars[var] - M * (1 - z[var]))
    @constraint(model, [var in eachindex(vars)], x <= vars[var])
    @constraint(model, sum(z) == 1)

    return x
end

# Boolean operators.

function Base.:!(a::AbstractVariableRef)::AffExpr
    _check_var_is_binary(a)
    return 1 - a
end

function Base.:~(a::AbstractVariableRef)::AffExpr
    _check_var_is_binary(a)
    return 1 - a
end

function Base.:&(a::AbstractVariableRef, b::AbstractVariableRef)::AbstractVariableRef
    model = owner_model(a)
    check_belongs_to_model(b, model)

    _check_var_is_binary(a)
    _check_var_is_binary(b)

    z = @variable(model, binary=true)
    @constraint(model, z <= a)
    @constraint(model, z <= b)
    @constraint(model, z >= a + b - 1)

    return z
end

function Base.all(vars::AbstractVector{<:AbstractVariableRef})::AbstractVariableRef
    model = owner_model(vars[1])
    for var in vars # TODO: Only check for first variable.
        check_belongs_to_model(var, model)
    end

    for var in vars
        _check_var_is_binary(var)
    end

    z = @variable(model, binary=true)
    @constraint(model, [var in eachindex(vars)], z <= vars[var])
    @constraint(model, z >= sum(vars) - length(vars) + 1)

    return z
end

function Base.:|(a::AbstractVariableRef, b::AbstractVariableRef)::AbstractVariableRef
    model = owner_model(a)
    check_belongs_to_model(b, model)

    _check_var_is_binary(a)
    _check_var_is_binary(b)

    z = @variable(model, binary=true)
    @constraint(model, z >= a)
    @constraint(model, z >= b)
    @constraint(model, z <= a + b)

    return z
end

function Base.any(vars::AbstractVector{<:AbstractVariableRef})::AbstractVariableRef
    model = owner_model(vars[1])
    for var in vars # TODO: Only check for first variable.
        check_belongs_to_model(var, model)
    end

    for var in vars
        _check_var_is_binary(var)
    end

    z = @variable(model, binary=true)
    @constraint(model, [var in eachindex(vars)], z >= vars[var])
    @constraint(model, z <= sum(vars))

    return z
end

function Base.:⊻(a::AbstractVariableRef, b::AbstractVariableRef)::AbstractVariableRef
    model = owner_model(a)
    check_belongs_to_model(b, model)

    _check_var_is_binary(a)
    _check_var_is_binary(b)

    z = @variable(model, binary=true)
    @constraint(model, z <= a + b)
    @constraint(model, z >= a - b)
    @constraint(model, z >= b - a)
    @constraint(model, z <= 2 - a - b)

    return z
end

# Don't override Pair or =>, as this might have unwanted side effects (could
# no more create a dictionary (variable, variable), for instance).
# function Base.:⇒(a::AbstractVariableRef, b::AbstractVariableRef)::AbstractVariableRef
#     model = owner_model(a)
#     check_belongs_to_model(b, model)
#
#     _check_var_is_binary(a)
#     _check_var_is_binary(b)
#
#     @constraint(model, a <= b)
#
#     # TODO: What should this function return?
# end

end


text:
using CompScienceMeshes, BEAST
using LinearAlgebra
using Profile
using TimerOutputs
using StaticArrays

function tau(x::SVector{U,T}) where {U,T}
    1.0-1.0/5.0
end

ntrc = X->ntrace(X,y)

T = tetmeshsphere(1.0,0.2)
X = nedelecd3d(T)
y = boundary(T)
@show numfunctions(X)

ϵ, μ, ω = 1.0, 1.0, 1.0; κ, η = ω * √(ϵ*μ), √(μ/ϵ)
ϵ_r =5.0
χ = tau
K, I, B = VIE.singlelayer(wavenumber=κ, tau=χ), Identity(), VIE.boundary(wavenumber=κ, tau=χ)
E = VIE.planewave(direction=ẑ, polarization=x̂, wavenumber=κ)
H = -1/(im*κ*η)*curl(E)

@hilbertspace j
@hilbertspace k
α = 1.0/ϵ_r
eq = @varform α*I[j,k]-K[j,k]-B[ntrc(j),k] == E[j]

dbvie = @discretise eq j∈X k∈X
u_m = solve(dbvie)
#=
#postprocessing
Φ, Θ = [0.0], range(0,stop=2π,length=100)
pts = [point(cos(ϕ)*sin(θ), sin(ϕ)*sin(θ), cos(θ)) for θ in Θ for ϕ in Φ]
ffd = potential(VIE.farfield(wavenumber=κ, tau=χ), pts, u, X)
ff = im*κ*ffd

using Plots

#Farfield
plot(xlabel="theta")
plot!(Θ, norm.(ff), label="far field", title="D-VIE")
=#
using Plots
#NearField
Z = range(-1,1,length=100)
Y = range(-1,1,length=100)
nfpoints = [point(0,y,z) for z in Z, y in Y]

Enear = BEAST.grideval(nfpoints,α.*u_m,X)
Enear = reshape(Enear,100,100)

contour(real.(getindex.(Enear,1)))
heatmap(Z, Y,  real.(getindex.(Enear,1)))

plot!(Y[2:99],real.(getindex.(Enear[2:99,50],1)),label="D-VIE (simplex)", linestyle=:dash, linecolor=:darkorange4)
#plot!(Y[2:99],real.(getindex.(Enear_simplex[2:99,50],1)),label="D-VIE (simplex)", linestyle=:solid, linecolor=:darkorange3)

#=
import Cairo
using DataFrames, Gadfly, RDatasets
D = dataset("datasets","HairEyeColor")
palette = ["skyblue","skyblue3"]

D = DataFrame(
   Type=["D-VIE \n (1,190 Tets.)","D-VIE \n (1,190 Tets.)","D-VIE \n (1,190 Tets.)","D-VIE \n (1,190 Tets.)","E-VIE \n (1,190 Tets.)","E-VIE \n (1,190 Tets.)","E-VIE \n (1,190 Tets.)","E-VIE \n (1,190 Tets.)","D-VIE \n (2,706 Tets.)","D-VIE \n (2,706 Tets.)","D-VIE \n (2,706 Tets.)","D-VIE \n (2,706 Tets.)","E-VIE \n (2,706 Tets.)","E-VIE \n (2,706 Tets.)","E-VIE \n (2,706 Tets.)","E-VIE \n (2,706 Tets.)"],
   Quad=["Simplex","Simplex","Classic","Classic","Simplex","Simplex","Classic","Classic","Simplex","Simplex","Classic","Classic","Simplex","Simplex","Classic","Classic"],
   Sing=["Regular","Singular","Regular","Singular","Regular","Singular","Regular","Singular","Regular","Singular","Regular","Singular","Regular","Singular","Regular","Singular"],
   Time=[24,4.77,23.1,27.9,41.9,8.49,41.5,45.7,125,13.1,125,68.6,222,20.5,222,115],
)

p = Gadfly.plot(D, x=:Quad, y=:Time, color=:Sing, xgroup=:Type,
    Geom.subplot_grid(Geom.bar(position=:stack)),
    Scale.color_discrete_manual(palette...),
    Guide.xlabel(""),Guide.ylabel("Assembly time [s]"), Guide.colorkey("Interaction"),Theme(major_label_font="Times",minor_label_font="Times",key_title_font="Times", key_label_font="Times",background_color=color("white")))

# draw(PNG("haireyecolor.png", 6.6inch, 4inch), p)
draw(PNG(6.6inch, 4inch), p)
=#


text:
# Get an unsigned position
function uposition(stream)
    UInt(position(stream))
end

# Function that generates a name based on the name of the cached function
function generate_cache_filename()
    _filename = "_" * string(hash(rand()), base=16) * "_.bin"
    return abspath(_filename)
end

# Function that retrieves one entry from a stream
function load_disk_cache_entry(io::IO, pos_start::UInt, pos_end::UInt,
                               decompressor=Noop)
    seek(io, pos_start)
    cbuf = read(io, pos_end-pos_start)
    buf = transcode(decompressor, cbuf)
    datum = deserialize(IOBuffer(buf))
    return datum
end

# Function that stores one entry to a stream
function store_disk_cache_entry(io::IO, pos_start::UInt, datum,
                                compressor=Noop)
    seek(io, pos_start)
    buf = IOBuffer()
    serialize(buf, datum)
    posbuf = uposition(buf)
    cbuf = transcode(compressor, buf.data[1:posbuf])
    write(io, cbuf)
    return uposition(io)
end

# Function that checks the consistency of the cache and `offset`
# field value of the Cache object
function check_disk_cache(filename::String,
                          offsets::D where D<:Dict,
                          output_type::Type)::Bool
    open(filename, "r") do fid
        _, decompressor = get_transcoders(filename)
        try
            O = deserialize(fid);
            @assert O == output_type
            for (_, (startpos, endpos)) in offsets
                load_disk_cache_entry(fid, startpos, endpos, decompressor)
            end
            return true
        catch excep
            return false
        end
    end
end

function get_transcoders(filename::T="") where T<:AbstractString
    ext = split(filename, ".")[end]
    if ext == "bz2" || ext == "bzip2"
        compressor = Bzip2Compressor
        decompressor = Bzip2Decompressor
    elseif ext == "gz" || ext == "gzip"
        compressor = GzipCompressor
        decompressor = GzipDecompressor
    else
        compressor = Noop  # no compression
        decompressor = Noop  # no compression
    end
    return compressor, decompressor
end


text:
using Cairo
using LSystems
using DataStructures

c = CairoRGBSurface(704, 1024);
cr = CairoContext(c);

save(cr);

set_source_rgb(cr, 1.0, 1.0, 1.0);
rectangle(cr, 0.0, 0.0, 704.0, 1024.0);
fill(cr);

restore(cr);
save(cr);

x, y, a, d = 352., 992., 0, 4
set_line_width(cr, 0.6);
move_to(cr, x, y);

kristall_start = "F"
kristall_trans = Dict([('F', "F+F--f+F-F++f-F"), ('f', "fff")])

function determine_new_position(x, y, a, d)
    x -= sind(a) * d
    y -= cosd(a) * d
    return x, y
end

s = Stack(Vector{Float64}) # [x, y, a]

level = 5
drawing_instructions = @task lindenmayer(kristall_start, kristall_trans, level)

for v in drawing_instructions
    if v == 'F' || v == 'f'
        x, y = determine_new_position(x, y, a, d)
        line_to(cr, x, y)
    elseif v == '+'
        a -= 90
    elseif v == '-'
        a += 90
    end
end

stroke(cr);
write_to_png(c, "kristall_$level.png");


text:
"""
A generally useful event scheduler module.

The `ExtensibleScheduler` package provides advanced and extensible [Julia](https://julialang.org/) events schedulers inspired by Python schedulers [APScheduler](https://apscheduler.readthedocs.io/), [schedule](https://github.com/dbader/schedule) and [sched](https://docs.python.org/3/library/sched.html).

It's also inspired by [Sched.jl](https://github.com/scls19fr/Sched.jl), a Julia event scheduler inspired by Python [sched](https://docs.python.org/3/library/sched.html).

Schedulers can use real time clock (system time) or simulated time (for simulation purpose).
"""
module ExtensibleScheduler

    export BlockingScheduler, run_pending
    export Action, add
    export Trigger
    export SimClock, set, NoSleep
    export iterate

    using Dates

    include("priority.jl")
    include("clock/base.jl")
    include("action.jl")
    include("triggers/base.jl")
    include("job.jl")
    include("events.jl")
    include("executors/base.jl")
    include("jobstores/base.jl")
    include("schedulers/base.jl")

end # module


text:
"""
    user-defined matrix generators
    to be populated with `include_generator` in user source code
"""
const USERMATRIXDICT = Dict{String,Function}()

"""
Associate names with matrix-generating functions
"""
const MATRIXDICT = Dict("hilb" => hilb, "hadamard" => hadamard,
                  "cauchy" => cauchy, "circul" => circul,
                  "dingdong" => dingdong, "frank" => frank,
                  "invhilb" => invhilb, "forsythe" => forsythe,
                  "magic" => magic, "grcar" => grcar,
                  "triw" => triw, "moler" => moler,
                  "pascal" => pascal, "kahan" => kahan,
                  "pei" => pei, "vand" => vand,
                  "invol" => invol, "chebspec" => chebspec,
                  "lotkin" => lotkin, "clement" => clement,
                  "fiedler" => fiedler, "minij" => minij,
                  "binomial" => binomialm, "tridiag" => tridiag,
                  "lehmer" => lehmer, "parter" => parter,
                  "chow" => chow, "randcorr" => randcorr,
                  "poisson" => poisson, "neumann" => neumann,
                  "rosser" => rosser, "sampling" => sampling,
                  "wilkinson" => wilkinson, "rando" => rando,
                  "randsvd" => randsvd, "rohess" => rohess,
                  "kms" => kms, "wathen" => wathen,
                  "oscillate" => oscillate, "toeplitz" => toeplitz,
                  "hankel" => hankel, "golub" => golub,
                  "companion" => companion,
                  "prolate" => prolate, "deriv2" => deriv2,
                  "shaw" => shaw, "wing" => wing,
                  "foxgood" => foxgood, "heat" => heat,
                  "baart" => baart, "phillips" => phillips,
                  "gravity" => gravity, "blur" => blur,
                  "spikes" => spikes, "ursell" => ursell,
                  "parallax" => parallax, "erdrey" => erdrey,
                  "gilbert" => gilbert, "smallworld" => smallworld
 )

"""
    predefined matrix classes (for the generated functions)
"""
const MATRIXCLASS = Dict(
            :symmetric => ["hilb", "cauchy", "circul", "dingdong",
                             "invhilb", "moler", "pascal", "pei",
                             "clement", "fiedler", "minij", "tridiag",
                             "lehmer", "randcorr", "poisson", "wilkinson",
                             "kms", "wathen", "oscillate", "prolate",
                             "hankel"],

            :inverse => ["hilb", "hadamard", "cauchy", "invhilb",
                           "forsythe", "magic", "triw", "moler", "pascal",
                           "kahan", "pei", "vand", "invol", "lotkin",
                           "clement", "fiedler", "minij", "tridiag",
                           "lehmer", "poisson", "kms" ],

            :illcond => ["hilb", "cauchy", "frank", "invhilb",
                            "forsythe", "triw", "moler", "pascal",
                            "kahan","pei", "vand", "invol", "lotkin",
                            "tridiag", "rosser", "randsvd", "kms",
                            "oscillate", "prolate", "golub"],

            :posdef => ["hilb", "cauchy", "circul", "invhilb",
                           "moler", "pascal", "pei", "minij", "tridiag",
                           "lehmer", "poisson", "kms", "wathen", "oscillate"],

            :eigen =>   ["hadamard", "circul", "dingdong", "frank",
                           "forsythe", "grcar", "pascal", "invol","chebspec",
                           "lotkin", "clement", "fiedler", "minij",
                           "tridiag", "parter", "chow", "poisson", "neumann",
                           "rosser", "sampling", "wilkinson","wathen",
                           "oscillate"],

             :sparse => ["poisson", "neumann", "wathen", "blur", "erdrey", "gilbert",
                          "smallworld"],

             :random => ["rosser", "rando", "randcorr", "randsvd", "rohess",
                          "wathen", "oscillate", "golub", "erdrey", "gilbert", "smallworld"],

             :regprob => ["deriv2", "shaw", "wing", "foxgood", "heat",
                           "baart", "phillips", "gravity", "blur",
                           "spikes", "ursell", "parallax"],

              :graph => ["erdrey", "gilbert", "smallworld"]
)

# remote parameters for several data sources
const TA_REMOTE = TURemoteType(RemoteParameters(
                    "https://sparse.tamu.edu/MM",
                    "https://sparse.tamu.edu/?per_page=All",
                    """<title>SuiteSparse Matrix Collection</title>""",
                    ("", """">Matrix Market""", 4, ".tar.gz", 2,
                     r"([^<]*)	"),
                    ".tar.gz"
                   ))

const UF_REMOTE = TURemoteType(RemoteParameters(
                    "https://www.cise.ufl.edu/research/sparse/MM",
                    "https://www.cise.ufl.edu/research/sparse/matrices/list_by_id.html",
                    """<title>UF Sparse Matrix Collection - sorted by id</title>""",
                    ("", """>MM""", 4, ".tar.gz", 2, r"([[:digit:]]*)	"),
                    ".tar.gz"
                   ))

const MM_REMOTE = MMRemoteType(RemoteParameters(
                    "ftp://math.nist.gov/pub/MatrixMarket2",
                    "http://math.nist.gov/MatrixMarket/matrices.html",
                    """<TITLE>The Matrix Market Matrices by Name</TITLE>""",
                    ("M", """<A HREF="/MatrixMarket/data/""", 2, ".html", 3, nothing),
                    ".mtx.gz"
                   ))

# preferred remote source for UF matrix collection (TAMU vs. UFl)
uf_remote = TA_REMOTE # may be altered
preferred(::Type{TURemoteType}) = uf_remote
preferred(::Type{MMRemoteType}) = MM_REMOTE
alternate(::Type{TURemoteType}) = uf_remote === TA_REMOTE ? UF_REMOTE : TA_REMOTE
alternate(::Type{MMRemoteType}) = MM_REMOTE
toggle_remote(T=TURemoteType) = begin global uf_remote = alternate(T) end
"""
    The place to store all matrix data in process
"""
const MATRIX_DB = MatrixDatabase()

# local storage directory
DATA_DIR = abspath(dirname(@__FILE__),"..", "data")
MY_DEPOT_DIR = abspath(dirname(@__FILE__), "..", "myMatrixDepot")


text:
#
# Stress test to check for memory leaks (#24)
#

using Mongoc

client = Mongoc.Client() # default locahost:27017
db = client["testDB"]
collection=db["testCollection"]

document = Mongoc.BSON()
document["name"] = "Felipe"
document["age"] = 35
document["preferences"] = [ "Music", "Computer", "Photography" ]

# generate some entries
for i in 1:100000
    push!(collection, document)
end

# read bson documents from collection
for i in 1:10000
    for entry in collection
    end
end


text:
__precompile__()
module ProgressiveHedgingSolvers

# Standard library
using LinearAlgebra
using SparseArrays
using Distributed
using Printf

# External libraries
using Parameters
using JuMP
using StochasticPrograms
using StochasticPrograms: _WS
using MathProgBase
using ProgressMeter

import Base: show, put!, wait, isready, take!, fetch
import StochasticPrograms: StructuredModel, internal_solver, optimize_structured!, fill_solution!, solverstr

const MPB = MathProgBase

export
    ProgressiveHedgingSolver,
    Fixed,
    Adaptive,
    Serial,
    Synchronous,
    Asynchronous,
    Crash,
    StructuredModel,
    optimsolver,
    optimize_structured!,
    fill_solution!,
    get_decision,
    get_objective_value

# Include files
include("types/types.jl")
include("penalties/penalization.jl")
include("execution/execution.jl")
include("ProgressiveHedging.jl")
include("spinterface.jl")

end # module


text:
function write_pccf_examples()
    T = 70; mu = 10.0; sigma = 2.0
    Pccf.write_pccf("./out/pccfExample1.txt", T, mu, sigma)
    mu = 15.0; sigma = 3
    Pccf.write_pccf("./out/pccfExample2.txt", T, mu, sigma)
end

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
export REMDQNLearner

mutable struct REMDQNLearner{Tq,Tt,Tf,R<:AbstractRNG} <: Any
    approximator::Tq
    target_approximator::Tt
    loss_func::Tf
    min_replay_history::Int
    update_freq::Int
    update_step::Int
    target_update_freq::Int
    sampler::NStepBatchSampler
    ensemble_num::Int
    ensemble_method::Symbol
    rng::R
    # for logging
    loss::Float32
end

"""
    REMDQNLearner(;kwargs...)

See paper: [An Optimistic Perspective on Offline Reinforcement Learning](https://arxiv.org/abs/1907.04543)

# Keywords

- `approximator`::[`AbstractApproximator`](@ref): used to get Q-values of a state.
- `target_approximator`::[`AbstractApproximator`](@ref): similar to `approximator`, but used to estimate the target (the next state).
- `loss_func`: the loss function.
- `γ::Float32=0.99f0`: discount rate.
- `batch_size::Int=32`
- `update_horizon::Int=1`: length of update ('n' in n-step update).
- `min_replay_history::Int=32`: number of transitions that should be experienced before updating the `approximator`.
- `update_freq::Int=4`: the frequency of updating the `approximator`.
- `ensemble_num::Int=1`: the number of ensemble approximators.
- `ensemble_method::Symbol=:rand`: the method of combining Q values. ':rand' represents random ensemble mixture, and ':mean' is the average.
- `target_update_freq::Int=100`: the frequency of syncing `target_approximator`.
- `stack_size::Union{Int, Nothing}=4`: use the recent `stack_size` frames to form a stacked state.
- `traces = SARTS`, set to `SLARTSL` if you are to apply to an environment of `FULL_ACTION_SET`.
- `rng = Random.GLOBAL_RNG`
"""
function REMDQNLearner(;
    approximator::Tq,
    target_approximator::Tt,
    loss_func::Tf,
    stack_size::Union{Int,Nothing} = nothing,
    γ::Float32 = 0.99f0,
    batch_size::Int = 32,
    update_horizon::Int = 1,
    min_replay_history::Int = 32,
    update_freq::Int = 1,
    ensemble_num::Int = 1,
    ensemble_method::Symbol = :rand,
    target_update_freq::Int = 100,
    traces = SARTS,
    update_step = 0,
    rng = Random.GLOBAL_RNG,
) where {Tq,Tt,Tf}
    copyto!(approximator, target_approximator)
    sampler = NStepBatchSampler{traces}(;
        γ = γ,
        n = update_horizon,
        stack_size = stack_size,
        batch_size = batch_size,
    )
    REMDQNLearner(
        approximator,
        target_approximator,
        loss_func,
        min_replay_history,
        update_freq,
        update_step,
        target_update_freq,
        sampler,
        ensemble_num,
        ensemble_method,
        rng,
        0.0f0,
    )
end

Functors.functor(x::REMDQNLearner) = (Q = x.approximator, Qₜ = x.target_approximator),
y -> begin
    x = @set x.approximator = y.Q
    x = @set x.target_approximator = y.Qₜ
    x
end

function (learner::REMDQNLearner)(env)
    s = send_to_device(device(learner.approximator), state(env))
    s = Flux.unsqueeze(s, ndims(s) + 1)
    q = reshape(learner.approximator(s), :, learner.ensemble_num)
    vec(mean(q, dims = 2)) |> send_to_host
end

function RLBase.update!(learner::REMDQNLearner, batch::NamedTuple)
    Q = learner.approximator
    Qₜ = learner.target_approximator
    γ = learner.sampler.γ
    loss_func = learner.loss_func
    n = learner.sampler.n
    batch_size = learner.sampler.batch_size
    ensemble_num = learner.ensemble_num
    D = device(Q)
    # Build a convex polygon to make a combination of multiple Q-value estimates as a Q-value estimate.
    if learner.ensemble_method == :rand
        convex_polygon = rand(Float32, (1, ensemble_num))
    else
        convex_polygon = ones(Float32, (1, ensemble_num))
    end
    convex_polygon ./= sum(convex_polygon)
    convex_polygon = send_to_device(D, convex_polygon)

    s, a, r, t, s′ = (send_to_device(D, batch[x]) for x in SARTS)
    a = CartesianIndex.(a, 1:batch_size)

    target_q = Qₜ(s′)
    target_q = convex_polygon .* reshape(target_q, :, ensemble_num, batch_size)
    target_q = dropdims(sum(target_q, dims = 2), dims = 2)

    if haskey(batch, :next_legal_actions_mask)
        l′ = send_to_device(D, batch[:next_legal_actions_mask])
        target_q .+= ifelse.(l′, 0.0f0, typemin(Float32))
    end

    q′ = dropdims(maximum(target_q; dims = 1), dims = 1)
    G = r .+ γ^n .* (1 .- t) .* q′

    gs = gradient(params(Q)) do
        q = Q(s)
        q = convex_polygon .* reshape(q, :, ensemble_num, batch_size)
        q = dropdims(sum(q, dims = 2), dims = 2)[a]

        loss = loss_func(G, q)
        ignore() do
            learner.loss = loss
        end
        loss
    end

    update!(Q, gs)
end


text:
module Localize

using Lints
using Fermi.Orbitals: LocalOrbitals

function boys!(ints::IntegralHelper; thresh=1E-4) where LO
    orbs = ints.orbs
    nao = length(orbs["*"][1].C)
end
end #module


text:
using Turkie
using Test

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