text:
import LaTeXStrings
import PGFPlotsX
import Statistics

"""
"""
function probability_calibration_scores_and_fractions(
        estimator::AbstractFittable,
        features_df::DataFrames.AbstractDataFrame,
        labels_df::DataFrames.AbstractDataFrame,
        single_label_name::Symbol,
        positive_class::AbstractString;
        window::Real = 0.1,
        )
    ytrue = Int.(
        singlelabelbinaryytrue(
            labels_df[single_label_name],
            positive_class,
            )
        )
    predictedprobabilitiesalllabels = predict_proba(estimator, features_df)
    yscore = Float64.(
        singlelabelbinaryyscore(
            predictedprobabilitiesalllabels[single_label_name],
            positive_class,
            )
        )
    scores, fractions = probability_calibration_scores_and_fractions(
        ytrue,
        yscore;
        window = window,
        )
    return scores, fractions
end

"""
"""
function probability_calibration_scores_and_fractions(
        ytrue::AbstractVector{<: Integer},
        yscore::AbstractVector{F};
        window::Real = 0.1,
        )::Tuple{Vector{F}, Vector{F}} where
        F <: AbstractFloat
    scores = sort(
        unique(
            vcat(
                yscore,
                zero(eltype(yscore)),
                one(eltype(yscore)),
                )
            )
        )
    fractions = fill(F(0), size(scores))
    num_rows = fill(Int(0), length(scores))
    for k = 1:length(scores)
        rows_that_have_approximately_the_kth_score = findall(
            ( scores[k] - window ) .<= ( yscore ) .<= ( scores[k] + window )
            )
        num_rows[k] = length(rows_that_have_approximately_the_kth_score)
        if length(rows_that_have_approximately_the_kth_score) == 0
            fractions[k] = NaN
        else
            fractions[k] = Statistics.mean(
                ytrue[rows_that_have_approximately_the_kth_score]
                )
        end
    end
    nonzero_indices = findall(num_rows .!= 0)
    scores_nonzero_indices::Vector{F} = scores[nonzero_indices]
    fractions_nonzero_indices::Vector{F} = fractions[nonzero_indices]
    return scores_nonzero_indices, fractions_nonzero_indices
end

"""
"""
function plot_probability_calibration_curve(
        estimator::AbstractFittable,
        features_df::DataFrames.AbstractDataFrame,
        labels_df::DataFrames.AbstractDataFrame,
        single_label_name::Symbol,
        positive_class::AbstractString;
        window::Real = 0.1,
        multiply_by::Real = 1.0,
        legend_pos::AbstractString = "outer north east",
        )
    scores, fractions = probability_calibration_scores_and_fractions(
        estimator,
        features_df,
        labels_df,
        single_label_name,
        positive_class;
        window = window,
        )
    result = plot_probability_calibration_curve(
        scores,
        fractions;
        multiply_by = multiply_by,
        legend_pos = legend_pos,
        )
    return result
end

"""
"""
function plot_probability_calibration_curve(
        scores::AbstractVector{F},
        fractions::AbstractVector{F};
        multiply_by::Real = 1.0,
        legend_pos::AbstractString = "outer north east",
        )::PGFPlotsXPlot where F <: AbstractFloat
    legend_pos::String = convert(String, legend_pos)
    scores = convert(Vector{F}, scores * multiply_by)
    fractions = convert(Vector{F}, fractions * multiply_by)
    score_versus_fraction_linearplotobject = PGFPlotsX.@pgf(
        PGFPlotsX.Plot(
            {
                only_marks,
                style = "black, fill = black",
                },
            PGFPlotsX.Coordinates(
                scores,
                fractions,
                )
            )
        )
    zero_coordinate = convert(F, zero(F))
    one_coordinate = convert(F, one(F) * multiply_by)
    perfectline_xs = [
        zero_coordinate,
        one_coordinate,
        ]
    perfectline_ys = [
        zero_coordinate,
        one_coordinate,
        ]
    perfectline_linearplotobject = PGFPlotsX.@pgf(
        PGFPlotsX.Plot(
            {
                no_marks,
                style = "dotted, red, color=red, fill=red",
                },
            PGFPlotsX.Coordinates(
                perfectline_xs,
                perfectline_ys,
                ),
            )
        )
    estimated_intercept, estimated_x_coefficient =
        simple_linear_regression(
            scores, # X
            fractions, # Y
            )
    bestfitline_xs = [
        zero_coordinate,
        one_coordinate,
        ]
    bestfitline_ys = [
        estimated_intercept,
        estimated_intercept + estimated_x_coefficient * one_coordinate,
        ]
    bestfitline_linearplotobject = PGFPlotsX.@pgf(
        PGFPlotsX.Plot(
            {
                no_marks,
                style = "dashed, blue, color=blue, fill=blue",
                },
            PGFPlotsX.Coordinates(
                bestfitline_xs,
                bestfitline_ys,
                ),
            )
        )
    p = PGFPlotsX.@pgf(
        PGFPlotsX.Axis(
            {
                xlabel = LaTeXStrings.LaTeXString(
                    "Probability of positive class"
                    ),
                ylabel = LaTeXStrings.LaTeXString(
                    "Fraction of positive class"
                    ),
                legend_pos = legend_pos,
                },
            score_versus_fraction_linearplotobject,
            perfectline_linearplotobject,
            bestfitline_linearplotobject,
            ),
        )
    wrapper = PGFPlotsXPlot(p)
    return wrapper
end

"""
"""
function probability_calibration_metrics(
        estimator::AbstractFittable,
        features_df::DataFrames.AbstractDataFrame,
        labels_df::DataFrames.AbstractDataFrame,
        single_label_name::Symbol,
        positive_class::AbstractString;
        window::Real = 0.1,
        )
    vectorofestimators = AbstractFittable[estimator]
    result = probability_calibration_metrics(
        vectorofestimators,
        features_df,
        labels_df,
        single_label_name,
        positive_class,
        )
    return result
end

"""
"""
function probability_calibration_metrics(
        vectorofestimators::AbstractVector{AbstractFittable},
        features_df::DataFrames.AbstractDataFrame,
        labels_df::DataFrames.AbstractDataFrame,
        single_label_name::Symbol,
        positive_class::AbstractString,
        window::Real = 0.1,
        )
    result = DataFrames.DataFrame()
    result[:metric] = [
        "R^2 (coefficient of determination)",
        "Brier score (binary formulation)",
        "Best fit line: estimated intercept",
        "Best fit line: estimated coefficient",
        ]
    for i = 1:length(vectorofestimators)
        ytrue = Int.(
            singlelabelbinaryytrue(
                labels_df[single_label_name],
                positive_class,
                )
            )
        predictedprobabilitiesalllabels = predict_proba(
                vectorofestimators[i],
                features_df,
                )
        yscore = Float64.(
            singlelabelbinaryyscore(
                predictedprobabilitiesalllabels[single_label_name],
                positive_class,
                )
        )
        binary_brier_score_value = binary_brier_score(ytrue, yscore)
        scores, fractions = probability_calibration_scores_and_fractions(
            ytrue,
            yscore;
            window = window,
            )
        r2_score_value = r2_score(scores, fractions)
        estimated_intercept, estimated_x_coefficient =
            simple_linear_regression(
                Float64.(scores), # X
                Float64.(fractions), # Y
                )
        result[Symbol(vectorofestimators[i].name)] = [
            r2_score_value,
            binary_brier_score_value,
            estimated_intercept,
            estimated_x_coefficient,
            ]
    end
    return result
end


text:
int main() {
    int x;
    int y;

    x = 10;
    y = 3;

    printInt(x + y);
    printInt(x - y);
    printInt(x * y);
    printInt(x / y);
    printInt(x % y);

    printInt(-x);

    x++;
    printInt(x);

    y--;
    printInt(y);

    return 0;
}

ext: jl

max_stars_repo_name: fredefox/javalette

max_issues_repo_name: fredefox/javalette

max_forks_repo_name: fredefox/javalette


text:
module CMUDict

export cmudict

using ..Data: deps, download_and_verify

const version = "0.7b"
const cache_prefix = "https://cache.julialang.org"

function load()
  suffixes_and_hashes = [(""       , "209a8b4cd265013e96f4658632a9878103b0c5abf62b50d4ef3ae1be226b29e4"),
                        (".phones" , "ffb588a5e55684723582c7256e1d2f9fadb130011392d9e59237c76e34c2cfd6"),
                        (".symbols", "408ccaae803641c6d7b626b6299949320c2dbca96b2220fd3fb17887b023b027")]
  if isdir(deps("cmudict"))
    if all(isfile(deps("cmudict", "cmudict$x")) for (x, _) in suffixes_and_hashes)
      return
    end
  end
  @info "Downloading CMUDict dataset"
  mkpath(deps("cmudict"))
  for (x, hash) in suffixes_and_hashes
    download_and_verify("$cache_prefix/https://svn.code.sf.net/p/cmusphinx/code/trunk/cmudict/cmudict-$version$x",
             deps("cmudict", "cmudict$x"), hash)
  end
end

function phones()
  load()
  Symbol.(first.(split.(split(read(deps("cmudict", "cmudict.phones"),String),
                        "\n", keepempty = false), "\t")))
end

function symbols()
  load()
  Symbol.(split(read(deps("cmudict", "cmudict.symbols"),String),
                "\n", keepempty = false))
end

function rawdict()
  load()
  Dict(String(xs[1]) => Symbol.(xs[2:end]) for xs in
       filter(!isempty, split.(split(read(deps("cmudict", "cmudict"),String), "\n"))))
end

validword(s) = isascii(s) && occursin(r"^[\w\-\.]+$", s)

cmudict() = filter(p -> validword(p.first), rawdict())

alphabet() = ['A':'Z'..., '0':'9'..., '_', '-', '.']

end


text:
import EnergyExpressions: showcoeff

"""
    LinearCombination(Ts::Vector{T}, coeffs::Vector{N})

Represents a general linear combination of objects of type `T`.
"""
struct LinearCombination{T,N<:Number}
    Ts::Vector{T}
    coeffs::Vector{N}
    function LinearCombination(Ts::Vector{T}, coeffs::Vector{N}) where {T,N}
        sel = .!iszero.(coeffs)
        new{T,N}(Ts[sel], coeffs[sel])
    end
end

Base.length(lc::LinearCombination) = length(lc.Ts)
Base.eltype(lc::LinearCombination{T,N}) where {T,N} = (T,N)

Base.:(==)(a::LinearCombination, b::LinearCombination) =
    a.Ts == b.Ts && a.coeffs == b.coeffs

Base.hash(lc::LinearCombination, h::UInt) =
    hash(lc.Ts, hash(lc.coeffs, h))

Base.zero(::LinearCombination{T,N}) where {T,N} =
    LinearCombination(Vector{T}(), Vector{N}())

Base.getindex(lc::LinearCombination, i) =
    (lc.Ts[i],lc.coeffs[i])

Base.iterate(lc::LinearCombination, i=1) =
    i > length(lc) ? nothing : (lc[i], i+=1)

function Base.show(io::IO, lc::LinearCombination)
    n = length(lc)
    for (i,(t,c)) in enumerate(lc)
        showcoeff(io, c, i > 1)
        write(io, " ")
        show(io, t)
        i < n && write(io, " ")
    end
end

Base.:(+)(A::LinearCombination, B::LinearCombination) =
    LinearCombination(vcat(A.Ts, B.Ts), vcat(A.coeffs, B.coeffs))

Base.:(-)(A::LinearCombination, B::LinearCombination) =
    LinearCombination(vcat(A.Ts, B.Ts), vcat(A.coeffs, -1 * B.coeffs))

Base.:(+)(A::LC, B::T) where {T,N,LC<:LinearCombination{<:T,N}} =
    LinearCombination(vcat(A.Ts, B), vcat(A.coeffs, one(N)))

Base.:(+)(A::T, B::LC) where {T,N,LC<:LinearCombination{<:T,N}} =
    LinearCombination(vcat(A, B.Ts), vcat(one(N), B.coeffs))

Base.:(-)(A::LC) where {LC<:LinearCombination} =
    LinearCombination(A.Ts, -1 * A.coeffs)

Base.:(*)(A::LinearCombination, b::Number) =
    LinearCombination(A.Ts, b * A.coeffs)

Base.:(*)(a::Number, B::LinearCombination) =
    LinearCombination(B.Ts, a * B.coeffs)

Base.:(/)(A::LinearCombination, b::Number) =
    LinearCombination(A.Ts, A.coeffs/b)

"""
    @linearly_combinable TT

Turns the type `TT` into a linearly combinable type, i.e. defines
arithmetic operators.

# Examples

```julia
julia> @linearly_combinable Symbol

julia> 4*(:x) - 5*(:y)
4 :x - 5 :y
```
"""
macro linearly_combinable(TT)
    quote
        Base.:(+)(A::T, B::U) where {T<:$(esc(TT)),U<:$(esc(TT))} =
            LinearCombination(promote_type(T,U)[A,B], [1,1])

        Base.:(-)(A::$(esc(TT))) =
            LinearCombination($(esc(TT))[A], [-1])

        Base.:(*)(a::Number, B::$(esc(TT))) =
            LinearCombination([B], [a])

        Base.:(*)(A::$(esc(TT)), b::Number) =
            LinearCombination([A], [b])
    end
end


text:
for (elty, ext) in ((:Float32, :s),
                    (:Float64, :d),
                    (:Complex64, :c),
                    (:Complex128, :z))

    @eval begin
        function print(A::DistMatrix{$elty}, title::String)
            err = ccall(($(string("ElPrintDist_", ext)), libEl), Cuint,
                (Ptr{Void}, Ptr{UInt8}),
                A.obj, title)
            err == 0 || throw(ElError(err))
        end
    end
end


