text:
# # Hydraulic functions

# This tutorial shows how to specify the hydraulic functions
# used in Richard's equation. In particular,
# we show how to choose the formalism for matric potential and hydraulic
# conductivity, and how to make the hydraulic conductivity account for
# the presence of ice as well as the temperature dependence of the
# viscosity of liquid water.

# ClimateMachine's Land model allows the user to pick between two hydraulics models,
# that of van Genuchten [vanGenuchten1980](@cite) or that of Brooks and Corey, see [BrooksCorey1964](@cite) or [Corey1977](@cite). The
# same model is consistently used for the matric potential
# and hydraulic conductivity.

# # Preliminary setup

# External modules
using Plots

# ClimateMachine modules
using ClimateMachine
using ClimateMachine.Land
using ClimateMachine.Land.SoilWaterParameterizations

const FT = Float32;
# # Specifying a  `hydraulics` model
# The van Genuchten model requires two free parameters, `α` and `n`.
# A third parameter, `m`, is computed from `n`. Of these, only `α` carries
# units, of inverse meters. The Brooks and Corey model also uses
# two free parameters, `ψ_b`, the magnitude of the matric potential at saturation,
#  and a constant `M`. `ψ_b` carries units of meters. The hydraulic conductivity
# requires an additional parameter, `Ksat` (m/s), which is the hydraulic conductivity
# in saturated soil. This parameter is the same between the two models for a given
# soil type, and is not stored in the `hydraulics` model, but rather as part of the
# [`SoilParamFunctions`](@ref ClimateMachine.Land.SoilParamFunctions).

# Below we show how to create two concrete examples of these hydraulics models,
# for sandy loam ([Bonan19a](@cite)). Note that the parameters chosen are a function of soil type.
vg_α = FT(7.5) # m^-1
vg_n = FT(1.89)
hydraulics = vanGenuchten{FT}(α = vg_α, n = vg_n)

ψ_sat = 0.09 # m
Mval = 0.228
hydraulics_bc = BrooksCorey{FT}(ψb = ψ_sat, m = Mval);
# # Matric Potential
# The matric potential `ψ` represents how much water clings to soil. Drier soil
# holds onto water more tightly, impeding the flow of water. As soil becomes
# wetter, the matric potential decreases in magnitude, enabling water flow.

# The van Genuchten expression for matric potential is
# ``
# ψ = -\frac{1}{α} S_l^{-1/(nm)}\times (1-S_l^{1/m})^{1/n},
# ``

# and the Brooks and Corey expression is
# ``
# ψ = -ψ_b S_l^{-M}.
# ``

# Here `S_l` is the effective saturation of liquid water, `θ_l/ν`, where `ν` is
# porosity of the soil. We neglect the residual pore space in the CliMA model.

# In the CliMA code, we use [multiple dispatch](https://en.wikipedia.org/wiki/Multiple_dispatch).
# With multiple dispatch, a function can have many
# ways of executing (called methods), depending on the *type* of the
# variables passed in. A simple example of multiple dispatch is the division operation.
# Integer division takes two numbers as input, and returns an integer - ignoring the decimal.
# Float division takes two numbers as input, and returns a floating point number, including the decimal.
# In Julia, we might write these as:

# ```julia
# function division(a::Int, b::Int)
#      return floor(Int, a/b)
# end
# ```
# ```julia
# function division(a::Float64, b::Float64)
#      return a/b
# end
# ```

# We can see that `division` is now a function with two methods.

# ```julia
# julia> division
# division (generic function with 2 methods)
# ```

# Now, using the same function signature, we can carry out integer
# division or floating point division, depending on the types of the
# arguments:

# ```julia
# julia> division(1,2)
# 0
#
# julia> division(1.0,2.0)
# 0.5
# ```

# Here is a more pertinent example:
# `hydraulics` is of type `vanGenuchten{Float32}` based on our choice of `FT`:

# ```julia
# julia> typeof(hydraulics)
# vanGenuchten{Float32}
# ```

# but meanwhile,

# ```julia
# julia> typeof(hydraulics_bc)
# BrooksCorey{Float32}
# ```

# The function `matric_potential` will execute different methods
# depending on if we pass a hydraulics model of type `vanGenuchten` or
# `BrooksCorey`. In both cases, it will return the correct value
# for `ψ`.

# Let's plot the matric potential as a function of the effective saturation `S_l = θ_l/ν`,
# which can range from zero to one.
S_l = FT.(0.01:0.01:0.99)
ψ = matric_potential.(Ref(hydraulics), S_l)
ψ_bc = matric_potential.(Ref(hydraulics_bc), S_l)
plot(
    S_l,
    log10.(-ψ),
    xlabel = "effective saturation",
    ylabel = "Log10(|ψ|)",
    label = "van Genuchten",
)
plot!(S_l, log10.(-ψ_bc), label = "Brooks and Corey")
savefig("bc_vg_matric_potential.png")
# ![](bc_vg_matric_potential.png)

# The steep slope in
# `ψ` near saturated and completely dry soil are part of the reason
# why Richard's equation is such a challenging numerical problem.

# # Hydraulic conductivity
# The hydraulic conductivity is a more complex function than the matric potential,
# as it depends on the temperature of the water, the volumetric ice fraction, and
# the volumetric liquid water fraction. It also depends on the `hydraulics` model
# chosen.

# We represent the hydraulic conductivity `K` as the product of four factors:
# `Ksat`, an `impedance_factor` (which accounts for the effect of ice on conductivity)
# a `viscosity_factor` (which accounts for the effect of temperature on the
# viscosity of liquid water, and how that in turn affects conductivity)
# and a `moisture_factor` (which accounts for the effect of liquid water, and is determined by the `hydraulics` model).

#  Let's start with ice and temperature independence, but moisture dependence.
# Below
# we choose additional parameters, consistent with the `hydraulics` parameters
# for sandy loam ([Bonan19a](@cite)).
ν = FT(0.41)
Ksat = FT(4.42 / (3600 * 100))
moisture_choice = MoistureDependent{FT}()
viscosity_choice = ConstantViscosity{FT}()
impedance_choice = NoImpedance{FT}();

# We are going to calculate `K = Ksat × viscosity_factor × impedance_factor × moisture_factor`.
# In the code, each of these factors is a function with multiple methods, except for `Ksat`.
#  Our function `hydraulic_conductivity` calls each of these functions in turn,
# and these functions use multiple dispatch to provide the correct value for `K`.

# Like we defined new type
# classes for `vanGenuchten{FT}` and `BrooksCorey{FT}`, we also created new type classes
# for the impedance choice, the viscosity choice, and the moisture choice.
#  For example, the function
# called `viscosity_factor`, when passed an object of type `ConstantViscosity{FT}`, executes a method that always
# returns 1. The same is true for the function `impedance_factor`, using the type
# `NoImpedance{FT}`, and for `moisture_factor`, using the type `MoistureIndependent{FT}`.

# In the case where the `moisture_factor = MoistureDependent{FT}()`, either the van
# Genuchten or Brooks and Corey expression is used based on the type of `hydraulics` model
# passed.

# The `moisture_factor` for the van Genuchten model is (denoting it as ``K_m``)

# ``
#  K_m = \sqrt{S_l}[1-(1-S_l^{1/m})^m]^2,
# ``

# for ``S_l < 1``,

# and for the Brooks and Corey model it is

# ``
# K_m = S_l^{2M+3},
# ``

# also for ``S_l<1``. When ``S_l\geq 1``, ``K_m = 1`` for each model.

# One side effect of this flexibility is that `hydraulic_conductivity`
# requires all the arguments it could possibly need passed to it, which is
# why here we must supply a value `T` and `θ_i`, even though they are not used
# in this particular example.
T = FT(0.0)
θ_i = FT(0.0)

K =
    Ksat .*
    hydraulic_conductivity.(
        Ref(impedance_choice),
        Ref(viscosity_choice),
        Ref(moisture_choice),
        Ref(hydraulics),
        Ref(θ_i),
        Ref(ν),
        Ref(T),
        S_l,
    );
# Let's also compute `K` when we include the effects of temperature
# and ice on the hydraulic conductivity.
# In the cases where a `TemperatureDependentViscosity{FT}` or
# `IceImpedance{FT}` type is passed, the correct factors are calculated,
# based on the temperature `T` and volumetric ice fraction `θ_i`.
# In these cases, the `viscosity_factor`, denoted here as ``K_v``,
# evaluates as:

# ``
# K_v = e^{\gamma(T-T_{ref})},
# ``

# where ``\gamma = 0.0264 \mbox{K}^{-1}`` is an empirical factor,
# and ``T_{ref} = 288`` K,
# and the `impedance_factor`, denoted ``K_i``, evaluates as:

# ``
# K_i = 10^{-\Omega f_i}
# ``

# where ``\Omega = 7`` is an empirical factor and
# ``f_i`` is the ratio of the volumetric
# ice fraction to total volumetric water fraction  ([Lundin1990](@cite)).

viscosity_choice_T = TemperatureDependentViscosity{FT}()
T = FT(300.0)
K_T =
    Ksat .*
    hydraulic_conductivity.(
        Ref(impedance_choice),
        Ref(viscosity_choice_T),
        Ref(moisture_choice),
        Ref(hydraulics),
        Ref(θ_i),
        Ref(ν),
        Ref(T),
        S_l,
    )
ice_impedance_I = IceImpedance{FT}()
θ_i = FT(0.1)
S_i = θ_i / ν;
# The total volumetric water fraction cannot
# exceed unity, so the effective liquid water saturation
# should have a max of 1-S_i.
S_l_accounting_for_ice = FT.(0.01:0.01:(0.99 - S_i))
K_i =
    Ksat .*
    hydraulic_conductivity.(
        Ref(ice_impedance_I),
        Ref(viscosity_choice),
        Ref(moisture_choice),
        Ref(hydraulics),
        Ref(θ_i),
        Ref(ν),
        Ref(T),
        S_l_accounting_for_ice,
    )
plot(
    S_l,
    log10.(K),
    xlabel = "total effective saturation, (θ_i+θ_l)/ν",
    ylabel = "Log10(K)",
    label = "Base case",
    legend = :bottomright,
)
plot!(S_l, log10.(K_T), label = "Temperature Dependent Viscosity; no ice")
plot!(
    S_l_accounting_for_ice .+ S_i,
    log10.(K_i),
    label = "Ice Impedance; S_i=0.24",
)
savefig("T_ice_K.png")
# ![](T_ice_K.png)

# If the user is not considering phase transitions
# and does not add in Freeze/Thaw source terms, the default is for zero
# ice in the model, for all time and space. In this case the ice impedance
# factor evaluates to 1 regardless of which type is passed.

# We can also look and see how the Brooks and Corey moisture factor differs from the
# van Genuchten moisture factor by changing the `hydraulics` model passed:
T = FT(0.0)
θ_i = FT(0.0)

K_bc =
    Ksat .*
    hydraulic_conductivity.(
        Ref(impedance_choice),
        Ref(viscosity_choice),
        Ref(moisture_choice),
        Ref(hydraulics_bc),
        Ref(θ_i),
        Ref(ν),
        Ref(T),
        S_l,
    )
plot(
    S_l,
    log10.(K),
    xlabel = "effective saturation",
    ylabel = "Log10(K)",
    label = "van Genuchten",
)
plot!(
    S_l,
    log10.(K_bc),
    xlabel = "effective saturation",
    ylabel = "Log10(K)",
    label = "Brooks and Corey",
)
savefig("bc_vg_k.png")
# ![](bc_vg_k.png)
# # Other features
# The user also has the choice of making the conductivity constant by choosing
# `MoistureIndependent{FT}()` along with `ConstantViscosity{FT}()` and
# `NoImpedance{FT}()` . This is useful for debugging!
no_moisture_dependence = MoistureIndependent{FT}()
K_constant =
    Ksat .*
    hydraulic_conductivity.(
        Ref(impedance_choice),
        Ref(viscosity_choice),
        Ref(no_moisture_dependence),
        Ref(hydraulics),
        Ref(θ_i),
        Ref(ν),
        Ref(T),
        S_l,
    );
# ```julia
# julia> unique(K_constant)
# 1-element Array{Float32,1}:
#  1.2277777f-5
# ```

# Note that choosing this option does not mean the matric potential
# is constant, as a hydraulics model is still required and employed.

# And, lastly, you might be wondering why we left `Ksat` out of the function
# for `hydraulic_conductivity`. It turns out it is also useful for debugging
# to be able to turn off the flow of water, by setting `Ksat = 0`.