text:
@coeff function K_T(reactor::Reactor)
  K_beta / K_G
end


text: # Dummy file just for CI to recognize a valid Julia package.


text:
# RELOG: Reverse Logistics Optimization
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using DataFrames
using CSV

function transportation_report(solution)::DataFrame
    df = DataFrame()
    df."source type" = String[]
    df."source location name" = String[]
    df."source latitude (deg)" = Float64[]
    df."source longitude (deg)" = Float64[]
    df."destination type" = String[]
    df."destination location name" = String[]
    df."destination latitude (deg)" = Float64[]
    df."destination longitude (deg)" = Float64[]
    df."product" = String[]
    df."year" = Int[]
    df."distance (km)" = Float64[]
    df."amount (tonne)" = Float64[]
    df."amount-distance (tonne-km)" = Float64[]
    df."transportation cost (\$)" = Float64[]
    df."transportation energy (GJ)" = Float64[]

    T = length(solution["Energy"]["Plants (GJ)"])
    for (dst_plant_name, dst_plant_dict) in solution["Plants"]
        for (dst_location_name, dst_location_dict) in dst_plant_dict
            for (src_plant_name, src_plant_dict) in dst_location_dict["Input"]
                for (src_location_name, src_location_dict) in src_plant_dict
                    for year = 1:T
                        push!(
                            df,
                            [
                                src_plant_name,
                                src_location_name,
                                round(src_location_dict["Latitude (deg)"], digits = 6),
                                round(src_location_dict["Longitude (deg)"], digits = 6),
                                dst_plant_name,
                                dst_location_name,
                                round(dst_location_dict["Latitude (deg)"], digits = 6),
                                round(dst_location_dict["Longitude (deg)"], digits = 6),
                                dst_location_dict["Input product"],
                                year,
                                round(src_location_dict["Distance (km)"], digits = 2),
                                round(
                                    src_location_dict["Amount (tonne)"][year],
                                    digits = 2,
                                ),
                                round(
                                    src_location_dict["Amount (tonne)"][year] *
                                    src_location_dict["Distance (km)"],
                                    digits = 2,
                                ),
                                round(
                                    src_location_dict["Transportation cost (\$)"][year],
                                    digits = 2,
                                ),
                                round(
                                    src_location_dict["Transportation energy (J)"][year] /
                                    1e9,
                                    digits = 2,
                                ),
                            ],
                        )
                    end
                end
            end
        end
    end
    return df
end

write_transportation_report(solution, filename) =
    CSV.write(filename, transportation_report(solution))

max_stars_repo_name: ANL-CEEESA/RELOG

max_issues_repo_name: ANL-CEEESA/RELOG

max_forks_repo_name: ANL-CEEESA/RELOG


text:
import Base.string
using Base.Test, MuKanren, FactCheck

# These tests are taken from interactive tutorial from miniKanren website (minikanren.org), currently not all of that tests are implemented here
facts("Tests for run/run* macro commands") do
  context("run") do
    #println(macroexpand(:(@fresh((x, y, z), equals(x, z), equals(3, y)))))
    #println(macroexpand(:(@run(1, (q), @fresh((x, y, z), equals(x, z), equals(3, y))))))
    @fact @run(1, (q), @fresh((x, y, z), equals(x, z), equals(3, y))) --> ["_.0"]
    #println(result)
    @fact @run(1, (q), @fresh((x, y), equals(x, q), equals(3, y))) --> ["_.0"]
    #(run 1 (y)
    #  (fresh (x z)
    #    (== x z)
    #    (== 3 y)))

    @fact @run(1, (y), @fresh((x, z), equals(x, z), equals(3, y))) --> [3]
    @fact @run(1, (q), @fresh((x, z), equals(x, z), equals(3, z), equals(q, x))) --> [3]
    #println(macroexpand(:( @run(1, (y), @fresh((x, y), equals(4, x), equals(x, y)), equals(3, y)) )))
    @fact @run(1, (y), @fresh((x, y), equals(4, x), equals(x, y)), equals(3, y)) --> [3]

    @fact @run(1, x, equals(4, 3)) --> []
    @fact @run(1, x, equals(x, 5), equals(x, 6)) --> []
  end
  context("complex run") do
    result = @run(2, q,
      @fresh((w, x, y),
        @conde(
            (
              equals(list(x, w, x), q),
              equals(y, w)
            ),
            (
              equals(list(w, x, w), q),
              equals(y, w)
              )
          )))
    @fact string(result) --> string([list("_.0", "_.1", "_.0"), list("_.0", "_.1", "_.0")])

  end
  context("infinite loop/stackoverflow") do
  #(run* (q)
#(let loop ()
#  (conde
#    ((== #f q))
#    ((== #t q))
#    ((loop)))))
  end

  context("anyo") do
    #println(macroexpand(:( @conde(g(), anyo(g)))))
    anyo = g -> @conde(g, anyo(g))

    @fact @run(10, q, anyo(@conde(
        equals(1, q),
        equals(2, q),
        equals(3, q)
        ))) --> [1, 2, 3, 1, 2, 3, 1, 2, 3, 1]
    #TODO: add example call with nevero and infinite loop

  end
end


text:
# import Particles package
using Particles

# import HDF5 package for storing solution
using HDF5

# parameters
Δt = 0.1    # time step size
nt = 100    # number of time steps
ns = 5      # number of samples

# output file
h5file = "pendulum.hdf5"

# random initial conditions
x₀ = [ [rand() .* 2π, rand() .* 2 .- 1] for i in 1:ns ]

# vector field
function pendulum!(ẋ, x)
    for i in eachindex(ẋ, x)
        ẋ[i][1] = x[i][2]
        ẋ[i][2] = - sin(x[i][1])
    end
end

# solution storage
function copy_to_hdf5(h5x, x, n)
    h5x[:,:,n+1] = reduce(hcat, x)
end

# create an Equation instance
equ = Equation(pendulum!, x₀)

# create HDF5 file and copy initial conditions
h5  = h5open(h5file, "w")
h5x = d_create(h5, "x", eltype(x₀[1]), ((2, ns, nt+1), (2, ns, -1)), "chunk", (2,ns,1))
copy_to_hdf5(h5x, x₀, 0)

# create Simulation instance
sim = Simulation(equ, ExplicitEuler, Δt, nt; f_postproc = (x, n) -> copy_to_hdf5(h5x, x, n))

# run simulation
try
    z = sim()
finally
    close(h5)
end

# load Plots package
using Plots

# read array from HDF5 file
x = h5read(h5file, "x")

# set plot ranges
xlim = (-π, +π)
ylim = (-3, +3)

# plot energy contour lines
plt = contour(
    LinRange(xlim..., 100),
    LinRange(ylim..., 100),
    (x,v) -> v^2 / 2  + 1 - cos(x),
    xlim = xlim,
    ylim = ylim,
    title = "Pendulum",
    legend = false,
    size = (800, 600)
)

# plot solutions
for i in axes(x,2)
    scatter!(plt, (x[1,i,:] .+ π) .% 2π .- π, x[2,i,:], marker=5)
end

# save figure to file
savefig("pendulum.png")

max_stars_repo_name: michakraus/JuliaAtHomeAssignments

max_issues_repo_name: michakraus/JuliaAtHomeAssignments

max_forks_repo_name: michakraus/JuliaAtHomeAssignments


text:
module TestGetDataframe

using Mimi
using Test

#------------------------------------------------------------------------------
#   1. Test with 1 dimension
#------------------------------------------------------------------------------

model1 = Model()

@defcomp testcomp1 begin
    var1 = Variable(index=[time])
    par1 = Parameter(index=[time])
    par_scalar = Parameter()

    function run_timestep(p, v, d, t)
        v.var1[t] = p.par1[t]
    end
end

late_first = 2030
early_last = 2100

@defcomp testcomp2 begin
    var2 = Variable(index=[time])
    par2 = Parameter(index=[time])

    function run_timestep(p, v, d, t)
        if TimestepValue(late_first) <= t <= TimestepValue(early_last)       # apply time constraints in the component
            v.var2[t] = p.par2[t]
        end
    end
end

years = collect(2015:5:2110)

set_dimension!(model1, :time, years)
add_comp!(model1, testcomp1)
set_param!(model1, :testcomp1, :par1, years)
set_param!(model1, :testcomp1, :par_scalar, 5.)

add_comp!(model1, testcomp2)

@test_throws ErrorException set_param!(model1, :testcomp2, :par2, late_first:5:early_last)
@test ! (:par2 in keys(model1.md.external_params))  # Test that after the previous error, the :par2 didn't stay in the model's parameter list
set_param!(model1, :testcomp2, :par2, years)

# Test running before model built
@test_throws ErrorException df = getdataframe(model1, :testcomp1, :var1)

# Now run model
run(model1)

# Test trying to load a dataframe for a scalar parameter
@test_throws ErrorException getdataframe(model1, :testcomp1, :par_scalar)

# Test trying to getdataframe from a variable that does not exist in the component
@test_throws ErrorException getdataframe(model1, :testcomp1, :var2)

# Regular getdataframe
df = getdataframe(model1, :testcomp1=>:var1, :testcomp1=>:par1, :testcomp2=>:var2, :testcomp2=>:par2)

dim = Mimi.dimension(model1, :time)
@test df.var1 == df.par1 == years
@test all(ismissing, df.var2[1 : dim[late_first]-1])
#@test all(ismissing, df.par2[1 : dim[late_first]-1])
@test df.var2[dim[late_first] : dim[early_last]] == df.par2[dim[late_first] : dim[early_last]] == late_first:5:early_last
@test all(ismissing, df.var2[dim[years[end]] : dim[early_last]])
@test all(ismissing, df.par2[dim[years[end]] : dim[early_last]])

# Test trying to load an item into an existing dataframe where that item key already exists
@test_throws UndefVarError _load_dataframe(model1, :testcomp1, :var1, df)

#------------------------------------------------------------------------------
#   2. Test with > 2 dimensions
#------------------------------------------------------------------------------
stepsize = 5
years = collect(2015:stepsize:2110)
regions = [:reg1, :reg2]
rates   = [0.025, 0.05]

nyears = length(years)
nregions = length(regions)
nrates = length(rates)

@defcomp testcomp3 begin
    par3 = Parameter(index=[time, regions, rates])
    var3 = Variable(index=[time])
end

# A. Simple case where component has same time length as model

model2 = Model()

set_dimension!(model2, :time, years)
set_dimension!(model2, :regions, regions)
set_dimension!(model2, :rates, rates)

data = Array{Int}(undef, nyears, nregions, nrates)
data[:] = 1:(nyears * nregions * nrates)

add_comp!(model2, testcomp3)
set_param!(model2, :testcomp3, :par3, data)

run(model2)

df2 = getdataframe(model2, :testcomp3, :par3)
@test size(df2) == (length(data), 4)

# Test trying to combine two items with different dimensions into one dataframe
@test_throws ErrorException getdataframe(model2, Pair(:testcomp3, :par3), Pair(:testcomp3, :var3))

# B. Test with shorter time than model

model3 = Model()
set_dimension!(model3, :time, years)
set_dimension!(model3, :regions, regions)
set_dimension!(model3, :rates, rates)

dim = Mimi.dimension(model3, :time)

late_first = 2030
early_last = 2100

add_comp!(model3, testcomp3)

indices = collect(late_first:stepsize:early_last)
nindices = length(indices)

valid_indices = collect(dim[late_first]:dim[early_last])
nvalid = length(valid_indices)

par3 = Array{Union{Missing,Float64}}(undef, nyears, nregions, nrates)
par3[:] .= missing

par3[valid_indices, :, :] = 1:(nindices * nregions * nrates)
set_param!(model3, :testcomp3, :par3, par3)
run(model3)

df3 = getdataframe(model3, :testcomp3 => :par3)
@test size(df3) == (nrates * nregions * nyears, 4)

# Test that times outside the component's time span are padded with `missing` values
@test all(ismissing, df3.par3[1 : (nrates * nregions * (dim[late_first] - 1))])

nmissing = (Int((years[end] - early_last) / stepsize) * nregions * nrates - 1)

@test all(ismissing, df3.par3[end - nmissing : end])

end #module


text:
const CLEXAMPLE = CoordinateList(:AC => (127.0, 119.0, 149.0),
                                 :PC => (128.0, 93.0, 141.0),
                                 :IH => (131.0, 114.0, 206.0))

img = NeuroMetaArray(rand(4,4))
props = NeuroMetadata()

for (property, default, valin, valout) in (
    (:acquisition_duration, 1.0u"s", 2.0u"s", 2.0u"s"),
    (:anatomical_landmark_coordinates, CoordinateList(), CLEXAMPLE, CLEXAMPLE),
    #:AnatomicalLandmarkCoordinateSystem,
    #:AnatomicalLandmarkCoordinateUnits,
    #:AnatomicalLandmarkCoordinateSystemDescription,
    #:CogAtlasID,
    #:CogPOID,
    (:coil_combination_method, "rSOS", "", ""),
    (:contrast_bolus_ingredient, UnkownContrast, "IODINE", IODINE),
    (:delay_after_trigger, 1.0u"s", 2.0u"s", 2.0u"s"),
    (:delay_time, 1.0u"s", 2.0u"s", 2.0u"s"),
    (:device_serial_number, "00000", "11111", "11111"),
    (:dewar_position, OneIntDeg, 2, 2u"°"),
    (:dwell_time, 1.0u"s", 2.0u"s", 2.0u"s"),
    (:echo_time, 1.0u"s", 2.0u"s", 2.0u"s"),
    (:effective_echo_spacing, 1.0u"s", 2.0u"s", 2.0u"s"),
    (:eeg_electrode_groups, Dict{Symbol,String}(), Dict(:grid1 => "10x8 grid on left temporal pole", :strip2 => "1x8 electrode strip on xxx"),  Dict(:grid1 => "10x8 grid on left temporal pole", :strip2 => "1x8 electrode strip on xxx")),
    (:eeg_ground, "", "mastoid", "mastoid"),
    (:eeg_placement_scheme, "", "surface strip and STN depth", "surface strip and STN depth"),
    (:electrical_stimulation, false, true, true),
    (:electrical_stimulation_parameters, "", "x", "x"),
    (:epoch_length, 1.0u"s", 2.0u"s", 2.0u"s"),
    (:fiducial_description, "", "pre-auricular","pre-auricular"),
    (:flip_angle,  OneIntDeg, 2, 2u"°"),
    (:gradient_set_type, "", "gradient", "gradient"),
    (:head_coil_coordinates, CoordinateList(), CLEXAMPLE, CLEXAMPLE),
    #:head_coil_coordinate_system,
    #:head_coil_coordinate_system_description,
    #:head_coil_coordinate_units,
    (:institution_address, "", "x", "x"),
    (:institution_name, "", "x", "x"),
    (:institutional_department_name, "", "x", "x"),
    (:instructions, "", "x", "x"),
    (:intended_for, String[], ["x"], ["x"]),
    (:inversion_time, 1.0u"s", 2.0u"s", 2.0u"s"),
    (:magnetic_field_strength, 3.0u"T", 1.5u"T", 1.5u"T"),
    (:manufacturer, "", "x", "x"),
    (:manufacturer_model_name, "", "x", "x"),
    (:matrix_coil_mode, "", "x", "x"),
    (:mr_transmit_coil_sequence, "", "x", "x"),
    (:multiband_acceleration_factor, "", "x", "x"),
    (:negative_contrast, false, true, true),
    (:number_of_volumes_discarded_by_scanner, 0, 1, 1),
    (:number_of_volumes_discarded_by_user, 0, 1, 1),
    (:number_shots, 0, 1, 1),
    (:nonlinear_gradient_correction, false, true, true),
    (:parallel_acquisition_technique, "", "x", "x"),
    (:parallel_reduction_factor, 0, 1, 1),
    (:partial_fourier, 1.0, 2.0, 2.0),
    (:partial_fourier_direction, "", "COMBINATION", "COMBINATION"),
    (:phase_encoding_direction, ipos, 2, jpos),
    (:power_line_frequency, 1u"Hz", 2u"Hz", 2u"Hz"),
    (:pulse_sequence, "", "x", "x"),
    (:pulse_sequence_details, "", "x", "x"),
    (:pulse_sequence_type, "", "x", "x"),
    (:receive_coil_active_elements, "", "x", "x"),
    (:receive_coil_name, "", "x", "x"),
    (:recording_duration, 1u"s", 2u"s", 2u"s"),
    (:recording_type, "", "x", "x"),
    (:repetition_time, 1u"s", 2u"s", 2u"s"),
    (:sampling_frequency, 1u"Hz", 2u"Hz", 2u"Hz"),
    (:scan_options, Dict{Symbol,Any}(), Dict{Symbol,Any}(:randmeta => 1), Dict{Symbol,Any}(:randmeta => 1)),
    (:scanning_sequence, "", "x", "x"),
    (:sequence_name, "", "x", "x"),
    (:sequence_varient, "", "x", "x"),
    (:subject_artefact_description, "", "x", "x"),
    (:slice_encoding_direction, ipos, 2, jpos),
    (:slice_timing, F64Sec[], [2u"s"], [2u"s"]),
    (:software_versions, "", "x", "x"),
    #(:start_time, 1u"s", 2u"s", 2u"s"),
    (:station_name, "", "x", "x"),
    (:task_description, "", "x", "x"),
    (:task_name, "", "x", "x"),
    (:total_readout_time, 1u"s", 2u"s", 2u"s"),
    (:volume_timing, F64Sec[], [2u"s"], [2u"s"]),
   )
    @testset "$property" begin
        eval(:(fxn(x) = getproperty(x, $(QuoteNode(property)))))
        eval(:(fxn!(x, val) = setproperty!(x, $(QuoteNode(property)), val)))

        @test @inferred(fxn(img)) == default
        @test @inferred(fxn(props)) == default
        if !in(property, NeuroCore.NO_SET_PROPERTIES)
            fxn!(img, valin)
            fxn!(props, valin)
            @test @inferred(fxn(img)) == valout
            @test @inferred(fxn(props)) == valout
        end
    end
end


text:
using NeuralPDE
using Flux
using DifferentialEquations
using LinearAlgebra
d = 100 # number of dimensions
X0 = fill(0.0f0, d) # initial value of stochastic control process
tspan = (0.0f0, 1.0f0)
λ = 1.0f0

g(X) = log(0.5f0 + 0.5f0 * sum(X .^ 2))
f(X, u, σᵀ∇u, p, t) = -λ * sum(σᵀ∇u .^ 2)
μ_f(X, p, t) = zero(X)  # Vector d x 1 λ
σ_f(X, p, t) = Diagonal(sqrt(2.0f0) * ones(Float32, d)) # Matrix d x d
prob = TerminalPDEProblem(g, f, μ_f, σ_f, X0, tspan)
hls = 10 + d # hidden layer size
opt = Flux.ADAM(0.01)  # optimizer
# sub-neural network approximating solutions at the desired point
u0 = Flux.Chain(Dense(d, hls, relu),
    Dense(hls, hls, relu),
    Dense(hls, 1))
# sub-neural network approximating the spatial gradients at time point
σᵀ∇u = Flux.Chain(Dense(d + 1, hls, relu),
    Dense(hls, hls, relu),
    Dense(hls, hls, relu),
    Dense(hls, d))
pdealg = NNPDENS(u0, σᵀ∇u, opt=opt)
@time ans = solve(prob, pdealg, verbose=true, maxiters=100, trajectories=100,
    alg=EM(), dt=1.2, pabstol=1.0f-2)
using OrdinaryDiffEq
function lorenz(du, u, p, t)
    du[1] = 10.0(u[2] - u[1])
    du[2] = u[1] * (28.0 - u[3]) - u[2]
    du[3] = u[1] * u[2] - (8 / 3) * u[3]
end
u0 = [1.0; 0.0; 0.0]
tspan = (0.0, 100.0)
prob = ODEProblem(lorenz, u0, tspan)
sol = solve(prob, Tsit5())
sol = solve(prob, Tsit5())

using Profile
using ProfileView
# Profile 1000 runs
ProfileView.@profview for i in 1:1000
    sol = solve(prob, Tsit5())
end
vscode.profiler()


text:
include("getdigit_n.jl")
include("myprintf.jl")
using Printf
using .GetDigit_n
using .MyPrintf

const DIGIT = 5
const MAXN = 500000

function basel(n)
    result = BigFloat(0)

    for i = 1:n
        result += BigFloat(1) / (BigFloat(i * i))
    end

    return sqrt(BigFloat(6) * result)
end

function main()
    n = GetDigit_n.getn(basel, MAXN, DIGIT, true)
    precision = GetDigit_n.getdigit(basel, n, DIGIT, true)
    precision = GetDigit_n.getdigit_2(basel, n, DIGIT, precision)

    setprecision(precision)
    MyPrintf.myprintf("π（計算値） = ", DIGIT, basel(n))
    MyPrintf.myprintf("π（厳密値） = ", DIGIT, BigFloat(π))
    @printf("n = %d, 計算%d桁で、%d桁求まりました\n", n, precision, DIGIT)
end

main()


text:
using Documenter
using CorrelationFunctions

makedocs(sitename = "CorrelationFunctions.jl documentation",
         format   = Documenter.HTML(prettyurls = false,
                                    assets     = [
                                        "assets/code.css"
                                    ]))

deploydocs(repo = "github.com/shamazmazum/CorrelationFunctions.jl.git")


text:
using Images

function cubicbezier!(xy::Matrix,
                      img::Matrix = fill(RGB(255.0, 255.0, 255.0), 17, 17),
                      col::ColorTypes.Color = convert(eltype(img), Gray(0.0)),
                      n::Int = 20)
    t = collect(0:n) ./ n
    M = hcat((1 .- t) .^ 3, # a
             3t .* (1 .- t) .^ 2, # b
             3t .^ 2 .* (1 .- t), # c
             t .^ 3) # d
    p = floor.(Int, M * xy)
    for i in 1:n
        drawline!(img, p[i, :]..., p[i+1, :]..., col)
    end
    return img
end

xy = [16 1; 1 4; 3 16; 15 11]
cubicbezier!(xy)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
export skewProductFun, skewpoints, skewtransform!, iskewtransform!

#
# These are bivariate constructors that avoid sampling the diagonal y = x
# in a TensorSpace of two of the same spaces.
#

function skewProductFun{D1,D2}(f::Function,sp::TensorSpace{Tuple{Chebyshev{D1},Chebyshev{D2}}};tol=100eps())
    for logn = 4:10
        X = coefficients(skewProductFun(f,sp,2^logn,2^logn+1;tol=tol))
        if size(X,1)<2^logn && size(X,2)<2^logn+1
            return ProductFun(X,sp;tol=tol,chopping=true)
        end
    end
    warn("Maximum grid size of ("*string(2^11)*","*string(2^11+1)*") reached")
    skewProductFun(f,sp,2^11,2^11+1;tol=tol)
end

function skewProductFun{D1,D2}(f::Function,sp::TensorSpace{Tuple{Laurent{D1},Laurent{D2}}};tol=100eps())
    for logn = 4:10
        X = coefficients(skewProductFun(f,sp,2^logn,2^logn;tol=tol))
        if size(X,1)<2^logn && size(X,2)<2^logn
            return ProductFun(X,sp;tol=tol,chopping=true)
        end
    end
    warn("Maximum grid size of ("*string(2^11)*","*string(2^11)*") reached")
    skewProductFun(f,sp,2^11,2^11;tol=tol)
end

function skewProductFun(f::Function,S::TensorSpace,M::Integer,N::Integer;tol=100eps())
    xy = ApproxFun.checkpoints(S)
    T = promote_type(eltype(f(first(xy)...)),eltype(S))
    ptsx,ptsy=skewpoints(S,M,N)
    vals=T[f(ptsx[k,j],ptsy[k,j]) for k=1:size(ptsx,1), j=1:size(ptsx,2)]
    ProductFun(skewtransform!(S,vals),S;tol=tol,chopping=true)
end

function skewtransform!{T,D1,D2}(S::TensorSpace{Tuple{Chebyshev{D1},Chebyshev{D2}}},M::Matrix{T};kindx::Int=1,kindy::Int=2)
    n=size(M,1)

    planc=plan_chebyshevtransform(M[:,1];kind=kindx)
    for k=1:size(M,2)
        M[:,k]=chebyshevtransform(M[:,k],planc;kind=kindx)
    end

    planr=plan_chebyshevtransform(vec(M[1,:]);kind=kindy)
    for k=1:n
        M[k,:]=chebyshevtransform(vec(M[k,:]),planr;kind=kindy)
    end
    M
end

function iskewtransform!{T,D1,D2}(S::TensorSpace{Tuple{Chebyshev{D1},Chebyshev{D2}}},M::Matrix{T};kindx::Int=1,kindy::Int=2)
    n=size(M,1)

    planc=plan_ichebyshevtransform(M[:,1];kind=kindx)
    for k=1:size(M,2)
        M[:,k]=ichebyshevtransform(M[:,k],planc;kind=kindx)
    end

    planr=plan_ichebyshevtransform(vec(M[1,:]);kind=kindy)
    for k=1:n
        M[k,:]=ichebyshevtransform(vec(M[k,:]),planr;kind=kindy)
    end
    M
end

function skewtransform!{T,D1,D2}(S::TensorSpace{Tuple{Laurent{D1},Laurent{D2}}},M::Matrix{T})
    n=size(M,1)

    planc=plan_transform(S[1],M[:,1])
    for k=1:size(M,2)
        M[:,k]=transform(S[1],M[:,k],planc)
    end

    planr=plan_transform(S[2],vec(M[1,:]))
    for k=1:n
        M[k,:]=process!(transform(S[2],vec(M[k,:]),planr))
    end
    M
end

function iskewtransform!{T,D1,D2}(S::TensorSpace{Tuple{Laurent{D1},Laurent{D2}}},M::Matrix{T})
    n=size(M,1)

    planc=plan_itransform(S[1],M[:,1])
    for k=1:size(M,2)
        M[:,k]=itransform(S[1],M[:,k],planc)
    end

    planr=plan_itransform(S[2],vec(M[1,:]))
    for k=1:n
        M[k,:]=itransform(S[2],iprocess!(vec(M[k,:])),planr)
    end
    M
end

skewpoints(d::TensorSpace,n,m;kwds...) = skewpoints(d,n,m,1;kwds...),skewpoints(d,n,m,2;kwds...)

function skewpoints{D1,D2}(d::TensorSpace{Tuple{Chebyshev{D1},Chebyshev{D2}}},n,m,k;kindx::Int=1,kindy::Int=2)
    ptsx=fromcanonical(d[1],chebyshevpoints(n;kind=kindx))
    ptsy=fromcanonical(d[2],chebyshevpoints(m;kind=kindy))
    promote_type(eltype(ptsx),eltype(ptsy))[fromcanonical(d,x,y)[k] for x in ptsx, y in ptsy]
end

function skewpoints{D1,D2}(d::TensorSpace{Tuple{Laurent{D1},Laurent{D2}}},n,m,k)
    ptsx=fromcanonical(d[1],fourierpoints(n))
    ptsy=fromcanonical(d[2],fourierpoints(m)+π/m)
    promote_type(eltype(ptsx),eltype(ptsy))[fromcanonical(d,x,y)[k] for x in ptsx, y in ptsy]