# # References
# - [vanGenuchten1980](@cite)
# - [BrooksCorey1964](@cite)
# - [Corey1977](@cite)
# - [Lundin1990](@cite)
# - [Bonan19a](@cite)


text:
"""

以下のプログラムを Julia に翻訳してみる。

星座グラフ
http://aoki2.si.gunma-u.ac.jp/R/Constellation.html

using RDatasets
iris = dataset("datasets", "iris")
dat = Matrix(iris[:, 1:4])
col = vcat(fill(1, 50), fill(2, 50), fill(3, 50))
constellationgraph(dat, col=col)

"""

using Statistics, Plots

function constellationgraph(dat; w=0, points=100, col=1, alpha=0.2)
    if w == 0
        nv = size(dat, 2)
        w = fill(1/nv, nv)
    end
    mx = maximum(dat, dims=1)
    mn = minimum(dat, dims=1)
    rg = mx-mn
    dat = π .* (dat .- mn) ./ rg
    y = sin.(dat) * w
    x = cos.(dat) * w
    pyplot()
    scatter(x, y, color=col, alpha=alpha, aspect_ratio=1,
            xlims=(-1.05, 1.05), ylims=(0, 1.05), grid=false,
            yticks=false, yshowaxis=false, label="")
    θ = range(0, π, length=points)
    plot!(cos.(θ), sin.(θ), color=:black, label="", linewidth=0.5)
    xticks!([-1, 0, 1], tick_direction=:out)
end

ext: jl

max_stars_repo_name: r-de-r/stats

max_issues_repo_name: r-de-r/stats

max_forks_repo_name: r-de-r/stats


text:
using Base.Test, DataFrames, EventHistory, StatsModels

function test_show(x)
    io = IOBuffer()
    show(io, x)
end

d = DataFrame(start = [1,2,5,2,1,7,3,4,8,8],
              stop  = [2,3,6,7,8,9,9,9,14,17],
              event = [1,1,1,1,1,1,1,0,0,0],
              x     = [1,0,0,1,0,1,1,1,0,0]);
d[:S] = EventHistory.Event([:start,:stop,:event], d);
mm = EventHistory.phreg(@formula(S~x), d)

@test EventHistory.coef(mm) ≈ [-0.02110521] atol=1e-8
@test mm.invhess ≈ [0.6323063] atol=1e-6

@test mm.grad[1]<1e-6

p = EventHistory.predict(mm, X=[0; 1], time=[6]; order=true);
@test isapprox(p[2:3], [0.3522254,0.3599854], atol=1e-6)

d = DataFrame(start = [1,2,5,2,1,7,3,4,8,8],
              stop  = [2,3,6,7,8,9,9.5,10,14,17],
              event = [1,1,1,1,1,1,1,0,0,0],
              x     = [1,0,0,1,0,1,1,1,0,0],
              z     = [1.0,0,2.0,0,3.0,0,4.0,0,5.0,0]);
d[:S] = EventHistory.Event([:start,:stop,:event], d);
mm = EventHistory.phreg(@formula(S~x+z), d)
@test length(coef(mm)) == 2


text:
# ------------------------------------------------------------------------------------------
# # 软件包
#
# Julia 有超过 2000 个登记在案的包，这使得第三方的软件包成为 Julia 语言生态的重要组成部分。
# 尽管如此，软件包生态依然有很大的成长空间。
#
# 值得一提的是，在 Julia 中调用其它语言的函数也是一等公民。Julia 提供了很棒的外部函数接口（FFI）。
# 我们可以很容易地通过 `PyCall` 或 `Rcall` 调用 Python 或 R。
# 这意味着你不必等到 Julia 的生态完全成熟就能迁移现有项目到 Julia，
# 而且迁移到 Julia 并不意味着放弃你在其他语言中惯用的包或者库！
#
# 查看全部可获得的包，请看
# https://pkg.julialang.org/
# 或者
# https://juliaobserver.com/
#
# 现在我们先学习一下该如何使用包。
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# 在安装好的 Julia 中第一次要使用一个包的时候，你需要用包管理器（package manager）明确添加（add）它：
# ------------------------------------------------------------------------------------------

using Pkg
Pkg.add("Example")

# ------------------------------------------------------------------------------------------
# 之后每次使用 Julia 时（比如打开一个新的 REPL 会话或者第一次打开一个 notebook），
# 你需要通过 `using` 关键字载入要用的包
# ------------------------------------------------------------------------------------------

using Example

# ------------------------------------------------------------------------------------------
# 在 `Example.jl` 的源码
# https://github.com/JuliaLang/Example.jl/blob/master/src/Example.jl
# 中我们可以看到声明了下面这个函数
#
# ```julia
# hello(who::String) = "Hello, $who"
# ```
#
# 载入了 `Example` 我们现在应该可以调用 `hello` 函数了
# ------------------------------------------------------------------------------------------

hello("it's me. I was wondering if after all these years you'd like to meet.")

# ------------------------------------------------------------------------------------------
# 现在让我们来玩 `Colors` 包吧
# ------------------------------------------------------------------------------------------

Pkg.add("Colors")

using Colors

# ------------------------------------------------------------------------------------------
# 首先创建一个 100 个不同颜色的调色板
# ------------------------------------------------------------------------------------------

palette = distinguishable_colors(100)

# ------------------------------------------------------------------------------------------
# 然后通过 `rand` 函数来创建随机矩阵
# ------------------------------------------------------------------------------------------

rand(3, 3)

rand(1:10, 3, 3)

rand(palette, 3, 3)

max_stars_repo_name: NahsiN/JuliaTutorials

max_issues_repo_name: NahsiN/JuliaTutorials

max_forks_repo_name: NahsiN/JuliaTutorials


text:
module FixedPointAcceleration
using LinearAlgebra: cond, pinv
using GLM: fit, LinearModel
using Dates: now
include("0_Structs_And_Enums.jl")
include("1_MainFunctions.jl")
export fixed_point, fixed_point_new_input # Main functionality.
export FunctionEvaluationResult, FixedPointResults # Structs.
export put_together_without_jumps, execute_function_safely # Auxillery functions. Mainly exported for testing.
# Exporting enums
export FixedPointAccelerationAlgorithm, Simple, Anderson, Aitken, Newton, VEA, SEA, MPE, RRE # Algorithms
export InvalidReplacement, NoAction, ReplaceElements, ReplaceVector # Strategies for handling invalids in proposed new inputs.
export FP_FunctionEvaluationError, NoError, ErrorExecutingFunction, LengthOfOutputNotSameAsInput, InputMissingsDetected
export InputNAsDetected, InputInfsDetected, OutputMissingsDetected, OutputNAsDetected, OutputInfsDetected # Errors in evaluating
export TerminationCondition, ReachedConvergenceThreshold, ReachedMaxIter, InvalidInputOrOutputOfIteration, FunctionIsNotTypeStable # Termination Conditions
end


text:
# push!(LOAD_PATH,"../src/")

using Documenter, FracFin

# makedocs()

makedocs(
    # options
    modules = FracFin,
    format = :html,
    sitename = "FracFin.jl",
    authors = "Yann Calec",
    pages = Any[
        "Home" => "index.md",
        "Getting started" => "getting_started.md",
        "Usage" => "usage.md"
    ]
)

# deploydocs(
#     deps = Deps.pip("mkdocs"),
#     repo = "github.com/yanncalec/FracFin.jl.git",
#     julia = "0.6"
# )

max_stars_repo_name: yanncalec/FracFin

max_issues_repo_name: yanncalec/FracFin

max_forks_repo_name: yanncalec/FracFin


text:
include("INCLUDEME.jl")

using WFLearning

function main(n; nbatch=20, epochs=200, learning_rate=0.1)
    pmap(nlayers->task(n, nlayers, nbatch, epochs, learning_rate), 2:20)
end

# main(6, 5, nbatch=10, nsamples=20, epochs=2, learning_rate=0.1)

max_stars_repo_name: Roger-luo/QCBM

max_issues_repo_name: Roger-luo/QCBM

max_forks_repo_name: Roger-luo/QCBM


text:
using Traitor
using Test

abstract type Size end

struct Big    <: Size end
struct Medium <: Size end
struct Small  <: Size end

Size(::Union{Type{Int32},Type{Int64}}) = Small
Size(::Type{BigInt}) = Big
Size(::Type{Int128}) = Medium

@traitor function howbig(x::Any::Big)
    "Huge!"
end

@traitor function howbig(x::::Medium)
    "So-so"
end
@test howbig(Int128(1)) == "So-so"
@test howbig(BigInt(1)) == "Huge!"

@traitor function howbig(x::Integer::Small)
    "Teensy..."
end

@test howbig(1) == "Teensy..."
# These next ones throw...
# @test howbig(Int128(1)) == "So-so"
# @test howbig(BigInt(1)) == "Huge!"
# ...until we define
@traitor howbig(x::Integer::Big) = "Huge!"
@traitor howbig(x::Integer::Medium) = "So-so"

@test howbig(1) == "Teensy..."
@test howbig(Int128(1)) == "So-so"
@test howbig(BigInt(1)) == "Huge!"

@test @inferred(supertrait(Small)) == Size
@test @inferred(supertrait(Union{Small,Medium})) == Size

abstract type Fooness end

struct FooA <: Fooness end
struct FooB <: Fooness end

Fooness(::Type{Any}) = FooA
Fooness(::Type{Int16}) = FooB

@test_throws ErrorException supertrait(Union{Small,FooB})

module SomeoneElsesCode
f(x) = x + 1
end

import .SomeoneElsesCode: f
betray!(f, Tuple{Any})

@traitor f(x::::Big) = x - 1

@test f(1) == 2
@test f(BigInt(1)) == 0

Size(::Type{Bool}) = Small
@test howbig(true) == "Teensy..."


text:
module Buffers

using MacroTools: @forward
using ZygoteRules
using ZygoteRules: AContext

export Buffer, bufferfrom

include("lib.jl")
include("buffer.jl")

end # module


text:
# ensure in right directory and environment
cd(dirname(@__FILE__))
using Pkg
Pkg.activate(".")

using Documenter, MuseInference, Turing

makedocs(
    sitename = "MuseInference",
    format = Documenter.HTML(
        assets = ["assets/muse.css"],
        disable_git = true,
    ),
    pages = [
        "index.md",
        "api.md",
    ],
    checkdocs = :none,
    strict = true,
    modules = [MuseInference]
)

deploydocs(
    repo = "https://github.com/marius311/MuseInference.jl",
    devbranch = "main",
    devurl = "latest",
    push_preview = true,
    forcepush = true
)


text:
using CBM
using Base.Test

function test_all()
    model = load_json(Pkg.dir() * "/CBM/Models/e_coli_core.json")

    println("test_add_reaction")
    test_add_reaction(model)
    println("test_change_objective")
    test_change_objective(model)
    println("test_change_reaction_bounds")
    test_change_reaction_bounds(model)
    println("test_fba")
    test_fba(model)
    println("test_find_blocked_reactions")
    test_find_blocked_reactions(model)
    println("test_find_deadend_metabolites")
    test_find_deadend_metabolites(model)
    println("test_find_exchange_reactions")
    test_find_exchange_reactions(model)
    println("test_find_reactions_from_gene")
    test_find_reactions_from_gene(model)
    println("test_find_reactions_from_metabolite")
    test_find_reactions_from_metabolite(model)
    println("test_fva")
    test_fva(model)
    println("test_gene_deletion")
    test_gene_deletion(model)
    println("test_remove_reaction")
    test_remove_reaction(model)
    println("test_robustness_analysis")
    test_robustness_analysis(model)
    println("test_solvers")
    test_solvers(model)
end

test_all()

max_stars_repo_name: isebarn/CBM

max_issues_repo_name: isebarn/CBM

max_forks_repo_name: isebarn/CBM


text:
export segy_scan