end

function process!{T}(c::Vector{T})
    N = length(c)
    twiddle = exp(im*π/N)
    twid = one(T)
    for n=2:2:N-1
        twid *= twiddle
        c[n] *= twid
        c[n+1] /= twid
    end
    if mod(N,2) == 0
        twid *= twiddle
        c[N] *= twid
    end
    c
end

function iprocess!{T}(c::Vector{T})
    N = length(c)
    twiddle = exp(-im*π/N)
    twid = one(T)
    for n=2:2:N-1
        twid *= twiddle
        c[n] *= twid
        c[n+1] /= twid
    end
    if mod(N,2) == 0
        twid *= twiddle
        c[N] *= twid
    end
    c
end


text:
function POMDPs.solve(solver::DeepQLearningSolver, env::AbstractEnvironment)
    train_graph = build_graph(solver, env)

    # init and populate replay buffer
    if solver.prioritized_replay
        replay = PrioritizedReplayBuffer(env, solver.buffer_size, solver.batch_size)
    else
        replay = ReplayBuffer(env, solver.buffer_size, solver.batch_size)
    end
    populate_replay_buffer!(replay, env, max_pop=solver.train_start)

    # init variables
    run(train_graph.sess, global_variables_initializer())
    # train
    dqn_train(solver, env, train_graph, replay)
    policy = DQNPolicy(train_graph.q, train_graph.s, env, train_graph.sess)
    return policy
end

function POMDPs.solve(solver::DeepQLearningSolver, problem::MDP)
    env = MDPEnvironment(problem, rng=solver.rng)
    # init session and build graph Create a TrainGraph object with all the tensors
    return solve(solver, env)
end

"""
    main training component
"""
function dqn_train(solver::DeepQLearningSolver,
                   env::AbstractEnvironment,
                   graph::TrainGraph,
                   replay::Union{ReplayBuffer, PrioritizedReplayBuffer})

    if solver.if_log == true
        summary_writer = tf.summary.FileWriter(solver.logdir)
    end

    # initialization
    obs = reset(env)
    done = false
    step = 0
    rtot = 0
    episode_rewards = Float64[0.0]
    episode_lengths = Int64[0]
    saved_best_reward = - Inf
    scores_eval = 0.
    non_corrected_scores_eval = 0.
    num_nmac = 0.
    non_corrected_num_nmac = 0.
    eval_episode_length = 0
    eval_action_stats = ones(n_actions(env)) / n_actions(env)
    eps = solver.eps_start
    weights = ones(solver.batch_size)
    model_saved = false
    policy = DQNPolicy(graph.q, graph.s, env, graph.sess)

    # train start:
    for t = 1:solver.max_steps
        # epsilon greedy policy
        if solver.train
            # if rand(solver.rng) > eps
            if rand() > eps
                action = get_action(graph, solver, env, obs) # feed into DQN
            else
                action = sample_action(env)
            end
        else
            action = get_action(graph, solver, env, obs)
        end
        # update epsilon
        if t < solver.eps_fraction * solver.max_steps
            eps = 1 - (1 - solver.eps_end) / (solver.eps_fraction * solver.max_steps) * t # decay
        else
            eps = solver.eps_end
        end
        ai = action_index(env.problem, action)

        # render_option = env.render
        # env.render = false
        op, rew, _,  done, info = step!(env, action)
        # restore render_option for visualiztion during eval
        # env.render = render_option

        exp = DQExperience(obs, ai, rew, op, done)
        add_exp!(replay, exp)
        obs = op
        step += 1
        episode_rewards[end] += rew
        episode_lengths[end] += 1

        if done || step >= solver.max_episode_length
            obs = reset(env)
            push!(episode_rewards, 0.0)
            push!(episode_lengths, 0)
            done = false
            step = 0
            rtot = 0
        end
        num_episodes = length(episode_rewards)
        avg100_reward = mean(episode_rewards[max(1, num_episodes - 101):end])

        # batch_train:
        if t % solver.train_freq == 0
            if solver.prioritized_replay
                s_batch, a_batch, r_batch, sp_batch, done_batch, indices, weights = sample(replay)
            else
                weights = ones(replay.batch_size)
                s_batch, a_batch, r_batch, sp_batch, done_batch = sample(replay)
            end

            # q_lo's are normalized
            q_lo_batch, q_lo_p_batch = get_q_lo_batch(graph, env, s_batch, sp_batch) # returns Array{Float64}

            feed_dict = Dict(graph.s => s_batch,
                             graph.a => a_batch,
                             graph.sp => sp_batch,
                             graph.r => r_batch,
                             graph.done_mask => done_batch,
                             graph.importance_weights => weights,
                             graph.q_lo => q_lo_batch,
                             graph.q_lo_p => q_lo_p_batch)

            loss_val, td_errors, grad_val, _ =
                run(graph.sess, [graph.loss, graph.td_errors, graph.grad_norm, graph.train_op], feed_dict)
            # loss_val, td_errors, grad_val, _, q =
            #     run(graph.sess, [graph.loss, graph.td_errors, graph.grad_norm, graph.train_op, graph.q], feed_dict)
        end

        if t % solver.target_update_freq == 0
            println("Target updated!")
            run(graph.sess, graph.update_op)
        end

        if t % solver.eval_freq == 0
            print(@sprintf("%5d / %5d  ", t, solver.max_steps))
            scores_eval, non_corrected_scores_eval, num_nmac, non_corrected_num_nmac, eval_episode_length, eval_action_stats =
                eval_q(graph, solver, env, t,
                    total_n_steps=solver.max_steps,
                    n_eval=solver.num_ep_eval,
                    max_episode_length=solver.max_episode_length,
                    verbose=solver.verbose)
        end

        if t % solver.log_freq == 0 && t >= solver.eval_freq
            if solver.if_log == true
                # log to tensorboard
                tb_avgr = logg_scalar(avg100_reward, "avg_reward")
                tb_evalr = logg_scalar(scores_eval, "eval_reward")
                tb_noncorrected_evalr = logg_scalar(non_corrected_scores_eval, "non_corrected_eval_reward")
                tb_num_nmac = logg_scalar(num_nmac, "num_nmac")
                tb_noncorrected_num_nmac = logg_scalar(non_corrected_num_nmac, "non_corrected_num_nmac")
                tb_loss = logg_scalar(loss_val, "loss")
                tb_tderr = logg_scalar(mean(td_errors), "mean_td_error")
                tb_grad = logg_scalar(grad_val, "grad_norm")
                if done
                    tb_epreward = logg_scalar(episode_rewards[end], "episode_reward")
                    tb_eplength = logg_scalar(episode_lengths[end], "episode_length")
                else # if the episode is not yet done
                    tb_epreward = logg_scalar(episode_rewards[end - 1], "episode_reward")
                    tb_eplength = logg_scalar(episode_lengths[end - 1], "episode_length")
                end
                tb_eps = logg_scalar(eps, "epsilon")

                tb_eval_episode_length = logg_scalar(eval_episode_length, "eval_episode_length")
                write(summary_writer, tb_eval_episode_length, t)

                for ai in 1:n_actions(env)
                    write(summary_writer,
                        logg_scalar(eval_action_stats[ai], "eval_action_freq_" * string(actions(env)[ai])), t)
                end

                write(summary_writer, tb_avgr, t)
                write(summary_writer, tb_evalr, t)
                write(summary_writer, tb_noncorrected_evalr, t)
                write(summary_writer, tb_num_nmac, t)
                write(summary_writer, tb_noncorrected_num_nmac, t)
                write(summary_writer, tb_loss, t)
                write(summary_writer, tb_tderr, t)
                write(summary_writer, tb_grad, t)
                write(summary_writer, tb_epreward, t)
                write(summary_writer, tb_eplength, t)
                write(summary_writer, tb_eps, t)
            end # if_log

            if solver.verbose
                logg = @sprintf("%5d / %5d eps %0.3f | avgR %1.3f | Loss %2.3f | Grad %2.6f | num_episodes %3d" ,
                                 t, solver.max_steps, eps, avg100_reward, loss_val, grad_val, num_episodes)
                println(logg)
            end
        end

        if t % solver.save_freq == 0 && solver.if_log == true
            if scores_eval >= saved_best_reward
                if solver.verbose
                    println("Saving new model with the best eval reward ", scores_eval)
                end
                policy = DQNPolicy(graph.q, graph.s, env, graph.sess)
                JLD.save(solver, policy,
                    weights_file=solver.logdir*"weights.jld",
                    problem_file=solver.logdir*"problem.jld")
                model_saved = true
                saved_best_reward = scores_eval
            end
        end

    end # for t = 1:solver.max_steps

    if model_saved
        if solver.verbose
            println("Restore model with eval reward ", saved_best_reward)
        end
    end

    return policy
end

"""
Evaluate a Q network
"""
function eval_q(graph::TrainGraph,
                solver::DeepQLearningSolver,
                env::AbstractEnvironment,
                t::Int64;
                total_n_steps::Int64=4000000,
                n_eval::Int64=100,
                max_episode_length::Int64=2000,
                verbose::Bool=true)
    # Evaluation
    dummy_problem = Dummy_Problem(0.99)
    avg_r = 0
    avg_num_nmac = 0
    eval_action_stats = zeros(n_actions(env))
    eval_episode_length = 0

    for i = 1 : n_eval
        done = false
        r_tot = 0.0
        step = 0
        obs = reset(env)
        while !done && step <= max_episode_length
            action = get_action(graph, solver, env, obs)
            ai = action_index(dummy_problem, action)
            eval_action_stats[ai] += 1
            obs, rew, num_nmac, done, info = step!(env, action)
            avg_num_nmac += num_nmac
            avg_r += rew
            step += 1
            eval_episode_length += 1
        end
    end
    avg_num_nmac /= n_eval
    avg_r /= n_eval
    eval_episode_length /= n_eval
    eval_action_stats /= sum(eval_action_stats)

    # Evaluate uncorrected env (using original VICAS on ego agent):
    # record the original correction option and nn_policy option
    correction_option = env.correction
    nn_policy_option = env.nn_policy
    stochastic_ego_policy_option = env.stochastic_ego_policy

    # temporarily deactivate settings
    env.correction = false
    env.nn_policy = false
    env.stochastic_ego_policy = false

    non_corrected_avg_r = 0
    non_corrected_avg_num_nmac = 0
    non_corrected_n_eval = n_eval
    for i = 1 : non_corrected_n_eval # less n_eval for non-corrected policy to save time
        done = false
        r_tot = 0.0
        step = 0
        obs = reset(env)
        while !done && step <= max_episode_length
            action = get_action(graph, solver, env, obs)
            obs, rew, num_nmac ,done, info = step!(env, action)
            non_corrected_avg_r += rew
            non_corrected_avg_num_nmac += num_nmac
            step += 1
        end
    end
    non_corrected_avg_num_nmac /= non_corrected_n_eval
    non_corrected_avg_r /= non_corrected_n_eval

    # restore correction option and nn_policy option
    env.correction = correction_option
    env.nn_policy = nn_policy_option
    env.stochastic_ego_policy = stochastic_ego_policy_option

    if verbose
        println("\nEvaluation ................. Avg Reward ", avg_r)
        println("Evaluation ... Non Corrected Avg Reward ", non_corrected_avg_r)
    end

    return  avg_r, non_corrected_avg_r, avg_num_nmac, non_corrected_avg_num_nmac,
            eval_episode_length, eval_action_stats
end

# unmodified code for recurrent net policy
##############################################################################################################

function POMDPs.solve(solver::DeepRecurrentQLearningSolver, problem::Union{MDP,POMDP})
    if !isa(problem, POMDP)
        env = MDPEnvironment(problem, rng=solver.rng)
    else
        env = POMDPEnvironment(problem, rng=solver.rng)
    end
    return solve(solver, env)
end

function POMDPs.solve(solver::DeepRecurrentQLearningSolver, env::AbstractEnvironment)
    #init session and build graph Create a TrainGraph object with all the tensors
    train_graph = build_graph(solver, env)
    # init and populate replay buffer
    replay = EpisodeReplayBuffer(env, solver.buffer_size, solver.batch_size, solver.trace_length)
    populate_replay_buffer!(replay, env, max_pop=solver.train_start)
    # init variables
    run(train_graph.sess, global_variables_initializer())
    #TODO save the training log somewhere
    drqn_train(solver, env, train_graph, replay)
    policy = train_graph.lstm_policy
    policy.sess = train_graph.sess
    return policy
end

function drqn_train(solver::DeepRecurrentQLearningSolver,
                   env::AbstractEnvironment,
                   graph::RecurrentTrainGraph,
                   replay::EpisodeReplayBuffer)
    summary_writer = tf.summary.FileWriter(solver.logdir)
    obs = reset(env)
    reset_hidden_state!(graph.lstm_policy)
    done = false
    step = 0
    rtot = 0
    episode = DQExperience[]
    sizehint!(episode, solver.max_episode_length)
    episode_rewards = Float64[0.0]
    saved_mean_reward = 0.
    model_saved = false
    scores_eval = 0.0
    eps = 1.0
    weights = ones(solver.batch_size*solver.trace_length)
    init_c = zeros(solver.batch_size, solver.arch.lstm_size)
    init_h = zeros(solver.batch_size, solver.arch.lstm_size)
    grad_val, loss_val = -1, -1 # sentinel value
    for t=1:solver.max_steps
        if rand(solver.rng) > eps
            action = get_action!(graph.lstm_policy, obs, graph.sess)
        else
            action = sample_action(env)
        end
        # update epsilon
        if t < solver.eps_fraction*solver.max_steps
            eps = 1 - (1 - solver.eps_end)/(solver.eps_fraction*solver.max_steps)*t # decay
        else
            eps = solver.eps_end
        end
        ai = action_index(env.problem, action)
        op, rew, _, done, info = step!(env, action)
        exp = DQExperience(obs, ai, rew, op, done)
        push!(episode, exp)
        obs = op
        step += 1
        episode_rewards[end] += rew
        if done || step >= solver.max_episode_length
            add_episode!(replay, episode)
            episode = DQExperience[] # empty episode
            obs = reset(env)
            reset_hidden_state!(graph.lstm_policy)
            push!(episode_rewards, 0.0)
            done = false
            step = 0
            rtot = 0
            # log episode reward

        end
        num_episodes = length(episode_rewards)
        avg100_reward = mean(episode_rewards[max(1, length(episode_rewards) - 101) : end])
        if t%solver.train_freq == 0
            s_batch, a_batch, r_batch, sp_batch, done_batch, trace_mask_batch = sample(replay)
            feed_dict = Dict(graph.s => s_batch,
                             graph.a => a_batch,
                             graph.sp => sp_batch,
                             graph.r => r_batch,
                             graph.done_mask => done_batch,
                             graph.trace_mask => trace_mask_batch,
                             graph.importance_weights => weights,
                             graph.hq_in.c => init_c,
                             graph.hq_in.h => init_h,
                             graph.hqp_in.c => init_c,
                             graph.hqp_in.h => init_h,
                             graph.target_hq_in.c => init_c,
                             graph.target_hq_in.h => init_h
                             )
            loss_val, td_errors_val, grad_val, _ = run(graph.sess,
                                                       [graph.loss, graph.td_errors, graph.grad_norm, graph.train_op],
                                                       feed_dict)
        end

        if t%solver.target_update_freq == 0
            run(graph.sess, graph.update_op)
        end

        if t%solver.eval_freq == 0
            # save hidden state before
            hidden_state = graph.lstm_policy.state_val
            scores_eval = eval_lstm(graph.lstm_policy,
                                     env,
                                     graph.sess,
                                     n_eval=solver.num_ep_eval,
                                     max_episode_length=solver.max_episode_length,
                                     verbose = solver.verbose)
            # reset hidden state
            graph.lstm_policy.state_val = hidden_state
        end

        if t%solver.log_freq == 0
            # log to tensorboard
            tb_avgr = logg_scalar(avg100_reward, "avg_reward")
            tb_evalr = logg_scalar(scores_eval[end], "eval_reward")
            tb_loss = logg_scalar(loss_val, "loss")
            tb_tderr = logg_scalar(mean(td_errors_val), "mean_td_error")
            tb_grad = logg_scalar(grad_val, "grad_norm")
            tb_epreward = logg_scalar(episode_rewards[end], "episode_reward")
            tb_eps = logg_scalar(eps, "epsilon")
            write(summary_writer, tb_avgr, t)
            write(summary_writer, tb_evalr, t)
            write(summary_writer, tb_loss, t)
            write(summary_writer, tb_tderr, t)
            write(summary_writer, tb_grad, t)
            write(summary_writer, tb_epreward, t)
            write(summary_writer, tb_eps, t)
            if  solver.verbose
                logg = @sprintf("%5d / %5d eps %0.3f |  avgR %1.3f | Loss %2.3f | Grad %2.3f | num_episodes %f",
                                 t, solver.max_steps, eps, avg100_reward, loss_val, grad_val, length(episode_rewards))
                println(logg)
            end
        end
        if t > solver.train_start && t%solver.save_freq == 0
            if scores_eval[end] >= saved_mean_reward
                if solver.verbose
                    println("Saving new model with eval reward ", scores_eval[end])
                end
                saver = tf.train.Saver()
                train.save(saver, graph.sess, solver.logdir*"weights.jld")
                model_saved = true
                saved_mean_reward = scores_eval[end]
            end
        end
    end
    if model_saved
        if solver.verbose
            println("Restore model with eval reward ", saved_mean_reward)
        end
    end
    return
end

function eval_lstm(policy::LSTMPolicy,
                env::AbstractEnvironment,
                sess;
                n_eval::Int64=100,
                max_episode_length::Int64=100,
                verbose::Bool=false)
    # Evaluation
    avg_r = 0
    for i=1:n_eval
        done = false
        r_tot = 0.0
        step = 0
        obs = reset(env)
        reset_hidden_state!(policy)
        # println("start at t=0 obs $obs")
        # println("Start state $(env.state)")
        while !done && step <= max_episode_length
            action = get_action!(policy, obs, sess)
            # println(action)
            obs, rew, _, done, info = step!(env, action)
            # println("state ", env.state, " action ", a)
            # println("Reward ", rew)
            # println(obs, " ", done, " ", info, " ", step)
            r_tot += rew
            step += 1
        end
        avg_r += r_tot
        # println(r_tot)
    end
    if verbose
        println("Evaluation ... Avg Reward ", avg_r/n_eval)
    end
    return  avg_r /= n_eval
end

max_stars_repo_name: parachutel/DeepCorrectionCAS


text:
using Reiwa
using Unicode
using Test

@testset "instantiate" begin

rstr = ReiwaString("\u337e\u337d\u337c\u337b\u32ff")
rstr2 = reiwa"㍾㍽㍼㍻㋿"
@test rstr == rstr2 == "\u337e\u337d\u337c\u337b\u32ff" == "㍾㍽㍼㍻㋿"

end

@testset "io" begin

rstr = reiwa"㍾㍽㍼㍻㋿"
@test sprint(print, rstr) == "㍾㍽㍼㍻㋿"
@test sprint(show, rstr) == "reiwa\"㍾㍽㍼㍻㋿\""

rstr3 = reiwa"\u337b31年度\t（\u32ff元年度）\n"
@test sprint(print, rstr3) == "㍻31年度\t（㋿元年度）\n"
@test sprint(show, rstr3) == "reiwa\"㍻31年度\\t（㋿元年度）\\n\""

end

@testset "iteration" begin

rstr = reiwa"㍾㍽㍼㍻㋿"
@test rstr[1] == '㍾' == '\u337e'
@test rstr[4] == '㍽' == '\u337d'
@test rstr[7] == '㍼' == '\u337c'
@test rstr[10] == '㍻' == '\u337b'
@test rstr[13] == '㋿' == '\u32ff'
@test collect(rstr) == Char[0x337e, 0x337d, 0x337c, 0x337b, 0x32ff]

end

@testset "isassigned" begin

@test Unicode.isassigned('\u337e')
@test Unicode.isassigned('\u337d')
@test Unicode.isassigned('\u337c')
@test Unicode.isassigned('\u337b')
@test Unicode.isassigned('\u32ff')
@test Unicode.isassigned(0x337e)
@test Unicode.isassigned(0x337d)
@test Unicode.isassigned(0x337c)
@test Unicode.isassigned(0x337b)
@test Unicode.isassigned(0x32ff)

end

@testset "normalize" begin

rstr = reiwa"\u337e\u337d\u337c\u337b\u32ff"
@test Unicode.normalize(rstr, :NFC) == "㍾㍽㍼㍻㋿"
@test Unicode.normalize(rstr, :NFD) == "㍾㍽㍼㍻㋿"
@test Unicode.normalize(rstr, :NFKC) == "明治大正昭和平成令和"
@test Unicode.normalize(rstr, :NFKD) == "明治大正昭和平成令和"
@test Unicode.normalize(rstr, compat=true) == "明治大正昭和平成令和"
@test Unicode.normalize(rstr) == "㍾㍽㍼㍻㋿"

rstr3 = reiwa"㍻31年度（\u32ff元年度）"
@test Unicode.normalize(rstr3, :NFC) == "㍻31年度（㋿元年度）"
@test Unicode.normalize(rstr3, :NFD) == "㍻31年度（㋿元年度）"
@test Unicode.normalize(rstr3, :NFKC) == "平成31年度(令和元年度)"
@test Unicode.normalize(rstr3, :NFKD) == "平成31年度(令和元年度)"
@test Unicode.normalize(rstr3, compat=true) == "平成31年度(令和元年度)"
@test Unicode.normalize(rstr3) == "㍻31年度（㋿元年度）"

end


text:
## Extend convert to handle JuMP expression types
# GenericAffExpr
function Base.convert(::Type{JuMP.GenericAffExpr{C, V}},
                      x::JuMP.GenericAffExpr{C, W}) where {C, V, W}
    if V != W
        return JuMP.GenericAffExpr{C,V}(x.constant, x.terms)
    else
        return x
    end
end

# GenericQuadExpr
function Base.convert(::Type{JuMP.GenericQuadExpr{C, V}},
                      x::JuMP.GenericQuadExpr{C, W}) where {C, V, W}
    if V != W
        return JuMP.GenericQuadExpr{C, V}(x.aff, x.terms)
    else
        return x
    end
end

# UnorderedPair
function Base.convert(::Type{JuMP.UnorderedPair{V}},
                      x::JuMP.UnorderedPair{W}) where {V, W}
    if V != W
        return JuMP.UnorderedPair{V}(x.a, x.b)
    else
        return x
    end
end

## Extend convert to handle JuMP containers
# Array -> SparseAxisArray
function Base.convert(::Type{JuMPC.SparseAxisArray}, arr::Array)
    data = Dict(Tuple(k) => arr[k] for k in CartesianIndices(arr))
    return JuMPC.SparseAxisArray(data)
end

# DenseAxisArray -> SparseAxisArray
function Base.convert(::Type{JuMPC.SparseAxisArray},
                      arr::JuMPC.DenseAxisArray)
    data = Dict(k.I => arr[k] for k in keys(arr))
    return JuMPC.SparseAxisArray(data)
end

## Extensions for ParameterBounds
# length
Base.length(bounds::ParameterBounds)::Int = length(bounds.intervals)

# equal to
function Base.:(==)(bounds1::ParameterBounds, bounds2::ParameterBounds)::Bool
    return bounds1.intervals == bounds2.intervals
end

# copy
Base.copy(bounds::ParameterBounds) = ParameterBounds(copy(bounds.intervals))

# function Base.convert(::Type{Array}, arr::JuMPC.SparseAxisArray)
#
# end
#
# function Base.convert(::Type{JuMPC.DenseAxisArray},
#                       arr::JuMPC.SparseAxisArray)
#
# end

# Extend to handle InfOptParameters correctly
function Base.:(==)(p1::InfOptParameter, p2::InfOptParameter)
    check1 = p1.set == p2.set
    check2 = isequal(p1.supports, p2.supports)
    check3 = p1.independent == p2.independent
    return (check1 && check2 && check3)
end

# Hack to make the keys function work for sparse arrays
Base.keys(d::JuMPC.SparseAxisArray) = keys(d.data)
# Hacky fix to compare SparseAxisArrays
function Base.isapprox(a::JuMPC.SparseAxisArray, b::JuMPC.SparseAxisArray)::Bool
    return all(isapprox.(a, b))
end

# Attempt to convert variable type of GenericAffExpr if possible
function _possible_convert(type::DataType,
                           aff::JuMP.GenericAffExpr{C, V}) where {C, V}
    valids = [k isa type for k in keys(aff.terms)]
    if all(valids)
        return JuMP.GenericAffExpr{C, type}(aff.constant, aff.terms)
    else
        return aff
    end
end

# Attempt to convert variable type of GenericQuadExpr if possible
function _possible_convert(type::DataType,
                           quad::JuMP.GenericQuadExpr{C, V}) where {C, V}
    valids_a = [k.a isa type for k in keys(quad.terms)]
    valids_b = [k.b isa type for k in keys(quad.terms)]
    valids_aff = [k isa type for k in keys(quad.aff.terms)]
    if all(valids_a) && all(valids_b) && all(valids_aff)
        aff = convert(JuMP.GenericAffExpr{C, type}, quad.aff)
        return JuMP.GenericQuadExpr{C, type}(aff, quad.terms)
    else
        return quad
    end
end

## Define functions to convert a JuMP array into a vector (need for @BDconstraint)
# AbstractArray
function _make_vector(arr::AbstractArray)
    return [arr[i] for i in keys(arr)]
end

# Array (no nothing)
function _make_vector(arr::Array)
    return arr
end

# Something else
function _make_vector(arr)
    return arr
end


text:
using MPI, MPIArrays

MPI.Init()
comm = MPI.COMM_WORLD
rank = MPI.Comm_rank(comm)
n_proc = MPI.Comm_size(comm)

@assert n_proc == 4

# A = MPIArray{Float64}(8, 8)
# b = MPIArray{Float64}(8)
# forlocalpart!(x -> fill!(x, rank), A)
# forlocalpart!(x -> fill!(x, rank), b)
# sync(A, b)
# B = A*b
# if rank == 0
#     x = localindices(A, 0)
#     println(typeof(x))
#     println("type of B ", typeof(B))
# end
# free(A)
# free(b)

# A = MPIArray{Float64}(8)
# forlocalpart!(x -> fill!(x, rank), A)
# sync(A)
# if rank == 0
#     x = localindices(A, 0)
#     println(typeof(x))
# end
# free(A)

# Cyclic mul!
# A = MPIArray{Float64}(comm, (2, 2), (2, 2), (8, 8))
# b = MPIArray{Float64}(comm, (4,), (2,), (8,))
# forlocalpart!(x -> fill!(x, rank), A)
# forlocalpart!(x -> fill!(x, rank), b)
# sync(A, b)
# if rank == 0
#     x = localindices(b, 0)
#     println(typeof(x))
#     println(length.(x))
#     len = map(x->length(localindices(b, x-1)), b.partitioning.ranks)
#     println(len)
#     # println("b get $(b[7])")
#     # println("b get $(b[localindices(b, 3)[1]])")
# end
# y = A*b
# if rank == 0
#     println("size $(size(A))")
#     show(stdout, "text/plain", A)
#     println()
#     println("y=A*x size $(size(y))")
#     show(stdout, "text/plain", y)
#     println()
# end
# free(A)
# free(b)