"""
    segy_scan(dir::String, filt::String, keys::Array{String,1}, blocksize::Int;
              chunksize::Int = CHUNKSIZE,
              pool::WorkerPool = WorkerPool(workers()),
              verbosity::Int = 1)

    returns: SeisCon

Scan header fields `keys` of files in `dir` matching the filter `filt` in blocks
containing `blocksize` continguous traces. The scanning of files is distributed to workers
in `pool`, the default pool is all workers.

`chunksize` determines how many MB of data will be loaded into memory at a time.

`CHUNKSIZE` defaults to 2048MB.

`verbosity` set to 0 silences updates on the current file being scanned.

"""
function segy_scan(dir::String, filt::Union{String, Regex}, keys::Array{String,1}, blocksize::Int;
                   chunksize::Int = CHUNKSIZE,
                   pool::WorkerPool = WorkerPool(workers()),
                   verbosity::Int = 1,
                   filter::Bool = true)

    endswith(dir, "/") ? nothing : dir *= "/"
    filter ? (filenames = searchdir(dir, filt)) : (filenames = [filt])
    files = map(x -> dir*x, filenames)
    files_sort = files[sortperm(filesize.(files), rev = true)]
    run_scan(f) = scan_file(f, keys, blocksize, chunksize=chunksize, verbosity=verbosity)
    s = pmap(run_scan, pool, files_sort)

    return merge(s)
end

"""
    segy_scan(dirs::Array{String,1}, filt::String, keys::Array{String,1}, blocksize::Int;
              chunksize::Int = CHUNKSIZE,
              pool::WorkerPool = WorkerPool(workers()),
              verbosity::Int = 1)

Scans all files whose name contains `filt` in each directory of `dirs` using `blocksize`.

"""
function segy_scan(dirs::Array{String,1}, filt::Union{String, Regex}, keys::Array{String,1}, blocksize::Int;
                   chunksize::Int = CHUNKSIZE,
                   pool::WorkerPool = WorkerPool(workers()),
                   verbosity::Int = 1,
                   filter::Bool = true)

    files = Array{String,1}()
    for dir in dirs
        endswith(dir, "/") ? nothing : dir *= "/"
        filter ? (filenames = searchdir(dir, filt)) : (filenames = [filt])
        append!(files, map(x -> dir*x, filenames))
    end
    files_sort = files[sortperm(filesize.(files), rev = true)]
    run_scan(f) = scan_file(f, keys, blocksize, chunksize=chunksize, verbosity=verbosity)
    s = pmap(run_scan, pool, files_sort)

    return merge(s)
end

"""
    segy_scan(dir::String, filt::String, keys::Array{String,1},
              chunksize::Int = CHUNKSIZE,
              pool::WorkerPool = WorkerPool(workers()),
              verbosity::Int = 1)

If no `blocksize` is specified, the scanner automatically detects source locations and returns
blocks of continguous traces for each source location, but each block no larger then CHUNKSIZE.

"""
function segy_scan(dir::String, filt::Union{String, Regex}, keys::Array{String,1};
                   chunksize::Int = CHUNKSIZE,
                   pool::WorkerPool = WorkerPool(workers()),
                   verbosity::Int = 1,
                   filter::Bool = true)

    endswith(dir, "/") ? nothing : dir *= "/"
    filter ? (filenames = searchdir(dir, filt)) : (filenames = [filt])
    files = map(x -> dir*x, filenames)
    files_sort = files[sortperm(filesize.(files), rev = true)]
    run_scan(f) = scan_file(f, keys, chunksize=chunksize, verbosity=verbosity)
    s = pmap(run_scan, pool, files_sort)

    return merge(s)
end

"""
    segy_scan(dirs::Array{String,1}, filt::String, keys::Array{String,1},
              chunksize::Int = CHUNKSIZE,
              pool::WorkerPool = WorkerPool(workers()),
              verbosity::Int = 1)

Scans all files whose name contains `filt` in each directory of `dirs`.

"""
function segy_scan(dirs::Array{String,1}, filt::Union{String, Regex}, keys::Array{String,1};
                   chunksize::Int = CHUNKSIZE,
                   pool::WorkerPool = WorkerPool(workers()),
                   verbosity::Int = 1,
                   filter::Bool = true)

    files = Array{String,1}()
    for dir in dirs
        endswith(dir, "/") ? nothing : dir *= "/"
        filter ? (filenames = searchdir(dir, filt)) : (filenames = [filt])
        append!(files, map(x -> dir*x, filenames))
    end
    files_sort = files[sortperm(filesize.(files), rev = true)]
    run_scan(f) = scan_file(f, keys, chunksize=chunksize, verbosity=verbosity)
    s = pmap(run_scan, pool, files_sort)

    return merge(s)
end

function searchdir(path, filt::String)
    filt = Regex(replace(filt, "*" => ".*."))
    filter!(x->occursin(filt, x),readdir(path))
end

function searchdir(path, filt::Regex)
    filter!(x->occursin(filt, x),readdir(path))
end


text:
## update inputs vars using calculated gradients

"Check if an object is of a struct type, i.e. not a number or array"
isstruct(::Type{T}) where T = !isempty(fieldnames(T))
isstruct(obj) = isstruct(typeof(obj))

function path_value_pairs(gx::Tangent; current_path=[])
    result = []
    ks = collect(keys(gx))
    vals = collect(gx)
    for (k, v) in zip(ks, vals)
        if v isa Tangent
            subresults = path_value_pairs(v; current_path=vcat(current_path, k))
            result = vcat(result, subresults)
        else
            pair = tuple(vcat(current_path, k)...) => v
            push!(result, pair)
        end
    end
    return result
end

function getfield_nested(s, path)
    v = s
    for name in path
        v = getfield(v, name)
    end
    return v
end

function setfield_nested!(s, path, val)
    container = s
    for name in path[1:end-1]
        container = getfield(container, name)
    end
    setfield!(container, path[end], val)
end

minus_updater(x::AbstractArray, gx) = (x .- gx)
minus_updater(x::Real, gx) = (x - gx)

function update!(x::AbstractArray{T,N}, gx::AbstractArray{T,N}, fn::Function=minus_updater) where {T,N}
    x .= fn(x, gx)
end

function update!(x::Real, gx::Real, fn::Function=minus_updater)
    error("Can't update value of $(typeof(x)) in-place!")
end

function update!(x, gx, fn::Function=minus_updater; ignore=Set())
    @assert isstruct(x) "Expected mutable struct as 1st argument"
    @assert gx isa Tangent "Gradient for structs should be Tangent"
    for (path, gv) in path_value_pairs(gx)
        if !in(path, ignore)
            v = getfield_nested(x, path)
            setfield_nested!(x, path, fn(v, gv))
        end
    end
end


text:
using ReachabilityAnalysis
using ExponentialUtilities
using Plots
using SparseArrays
using ReachabilityAnalysis: dim

_homogeneize(p::SecondOrderAffineContinuousSystem) = _homogeneize(p.M, p.K, p.b)