A = MPIArray{Float64}(comm, (2, 2), (2, 2), (9, 7))
# A = MPIArray{Float64}(9, 7)
forlocalpart!(x -> fill!(x, rank), A)
sync(A)
if rank == 0
    println("size $(size(A))")
    show(stdout, "text/plain", A)
    println()
end
free(A)

A = MPIArray{Float64}(comm, (2, 2), (2, 2), (5, 7))
# A = MPIArray{Float64}(9, 7)
forlocalpart!(x -> fill!(x, rank), A)
sync(A)
if rank == 0
    println("size $(size(A))")
    show(stdout, "text/plain", A)
    println()
end
free(A)

A = MPIArray{Float64}(comm, (4, ), (2, ), (16, ))
forlocalpart!(x -> fill!(x, rank), A)
sync(A)
if rank == 0
    println("size $(size(A))")
    show(stdout, "text/plain", A)
    println()
end
free(A)

# useful
A = CyclicMPIArray(Float64, 7, 5)
forlocalpart!(x -> fill!(x, rank), A)
sync(A)
if rank == 0
    println("size $(size(A))")
    show(stdout, "text/plain", A)
    println()
end
free(A)

A = CyclicMPIArray(Float64, 7, 5)
forlocalpart!(x -> fill!(x, rank), A)
sync(A)
if rank == 0
    println("size $(size(A))")
    show(stdout, "text/plain", A)
    println()
end
free(A)

MPI.Finalize()


text:
glb(x::Vnum{N,ES} where {N,ES}) = isulp(x) ? prev(x) : x
lub(x::Vnum{N,ES}) where {N,ES} = isulp(x) ? next(x) : x
upper_ulp(x::Vnum{N,ES}) where {N,ES} = isulp(x) ? x : next(x)
lower_ulp(x::Vnum{N,ES}) where {N,ES} = isulp(x) ? x : prev(x)

outer_exact(x::Vnum{N,ES}) where {N,ES} = reinterpret((@Int), x) < 0 ? lub(x) : glb(x)
inner_exact(x::Vnum{N,ES}) where {N,ES} = reinterpret((@Int), x) < 0 ? glb(x) : lub(x)

function outer_ulp(x::Vnum{N,ES}) where {N,ES}
  xint = reinterpret(@Int, x)
  if xint == 0
    throw(NaNError(:inner_exact, [x]))
  elseif xint < 0
    lower_ulp(x)
  else
    upper_ulp(x)
  end
end

function inner_ulp(x::Vnum{N,ES}) where {N,ES}
  xint = reinterpret(@Int, x)
  if xint == 0
    throw(NaNError(:inner_exact, [x]))
  elseif xint < 0
    upper_ulp(x)
  else
    lower_ulp(x)
  end
end

max_issues_repo_name: MohHizzani/SigmoidNumbers

max_forks_repo_name: MohHizzani/SigmoidNumbers


text:
using TreeArrays
using OffsetArrays
using Test

using TreeArrays: isactive, isallocated, rootnode

include("utils.jl")
include("MaskedArray.jl")
include("LeafNode.jl")
include("inner_nodes.jl")
include("TreeView.jl")
include("TreeArray.jl")
include("ContinuousView.jl")
include("iterators.jl")


text:
#2020 Klaus Stricker - Tests for ReadWriteDlm2
#License is MIT: http://julialang.org/license

# rwd2tests_4.jl

# Tests for readcsv2 and writecsv2
# ================================

@testset "4_readwritecsv2" begin
    # test comments with readcsv2
    @test isequaldlm(readcsv2(IOBuffer("#this is comment\n1,2,3\n#one more comment\n4,5,6"), comments=true), Any[1 2 3;4 5 6], Any)
    @test isequaldlm(readcsv2(IOBuffer("#this is \n#comment\n1,2,3\n#one more \n#comment\n4,5,6"), comments=true), Any[1 2 3;4 5 6], Any)

    # test readcsv2 and writecsv2 with alternative decimal
    a = Float64[1.1 1.2;2.1 2.2]
    writecsv2("test.csv", a, decimal='€')
    @test read("test.csv", String) == "1€1,1€2\n2€1,2€2\n"
    b = readcsv2("test.csv", Any, rs=(r"(\d)€(\d)", s"\1.\2"), decimal='n')
    rm("test.csv")
    @test a == b

    #  Test different types with header for readcsv2 and writecsv2
    a = Any["Nr" "Value";1 Date(2017);2 DateTime(2018);3 Dates.Time(23,54,45,123,456,78);4 1.5e10+5im;5 1500//5;6 1.5e10]
    writecsv2("test.csv", a)
    @test read("test.csv", String) ==
    """
    Nr,Value
    1,2017-01-01
    2,2018-01-01T00:00:00.0
    3,23:54:45.123456078
    4,1.5e10+5.0im
    5,300//1
    6,1.5e10
    """
    b = readcsv2("test.csv", header=true)
    rm("test.csv")
    @test a[2:end,:] == b[1]
    @test a[1:1,:] == b[2]

    # Test size for readcsv2
    @test size(readcsv2(IOBuffer("1,2,3,4"))) == (1,4)
    @test size(readcsv2(IOBuffer("1,2,3,"))) == (1,4)
    @test size(readcsv2(IOBuffer("1,2,3,4\n"))) == (1,4)
    @test size(readcsv2(IOBuffer("1,2,3,\n"))) == (1,4)
    @test size(readcsv2(IOBuffer("1,2,3,4\n1,2,3,4"))) == (2,4)
    @test size(readcsv2(IOBuffer("1,2,3,4\n1,2,3,"))) == (2,4)
    @test size(readcsv2(IOBuffer("1,2,3,4\n1,2,3"))) == (2,4)

    @test size(readcsv2(IOBuffer("1,2,3,4\r\n"))) == (1,4)
    @test size(readcsv2(IOBuffer("1,2,3,4\r\n1,2,3\r\n"))) == (2,4)
    @test size(readcsv2(IOBuffer("1,2,3,4\r\n1,2,3,4\r\n"))) == (2,4)
    @test size(readcsv2(IOBuffer("1,2,3,\"4\"\r\n1,2,3,4\r\n"))) == (2,4)

    #Time types for readcsv2 and writecsv2
    a = [Dates.Time(23,55,56,123,456,789) Dates.Time(23,55,56,123,456) Dates.Time(23,55,56,123) Dates.Time(12,45) Dates.Time(11,23,11)]
    writecsv2("test.csv", a)
    @test read("test.csv", String) == "23:55:56.123456789,23:55:56.123456,23:55:56.123,12:45:00,11:23:11\n"
    b = readcsv2("test.csv")
    rm("test.csv")
    @test b == a

    # Test readcsv2/writecsv2 with Complex - Rationals
    a = Complex[complex(-1//3,-7//5) complex(1,-1//3) complex(-1//2,3e-19)]
    writecsv2("test.csv", a)
    @test read("test.csv", String) == "-1//3-7//5*im,1//1-1//3*im,-0.5+3.0e-19im\n"
    b = readcsv2("test.csv", Complex)
    rm("test.csv")
    @test a == b
    @test typeof(b) == Array{Complex,2}
end

@testset "4_empty IO-data" begin
    # Tests for empty IO-Data
    a = ""
    writedlm2("test.csv", a)
    @test read("test.csv", String) == ""
    b = readdlm2("test.csv")
    @test typeof(b) == Array{Any,2}
    @test typeof(readdlm2("test.csv", Any)) == Array{Any,2}
    rm("test.csv")
    @test isempty(b)
    writecsv2("test.csv", a)
    @test read("test.csv", String) == ""
    b = readcsv2("test.csv")
    @test typeof(b) == Array{Any,2}
    @test typeof(readcsv2("test.csv", Float64)) == Array{Float64,2}
    rm("test.csv")
    @test isempty(b)

    a = [""]
    writedlm2("test.csv", a)
    @test read("test.csv", String) == "\n"
    b = readdlm2("test.csv")
    @test typeof(b) == Array{Any,2}
    @test typeof(readdlm2("test.csv", Any)) == Array{Any,2}
    rm("test.csv")
    @test isempty(b)
    writecsv2("test.csv", a)
    @test read("test.csv", String) == "\n"
    b = readcsv2("test.csv")
    @test typeof(b) == Array{Any,2}
    @test typeof(readcsv2("test.csv", Float64)) == Array{Float64,2}
    rm("test.csv")
    @test isempty(b)

    a = missing
    writedlm2("test.csv", a)
    @test read("test.csv", String) == "na\n"
    b = readcsv2("test.csv")
    rm("test.csv")
    @test typeof(b) == Array{Any,2}
    @test isequal(a, b[1])
    writecsv2("test.csv", a)
    @test read("test.csv", String) == "na\n"
    b = readcsv2("test.csv")
    rm("test.csv")
    @test typeof(b) == Array{Any,2}
    @test isequal(a, b[1])

    a = [nothing]
    a = reshape(a, 1,1)
    writedlm2("test.csv", a)
    @test read("test.csv", String) == "nothing\n"
    b = readdlm2("test.csv")
    rm("test.csv")
    @test typeof(b) == Array{Any,2}
    @test isequal(a, b)
    writecsv2("test.csv", a)
    @test read("test.csv", String) == "nothing\n"
    b = readdlm2("test.csv")
    rm("test.csv")
    @test typeof(b) == Array{Any,2}
    @test isequal(a, b)

    a = [1.2 NaN "" nothing missing]
    writedlm2("test.csv", a)
    @test read("test.csv", String) == "1,2;NaN;;nothing;na\n"
    b = readdlm2("test.csv")
    rm("test.csv")
    @test isequal(a, b)

    a = [1.2 NaN "" nothing missing]
    writecsv2("test.csv", a)
    @test read("test.csv", String) == "1.2,NaN,,nothing,na\n"
    b = readcsv2("test.csv")
    rm("test.csv")
    @test isequal(a, b)

end

@testset "4_  random data" begin
    # Test write and read of random Any array with 10000 rows
    n = 1000
    A = Array{Any}(undef, n, 9)
    for i = 1:n
        A[i,:] = Any[randn() rand(Int) rand(Bool) rand(Date(1980,1,1):Day(1):Date(2017,12,31)) rand(Dates.Time(0,0,0,0,0,0):Dates.Nanosecond(1):Dates.Time(23,59,59,999,999,999)) rand(DateTime(1980,1,1,0,0,0,0):Dates.Millisecond(1):DateTime(2017,12,31,23,59,59,999)) randstring(24) complex(randn(), randn()) (rand(Int)//rand(Int))]
    end

    writedlm2("test.csv", A)
    B = readdlm2("test.csv")
    rm("test.csv")
    @test isequaldlm(A, B, Any)

    writecsv2("test.csv", A)
    B = readcsv2("test.csv")
    rm("test.csv")
    @test isequaldlm(A, B, Any)
end

@testset "4_abstract type" begin

    #Abstract Time Types
    a = TimeType[Date(2017, 1, 1) DateTime(2017, 2, 15, 23, 0, 0)]
    writedlm2("test.csv", a, dfs="mm/yyyy", dtfs="dd.mm.yyyy/HH.h")
    b = readdlm2("test.csv", TimeType, dfs="mm/yyyy", dtfs="dd.mm.yyyy/HH.h")
    @test isequaldlm(a, b, TimeType)
    b = readdlm2("test.csv", TimeType, rs=(), dfs="mm/yyyy", dtfs="dd.mm.yyyy/HH.h")
    rm("test.csv")
    @test isequaldlm(a, b, TimeType)

    # Abstract Number
    a = Number[1 1.1 1//3 complex(-1,-2) complex(1.2,-2) complex(-1e9,3e-19)]
    writedlm2("test.csv", a)
    @test read("test.csv", String) == "1;1,1;1//3;-1-2im;1,2-2,0im;-1,0e9+3,0e-19im\n"
    b = readdlm2("test.csv", Number)
    rm("test.csv")
    @test a == b
    @test typeof(b) == Array{Number,2}

    # Abstract real
    a = readdlm2(IOBuffer("1,1;2\n3;4,5\n5;1//3\n"), Real)
    b = Real[1.1 2; 3 4.5; 5 1//3]
    @test a == b
    @test typeof(a) == typeof(b)

end

@testset "4_    1 0 1 0 0" begin

    # Different Types for "1 0 1 0 0"
    a = readdlm2(IOBuffer("1;0;1;0;0\n"))
    b = Any[1 0 1 0 0]
    @test a == b
    @test typeof(a) == typeof(b)

    a = readdlm2(IOBuffer("1;0;1;0;0\n"), Float64)
    b = [1.0 0.0 1.0 0.0 0.0]
    @test a == b
    @test typeof(a) == typeof(b)

    a = readdlm2(IOBuffer("1;0;1;0;0\n"), Float16)
    b = Float16[1.0 0.0 1.0 0.0 0.0]
    @test a == b
    @test typeof(a) == typeof(b)

    a = readdlm2(IOBuffer("1;0;1;0;0\n"), Int)
    b = [1 0 1 0 0]
    @test a == b
    @test typeof(a) == typeof(b)

    a = readdlm2(IOBuffer("1;0;1;0;0\n"), Char)
    b = ['1' '0' '1' '0' '0']
    @test a == b
    @test typeof(a) == typeof(b)

    a = readdlm2(IOBuffer("1;0;1;0;0\n"), AbstractString)
    b = AbstractString["1" "0" "1" "0" "0"]
    @test a == b
    @test typeof(a) == typeof(b)

    a = readdlm2(IOBuffer("1;0;1;0;0\n"), String)
    b = ["1" "0" "1" "0" "0"]
    @test a == b
    @test typeof(a) == typeof(b)

    a = readdlm2(IOBuffer("1;0;1;0;0\n"), Bool)
    b = [true false true false false]
    @test a == b
    @test typeof(a) == typeof(b)

    a = readdlm2(IOBuffer("1;0;1;0;0\n"), Complex{Int})
    b = [1+0im 0+0im 1+0im 0+0im 0+0im]
    @test a == b
    @test typeof(a) == typeof(b)

    a = readdlm2(IOBuffer("1;0;1;0;0\n"), Complex{Float64})
    b = [1.0+0.0im 0.0+0.0im 1.0+0.0im 0.0+0.0im 0.0+0.0im]
    @test a == b
    @test typeof(a) == typeof(b)

    a = readdlm2(IOBuffer("1;0;1;0;0\n"), Complex{Rational{Int}})
    b = [1//1+0//1*im 0//1+0//1*im 1//1+0//1*im 0//1+0//1*im 0//1+0//1*im]
    @test a == b
    @test typeof(a) == typeof(b)

    a = readdlm2(IOBuffer("1;0;1;0;0\n"), Rational)
    b = Rational[1//1 0//1 1//1 0//1 0//1]
    @test a == b
    @test typeof(a) == typeof(b)

end


text:
__precompile__()
module SuperStiff

include("PeriodizeSC.jl")
include("Stiffness.jl")

end # module

max_stars_repo_name: oliviersimard/SuperStiff

max_issues_repo_name: oliviersimard/SuperStiff

max_forks_repo_name: oliviersimard/SuperStiff


text:
module TestPrograms
    using Base.Test
    using Brainfuck

    text = readall(Pkg.dir("Brainfuck", "test", "programs", "hello2.bfk"))
    ops = Brainfuck.lexandparse(text)
    io = IOBuffer()
    Brainfuck.interpret(ops, io_out = io)
    @test takebuf_string(io) == "Hello world!\n"

    text = readall(Pkg.dir("Brainfuck", "test", "programs", "99bottles.bfk"))
    io = IOBuffer()
    Brainfuck.brainfuck(text, io_out = io)
end


text:
module MixtureModels
	using NumericExtensions
	using MLBase
	using Distributions

	export

	# types
	Mixture, ncomponents,

	# utils
	qmatrix, qmatrix!,

	# fmm
	FiniteMixtureEM, FiniteMixtureEMResults, fmm_em,
	fit_fmm!, fit_fmm

	include("types.jl")
	include("utils.jl")
	include("fmm_em.jl")
end


text:
# ---
# title: 207. Course Schedule
# id: problem207
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Depth-first Search, Breadth-first Search, Graph, Topological Sort
# link: <https://leetcode.com/problems/course-schedule/description/>
# hidden: true
# ---
#
# There are a total of `numCourses` courses you have to take, labeled from `0`
# to `numCourses-1`.
#
# Some courses may have prerequisites, for example to take course 0 you have to
# first take course 1, which is expressed as a pair: `[0,1]`
#
# Given the total number of courses and a list of prerequisite **pairs** , is it
# possible for you to finish all courses?
#
#
#
# **Example 1:**
#
#
#
#     Input: numCourses = 2, prerequisites = [[1,0]]
#     Output: true
#     Explanation:  There are a total of 2 courses to take.
#                  To take course 1 you should have finished course 0. So it is possible.
#
#
# **Example 2:**
#
#
#
#     Input: numCourses = 2, prerequisites = [[1,0],[0,1]]
#     Output: false
#     Explanation:  There are a total of 2 courses to take.
#                  To take course 1 you should have finished course 0, and to take course 0 you should
#                  also have finished course 1. So it is impossible.
#
#
#
#
# **Constraints:**
#
#   * The input prerequisites is a graph represented by **a list of edges** , not adjacency matrices. Read more about [how a graph is represented](https://www.khanacademy.org/computing/computer-science/algorithms/graph-representation/a/representing-graphs).
#   * You may assume that there are no duplicate edges in the input prerequisites.
#   * `1 <= numCourses <= 10^5`
#
#
## @lc code=start
using LeetCode

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


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Atomsk_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Atomsk")
JLLWrappers.@generate_main_file("Atomsk", UUID("eddcc1f7-2d42-55ad-abc0-463b4ff6875a"))
end  # module Atomsk_jll


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

name = "Atomsk"
version = v"0.11.2"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/pierrehirel/atomsk.git", "3333858281e0ebd6279825b83ab871cf4d050a8d")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/atomsk/src/
mkdir -p ${bindir}
# The makefile doesn't handle parallel builds
make atomsk BIN="atomsk${exeext}"
make install INSTPATH=${prefix} BIN="atomsk${exeext}"
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)
platforms = expand_gfortran_versions(platforms)

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

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae")),
    Dependency(PackageSpec(name="LAPACK_jll", uuid="51474c39-65e3-53ba-86ba-03b1b862ec14")),
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
"""
Compute Rayleigh scattering cross-sections.

Most of the expressions are obtained from
Bodhaine et al. (1999) J. Atmosph. Ocean. Tech., 16 1854.

If invoked as julia Rayleigh.jl, computes the cross sections for 337 and 777 nm.
From the julia REPL you can compute cross sections for any wavelength as e.g.
```julia-repl
julia> include("rayleigh.jl")
julia> using Main.Ray
julia> rayleigh(337 * co.nano)
3.4371080894577986e-30
```
Alejandro Luque, IAA-CSIC, 2019
"""
module Ray

export AIR_STP, rayleigh, co

include("constants.jl")
const co = constants

# Depolarization function for different gases.  λ here is μm.
const F = Dict("N2" => λ -> 1.034 + 3.17e-4 / λ^2,
               "O2" => λ -> 1.096 + 1.385e-3 / λ^2 + 1.448e-4 / λ^4,
               "Ar" => λ -> 1.00,
               "CO2" => λ -> 1.15)

# These are densities from ICAO, extracted from Wikipedia.  Note that
# the expressions in Bodhaine assume 300 ppmv CO2 but we consider this
# close enough for our purposes.
const AIR_STP = Dict("N2" => 780840,
                     "O2" => 209476,
                     "Ar" => 9340,
                     "CO2" => 314)

# The density of air in constants.jl is defined at (273.15 K, 101.325 kPa)
# but the functions in Bodhaine et al. 1999 are valid for 288.15 K so we have to
# rescale
nair = co.nair * 273.15 / 288.15

"""
Compute the refractive index of air.

Takes the wavelength λ in m and returns (n - 1).
Uses the formula from Edlén (1953) referenced from Bodhaine et al. 1999.

"""
function edlen(λ)
    a, b, c = 6432.8, 2949810, 25540
    p, q = 146, 41

    λ_μm = λ / co.micro
    λ2 = 1 / λ_μm^2

    1.0 + 1e-8 * (a + b / (p - λ2) + c / (q - λ2))
end

"""
Compute the refractive index of air.

Takes the wavelength λ in m and returns (n - 1).
Uses the formula from Peck and Reeder (1953) referenced from Bodhaine et al. 1999.

"""
function peck(λ)
    a, b, c = 8060.51, 2480990, 17455.7
    p, q = 132.274, 39.32957

    λ_μm = λ / co.micro
    λ2 = 1 / λ_μm^2

    1.0 + 1e-8 * (a + b / (p - λ2) + c / (q - λ2))
end

"""
Compute the Rayleigh scattering cross-section.
"""
function σ(λ; composition=AIR_STP)
    n = peck(λ)
    P = sum(ppmv * F[gas](λ / co.micro) for (gas, ppmv) in composition) / 1e6

    N = (n^2 - 1)^2 / (n^2 + 2)^2
    24π^3 * N * P / (λ^4 * nair^2)
end

function main()
    println("# Rayleigh scattering cross-section for 337 nm")
    σ = rayleigh(337 * co.nano)
    println("σ: $σ")

    println()
    println("# Rayleigh scattering cross-section for 777 nm")
    σ = rayleigh(777 * co.nano)
    println("σ: $σ")
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end

# Depolarization functions
end


text:
using StatsBase: sample
using UnicodePlots: spy
using Printf

abstract type SEM end

struct GaussianSEM <: SEM
    p       ::Int64
    B       ::Matrix{Float64}
    err_var ::Vector{Float64}
    function GaussianSEM(B, err_var)
        @assert size(B, 1) == size(B, 2) == length(err_var)
        @assert all(err_var .> 0)
        p = length(err_var)
        new(p, B, err_var)
    end
end

function Base.show(io::IO, sem::GaussianSEM)
    print(io, "Gaussian SEM with $(sem.p) variables:\n")
    print(io, "B = \n")
    print(io, spy(sem.B))
    print(io, "σ² = $(sem.err_var)")
end

"""
    simulate(sem, [n])
    simulate(sem, [do_variables, do_values], [n])

Simulate from a Gaussian SEM `sem`. `n` is the sample size.
do-interventions can be performed by specifying vectors of `do_variables` and `do_values`.
"""
function simulate(sem::GaussianSEM)
    p = sem.p
    ϵ = randn(p) .* sqrt.(sem.err_var)
    return (I - sem.B) \ ϵ
end

function simulate(sem::GaussianSEM, do_variables::Vector{Int64}, do_values::Vector{Float64})
    @assert length(do_variables) == length(do_values)
    p = sem.p
    ϵ = randn(p) .* sqrt.(sem.err_var)
    ϵ[do_variables] .= do_values
    B = copy(sem.B)
    B[do_variables, :] .= 0
    return (I - B) \ ϵ
end

function simulate(sem::GaussianSEM, n::Int64)
    return vcat(map(i -> simulate(sem), 1:n)'...)
end

function simulate(sem::GaussianSEM, do_variables::Vector{Int64}, do_values::Vector{Float64}, n::Int64)
    return vcat(map(i -> simulate(sem, do_variables, do_values), 1:n)'...)
end

function causes(sem::SEM, i::Int64)
    @assert 1 <= i <= sem.p
    return (1:sem.p)[sem.B[i, :].!=0]
end

function cov(sem::GaussianSEM)
    S = inv(I - sem.B)
    return S * diagm(sem.err_var) * S'
end

"""
    random_gaussian_SEM(p, k; [lb=-2, ub=2, var_min=0.5, var_max=2])

Generate a random-graph acyclic SEM with `p` variables and `k` average degree, and random coefficients.
* `lb`, `ub`: coeff  ~ unif[`lb`, `ub`] with random sign
* `var_min`, `var_max`: var of error ~ unif[`var.min`, `var.max`]
"""
function random_gaussian_SEM(p::Int64, k::Int64; lb=-2, ub=2, var_min=0.5, var_max=2)
    B = zeros(p, p)
    B[rand(p, p) .< 2k / (p-1)] .= 1
    B[UpperTriangular(B).!=0] .= 0
    m = sum(B.==1)
    B[B.==1] .= (rand(m) * (ub - lb) .+ lb) .* sign.(randn(m))
    err_var = rand(p) * (var_max - var_min) .+ var_min
    _order = sample(1:p, p, replace=false)
    B = B[_order, _order]
    return GaussianSEM(B, err_var)
end

"""
    random_noise_intervened_SEM(sem::GaussianSEM, [p_intervened=2, noise_multiplier_min=0.5, noise_multiplier_max=2., avoid=[],
                                                   prob_coeff_unchanged=2/3, lb=-2, ub=2])

Produce a new SEM based on original SEM by changing coefficients and noise variances.
* `p_intervened`: randomly choose `p_intervened` variables to intervene; will avoid those specified in `avoid`
* [`noise_multiplier_min`, `noise_multiplier_max`]: interval that noise multiplier is uniformly sampled from
* `prob.coeff.unchanged`: probability that coefficient is not changed
* `[lb, ub]`: if to change, coefficient is drawn uniformly from this interval with random sign

Return: `sem_new`, `intervened_variables`
"""
function random_noise_intervened_SEM(sem::GaussianSEM;
                                     p_intervened=2, noise_multiplier_min=0.5, noise_multiplier_max=2., avoid=[],
                                     prob_coeff_unchanged=2/3, lb=-2, ub=2)
    B = copy(sem.B)
    p = sem.p
    err_var = copy(sem.err_var)
    vars = sample(setdiff(collect(1:p), avoid), p_intervened, replace=false)
    for i in vars
        noise_multiplier = rand() * (noise_multiplier_max - noise_multiplier_min) + noise_multiplier_min
        err_var[i] = err_var[i] * noise_multiplier
        if rand() > prob_coeff_unchanged
            _J = (1:p)[B[i, :].!=0]
            B[i, _J] .= rand(length(_J)) * (ub - lb) .+ lb
        end
    end
    return GaussianSEM(B, err_var), vars
end

max_issues_repo_name: richardkwo/InvariantCausal

max_forks_repo_name: richardkwo/InvariantCausal


text:
abstract type PseudoAbsenceGenerator end

struct WithinRadius <: PseudoAbsenceGenerator
end

struct SurfaceRangeEnvelope <: PseudoAbsenceGenerator
end

struct RandomSelection <: PseudoAbsenceGenerator
end


text:
using Documenter, SPPL

makedocs(modules=[SPPL],
         sitename="SPPL",
         authors="McCoy R. Becker and other contributors",
         pages=["API Documentation" => "index.md"])

deploydocs(repo = "github.com/probcomp/SPPL.jl.git")


text:
"""
Return the Inorder traversal of a complete binary tree with given `height`.
(Root node will be `0` and placed in the middle of the array
"""
function binary_inorder(height::Integer)::Vector{Int}
    n = (1 << Int(height)) - 1
    out = Vector{Int}(undef, n)

    s = 2
    while n >= 1
        t = s >> 1
        b, a = n-1, n >> 1
        for (i, j) in enumerate(a:b)
            out[i*s - t] = j
        end
        n >>= 1
        s <<= 1
    end

    return out
end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
"""
    frft(signal, α)

By entering the signal and the order, we can obtain the Fractional Fourier transform value.

# Example

```julia-repl
julia> frft([1,2,3], 0.5)
0.10046461358821344 - 0.25940948097727745im
3.0746001042331557 + 0.19934938154063286im
1.3156643288728376 - 0.9364535656119107im
```
"""
function frft(signal, α)::Vector{ComplexF64}

    #Ensure the input signal is proper
    isa(signal, AbstractArray) ? nothing : error("Please input a proper signal")

    N = length(signal)
    M = Int64((N-1)/2)

    if rem(N, 2)==0
        error("Signal must be odd")
    end

    ##For α special cases
    #If the order is an integer, the fourier transform would be its n-th order Fourier transform.
    α = mod(α, 4)
    if α == 0
        return signal
    elseif α == 1
        y = fft([signal[M+1:N]; signal[1:M]])/sqrt(N)
        y = [y[M+2 : N]; y[1:M+1]]
        return y
    elseif α == 2
        y = reverse(signal, dims=1)
        return y
    elseif α == 3
        y = fft([signal[M+1:N]; signal[1:M]])*sqrt(N)
        y = [y[M+2 : N]; y[1:M+1]]
        return y
    end

    if α > 2
        signal = reverse(signal, dims=1)
        α = α-2
    end
    if 1.5 < α < 2
        y = fft([signal[M+1:N]; signal[1:M]]) ./sqrt(N)
        y = [y[M+2 : N]; y[1:M+1]]
        α = α-1
    end
    if α < 0.5
        y = fft([signal[M+1:N]; signal[1:M]]) .*sqrt(N)
        y = [y[M+2 : N]; y[1:M+1]]
        α = α+1
    end

    # We need to distinguish **order** and **angle**
    ϕ = α*π/2

    # Set the sampling rate as 3.
    signal = sinc_interp(signal, 3)
    signal = [zeros(ComplexF64, 2*M); signal; zeros(ComplexF64, 2*M)]

    x2 = freq_shear(signal, 2*pi/N*(cot(ϕ)-csc(ϕ))/3^2)
    x3 = time_shear(x2, 2*pi/N*csc(ϕ)/3^2)
    y = freq_shear(x3, 2*pi/N*(cot(ϕ)-csc(ϕ))/3^2)

    y = y[2*M+1:end-2*M]
    y = y[1:3:end]
    y = @. exp(-im*(pi*sign(sin(ϕ))/4-ϕ/2))*y

    return y
end

function freq_shear(x, c)
    x = x[:]
    N = length(x)
    if rem(N, 2) == 0
        error("Signal must be odd")
    end

    M = (N-1)/2
    N = collect(-M:M)

    y = @. x*exp(im*c/2*N^2)
    return y
end

function time_shear(x, c)
    x = x[:]
    N = length(x)

    if rem(N, 2) == 0
        error("Signal must be odd")
    end
    M = Int64((N-1)/2)

    interp = ceil(Int, 2*abs(c)/(2*pi/N))
    xx = sinc_interp(x, interp) ./interp
    n = collect(-2*M:1/interp:2*M)
    z = @. exp(im*c/2*n^2)
    y = conv(xx, z)
    center = (length(y)+1)/2
    y = y[Int64(center-interp*M):Int64(interp):Int64(center+interp*M)]
    y = @. y*sqrt(c/2/pi)
end

"""
    sinc_interp(signal, rate)

```Sinc interpolation``` of **signal** at rate **rate**.

For more details, please refer to [Whittaker-Shannon interpolation](https://en.wikipedia.org/wiki/Whittaker%E2%80%93Shannon_interpolation_formula).
"""
function sinc_interp(x, rate)
    x = x[:]
    N = length(x)
    M = rate*N - rate + 1

    y = zeros(ComplexF64, M)
    y[1:rate:M] = x

    h =  sinc.(collect(-(N-1-1/rate):1/rate:(N-1-1/rate)))
    out = conv(y, h)
    out = out[(rate*N-rate):(end-rate*N+rate+1)]
end


text:
# This file is a part of LegendDataTypes.jl, licensed under the MIT License (MIT).

function _match_or_throw(expr::Regex, s::AbstractString)
    m = match(expr, s)
    m != nothing || throw(ErrorException("Invalid input, line didnt match $expr"))
    m
end


text:
"""
    StaticUnitRange{T,B,E,L}

Supertype for static ranges with a step size of oneunit(T) with elements of type T.
UnitRange and other types are subtypes of this

# Examples
```jldoctest
```
"""
abstract type StaticUnitRange{T,B,E,L} <: StaticStartRange{T,B,E,L} end

# TODO: these should be static
@inline step(::StaticUnitRange{T}) where T = one(T)::T
@inline step(::Type{<:StaticUnitRange{T}}) where T = one(T)::T

# if length is dynamic
length(r::StaticUnitRange{T,B,E,Dynamic}) where {T,B,E,Dynamic} = Integer(last(r) - first(r) + step(r))
function length(r::StaticUnitRange{T,B,E,Dynamic}) where {T<:Union{Int,Int64,Int128},B,E}
    Base.@_inline_meta
    (last(r) - first(r)) + one(T)
end

length(r::StaticUnitRange{T,B,E,Dynamic}) where {T<:Union{UInt,UInt64,UInt128},B,E} =
    last(r) < first(r) ? zero(T) : (last(r) - first(r) + one(T))

"UniteSRange - Static parametric type  variant of StaticUnitRange"
struct UnitSRange{T,B,E,L} <: StaticUnitRange{T,B,E,L} end

"UniteMRange - Mutable variant of StaticUnitRange"
mutable struct UnitMRange{T,B,E} <: StaticUnitRange{T,B,E,Dynamic}
    start::B
    stop::E
end

UnitMRange{Bool}(start::Bool, stop::SBool{E}) where E =
   UnitMRange{Bool,Bool,SBool{E}}(start, stop)

UnitMRange{T}(start::B, stop::T) where {T<:BaseInteger,B<:SInteger} =
   UnitMRange{T,B,T}(start, ifelse(stop >= start, stop, ofeltype(T, start - oneunit(stop - start))))

UnitMRange{T}(start::T, stop::E) where {T<:BaseInteger,E<:SInteger} =
   UnitMRange{T,T,E}(start, ifelse(stop >= start, stop, ofeltype(T,start-oneunit(stop - start))))

UnitMRange{T}(start::T, stop::T) where T =
    UnitMRange{T,T,T}(start, ifelse(stop >= start, stop, ofeltype(T,start-oneunit(stop - start))))

function UnitMRange(start::B, stop::E) where {B,E}
    T = Base.promote_eltype(B, E)
    UnitMRange{T}(ofeltype(T, start), ofeltype(T, stop))
end

function UnitMRange{T}(start::B, stop::T) where {T,B<:SReal}
    if stop >= start
        UnitMRange{T,B,T}(start, ofeltype(T, start + floor(stop - start)))
    else
        UnitMRange{T,B,T}(start, ofeltype(T, oneunit(stop - start)))
    end
end

function UnitMRange{T}(start::T, stop::SReal{E}) where {T,E}
    if stop >= start
        last = ofeltype(T, start + floor(E::T - start))
    else
        last = ofeltype(T, start - floor(E::T - start))
    end
    UnitMRange{T,T,typeof(last)}(start, last)
end

UnitSRange{Bool}(start::SBool{B}, stop::SBool{E}) where {B,E} =
    UnitSRange{Bool,SBool{B},SBool{E}}(start, stop)

function UnitSRange{T}(start::B, stop::E) where {B<:SReal,E<:SReal,T<:BaseInteger}
    if stop >= start
        return UnitSRange{T,B}(start, stop)
    else
        return UnitSRange{T,B}(start, ofeltype(T, start-one(start)))
    end
end

function UnitSRange{T}(start::B, stop::E) where {B<:SReal,E<:SReal,T}
    if stop >= start
        UnitSRange{T,B}(start, ofeltype(T, start+floor(stop-start)))
    else
        UnitSRange{T,B}(start, ofeltype(T, start-one(start - stop)))
    end
end

UnitSRange{T,B}(start::B, stop::E) where {T<:Union{Int,Int64,Int128},B,E} =
    UnitSRange{T,B,E}(start, stop, ofeltype(Int, (stop - start) + one(B)))

UnitSRange{T,B}(start::B, stop::E) where {T<:Union{UInt,UInt64,UInt128},B,E} =
    UnitSRange{T,B,E}(start, stop, stop < stop ? zero(B) : (stop - start) + one(start))

UnitSRange{T,B}(start::B, stop::E) where {T<:Real,B,E} =
    UnitSRange{T,B,E}(start, stop, integer(start - stop + oneunit(B)))

UnitSRange{T,B,E}(start::B, stop::E, len::L) where {T,B,E,L} = UnitSRange{T,B,E,L}()

isstatic(::UnitSRange) = true
isstatic(::Type{<:UnitSRange}) = true
isstatic(::UnitMRange) = true
isstatic(::Type{<:UnitMRange}) = false

#promote_rule(a::Type{UnitRange{T1}}, b::Type{UnitRange{T2}}) where {T1,T2} = el_same(promote_type(T1,T2), a, b)
Base.UnitRange(r::StaticUnitRange) = UnitRange(values(first(r)), values(last(r)))
Base.UnitRange{T}(r::StaticUnitRange) where {T<:Real} = UnitRange{T}(values(first(r)), values(last(r)))

UnitSRange(r::Union{AbstractUnitRange,UnitMRange}) = UnitSRange(tostatic(first(r)), tostatic(last(r)))
UnitMRange(r::Union{AbstractUnitRange,UnitSRange}) = UnitMRange(first(r), last(r))
UnitMRange(r::UnitMRange) = r
UnitSRange(r::UnitSRange) = r

UnitSRange{T}(r::Union{AbstractUnitRange,UnitMRange}) where T = UnitSRange(tostatic(first(r)), tostatic(last(r)))
UnitMRange{T}(r::Union{AbstractUnitRange,UnitSRange}) where T = UnitMRange(first(r), last(r))
UnitMRange{T}(r::UnitMRange{T}) where T = r
UnitSRange{T}(r::UnitSRange{T}) where T = r
UnitMRange{T}(r::UnitMRange) where T = UnitSRange(tostatic(first(r)), tostatic(last(r)))
UnitSRange{T}(r::UnitSRange) where T = UnitMRange(first(r), last(r))

#promote_rule(a::Type{UnitRange{T1}}, ::Type{UR}) where {T1,UR<:AbstractUnitRange} =
#    promote_rule(a, UnitRange{eltype(UR)})

StaticUnitRange{T}(start::SVal, stop::SVal) where T = UnitSRange{T}(start, stop)
StaticUnitRange{T}(start::SVal, stop::Any) where T = UnitMRange{T}(start, stop)
StaticUnitRange{T}(start::Any, stop::SVal) where T = UnitMRange{T}(start, stop)

const OverflowSafe = Union{Bool,Int8,Int16,Int32,Int64,Int128,
                           UInt8,UInt16,UInt32,UInt64,UInt128}

_in_unit_range(v::StaticUnitRange, val, i::Integer) =
    i > 0 && val <= last(v) && val >= first(v)

function getindex(v::UnitMRange{T}, i::Integer) where {T<:OverflowSafe}
    Base.@_inline_meta
    val = first(v) + (i - one(i))
    @boundscheck _in_unit_range(v, val, i) || Base.throw_boundserror(v, i)
    val % T
end

function getindex(r::UnitSRange{T,B,E,L}, i::Integer) where {T,B,E,L}
    Base.@_inline_meta
    @boundscheck if i < 1 || i > values(L)
        throw(BoundsError(r, i))
    end
    unsafe_getindex(r, T(i))
end

@pure unsafe_getindex(::UnitSRange{T,B,E,L}, i::T) where {T,B,E,L} = values(B) + i - one(T)

function getindex(r::StaticUnitRange, s::Union{AbstractUnitRange{<:Integer},StaticUnitRange{<:Integer}})
    Base.@_inline_meta
    @boundscheck checkbounds(r, s)
    f = first(r)
    st = ofeltype(f, f + first(s)-SOne)
    if st isa SReal
        return range(st, length=SInteger(length(s)))
    else
        return range(st, length=length(s))
    end
end

@pure Base.checkindex(::Type{Bool}, r::UnitSRange{T,B,E,L}, i::Real) where {T,B,E,L} =
    SOne <= i && i <= L()::L

@inline Base.checkindex(::Type{Bool}, r::UnitSRange{T,B,E,L}, i::AbstractRange) where {T,B,E,L} =
    checkindex(Bool, r, first(i)) && checkindex(Bool, r, last(i))

function Base.checkindex(::Type{Bool}, inds::StaticUnitRange{T,B,E,L}, i::Real) where {T,B,E,L}
    Base.@_inline_meta
    _checkfirst(first(inds)::B, i) & _checklast(last(inds)::E, i)
end

_checkfirst(b::BaseReal, i::Real) = b <= i
@pure _checkfirst(::B, i::Real) where B<:SReal = values(B) <= i

_checklast(l::BaseReal, i::Real) = i <= l
@pure _checklast(::L, i::Real) where L<:SReal = i <= values(L)

@pure function Base.iterate(r::UnitSRange{T,B,E,L}) where {T,B,E,L}
    if values(L) < 1
        return nothing
    else
        return (values(B)::T, values(B)::T)::Tuple{T,T}
    end
end

@inline function Base.iterate(r::UnitSRange{T,B,E,L}, i::T) where {T,B,E,L}
    if i < values(E)::T
        return (i + one(T), i + one(T))::Tuple{T,T}
    else
        return nothing
    end
end

showrange(io::IO, r::StaticUnitRange) = print(io, "$(first(r)):$(last(r))")