function _homogeneize(M, K, F)
    m = size(K, 1)
    Zm = spzeros(m, m)
    Im = sparse(1.0I, m, m)
    M = Diagonal(M)
    invM = inv(M)

    u0 = invM * F
    Aext = [Zm           Im      spzeros(m) ;
           -invM*K       Zm         u0      ;
            spzeros(1, 2m+1)                ]

    S = @system(x' = Aext*x) # sistema homogeneizado
end

function _reach_homog_krylov_LGG09_modif!(out, Ω₀::LazySet, Aᵀδ::AbstractMatrix,
                                          ℓ::AbstractVector, NSTEPS;
                                          hermitian=false, m=min(30, size(Aᵀδ, 1)), tol=1e-7)

    @assert size(out, 1) == 2
    @assert size(out, 2) == NSTEPS

    # initialization of the krylov subspace
    # esta parte genera una base de
    # \{Aᵀδ * ℓ, (Aᵀδ)^2 * ℓ, (Aᵀδ)^3 * ℓ, ..., (Aᵀδ)^m * ℓ}
    TA, Tb = eltype(Aᵀδ), eltype(ℓ)
    T = promote_type(TA, Tb)
    Ks = KrylovSubspace{T, real(T)}(length(ℓ), m)
    arnoldi!(Ks, Aᵀδ, ℓ; m=m, ishermitian=hermitian, tol=tol)

    # rᵢ stores is the cache for each vector: (Φᵀ)^i ℓ
    rᵢ = deepcopy(ℓ)

    @inbounds for i in 1:NSTEPS
        out[1, i] = ρ(rᵢ, Ω₀)
        out[2, i] = ρ(-rᵢ, Ω₀)
        expv!(rᵢ, i*1.0, Ks) # rᵢ <- exp(A^T δ * i) * ℓ
    end
    return out
end

# ------------------------------------------------------------------------
# Auxiliary methods that are not yet available in ReachabilityAnalysis.jl
# ------------------------------------------------------------------------

# calcula: out <- exp(A * NSTEPS * dt) * b
function _expv(A, b, NSTEPS, dt; hermitian=false, m=min(30, size(A, 1)), tol=1e-7)

    # initialization of the krylov subspace
    TA, Tb = eltype(A), eltype(b)
    T = promote_type(TA, Tb)
    Ks = KrylovSubspace{T, real(T)}(length(b), m)
    arnoldi!(Ks, A, b; m=m, ishermitian=hermitian, tol=tol)

    out = similar(b)
    expv!(out, NSTEPS*dt, Ks)

    return out
end

# calcula: out <- exp(A * NSTEPS * dt) * b
function _phiv(A, b, NSTEPS, dt; hermitian=false, m=min(30, size(A, 1)), tol=1e-7)

    # initialization of the krylov subspace
    TA, Tb = eltype(A), eltype(b)
    T = promote_type(TA, Tb)
    Ks = KrylovSubspace{T, real(T)}(length(b), m)
    arnoldi!(Ks, A, b; m=m, ishermitian=hermitian, tol=tol)

    out = Matrix{Float64}(undef, size(A, 1), 3)
    phiv!(out, NSTEPS*dt, Ks, 2)

    return view(out, :, 3) .* dt^2
end

#
# E₊ = ⊡(Φ₂ ⊡(A^2 X₀))
#
function _build_Eplus(A, X0, δt; m=min(30, size(A, 1)), tol=1e-7, Φ₂can=nothing)
    if isa(X0, Vector)
        X0 = Singleton(X0)
    end
    n = dim(X0)

    # build box E+
    H = symmetric_interval_hull(A^2 * X0)
    E₊_high = Vector{Float64}(undef, n)

    if isnothing(Φ₂can)
        Φ₂can = Φ₂_canonical(A, δt, m=m, tol=tol)
    end

    for i in 1:n
        Φ₂col = view(Φ₂can, :, i)
        E₊_high[i] = max(ρ(Φ₂col, H), abs(ρ(-Φ₂col, H)))
    end

    E₊_low = -E₊_high # since it is a SIH
    E₊ = Hyperrectangle(low=E₊_low, high=E₊_high)
end

function Φ₂_canonical(A, δt; m=min(size(A, 1)), tol=1e-7)
    n = size(A, 1)
    P3n = ReachabilityAnalysis._P_3n(abs.(A), δt, n)
    P3n_tr = copy(transpose(P3n))

    Φ₂can = Matrix{Float64}(undef, n, n)
    b = zeros(3n)
    for i in 1:n
        b[i] = 1.0
        out = _expv(P3n_tr, b, 1, 1.0, m=m, tol=tol)
        Φ₂can[:, i] .= view(out, 2n+1:3n)
        b[i] = 0.0
    end
    return Φ₂can
end

function _support_function_Eplus_baseline(d, A, δ, V)
    P = ReachabilityAnalysis.Φ₂(abs.(A), δ)
    E = symmetric_interval_hull(P * V)
    return ρ(d, E)
end

function _support_function_Eplus_krylov(d, A::AbstractMatrix{N}, δ,
                           V::AbstractHyperrectangle;
                           cutoff=eps(N), mk=min(30, size(A, 1)), tol=1e-7) where {N}

    # use only the indices in which V is not flat
    idx = findall(ri -> abs(ri) > cutoff, radius_hyperrectangle(V))
    Aabsᵀ = copy(transpose(abs.(A))) # sparse

    out = zero(N)
    n = size(A, 1)
    ei = zeros(N, n)
    @inbounds for i in idx
        ei[i] = one(N)

        R1 = phiv(δ, Aabsᵀ, ei, 2, m=mk, tol=tol) .* δ^2
        α = ρ(view(R1, :, 3), V)

        R2 = phiv(δ, Aabsᵀ, -ei, 2, m=mk, tol=tol) .* δ^2
        β = ρ(view(R2, :, 3), V)

        out += abs(d[i]) * max(α, abs(β))
        ei[i] = zero(N)
    end
    return out
end

function _support_function_Eplus_krylov_optim(d::AbstractVector{N}, A::AbstractMatrix{N}, δ,
                                              V::AbstractHyperrectangle;
                                              cutoff=eps(N), mk=min(30, size(A, 1)), tol=1e-7) where {N}

    # use only the indices in which V is not flat
    idx = findall(ri -> abs(ri) > cutoff, radius_hyperrectangle(V))
    Aabs = copy(abs.(A)) # sparse

    v = V.radius[idx]
    Maux = phiv(δ, Aabs, v, 2, m=mk, tol=tol)
    Pv = view(Maux, :, 3) .* δ^2 # rescale

    return dot(abs.(d), Pv)
end

function _build_Eplus_optim(A::SparseMatrixCSC{N, D}, X0::AbstractVector{N}, δt; m=min(30, size(A, 1)), tol=1e-7, cutoff=eps(N)) where {N, D}
    _build_Eplus_optim(A, Singleton(X0), δt; m=m, tol=tol, cutoff=cutoff)
end

# optimized version of the constructor of the Eplus set using a single Krylov subspace evaluation
# uses the idea in _support_function_Eplus_krylov_optim
function _build_Eplus_optim(A::SparseMatrixCSC{N, D}, X0::AbstractSingleton{N}, δt; m=min(30, size(A, 1)), tol=1e-7, cutoff=eps(N)) where {N, D}
    n = dim(X0)

    # compute V
    AX0 = linear_map(A, X0)
    V = symmetric_interval_hull(A * AX0)
    v = V.radius

    # compute Pv using Krylov
    Aabs = copy(abs.(A))
    Pv = _phiv(Aabs, v, 1, δt; m=m, tol=tol)

    # return the symmetric interval hull whose radius is Pv
    E⁺ = Hyperrectangle(zeros(n), Pv)

    return E⁺
end

function _build_Eplus_optim(A::SparseMatrixCSC{N, D}, X0::Hyperrectangle{N}, δt; m=min(30, size(A, 1)), tol=1e-7, cutoff=eps(N)) where {N, D}
    n = dim(X0)
    A2 = A * A; # fast if A sparse
    V = symmetric_interval_hull(A2 * X0)
    v = V.radius
    Aabs = copy(abs.(A))
    Pv = _phiv(Aabs, v, 1, δt; m=m, tol=tol)
    E⁺ = Hyperrectangle(zeros(n), Pv)
end

max_stars_repo_name: JuliaReach/SetPropagation-FEM-Examples

max_issues_repo_name: JuliaReach/SetPropagation-FEM-Examples

max_forks_repo_name: JuliaReach/SetPropagation-FEM-Examples


text:
#=
Functions for the McEliece cryptosystem using Goppa codes.

Author: Ivan A. Moreno Soto
Last updated: 2019/July/15
=#

module McEliece

export
    generateKeys,
    saveKeys,
    encrypt,
    decrypt,
    haltingCheck

include("./Goppa.jl")
using .GoppaCode

include("./FiniteFieldMatrix.jl")
using .FiniteFieldMatrix

include("./Misc.jl")
using .Misc

using Random
using StatsBase
using DelimitedFiles

"""
    haltingCheck(publicKey, privateKey, samples=2000)

Checks that a given pair of keys correctly finishes, i.e. always determines whether or not a cyphertext
is valid or not, by trying to decrypt random vectors.
"""
function haltingCheck(publicKey, privateKey, samples=2000)
    t = publicKey[2]
    m = privateKey[5]
    k = size(publicKey[1], 1)
    n = size(publicKey[1], 2)

    failures = 0
    for i in 1:samples
        test = [Array{Bool, 1}([rand([0, 1]) for _ in 1:n])]
        try
            decrypt(test, publicKey[1], privateKey[1], privateKey[2], privateKey[3], privateKey[4], m, t)
        catch exc
            if isa(exc, InvalidCodewordError)
                failures += 1
                if i % 100 == 0
                    println("Already $(i) steps.")
                end
                continue
            else
                throw(exc)
            end
        end
    end

    return failures
end

"""
    generateKeys(n, t, m)

Generates the public and private keys for a McEliece cryptosystem with length n
and with capacity to correct up to t errors.
"""
function generateKeys(n, t, m)
    if n > 2^m
        error("Extension field can't support $(n) elements.")
    end
    if n - m*t < 1
        error("Incompatible parameters: $(n) - $(m)*$(t) < 0.")
    end

    G = GInverse = goppaPolynomial = codeSupport = k = [] # Just to make sure these variables exist in this scope.
    keysCreated = false
    while !keysCreated
        try
            # Generating keys until we have a G with full rank.
            G, goppaPolynomial, codeSupport, k = generateGoppaCode(n, t, m)
            GInverse = GeneratorInverse!(G)
            keysCreated = true
        # To save time, it always tries to invert G right away. If G is not invertible,
        # the function throws an exception that gets catched and only indicates that another
        # goppa code is needed.
        catch exc
            # throw(exc)
            continue
        end
    end

    P = getPermutation(n)
    S, SInverse = getScrambler(k)

    GPub = (S * G) * P
    GPub = GPub .% 2 # S * GP can break out of GF(2)

    PInverse = Array{UInt32, 2}(inv(P))
    GSInverse = (GInverse * SInverse) .% 2

    return (GPub, t), (GSInverse, PInverse, goppaPolynomial, codeSupport, m)
end

"""
    saveKeys(publicKey, privateKey, publicPrefix="public", privatePrefix="private")

Saves the given keys to text files in a format friendly for Julia.
"""
function saveKeys(publicKey, privateKey, publicPrefix="public", privatePrefix="private")
    writedlm(publicPrefix * "_g.txt", publicKey[1])
    open(publicPrefix * "_t.txt", "w") do file
        write(file, string(publicKey[2]))
    end

    writedlm(privatePrefix * "_gs.txt", privateKey[1])
    writedlm(privatePrefix * "_p.txt", privateKey[2])
    writedlm(privatePrefix * "_polynomial.txt", privateKey[3])
    writedlm(privatePrefix * "_support.txt", privateKey[4])

    open(privatePrefix * "_m.txt", "w") do file
        write(file, string(privateKey[5]))
    end
end

"""
    getErrorVector(n, t)

Returns a random binary vector of n components with weigth t.
"""
function getErrorVector(n, t)
    positions = sample(1:n, t, replace=false)
    error = falses(n)

    for pos in positions
        error[pos] = true
    end

    return error
end

"""
    encrypt(plaintext, GPub, t)

Takes an arbitrarily long plaintext and encrypts with a public key from a McEliece
cryptosystem. Returns size(plaintext)/k codewords with t errors.
"""
function encrypt(plaintext, GPub, t)
    k = size(GPub)[1]
    n = size(GPub)[2]
    m = string2bits(plaintext, k)

    cipher = []
    for i in 1:div(size(m)[1], k)
        e = getErrorVector(n, t)
        codeword = (transpose(m[((i-1)*k + 1):(i*k)]) * GPub) .% 2

        subc = falses(n)
        for i in 1:n
            subc[i] = (codeword[i] + e[i]) % 2
        end

        push!(cipher, subc)
    end

    return cipher
end

"""
    isPlaintextValid(plaintext, cyphertext, G, t)

Returns true if (plaintext * G) + cyphertext has t ones, false otherwise.
"""
function isPlaintextValid(plaintext, cyphertext, G, t)
    codeword = (plaintext * G) .% 2
    errorVector = falses(length(codeword))
    for i in 1:length(codeword)
        errorVector[i] = (codeword[i] + cyphertext[i]) % 2
    end

    errors = 0
    for e in errorVector
        if e == 1
            errors += 1
        end
    end

    return errors == t
end

"""
    decryptWord!(cyphertext, G, GInverse, P, goppaPolynomial, codeSupport, m, t)

Decodes and decrypts a 1xn array of bits. It expects that P is the
invers of its respective matrix.
"""
function decryptWord!(cyphertext, G, GInverse, P, goppaPolynomial, codeSupport, m, t)
    permutedcypher = (transpose(cyphertext) * P) .% 2
    codeword = decode!(permutedcypher, goppaPolynomial, codeSupport, m, t)
    plaintext = (codeword * GInverse) .% 2

    if !isPlaintextValid(plaintext, cyphertext, G, t)
        error("The message is invalid.")
    end

    return plaintext
end

"""
    decrypt(msg, G, GInverse, P, goppaPolynomial, codeSupport, m, t)

Decrypts an array of codewords (msg) using a private key, and
returns the plaintext as a string.
"""
function decrypt(msg, G, GInverse, P, goppaPolynomial, codeSupport, m, t)
    plainbits = [decryptWord!(word, G, GInverse, P, goppaPolynomial, codeSupport, m, t) for word in msg]
    plainbits = [bit for bit in Iterators.flatten(plainbits)]
    return bits2string(plainbits)
end

end # Module McEliece

max_stars_repo_name: rexemin/McEliece

max_issues_repo_name: rexemin/McEliece

max_forks_repo_name: rexemin/McEliece


text:
# This file is a part of DataDeps.jl. License is MIT.

"""
    `datadep"Name"` or `datadep"Name/file"`

Use this just like you would a file path, except that you can refer by name to the datadep.
The name alone will resolve to the corresponding folder.
Even if that means it has to be downloaded first.
Adding a path within it functions as expected.
"""
macro datadep_str(namepath)
    :(resolve($(esc(namepath)), @__FILE__))
end

"""
    resolve(datadep, inner_filepath, calling_filepath)

Returns a path to the folder containing the datadep.
Even if that means downloading the dependancy and putting it in there.

     - `inner_filepath` is the path to the file within the data dir
     - `calling_filepath` is a path to the file where this is being invoked from

This is basically the function the lives behind the string macro `datadep"DepName/inner_filepath"`.
"""
function resolve(datadep::AbstractDataDep, inner_filepath, calling_filepath)::String
    while true
        dirpath = _resolve(datadep, calling_filepath)
        filepath = joinpath(dirpath, inner_filepath)

        if can_read_file(filepath)
            return realpath(filepath) # resolve any symlinks for maximum compatibility with external applications
        else # Something has gone wrong
            @warn("DataDep $(datadep.name) found at \"$(dirpath)\". But could not read file at \"$(filepath)\".")
            println("Something has gone wrong. What would you like to do?")
            input_choice(
                ('A', "Abort -- this will error out",
                    ()->abort("Aborted resolving data dependency, program could not continue.")
                ),
                ('R', "Retry -- do this after fixing the problem outside of this script",
                    ()->nothing  # nothing to do
                ),
                ('X', "Remove directory and retry  -- will retrigger download if there isn't another copy elsewhere",
                    ()->rm(dirpath, force=true, recursive=true)
                )
            )
        end
    end
end

function resolve(datadep_name::AbstractString, inner_filepath, calling_filepath)::String
    resolve(registry[datadep_name], inner_filepath, calling_filepath)
end

"""
    resolve("name/path", @__FILE__)

Is the function that lives directly behind the `datadep"name/path"` macro.
If you are working the the names of the datadeps programatically,
and don't want to download them by mistake;
it can be easier to work with this function.

Note though that you must include `@__FILE__` as the second argument,
as DataDeps.jl uses this to allow reading the package specific `deps/data` directory.
Advanced usage could specify a different file or `nothing`, but at that point you are on your own.
"""
function resolve(namepath::AbstractString, calling_filepath)
    parts = splitpath(namepath)
    name = first(parts)
    inner_path = length(parts) > 1 ? joinpath(Iterators.drop(parts, 1)...) : ""
    resolve(name, inner_path, calling_filepath)
end

"The core of the resolve function without any user friendly file stuff, returns the directory"
function _resolve(datadep::AbstractDataDep, calling_filepath)::String
    lp = try_determine_load_path(datadep.name, calling_filepath)
    if lp != nothing
        lp
    else
        handle_missing(datadep, calling_filepath)
    end
end


text:
"""
    MonteCarloZoo

A broad collection of Monte Carlo algorithms implemented in Julia.
"""
module MonteCarloZoo

abstract type Sampler end
export sample

include("Diagnostics.jl")
include("Filtering.jl")
include("Generators.jl")
include("MarkovChain.jl")
include("Simple.jl")
include("Transformation.jl")
include("Utils.jl")

end  # module


text:
# --------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120--------130
function get_fluence(time::Time, pulse::Pulse)
	env = get_env(time.t, pulse)
	dt = au2as(time.dt) * 1e-18 # dt in second
	flux = pulse.I * env.^2
	fluence = copy(flux)
	fluence[1] = flux[1]
	fluence[2] = (flux[1] + flux[2])
	for it in 3:length(fluence)
		fluence[it] = fluence[it - 1] + (flux[it - 1] + flux[it]) / 2
	end
	return fluence .* dt
end


text:
"""
Generate a wavelet iteratively adjusted amplitude adjusted Fourier transform (WIAAFT) [1]
surrogate series for `ts`. Preserves mean and variance structure of the signal, but
randomises nonlinear properties of the signal (i.e. Hurst exponents) [1].

[1] C.J. Keylock (2006). "Constrained surrogate time series with preservation of the mean
and variance structure". Phys. Rev. E. 73: 036707. doi:10.1103/PhysRevE.73.036707.

"""
function wiaaft(ts; n_maxiter = 100, tol = 1e-5, n_windows = 50)

    # Find the maximum number of levels possible with these data
    L = maxtransformlevels(ts)

    # Create a Daubechies wavelet with 16 vanishing moments and periodic boundary conditions
    wt = wavelet(WT.Daubechies{16}(), WT.Periodic)

    # Do the discrete wavelet transform across levels 1:L
    t = [dwt(ts, wt, i) for i = 1:L]

    # Constrained realization of the coefficients at each level by
    # iAAFT algorithm (treating the coefficients as a time series)
    coeff_surrogates = [iaaft(t[i]) for i in 1:L]

    # Mirror the surrogates
    mirror_surrogates = [reverse(x) for x in coeff_surrogates]

    n = length(ts)
    surrogate_coeffs = Vector{Float64}(undef, n)
    mirror_coeffs = Vector{Float64}(undef, n)
    selected_coeffs = Vector{Float64}(undef, n)
    surrogate = Vector{Float64}(undef, n)
    surrogates = Vector{Tuple{Int, Vector{Float64}, Vector{Float64}, Vector{Float64}}}(undef, 0)
    for dyadic_scale = 1:L
        surrogate_coeffs[:] = iaaft(coeff_surrogates[dyadic_scale])
        mirror_coeffs[:] = reverse(surrogate_coeffs)

        # Circularly rotate to minimize error function
        min_errors = Vector{Float64}(undef, 2)
        [min_errors[i] = Inf for i in 1:2]

        errors = Vector{Float64}(undef, 2)
        shift_minimising_errfunc = Vector{Int}(undef, 2)
        for i in 0:n-1
            errors[1] = rmsd(coeff_surrogates[dyadic_scale],
                            circshift(surrogate_coeffs, i))
            errors[2] = rmsd(coeff_surrogates[dyadic_scale],
                            circshift(mirror_coeffs, i))

            if errors[1] < min_errors[1]
                min_errors[1] = errors[1]
                shift_minimising_errfunc[1] = i
            end

            if errors[2] < min_errors[2]
                min_errors[2] = errors[2]
                shift_minimising_errfunc[2] = i
            end
        end
        if min_errors[1] < min_errors[2]
            selected_coeffs[:] = circshift(surrogate_coeffs, shift_minimising_errfunc[1])
        else
            selected_coeffs[:] = circshift(mirror_coeffs, shift_minimising_errfunc[2])
        end

        # Invert the DWT on original
        inverted =  idwt(coeff_surrogates[dyadic_scale], wt, dyadic_scale)
        inverted_surrogate = idwt(selected_coeffs, wt, dyadic_scale)

        surrogate[sortperm(inverted_surrogate)] = sort(inverted)
        surrogate[sortperm(inverted)] = sort(ts)

        push!(surrogates, (dyadic_scale, inverted, inverted_surrogate, surrogate[:]))
    end

    return surrogates
end

export wiaaft


text:
# parameters

# wrapper of C types defined in include/bayesopt/parameters.h

@enum learning_type L_FIXED L_EMPIRICAL L_DISCRETE L_MCMC L_ERROR=-1
@enum score_type SC_MTL SC_ML SC_MAP SC_LOOCV SC_ERROR=-1

immutable kernel_parameters
    name::Cstring
    hp_mean::NTuple{128,Cdouble}
    hp_std::NTuple{128,Cdouble}
    n_hp::Csize_t
end

immutable mean_parameters
    name::Cstring
    coef_mean::NTuple{128,Cdouble}
    coef_std::NTuple{128,Cdouble}
    n_coef::Csize_t
end

type bopt_params
    n_iterations::Csize_t
    n_inner_iterations::Csize_t
    n_init_samples::Csize_t
    n_iter_relearn::Csize_t
    init_method::Csize_t
    random_seed::Cint
    verbose_level::Cint
    log_filename::Cstring
    load_save_flag::Csize_t
    load_filename::Cstring
    save_filename::Cstring
    surr_name::Cstring
    sigma_s::Cdouble
    noise::Cdouble
    alpha::Cdouble
    beta::Cdouble
    sc_type::score_type
    l_type::learning_type
    l_all::Cint
    epsilon::Cdouble
    force_jump::Csize_t
    kernel::kernel_parameters
    mean::mean_parameters
    crit_name::Cstring
    crit_params::NTuple{128,Cdouble}
    n_crit_params::Csize_t
end

function Base.show(io::IO, params::bopt_params)
    print(io, "BayesOpt parameters:")
    width = maximum([length(string(f)) for f in fieldnames(bopt_params)])
    for f in fieldnames(bopt_params)
        val = params.(f)
        if f == :sc_type
            val = score2str(val)
        elseif f == :l_type
            val = learn2str(val)
        elseif f in (:kernel, :mean)
            val = val.name
        end
        if isa(val, Number)
            val = string(val)
        elseif isa(val, Cstring)
            val = bytestring(val)
        else
            # skip
            continue
        end
        println(io)
        print(io, "  ", lpad(f, width), " = ", val)
    end
end

function learn2str(name::learning_type)
    return ccall((:learn2str, libbayesopt), Cstring, (learning_type,), name)
end

function score2str(name::score_type)
    return ccall((:score2str, libbayesopt), Cstring, (score_type,), name)
end

# setters
for f in [
      :set_kernel,
      :set_mean,
      :set_criteria,
      :set_surrogate,
      :set_log_file,
      :set_load_file,
      :set_save_file,
      :set_learning,
      :set_score]
    @eval function $f(params::bopt_params, name::String)
        ccall(
            ($(string(f)), libbayesopt),
            Void,
            (Ref{bopt_params}, Cstring),
            params, name
        )
    end
end

function initialize_parameters_to_default()
    ccall(
        (:initialize_parameters_to_default, libbayesopt),
        bopt_params, ())
end


text:
#include("setup_1d.jl")
using Wavelets

wf = waveletfilter(WT.db2)
wl = waveletls(WT.db2)
N = 8*1024;
x0 = rand(N);
L = nscales(N)  # int(log2(N)-2)
tn = 30

wt = wl
TI = false

println("denoise (N=",N,"), ", L, " levels")
f(x0, L, wt,TI) = for i = 1:tn; denoise(x0, wt=wt, TI=TI); end
f(x0, L, wt,TI);
@time f(x0, L, wt,TI);

#@profile f(x0,L,wt,TI);
#Profile.print()


text:
using PyPlot
using PyCall

# need a separate Python helper file for
# centered lower figure as colon operator
# not implemented fully implemented in PyCall
gs = pyimport("matplotlib.gridspec")
pushfirst!(PyVector(pyimport("sys")."path"), joinpath("assets"))
plothelper = pyimport("plothelper")

function plotmodel(modelvect; ymaxG = nothing, ymaxJ = nothing)

    dt = 0.025
    # Create a time only dataset
    dϵ = timeline(t_start = 0.0, t_end = 10.0, step = dt)
    dσ = timeline(t_start = 0.0, t_end = 10.0, step = dt)

    # calculate a strain/stress data by appling a function of time (by defalut a unit step otherwise substitute hstep(amp = 2.))
    dϵ = strainfunction(dϵ, hstep())
    dσ = stressfunction(dσ, hstep()) - stressfunction(dσ, hstep(offset = 5.0))

    colplot = ["red", "#dc7633", "#229954", "#d4ac0d", "blue"]

    fig, spec, ax1, ax2, ax3 = plothelper.get_fig_axes()
    spec.update(wspace = 0.25,
                hspace = 0.25,
                top = 0.96,
                bottom = 0.07,
                left = 0.07,
                right = 0.98)

    for i = 1:1:length(modelvect)
        #Relaxation modulus
        dG = modelpredict(dϵ, modelvect[i])
        ax1.plot(dG.t, dG.σ, color=colplot[i])
        ax1.plot([-5*dt,0,dG.t[1]],[0.0, 0.0,dG.σ[1]], color=colplot[i])

        # Creep modulus
        dJ = modelpredict(dσ, modelvect[i,1])
        ax2.plot(dJ.t, dJ.ϵ, color=colplot[i])
        ax2.plot([-5*dt,0,dJ.t[1]],[0.0, 0.0,dJ.ϵ[1]], color=colplot[i])
    end

    # Plot settings

    ax1.set_xlabel("Time", fontsize = 10)
    ax1.set_ylabel("Stress", fontsize = 10)
    ax1.tick_params("both", labelsize = 9)
    ax1.set_title("Relaxation response", fontsize = 10)
    isnothing(ymaxJ) ? ax1.set_ylim(bottom = -dt, top = ymaxG) : ax1.set_ylim(bottom = -dt)
    ax1.tick_params("both", labelsize = 9)
    ax1.tick_params(axis="x", which="minor", bottom=true)
    ax1.grid(b=true, alpha = 0.2, axis="y", which="major")
    ax1.grid(b=true, alpha = 0.2, axis="x", which="major")

    ax2.set_xlabel("Time", fontsize = 10)
    ax2.set_ylabel("Strain", fontsize = 10)
    ax2.tick_params("both", labelsize = 9)
    ax2.set_title("Creep response", fontsize = 10)
    isnothing(ymaxJ) ? ax2.set_ylim(bottom = -dt, top = ymaxJ) : ax2.set_ylim(bottom = -dt)
    ax2.tick_params("both", labelsize = 9)
    ax2.tick_params(axis="x", which="minor", bottom=true)
    ax2.grid(b=true, alpha = 0.2, axis="y", which="major")
    ax2.grid(b=true, alpha = 0.2, axis="x", which="major")

    dω = frequencyspec(ω_start = 1e-2, ω_end = 1e2, logstep = log10(1e2/1e-2)/30)

    for i = 1:1:length(modelvect)
        # Storage and Loss moduli
        d_Gp = dynamicmodelpredict(dω, modelvect[i])
        ax3.loglog(d_Gp.ω, d_Gp.Gp, color=colplot[i], "-")
        ax3.loglog(d_Gp.ω, d_Gp.Gpp, color=colplot[i], "--")
    end

    ax3.set_xlabel("Frequency", fontsize = 10)
    ax3.set_ylabel("Storage (—) and Loss (- -) moduli", fontsize = 10)
    ax3.tick_params("both", labelsize = 9)
    ax3.grid(b=true, alpha = 0.2, axis="y", which="major")
    ax3.grid(b=true, alpha = 0.2, axis="x", which="major")
    ax3.set_title("Frequency response", fontsize = 10)

    return fig

end


text:
module DirtyReport

using StatsBase

export Report

const uid = ((n)->()->n+=1)(0)

include("JsonBuilder.jl")
include("Data.jl")
include("Layout.jl")
include("Section.jl")
include("Report.jl")
include("render.jl")

end


text:
using ITensors

function main(; Nx = 6,
                Ny = 3,
                U = 4.0,
                t = 1.0)
  N = Nx * Ny

  sweeps = Sweeps(10)
  maxdim!(sweeps, 100, 200, 400, 800, 1600)
  cutoff!(sweeps, 1e-6)
  noise!(sweeps, 1e-6, 1e-7, 1e-8, 0.0)
  @show sweeps

  sites = siteinds("Electron", N;
                   conserve_qns = true)

  lattice = square_lattice(Nx, Ny;
                           yperiodic = true)

  ampo = AutoMPO()
  for b in lattice
    ampo += -t, "Cdagup", b.s1, "Cup", b.s2
    ampo += -t, "Cdagup", b.s2, "Cup", b.s1
    ampo += -t, "Cdagdn", b.s1, "Cdn", b.s2
    ampo += -t, "Cdagdn", b.s2, "Cdn", b.s1
  end
  for n in 1:N
    ampo += U, "Nupdn", n
  end
  H = MPO(ampo,sites)

  # Half filling
  state = [isodd(n) ? "Up" : "Dn" for n in 1:N]

  # Initialize wavefunction to a random MPS
  # of bond-dimension 10 with same quantum
  # numbers as `state`
  psi0 = randomMPS(sites, state)

  energy,psi = dmrg(H, psi0, sweeps)
  @show t, U
  @show flux(psi)
  @show maxlinkdim(psi)
  @show energy

  return
end

main()


text:
using SimpleHypergraphs
"""
  h2txt(h::Hypergraph, fname::AbstractString)

Translate hypergraph to text file. This is in the following form.
```
N M
```
N is number of hypergraph vertices.
M is number of hypereedges.

"""
function h2txt(h::Hypergraph, fname)
  out = open(fname, "w")
  print(out, nhv(h), ' ', nhe(h), '\n')
  for he in 1:nhe(h)
    print(out, length(getvertices(h, he)), ' ')
    for node in keys(getvertices(h, he))
      print(out, node, ' ')
    end
    print(out, '\n')
  end
  close(out)
end

function txt2h(fname)
  f = open(fname)
  (n, m) = parse.(Int, split(readline(f), ' '))
  h = Hypergraph{Int}(n, m)
  he_i = 1
  while(!eof(f))
    line = split(readline(f), ' ')
    if length(line) == 0 continue end
    vᵢ= line[1]
    he = parse.(Int, line[2:end-1])
    h[he, he_i] .= 1
    he_i += 1
  end
  close(f)

  return h
end


text:
abstract type ModelParameter end

mutable struct GammaParameter <: ModelParameter
	prior_a::Float64
	prior_b::Float64
	prior_distribution::Distributions.Gamma

	variational_a::Float64
	variational_b::Float64

	expectation_trace::Vector{Float64}

	var_name::String
	function GammaParameter(a::Float64, b::Float64; variable_name = "GammaParameter")
		gp = new()
		gp.prior_a = a
		gp.prior_b = b
		gp.prior_distribution = Distributions.Gamma(a, 1. / b)

		gp.variational_a = a
		gp.variational_b = b
		gp.expectation_trace = Vector{Float64}(undef, 0)
		gp.var_name = variable_name
		gp
	end
end

const VectorGammaParameter = Vector{GammaParameter}

mutable struct NormalParameter <: ModelParameter
	prior_mean::Float64
	prior_variance::Float64
	prior_distribution::Distributions.Normal

	variational_mean::Float64
	variational_variance::Float64

	expectation_trace::Vector{Float64}

	var_name::String
	function NormalParameter(m_0::Float64, s_0::Float64; variable_name = "NormalParameter")
		np = new()
		np.prior_mean = m_0
		np.prior_variance = s_0
		np.prior_distribution = Distributions.Normal(m_0, s_0)

		np.variational_mean = rand(np.prior_distribution)
		np.variational_variance = s_0
		np.expectation_trace = Vector{Float64}(undef, 0)
		np.var_name = variable_name
		np
	end
end

mutable struct MvNormalParameter <: ModelParameter
	prior_mean::Vector{Float64}
	prior_covariance::Matrix{Float64}
	prior_distribution::Distributions.MvNormal

	variational_mean::Vector{Float64}
	variational_covariance::Matrix{Float64}

	expectation_trace::Vector{Vector{Float64}}

	var_name::String
	function MvNormalParameter(m_0::Float64, s_0::T, D::Int64; variable_name = "MvNormalParameter") where {T<:Union{Float64, Vector{Float64}}}
		np = new()
		np.prior_mean = m_0.*ones(D)
		np.prior_covariance = s_0 .* Matrix(1.0I, D, D)
		np.prior_distribution = Distributions.MvNormal(np.prior_mean, np.prior_covariance)

		np.variational_mean = rand(np.prior_distribution)
		np.variational_covariance = s_0 .* Matrix(1.0I, D, D)
		np.expectation_trace = Vector{Vector{Float64}}()
		np.var_name = variable_name
		np
	end
end

set_variable_name(p::ModelParameter, name::String) = p.var_name = name

function check_params(p::GammaParameter)
	if p.variational_a < 0
		# @warn "paramater negative:" p.var_name  p.variational_a
		p.variational_a = 1e-7
	end
	if p.variational_b < 0
		# @warn "paramater negative:" p.var_name  p.variational_b
		p.variational_b = 1e-7
	end
end

function check_params(p::NormalParameter)
	if p.variational_variance ≤ 0
		# @warn "variance non-positive" p.var_name p.variational_variance
	end
	p.variational_variance = 1e-7
end

function check_params(p::MvNormalParameter)
	if !issymmetric(p.variational_covariance) || !isposdef(p.variational_covariance)
		p.variational_covariance = Array(Hermitian(p.variational_covariance))
	end
end

expected_value(p::GammaParameter) = p.variational_a/p.variational_b
expected_log_value(p::GammaParameter) = digamma(p.variational_a) - log(p.variational_b)

expected_value(p::Union{NormalParameter, MvNormalParameter}) = p.variational_mean
expected_squared_value(p::NormalParameter) = p.variational_mean^2 + p.variational_variance
expected_squared_value(p::MvNormalParameter) = p.variational_mean*p.variational_mean' + p.variational_covariance

function elbo(gp::GammaParameter)
	check_params(gp)
	q_dist = Distributions.Gamma(gp.variational_a, 1. / gp.variational_b)
	compute_elbo(q_dist, gp)
end

function elbo(np::NormalParameter)
	check_params(np)
	# p_dist = Distributions.Normal(np.prior_mean, np.prior_variance)
	q_dist = Distributions.Normal(np.variational_mean, np.variational_variance)
	compute_elbo(q_dist, np)
end

function elbo(np::MvNormalParameter)
	check_params(np)
	# p_dist = Distributions.MvNormal(np.prior_mean, np.prior_covariance)
	q_dist = Distributions.MvNormal(np.variational_mean, np.variational_covariance)
	compute_elbo(q_dist, np)
end

function compute_elbo(q_dist::Distributions.Distribution, param::ModelParameter)
	expected_param = expected_value(param)
	push!(param.expectation_trace, expected_param)
	Distributions.logpdf(param.prior_distribution, expected_param) - Distributions.logpdf(q_dist, expected_param)
end


text:
using Test, bad

# construct the usual prior
pnull  = .2 # response under TAU
pmcr   = pnull + .1 # minimal clinically relevant
prior1 = Beta(mean = .35, sd = .1)  # start with an subjective prior
prior2 = update(prior1, 4, 10)      # update with phase I data
prior3 = .8*prior2 + .2*Beta(1, 1)  # robustify
prior  = prior3 <= min(2*pmcr, 1.0) # restrict to plausible range

p = 0:.01:1

# probability of succes is simply power * Pr[ p >= pmcr ]
pos   =  (1 - cdf(pmcr, prior)) * Power(prior >= pmcr)
# probability of an type one error is avtoer * Pr [p <= pnull ]
potoe = cdf(pnull, prior) * Power(prior <= pnull)
# finally, expected sample size is the usual
ess   = SampleSize(prior)
# overall utility is then goverend by two factors, both on the scale of the
# average per-patient costs in phase II
var"cost failed phase III" = 12000.
var"risk weighted profit of successful phaseIII" = 200.
utility = -var"cost failed phase III"*potoe +
    var"risk weighted profit of successful phaseIII"*pos -
    ess

# first consider the 'standard design'
α, β = .05, .2
problem = Problem(
    minimise(ess),
    Power(prior  | pnull) <= α,
    Power(prior >= pmcr)  >= 1 - β
)

design = optimise(problem; verbosity = 0)
Power(prior >= pmcr)(design), Power(prior | pnull)(design), ess(design), utility(design)

# now, lets relax the power and type one error rate constraints and maximise
# utility instead! Note that we need to manually increase the marginal
# feasible space since the heuristics based on the given constraints are
# now favouring small designs
uproblem = Problem(
    maximise(utility),
    Power(prior  | pnull) <= 4*α,
    Power(prior >= pmcr)  >= 1 - 2*β,
    n1values = 5:35,
    nmax     = 100
)
udesign = optimise(uproblem; verbosity = 0)
Power(prior >= pmcr)(udesign), Power(prior | pnull)(udesign), ess(udesign), utility(udesign)

@test utility(udesign) > utility(design)


text:
int main() {
    if (true) return 0;
}

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
include(joinpath(pwd(), "src/direct_policy_optimization/dpo.jl"))
include(joinpath(@__DIR__, "cartpole_friction.jl"))

# Additive noise model
model_friction = CartpoleFriction(4, 7, 4, 1.0, 0.2, 0.5, 9.81, μ0)

function fd(model::CartpoleFriction, x⁺, x, u, w, h, t)
    midpoint_implicit(model, x⁺, x, u, w, h) - w
end

# Nominal solution
X̄, Ū = unpack(Z̄_friction, prob_friction)
prob_nom = prob_friction.prob

# DPO
N = 2 * model_friction.n
D = 2 * model_friction.d

β = 1.0
δ = 1.0e-2

# initial samples
x1_sample = resample(x1, Diagonal(ones(model_friction.n)), 1.0)

# mean problem
prob_mean = trajectory_optimization(
				model_friction,
				EmptyObjective(),
				T,
				h = h,
				dynamics = false
				)

# sample problems
prob_sample = [trajectory_optimization(
				model_friction,
				EmptyObjective(),
				T,
				h = h,
				xl = state_bounds(model_friction, T, x1 = x1_sample[i])[1],
				xu = state_bounds(model_friction, T, x1 = x1_sample[i])[2],
				ul = ul,
				uu = uu,
				dynamics = false,
				con = con_friction
				) for i = 1:N]

# sample objective
Q = [(t < T ? Diagonal(10.0 * ones(model_friction.n))
    : Diagonal(100.0 * ones(model_friction.n))) for t = 1:T]
R = [Diagonal([1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]) for t = 1:T-1]

obj_sample = sample_objective(Q, R)
policy = linear_feedback(model.n, 1)
dist = disturbances([Diagonal(δ * ones(model_friction.n)) for t = 1:T-1])
sample = sample_params(β, T)

prob_dpo = dpo_problem(
	prob_nom, prob_mean, prob_sample,
	obj_sample,
	policy,
	dist,
	sample)

# TVLQR policy
K = tvlqr(model, X̄, [Ū[t][1:1] for t = 1:T-1],
 	Q, [R[t][1:1, 1:1] for t = 1:T-1], h)

z0_dpo = zeros(prob_dpo.num_var)
z0_dpo[prob_dpo.prob.idx.nom] = pack(X̄, Ū, prob_nom)
z0_dpo[prob_dpo.prob.idx.mean] = pack(X̄, Ū, prob_nom)
for i = 1:N
	z0_dpo[prob_dpo.prob.idx.sample[i]] = pack(X̄, Ū, prob_nom)
end
for j = 1:(N + D)
	z0_dpo[prob_dpo.prob.idx.slack[j]] = vcat(X̄[2:end]...)
end
for t = 1:T-1
	z0_dpo[prob_dpo.prob.idx.policy[prob_dpo.prob.idx.θ[t]]] = vec(copy(K[t]))
end

include("/home/taylor/.julia/dev/SNOPT7/src/SNOPT7.jl")

# Solve
Z = solve(prob_dpo, copy(z0_dpo),
	nlp = :SNOPT7,
	tol = 1.0e-3, c_tol = 1.0e-3,
	# max_iter = 1000)
	time_limit = 600)

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
using DFTShims: ConstantArrays
using AxisArrays

@test size(ConstantArray(1.0, (3, 4))) == (3, 4)
array = AxisArray(ConstantArray(1, (2, 2)), Axis{:spin}((:u, :d)))
@test size(array) == (2, 2)
@test axes(array, 1) == Axis{:spin}((:u, :d))

max_stars_repo_name: mdavezac/DFTShims

max_issues_repo_name: mdavezac/DFTShims

max_forks_repo_name: mdavezac/DFTShims


text:
module GenomicIntersections

using DataFrames, GenomicFeatures, Statistics

export countoverlap, intersecttable, overlappinglocations, peakintersect, markintersection!, annotatecol!, countintersection!

include("intersections.jl")

end


text:
module TestRDS
using Test
using Dates
using DataFrames
using CategoricalArrays
using RData
using TimeZones

@testset "Loading RDS files (version=$ver)" for ver in (2, 3)
    rdata_path = joinpath(dirname(@__FILE__), "data_v$ver")

    @testset "loading/converting basic R types" begin
        df = DataFrame(num = [1.1, 2.2],
                       int = Int32[1, 2],
                       logi = [true, false],
                       chr = ["ab", "c"],
                       factor = categorical(["ab", "c"], compress=true),
                       cplx = ComplexF64[1.1+0.5im, 1.0im])
        rdf = sexp2julia(load(joinpath(rdata_path, "types.rds"), convert=false))
        @test rdf isa DataFrame
        @test eltype.(eachcol(rdf)) == eltype.(eachcol(df))
        @test isequal(rdf, df)

        rdf_ascii = sexp2julia(load(joinpath(rdata_path, "types_ascii.rds"), convert=false))
        @test rdf_ascii isa DataFrame
        @test eltype.(eachcol(rdf_ascii)) == eltype.(eachcol(df))
        @test isequal(rdf_ascii, df)

        rdf_decomp = sexp2julia(load(joinpath(rdata_path, "types_decomp.rds"), convert=false))
        @test rdf_decomp isa DataFrame
        @test eltype.(eachcol(rdf_decomp)) == eltype.(eachcol(df))
        @test isequal(rdf_decomp, df)

        rdf = load(joinpath(rdata_path, "types.rds"))
        @test rdf isa DataFrame
        @test eltype.(eachcol(rdf)) == eltype.(eachcol(df))
        @test isequal(rdf, df)

        rdf_ascii = load(joinpath(rdata_path, "types_ascii.rds"))
        @test rdf_ascii isa DataFrame
        @test eltype.(eachcol(rdf_ascii)) == eltype.(eachcol(df))
        @test isequal(rdf_ascii, df)

        rdf_decomp = load(joinpath(rdata_path, "types_decomp.rds"))
        @test rdf_decomp isa DataFrame
        @test eltype.(eachcol(rdf_decomp)) == eltype.(eachcol(df))
        @test isequal(rdf_decomp, df)
    end

    @testset "Date conversion" begin
        numdates = load(joinpath(rdata_path, "numdates.rds"))
        @test numdates == Date("2017-01-01") + Dates.Day.(1:4)

        numdates_ascii = load(joinpath(rdata_path, "numdates_ascii.rds"))
        @test numdates_ascii == Date("2017-01-01") + Dates.Day.(1:4)

        intdates = load(joinpath(rdata_path, "intdates.rds"))
        @test intdates == Date("2017-01-01") + Dates.Day.(1:4)

        intdates_ascii = load(joinpath(rdata_path, "intdates_ascii.rds"))
        @test intdates_ascii == Date("2017-01-01") + Dates.Day.(1:4)
    end

    @testset "DateTime conversion" begin
        datetimes = load(joinpath(rdata_path, "datetimes.rds"))
        testdts = map(i -> ZonedDateTime(DateTime("2017-01-01T13:23") + Dates.Second(i),
                                 TimeZone("UTC")), 1:4)
        @test datetimes[1] == testdts
        @test datetimes[2] == testdts[1]
        @test datetimes[3] isa DictoVec
        @test datetimes[3].data == testdts
        @test [datetimes[3].index2name[i] for i in 1:length(datetimes[3])] == ["A", "B", "C", "D"]
        @test datetimes[4] isa DictoVec
        @test datetimes[4].data == [testdts[1]]
        @test datetimes[4].index2name[1] == "A"
    end

    @testset "Date and DateTime in a DataFrame" begin
        rdfs = load(joinpath(rdata_path, "datedfs.rds"))
        df = DataFrame(date=map(i -> Date("2017-01-01") + Dates.Day(i), 1:4),
                       datetime=map(i -> ZonedDateTime(DateTime("2017-01-01T13:23") + Dates.Second(i), tz"UTC"), 1:4))
        @test length(rdfs) == 2
        @test rdfs[1] isa DataFrame
        @test rdfs[2] isa DataFrame
        @test eltype.(eachcol(df)) == eltype.(eachcol(rdfs[1]))
        @test eltype.(eachcol(df)) == eltype.(eachcol(rdfs[2]))
        @test isequal(df[1:1, :], rdfs[1])
        @test isequal(df, rdfs[2])
    end

    @testset "NA Date and DateTime conversion" begin
        dates = load(joinpath(rdata_path, "datesNA.rds"))

        testdates = [Date("2017-01-01") + Dates.Day.(1:4); missing]
        @test all(dates[1] .=== testdates)

        testdts = [map(i -> ZonedDateTime(DateTime("2017-01-01T13:23") + Dates.Second(i), tz"UTC"), 1:4);
                   missing]
        @test all(dates[2] .=== testdts)
    end

    @testset "DateTime timezones" begin
        # tz"CST" is not supported by TimeZones.jl
        datetimes = @test_logs (:warn, "Could not determine the timezone of 'CST', treating as 'UTC'") begin
            load(joinpath(rdata_path, "datetimes_tz.rds"))
        end
        # assumes generate_rda.R was generated on system set to PST!
        @test datetimes[1] == ZonedDateTime(DateTime("2017-01-01T21:23"), tz"UTC")
        # should be tz"CST", but gets substituted to tz"UTC"
        # FIXME update the test when CST is supported
        @test datetimes[2] == ZonedDateTime(DateTime("2017-01-01T13:23"), tz"UTC")
        @test datetimes[3] == ZonedDateTime(DateTime("2017-01-01T13:23"), tz"America/Chicago")
    end

end # for ver in ...

end # module TestRDS


text:
using AlgebraicRelations.Presentations
using Catlab

# Initialize workflow object
wf = Presentation()

# Add Products to workflow
Files, Images, NeuralNet,
       Accuracy, Metadata = add_types!(wf, [(:Files, String),
                                               (:Images, String),
                                               (:NeuralNet, String),
                                               (:Accuracy, Real),
                                               (:Metadata, String)]);

# Add Processes to workflow
extract, split_im, train, evaluate = add_processes!(wf, [(:extract, Files, Images),
                                                         (:split_im, Images, Images⊗Images),
                                                         (:train, NeuralNet⊗Images, NeuralNet⊗Metadata),
                                                         (:evaluate, NeuralNet⊗Images, Accuracy⊗Metadata)]);
# Convert to Schema
@present_to_schema TrainDB(wf);
g = draw_schema(wf)

@test wf isa Catlab.Present.Presentation
@test TrainDB() isa Catlab.CategoricalAlgebra.ACSet
@test g isa Catlab.Graphics.Graphviz.Graph


text:
# Add docs/lib to LOAD_PATH so that we could load the Showcase module like a library.
let doclib = abspath(joinpath(@__DIR__, "Showcase", "src"))
    doclib in LOAD_PATH || pushfirst!(LOAD_PATH, doclib)
end

using Documenter, DocStringExtensions
import Showcase

makedocs(
    sitename = "DocStringExtensions.jl",
    modules = [DocStringExtensions, Showcase],
    clean = false,
    pages = Any[
        "Home" => "index.md",
        "Showcase" => "showcase.md",
    ],
    format = Documenter.HTML(
        assets = ["assets/favicon.ico"],
        analytics = "UA-136089579-3",
    ),
)

deploydocs(
    repo = "github.com/JuliaDocs/DocStringExtensions.jl.git",
)


text:
isdefined(Base, :__precompile__) && __precompile__()
module PowerSimulations

#################################################################################
# Exports

# Base Models
export Simulation
export DecisionModel
export EmulationModel
export ProblemResults
export ProblemTemplate
export InitialCondition
export SimulationModels
export SimulationSequence
export SimulationResults

# Network Relevant Exports
export NetworkModel
export StandardPTDFModel
export PTDFPowerModel
export CopperPlatePowerModel
export AreaBalancePowerModel

######## Device Models ########
export DeviceModel
export FixedOutput
######## Service Models ########
export ServiceModel
export RangeReserve
export RampReserve
export StepwiseCostReserve
export NonSpinningReserve
export PIDSmoothACE
export GroupReserve
######## Branch Models ########
export StaticBranch
export StaticBranchBounds
export StaticBranchUnbounded
export HVDCLossless
export HVDCDispatch
export HVDCUnbounded
# export VoltageSourceDC
######## Load Models ########
export StaticPowerLoad
export InterruptiblePowerLoad
export DispatchablePowerLoad
######## Renewable Formulations ########
export RenewableFullDispatch
export RenewableConstantPowerFactor
######## Hydro Formulations ########
export HydroDispatchRunOfRiver
export HydroDispatchReservoirBudget
export HydroDispatchReservoirStorage
export HydroCommitmentRunOfRiver
export HydroCommitmentReservoirBudget
export HydroCommitmentReservoirStorage
export HydroDispatchPumpedStorage

######## Renewable Formulations ########
export BookKeeping
export BatteryAncillaryServices
export EnergyTarget

######## Thermal Formulations ########
export ThermalStandardUnitCommitment
export ThermalBasicUnitCommitment
export ThermalBasicCompactUnitCommitment
export ThermalBasicDispatch
export ThermalStandardDispatch
export ThermalDispatchNoMin
export ThermalMultiStartUnitCommitment
export ThermalCompactUnitCommitment
export ThermalCompactDispatch

###### Regulation Device Formulation #######
export DeviceLimitedRegulation
export ReserveLimitedRegulation

######## Hybrid Formulations ########
export BasicHybridDispatch
export StandardHybridDispatch

# feedforward models
export UpperBoundFeedforward
export LowerBoundFeedforward
export SemiContinuousFeedforward
export EnergyLimitFeedforward
export FixValueFeedforward
export EnergyTargetFeedforward

# InitialConditions chrons
export InterProblemChronology
export IntraProblemChronology

# Initial Conditions Quantities
export DevicePower
export DeviceStatus
export InitialTimeDurationOn
export InitialTimeDurationOff
export InitialEnergyLevel

# operation_models
export GenericOpProblem
export UnitCommitmentProblem
export EconomicDispatchProblem
# export OptimalPowerFlow

# Functions
export build!
## Op Model Exports
export get_initial_conditions
export serialize_problem
export serialize_results
export serialize_optimization_model
## Decision Model Export
export solve!
## Emulation Model Exports
export run!
## Sim Model Exports
export execute!
export get_simulation_model
## Template Exports
export template_economic_dispatch
export template_unit_commitment
export template_agc_reserve_deployment
export EconomicDispatchProblem
export UnitCommitmentProblem
export AGCReserveDeployment
export run_economic_dispatch
export run_unit_commitment
export set_device_model!
export set_service_model!
export set_network_model!
export get_network_formulation
## Results interfaces
export SimulationResultsExport
export ProblemResultsExport
export export_results
export export_realized_results
export export_optimizer_stats
export get_variable_values
export get_dual_values
export get_parameter_values
export get_aux_variable_values
export get_expression_values
export get_timestamps
export get_model_name
export get_decision_problem_results
export get_emulation_problem_results
export get_system
export get_system!
export set_system!
export list_variable_keys
export list_dual_keys
export list_parameter_keys
export list_aux_variable_keys
export list_expression_keys
export list_variable_names
export list_dual_names
export list_parameter_names
export list_aux_variable_names
export list_expression_names
export list_decision_problems
export list_supported_formats
export load_results!
export read_variable
export read_dual
export read_parameter
export read_aux_variable
export read_expression
export read_variables
export read_duals
export read_parameters
export read_aux_variables
export read_expressions
export read_realized_variable
export read_realized_dual
export read_realized_parameter
export read_realized_aux_variable
export read_realized_expression
export read_realized_variables
export read_realized_duals
export read_realized_parameters
export read_realized_aux_variables
export read_realized_expressions
export get_realized_timestamps
export get_problem_base_power
export get_objective_value
export read_optimizer_stats

## Utils Exports
export get_all_constraint_index
export get_all_variable_index
export get_constraint_index
export get_variable_index
export list_recorder_events
export show_recorder_events
export list_simulation_events
export show_simulation_events
export export_realized_results

## Enums
export BuildStatus
export RunStatus

# Variables
export ActivePowerVariable
export ActivePowerInVariable
export ActivePowerOutVariable
export HotStartVariable
export WarmStartVariable
export ColdStartVariable
export EnergyVariable
export EnergyVariableUp
export EnergyVariableDown
export EnergyShortageVariable
export EnergySurplusVariable
export LiftVariable
export OnVariable
export ReactivePowerVariable
export ReservationVariable
export ActivePowerReserveVariable
export ServiceRequirementVariable
export WaterSpillageVariable
export StartVariable
export StopVariable
export SteadyStateFrequencyDeviation
export AreaMismatchVariable
export DeltaActivePowerUpVariable
export DeltaActivePowerDownVariable
export AdditionalDeltaActivePowerUpVariable
export AdditionalDeltaActivePowerDownVariable
export SmoothACE
export SystemBalanceSlackUp
export SystemBalanceSlackDown
export ReserveRequirementSlack
export VoltageMagnitude
export VoltageAngle
export FlowActivePowerVariable
export FlowActivePowerFromToVariable
export FlowActivePowerToFromVariable
export FlowReactivePowerFromToVariable
export FlowReactivePowerToFromVariable
export PowerAboveMinimumVariable

# Auxiliary variables
export TimeDurationOn
export TimeDurationOff
export PowerOutput
export EnergyOutput

# Constraints
export AbsoluteValueConstraint
export ActiveConstraint
export ActivePowerVariableLimitsConstraint
export ActivePowerVariableTimeSeriesLimitsConstraint
export ActiveRangeConstraint
export ActiveRangeICConstraint
export AreaDispatchBalanceConstraint
export AreaParticipationAssignmentConstraint
export BalanceAuxConstraint
export CommitmentConstraint
export CopperPlateBalanceConstraint
export DurationConstraint
export EnergyBalanceConstraint
export EnergyBudgetConstraint
export EnergyCapacityConstraint
export EnergyCapacityDownConstraint
export EnergyCapacityUpConstraint
export EnergyLimitConstraint
export EnergyShortageVariableLimitsConstraint
export EnergyTargetConstraint
export EqualityConstraint
export FeedforwardSemiContinousConstraint
export FeedforwardUpperBoundConstraint
export FeedforwardLowerBoundConstraint
export FeedforwardIntegralLimitConstraint
export FlowActivePowerConstraint
export FlowActivePowerFromToConstraint
export FlowActivePowerToFromConstraint
export FlowLimitConstraint
export FlowLimitFromToConstraint
export FlowLimitToFromConstraint
export FlowRateConstraint
export FlowRateConstraintFromTo
export FlowRateConstraintToFrom
export FlowReactivePowerConstraint
export FlowReactivePowerFromToConstraint
export FlowReactivePowerToFromConstraint
export FrequencyResponseConstraint
export HVDCPowerBalance
export HVDCTotalPowerDeliveredVariable
export InflowRangeConstraint
export InputActivePowerVariableLimitsConstraint
export InputPowerRangeConstraint
export MustRunConstraint
export NetworkFlowConstraint
export NodalBalanceActiveConstraint
export NodalBalanceReactiveConstraint
export OutputActivePowerVariableLimitsConstraint
export PieceWiseLinearCostConstraint
export PowerOutputRangeConstraint
export ParticipationAssignmentConstraint
export RampConstraint
export RampLimitConstraint
export RangeLimitConstraint
export RateLimitConstraint
export RateLimitConstraintFromTo
export RateLimitConstraintToFrom
export ReactivePowerVariableLimitsConstraint
export ReactiveRangeConstraint
export RegulationLimitsConstraint
export RequirementConstraint
export ReserveEnergyConstraint
export ReservePowerConstraint
export SACEPIDAreaConstraint
export StartTypeConstraint
export StartupInitialConditionConstraint
export StartupTimeLimitTemperatureConstraint

# Parameters
# Time Series Parameters
export ActivePowerTimeSeriesParameter
export ReactivePowerTimeSeriesParameter
export RequirementTimeSeriesParameter
export EnergyTargetTimeSeriesParameter
export EnergyBudgetTimeSeriesParameter

# Feedforward Parameters
export OnStatusParameter
export UpperBoundValueParameter

# Expressions
export SystemBalanceExpressions
export RangeConstraintLBExpressions
export RangeConstraintUBExpressions
export CostExpressions
export ActivePowerBalance
export ReactivePowerBalance
export EmergencyUp
export EmergencyDown
export RawACE
export ProductionCostExpression
export ActivePowerRangeExpressionLB
export ComponentActivePowerRangeExpressionLB
export ReserveRangeExpressionLB
export ActivePowerRangeExpressionUB
export ReserveRangeExpressionUB
export ComponentActivePowerRangeExpressionUB
export ComponentReserveUpBalanceExpression
export ComponentReserveDownBalanceExpression

#################################################################################
# Imports
import DataStructures: OrderedDict, Deque, SortedDict
import Logging
import Serialization
# Modeling Imports
import JuMP
# so that users do not need to import JuMP to use a solver with PowerModels
import JuMP: optimizer_with_attributes
import JuMP.Containers: DenseAxisArray, SparseAxisArray
export optimizer_with_attributes
import MathOptInterface
import ParameterJuMP
import LinearAlgebra
import JSON3
import PowerSystems
import InfrastructureSystems
import InfrastructureSystems: @assert_op, list_recorder_events, get_name
export get_name
export get_model_base_power
export get_optimizer_stats
export get_timestamps
export get_resolution
import PowerModels
import TimerOutputs
import ProgressMeter

# Base Imports
import Base.getindex
import Base.length
import Base.first

# TimeStamp Management Imports
import Dates
import TimeSeries

# I/O Imports
import DataFrames
import JSON
import CSV
import SHA
import HDF5
import PrettyTables

# PowerModels exports
export ACPPowerModel
export ACRPowerModel
export ACTPowerModel
export DCPPowerModel
export NFAPowerModel
export DCPLLPowerModel
export LPACCPowerModel
export SOCWRPowerModel
export SOCWRConicPowerModel
export QCRMPowerModel
export QCLSPowerModel

################################################################################

# Type Alias From other Packages
const PM = PowerModels
const PSY = PowerSystems
const PSI = PowerSimulations
const IS = InfrastructureSystems
const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities
const PJ = ParameterJuMP
const MOPFM = MOI.FileFormats.Model
const TS = TimeSeries

################################################################################

function progress_meter_enabled()
    return isa(stderr, Base.TTY) &&
           (get(ENV, "CI", nothing) != "true") &&
           (get(ENV, "RUNNING_PSI_TESTS", nothing) != "true")
end

using DocStringExtensions

@template DEFAULT = """
                    $(TYPEDSIGNATURES)
                    $(DOCSTRING)
                    $(METHODLIST)
                    """
# Includes

include("core/definitions.jl")

# Core components
include("core/formulations.jl")
include("core/abstract_simulation_store.jl")
include("core/operation_model_abstract_types.jl")
include("core/optimization_container_types.jl")
include("core/abstract_feedforward.jl")
include("core/optimization_container_keys.jl")
include("core/network_model.jl")
include("core/parameters.jl")
include("core/service_model.jl")
include("core/device_model.jl")
include("core/variables.jl")
include("core/auxiliary_variables.jl")
include("core/constraints.jl")
include("core/expressions.jl")
include("core/initial_conditions.jl")
include("core/settings.jl")
include("core/cache_utils.jl")
include("core/optimizer_stats.jl")
include("core/dataset.jl")
include("core/dataset_container.jl")

include("core/optimization_container.jl")
include("core/store_common.jl")

# Order Required
include("initial_conditions/initial_condition_chronologies.jl")

include("operation/problem_template.jl")
include("operation/operation_model_interface.jl")
include("operation/model_store_params.jl")
include("operation/abstract_model_store.jl")
include("operation/decision_model_store.jl")
include("operation/emulation_model_store.jl")
include("operation/initial_conditions_update_in_memory_store.jl")
include("operation/model_internal.jl")
include("operation/decision_model.jl")
include("operation/emulation_model.jl")
include("operation/problem_results_export.jl")
include("operation/problem_results.jl")
include("operation/operation_model_serialization.jl")
include("operation/time_series_interface.jl")
include("operation/optimization_debugging.jl")
include("operation/model_numerical_analysis_utils.jl")

include("initial_conditions/add_initial_condition.jl")
include("initial_conditions/update_initial_conditions.jl")
include("initial_conditions/calculate_initial_condition.jl")

include("parameters/add_parameters.jl")
include("parameters/update_parameters.jl")

include("feedforward/feedforwards.jl")
include("feedforward/feedforward_arguments.jl")
include("feedforward/feedforward_constraints.jl")

include("simulation/optimization_output_cache.jl")
include("simulation/optimization_output_caches.jl")
include("simulation/simulation_models.jl")
include("simulation/simulation_state.jl")
include("simulation/initial_condition_update_simulation.jl")
include("simulation/simulation_store_params.jl")
include("simulation/hdf_simulation_store.jl")
include("simulation/in_memory_simulation_store.jl")
include("simulation/simulation_problem_results.jl")
include("simulation/realized_meta.jl")
include("simulation/decision_model_simulation_results.jl")
include("simulation/emulation_model_simulation_results.jl")
include("simulation/simulation_sequence.jl")
include("simulation/simulation_internal.jl")
include("simulation/simulation.jl")
include("simulation/simulation_results_export.jl")
include("simulation/simulation_results.jl")

include("devices_models/devices/common/objective_functions.jl")
include("devices_models/devices/common/range_constraint.jl")
include("devices_models/devices/common/add_variable.jl")
include("devices_models/devices/common/add_auxiliary_variable.jl")
include("devices_models/devices/common/add_constraint_dual.jl")
include("devices_models/devices/common/rateofchange_constraints.jl")
include("devices_models/devices/common/duration_constraints.jl")
include("devices_models/devices/common/get_time_series.jl")

# Device Modeling components
include("devices_models/devices/interfaces.jl")
include("devices_models/devices/common/add_to_expression.jl")
include("devices_models/devices/common/set_expression.jl")
include("devices_models/devices/renewable_generation.jl")
include("devices_models/devices/thermal_generation.jl")
include("devices_models/devices/electric_loads.jl")
include("devices_models/devices/AC_branches.jl")
include("devices_models/devices/DC_branches.jl")
include("devices_models/devices/storage.jl")
include("devices_models/devices/hydro_generation.jl")
include("devices_models/devices/regulation_device.jl")
include("devices_models/devices/hybrid_generation.jl")

# Services Models
include("services_models/agc.jl")
include("services_models/reserves.jl")
include("services_models/group_reserve.jl")
include("services_models/service_slacks.jl")
include("services_models/services_constructor.jl")

# Network models
include("network_models/copperplate_model.jl")
include("network_models/powermodels_interface.jl")
include("network_models/pm_translator.jl")
include("network_models/network_slack_variables.jl")
include("network_models/area_balance_model.jl")

include("initial_conditions/initialization.jl")

# Device constructors
include("devices_models/device_constructors/constructor_validations.jl")
include("devices_models/device_constructors/thermalgeneration_constructor.jl")
include("devices_models/device_constructors/hydrogeneration_constructor.jl")
include("devices_models/device_constructors/branch_constructor.jl")
include("devices_models/device_constructors/renewablegeneration_constructor.jl")
include("devices_models/device_constructors/load_constructor.jl")
include("devices_models/device_constructors/storage_constructor.jl")
include("devices_models/device_constructors/regulationdevice_constructor.jl")
include("devices_models/device_constructors/hybridgeneration_constructor.jl")
# Network constructors
include("network_models/network_constructor.jl")

# Templates for Operation Problems
include("operation/operation_problem_templates.jl")

# Operations Decision Problems
include("operation/decision_problems.jl")

# Utils
include("utils/printing.jl")
include("utils/file_utils.jl")
include("utils/logging.jl")
include("utils/dataframes_utils.jl")
include("utils/jump_utils.jl")
include("utils/powersystems_utils.jl")
include("utils/recorder_events.jl")
include("utils/datetime_utils.jl")

end


text:
# phylo/metadata.jl
# =================
#
# Types and methods for metadata attached to phylogenies.
#
# Part of the Bio.Phylo module.
#
# This file is a part of BioJulia. License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

abstract BranchMetaData

# Every Metadata type should be immutable, and should have the following
# methods defined:
#
# A no-arg constructor. e.g. BasicBranch()
# A method called new_branchlength, which constructs a new metadata variable,
# from the old metadata variable, and a new branchlength. The new branchlength
# should be of type Nullable{Float64}.
# A method branchlength which returns the branchlength

immutable BasicBranch <: BranchMetaData
    len::Nullable{Float64}
    conf::Nullable{Float64}
end

@inline function BasicBranch()
    return BasicBranch(Nullable{Float64}(), Nullable{Float64}())
end

@inline function new_branchlength(x::BasicBranch, len::Nullable{Float64})
    return BasicBranch(len, x.conf)
end

branchlength(x::BasicBranch) = x.len


text:
using MLStyle
using YaoLocations
using YaoCompiler
using YaoCompiler.Intrinsics
using OpenQASM
using YaoTargetQASM
using YaoTargetQASM.QASM
using Test

@testset "reginfo" begin
    include("reginfo.jl")
end

@testset "codegen" begin
    include("codegen.jl")
end

@testset "frontend" begin
    include("frontend.jl")
end