text:
"""
    AbstractEmpiricalCDF

Concrete types are `EmpiricalCDF` and `EmpiricalCDFHi`.
"""
abstract type AbstractEmpiricalCDF end

struct EmpiricalCDF{T <: Real} <: AbstractEmpiricalCDF
    xdata::Vector{T}
end

Base.:(==)(cdf1::EmpiricalCDF, cdf2::EmpiricalCDF) = cdf1.xdata == cdf2.xdata
Base.copy(cdf::EmpiricalCDF) = EmpiricalCDF(copy(cdf.xdata))

"""
    data(cdf::AbstractEmpiricalCDF)

return the array holding samples for `cdf`.
"""
data(cdf::AbstractEmpiricalCDF) = cdf.xdata

# Extend base functions
for f in (:length, :size, :minimum, :maximum, :extrema, :issorted, :iterate, :getindex, :lastindex, :firstindex)
    @eval begin
        Base.$(f)(cdf::AbstractEmpiricalCDF, args...) = $(f)(cdf.xdata, args...)
    end
end

# Extend Statistics functions
for f in (:mean, :median, :middle, :std, :stdm, :var, :varm, :quantile)
    @eval begin
        Statistics.$(f)(cdf::AbstractEmpiricalCDF, args...; kws...) = Statistics.$(f)(cdf.xdata, args...; kws...)
    end
end

# Same as above, but the return value is the cdf

@doc """
    sort!(cdf::AbstractEmpiricalCDF)

Sort the data collected in `cdf`. You must call `sort!` before using `cdf`.
""" sort!

for f in (:sort!, :empty!)
    @eval begin
        Base.$(f)(cdf::AbstractEmpiricalCDF, args...) = ($(f)(cdf.xdata,args...); cdf)
    end
end

getcdfindex(cdf::AbstractEmpiricalCDF, x::Real) = searchsortedlast(cdf.xdata, x)

# With several tests, this is about the same speed or faster than the routine borrowed from StatsBase
function _getinds(cdf::AbstractEmpiricalCDF, v::AbstractArray{T}) where T <: Real
    r = Array{eltype(cdf.xdata)}(undef, size(v)...)
    for (i, x) in enumerate(v)
        @inbounds r[i] = cdf(x)
    end
    return r
end

(cdf::EmpiricalCDF)(x::Real) = _val_at_index(cdf, getcdfindex(cdf, x))

(cdf::EmpiricalCDF)(v::AbstractArray) = _getinds(cdf, v)

# This is not done by getindex, but rather by calling the object.
#Base.getindex(cdf::AbstractEmpiricalCDF, v::AbstractArray{T}) where {T <: Real} = _getinds(cdf,v)

"""
    EmpiricalCDF{T=Float64}()

Construct an empirical CDF. After inserting elements with `push!` or
`append!`, and before using most of the functions below, the CDF must be sorted with `sort!`.

`EmpiricalCDF` and `EmpiricalCDFHi` are callable objects. For `cdf::AbstractEmpiricalCDF`,
`cdf(x)` returns the estimate of the CDF at `x`. By contrast, `cdf[inds]` indexes into the underlying
data array.

```julia-repl
julia> cdf = EmpiricalCDF();
julia> append!(cdf,randn(10^6));
julia> sort!(cdf);
julia> cdf(0.0)
0.499876
julia> cdf(1.0)
0.840944
julia> cdf(-1.0)
0.158494
```
In this example, we collected \$10^6\$ samples from the unit normal distribution. About half of the
samples are greater than zero. Approximately the same mass is between zero and one as is between
zero and minus one.
"""
EmpiricalCDF{T}() where T = EmpiricalCDF{T}(Array{T}(undef,0))
EmpiricalCDF() = EmpiricalCDF{Float64}()

"""
    EmpiricalCDFHi{T <: Real} <: AbstractEmpiricalCDF

Empirical CDF with lower cutoff. That is, keep only the tail.
"""
struct EmpiricalCDFHi{T <: Real} <: AbstractEmpiricalCDF
    lowreject::T  # reject counts smaller than this
    rejectcounts::Array{Int,1}  # how many rejections have we done
    xdata::Array{T,1}  # death times
end

(cdf::EmpiricalCDFHi)(x::Real) = _val_at_index(cdf, getcdfindex(cdf, x))

Base.:(==)(cdf1::EmpiricalCDFHi, cdf2::EmpiricalCDFHi) = cdf1.lowreject == cdf2.lowreject &&
    cdf1.rejectcounts == cdf2.rejectcounts && cdf1.xdata == cdf2.xdata

Base.copy(cdf::EmpiricalCDFHi) = EmpiricalCDFHi(cdf.lowreject, copy(cdf.rejectcounts), copy(cdf.xdata))

"""
    EmpiricalCDFHi(lowreject::Real)

Return a CDF that does not store samples whose values are less than `lowreject`. The
sorted CDF will still be properly normalized.
"""
function EmpiricalCDFHi(lowreject::Real)
    cdf = EmpiricalCDFHi(lowreject, Array{Int}(undef,1), Array{typeof(lowreject)}(undef,0))
    cdf.rejectcounts[1] = 0
    cdf
end

function Base.empty!(cdfhi::EmpiricalCDFHi)
    empty!(cdfhi.xdata)
    cdf.rejectcounts[1] = 0
    return cdfhi
end

"""
    EmpiricalCDF(lowreject::Real)

If `lowreject` is finite return `EmpiricalCDFHi(lowreject)`. Otherwise
return `EmpiricalCDF()`.
"""
function EmpiricalCDF(lowreject::Real)
    if isfinite(lowreject)
        EmpiricalCDFHi(lowreject)
    else
        EmpiricalCDF()
    end
end

"""
    counts(cdf::AbstractEmpiricalCDF)

Return the number of counts added to `cdf`. This includes counts
that may have been discarded because they are below of the cutoff.
"""
counts(cdf::AbstractEmpiricalCDF) = _total_counts(cdf)

# Number of points accepted plus number of points rejected
_total_counts(cdf::EmpiricalCDF) = length(cdf.xdata)
_total_counts(cdf::EmpiricalCDFHi) = length(cdf.xdata) + _rejectcounts(cdf)

# Value of the cdf corresponding to index i
_val_at_index(cdf::EmpiricalCDF, i) =  i/_total_counts(cdf)
_val_at_index(cdf::EmpiricalCDFHi, i) =  (i+ _rejectcounts(cdf) )/_total_counts(cdf)

_increment_rejectcounts(cdf::EmpiricalCDFHi) = cdf.rejectcounts[1] += 1

_addto_rejectcounts(cdf::EmpiricalCDFHi, n) = cdf.rejectcounts[1] += n

_rejectcounts(cdf::EmpiricalCDFHi) = cdf.rejectcounts[1]

"""
    push!(cdf::EmpiricalCDF,x::Real)

add the sample `x` to `cdf`.
"""
Base.push!(cdf::EmpiricalCDF,x::Real) = (push!(cdf.xdata,x) ; cdf)

function _push!(cdf::EmpiricalCDFHi,x::Real)
    if  x >= cdf.lowreject
        push!(cdf.xdata,x)
    else
        _increment_rejectcounts(cdf)
    end
end

function Base.push!(cdf::EmpiricalCDFHi, x::Real)
    if isinf(cdf.lowreject)
        push!(cdf.xdata,x)
    else
        _push!(cdf,x)
    end
    cdf
end

Base.append!(cdf_to::EmpiricalCDF, cdf_from::AbstractEmpiricalCDF) = append!(cdf_to, cdf_from.xdata)

function Base.append!(cdf_to::EmpiricalCDFHi, cdf_from::AbstractEmpiricalCDF)
    append!(cdf_to, cdf_from.xdata)
    _addto_rejectcounts(cdf_to, _rejectcounts(cdf_from))
    return cdf_to
end

"""
    append!(cdf::EmpiricalCDF, a::AbstractArray)
    append!(cdf_to::AbstractEmpiricalCDF, cdf_from::AbstractEmpiricalCDF)

Add samples in `a` to `cdf`. Add samples in `cdf_from` to `cdf_to`.
"""
Base.append!(cdf::EmpiricalCDF, a::AbstractArray) =  (append!(cdf.xdata, a); cdf)

function Base.append!(cdf::EmpiricalCDFHi, times)
    if  isinf(cdf.lowreject)
        append!(cdf.xdata, times)
    else
        for x in times
            _push!(cdf,x)
        end
    end
    cdf
end

function _inverse(cdf::EmpiricalCDF, x)
    ind = floor(Int, _total_counts(cdf)*x) + 1
    cdf.xdata[ind]
end

function _inverse(cdf::EmpiricalCDFHi, x)
    ind = floor(Int, _total_counts(cdf)*x) + 1 - _rejectcounts(cdf)
    cdf.xdata[ind]
end

"""
    rand(cdf::EmpiricalCDF)

Pick a random sample from the distribution represented by `cdf`.
"""
Base.rand(cdf::AbstractEmpiricalCDF) = _inverse(cdf,rand())

"""
    finv(cdf::AbstractEmpiricalCDF) --> Function

Return the quantile function, that is, the functional inverse of `cdf`.
`cdf` is a callable object. Note that finv differs slightly from `quantile`.

### Examples
Here, `cdf` contains \$10^6\$ samples from the unit normal distribution.
```julia-repl
julia> icdf = finv(cdf);
julia> icdf(.5)
-0.00037235611091389375
julia> icdf(1.0-eps())
4.601393290425543
julia> maximum(cdf)
4.601393290425543
```
"""
function finv(cdf::EmpiricalCDF)
    function (c::Real)
        (c < 0 || c >= 1) && throw(DomainError())
        _inverse(cdf,c)
    end
end

function finv(cdf::EmpiricalCDFHi)
    function (c::Real)
        (c < cdf.lowreject || c >= 1) && throw(DomainError())
        _inverse(cdf,c)
    end
end

function Base.show(io::IO, cdf::EmpiricalCDF)
    print(io, string(typeof(cdf)))
    print(io, "(n=")
    print(io,length(cdf))
    print(io,')')
end

function Base.show(io::IO, cdf::EmpiricalCDFHi)
    print(io, string(typeof(cdf)))
    print(io, "(n=")
    print(io,length(cdf))
    print(io,",lowreject=")
    print(io,cdf.lowreject)
    print(io,')')
end

function Base.print(cdf::AbstractEmpiricalCDF,fn::String)
    io = open(fn,"w")
    print(io,cdf)
    close(io)
end

@doc """
    print(io::IO, cdf::AbstractEmpiricalCDF)

Call logprint(io,cdf)
""" print

Base.print(io::IO, cdf::AbstractEmpiricalCDF) = logprint(io,cdf)

"`linprint(io::IO ,cdf::AbstractEmpiricalCDF, n=2000)` print (not more than) `n` linearly spaced points after sorting the data."
linprint(io::IO, cdf::AbstractEmpiricalCDF,nprint_pts=2000; lastpt=false) = _printcdf(io,cdf,false,nprint_pts, lastpt=lastpt)
#linprint(io::IO, cdf::AbstractEmpiricalCDF, nprint_pts::Integer) = _printcdf(io,cdf,false, nprint_pts)

"`logprint(io::IO, cdf::EmpiricalCDF, n=2000)` print (not more than) `n` log spaced points after sorting the data."
logprint(io::IO, cdf::AbstractEmpiricalCDF, nprint_pts::Integer) = _printcdf(io,cdf,true, nprint_pts)
logprint(io::IO, cdf::AbstractEmpiricalCDF) = _printcdf(io,cdf,true, 2000)

Base.print(io::IO, cdf::AbstractEmpiricalCDF, nprint_pts::Integer) = logprint(io,cdf,nprint_pts)

Base.print(io::IO, cdf::AbstractEmpiricalCDF, prpts::AbstractArray) = _printcdf(io,cdf,prpts)

function Base.print(cdf::AbstractEmpiricalCDF, prpts::AbstractArray, fname::String)
    io = open(fname,"w")
    print(io,cdf,prpts)
    close(io)
end

"""
    linprint(fn::String, cdf::AbstractEmpiricalCDF, n=2000)

print `cdf` to file `fn`. Print no more than `n` linearly spaced points.
"""
function linprint(fn::String, cdf::AbstractEmpiricalCDF, n=2000; lastpt=false)
    io = open(fn,"w")
    linprint(io,cdf,n; lastpt=lastpt)
    close(io)
end

function logprint(fn::String, cdf::AbstractEmpiricalCDF, n=2000)
    io = open(fn,"w")
    logprint(io,cdf,n)
    close(io)
end

# print many fields, prpts is iterable
function _printfull(io, cdf::AbstractEmpiricalCDF, prpts; lastpt=false)
    n = length(cdf)
    println(io, "# log10(t)  log10(P(t))  log10(1-P(t))  t  1-P(t)   P(t)")
    (p, state) = iterate(prpts)
    ulim = lastpt ? n : n - 1 # usually don't print ordinate value of 1
    for i in 1:ulim
        @inbounds xp = cdf.xdata[i]
        next = iterate(prpts, state)
        if next == nothing
            break
        end
        if xp >= p
            (p, state) = next
            cdf_val = _val_at_index(cdf,i)
            Printf.@printf(io,"%e\t%e\t%e\t%e\t%e\t%e\n", log10(abs(xp)),  log10(abs(1-cdf_val)), log10(abs(cdf_val)), xp, cdf_val, 1-cdf_val)
        end
    end
end

function _printextraheader(io,cdf::AbstractEmpiricalCDF)
end
function _printextraheader(io,cdf::EmpiricalCDFHi)
    nreject = _rejectcounts(cdf)
    n = length(cdf)
    println(io, "# cdf: lowreject = ", cdf.lowreject)
    println(io, "# cdf: lowreject counts = ", nreject)
    Printf.@printf(io, "# cdf: fraction kept = %f\n", n/(n+nreject))
end

# Note that we sort in place
function _printcdf(io::IO, cdf::AbstractEmpiricalCDF, logprint::Bool, nprint_pts::Integer; lastpt=false)
    if length(cdf) == 0
        error("Trying to print empty cdf")
    end
    x = cdf.xdata
    sort!(x)
    xmin = x[1]
    xmax = x[end]
    local prpts
    if logprint && xmin > 0
        println(io, "# cdf: log spacing of coordinate")
        prpts = 10 .^ range(log10(xmin), stop=log10(xmax), length=nprint_pts)
    else
        println(io, "# cdf: linear spacing of coordinate")
        prpts = range(xmin, stop=xmax, length=nprint_pts)
    end
    _printcdf(io, cdf, prpts; lastpt=lastpt)
end

# FIXME: We need to avoid resorting
function _printcdf(io::IO, cdf::AbstractEmpiricalCDF, prpts::AbstractArray; lastpt=false)
    x = cdf.xdata
    sort!(x)
    n = length(x)
    if n == 0
        error("Trying to print empty cdf")
    end
    println(io, "# cdf of survival times")
    println(io, "# cdf: total     counts = ", _total_counts(cdf))
    _printextraheader(io,cdf)
    println(io, "# cdf: number of points in cdf: $n")
    _printfull(io,cdf,prpts; lastpt=lastpt)
end


text:
# Create example data
x = repeat(LinRange(-2, 2, 20), inner=10)
y = repeat(LinRange(0, pi, 10), outer=20)
z = sin.(x) .+ cos.(y)
u = 0.1ones(200)
v = zeros(200)
w = 0.5z
# Plot vectors
quiver3(x, y, z, u, v, w, "o", markersize=0.5)

ext: jl


text:
# Generate julialatex.el file
using LaTeX_Entities
const LE = LaTeX_Entities

open("julialatex.el", "w") do io
    def = LE.default
    for nam in def.nam
        println(io, "(puthash \"", word, "\" \"", LE.lookupname(def, nam), "\" julia-latexsubs)")
    end
end


text:
@testset "Smoothing" begin
  @testset "TaubinSmoothing" begin
    # TODO
  end
end


text:
using YaoBase, YaoArrayRegister
export GeneralMatrixBlock, matblock

"""
    GeneralMatrixBlock{D, MT} <: PrimitiveBlock{D}

General matrix gate wraps a matrix operator to quantum gates. This is the most
general form of a quantum gate.
"""
struct GeneralMatrixBlock{D,T,MT<:AbstractMatrix{T}} <: PrimitiveBlock{D}
    m::Int
    n::Int
    mat::MT

    function GeneralMatrixBlock{D}(m::Int, n::Int, A::MT) where {D,T,MT<:AbstractMatrix{T}}
        (D^m, D^n) == size(A) ||
            throw(DimensionMismatch("expect a $(D^m) x $(D^n) matrix, got $(size(A))"))

        return new{D,T,MT}(m, n, A)
    end
end
YaoBase.nqudits(m::GeneralMatrixBlock) = m.n

GeneralMatrixBlock(m::AbstractMatrix; nlevel=2) = GeneralMatrixBlock{nlevel}(logdi.(size(m), nlevel)..., m)

"""
    matblock(m::AbstractMatrix; nlevel=2)

Create a [`GeneralMatrixBlock`](@ref) with a matrix `m`.

# Example

```jldoctest; setup=:(using YaoBlocks)
julia> matblock(ComplexF64[0 1;1 0])
matblock(...)
```

!!!warn

    Instead of converting it to the default data type `ComplexF64`,
    this will return its contained matrix when calling `mat`.
"""
matblock(m::AbstractMatrix; nlevel=2) = GeneralMatrixBlock(m; nlevel=nlevel)

"""
    matblock(m::AbstractMatrix; nlevel=2)

Create a [`GeneralMatrixBlock`](@ref) with a matrix `m`.
"""
matblock(m::AbstractBlock{D}) where {D} = GeneralMatrixBlock(mat(m); nlevel=D)

cache_key(m::GeneralMatrixBlock) = hash(m.mat)

"""
    mat(A::GeneralMatrixBlock)

Return the matrix of general matrix block.

!!!warn

    Instead of converting it to the default data type `ComplexF64`,
    this will return its contained matrix.
"""
mat(A::GeneralMatrixBlock) = A.mat

function mat(::Type{T}, A::GeneralMatrixBlock) where {T}
    if eltype(A.mat) == T
        return A.mat
    else
        # this errors before, but since we allow one to specify T in mat
        # this should be allowed but with a suggestion
        @warn "converting $(eltype(A.mat)) to eltype $T, consider create another matblock with eltype $T"
        return copyto!(similar(A.mat, T), A.mat)
    end
end

Base.:(==)(A::GeneralMatrixBlock, B::GeneralMatrixBlock) = A.mat == B.mat
Base.copy(A::GeneralMatrixBlock) = GeneralMatrixBlock(copy(A.mat))
Base.adjoint(x::GeneralMatrixBlock) = Daggered(x)


text:
##############
# Issue type #
##############

@ghdef mutable struct Issue
    id::Union{Int, Nothing}
    number::Union{Int, Nothing}
    comments::Union{Int, Nothing}
    title::Union{String, Nothing}
    state::Union{String, Nothing}
    body::Union{String, Nothing}
    user::Union{Owner, Nothing}
    assignee::Union{Owner, Nothing}
    closed_by::Union{Owner, Nothing}
    created_at::Union{Dates.DateTime, Nothing}
    updated_at::Union{Dates.DateTime, Nothing}
    closed_at::Union{Dates.DateTime, Nothing}
    labels::Union{Vector{Dict}, Nothing}
    milestone::Union{Dict, Nothing}
    pull_request::Union{PullRequest, Nothing}
    url::Union{HTTP.URI, Nothing}
    html_url::Union{HTTP.URI, Nothing}
    labels_url::Union{HTTP.URI, Nothing}
    comments_url::Union{HTTP.URI, Nothing}
    events_url::Union{HTTP.URI, Nothing}
    locked::Union{Bool, Nothing}
end

Issue(number::Real) = Issue(Dict("number" => number))

namefield(issue::Issue) = issue.number

###############
# API Methods #
###############

@api_default function issue(api::GitHubAPI, repo, issue_obj; options...)
    result = gh_get_json(api, "/repos/$(name(repo))/issues/$(name(issue_obj))"; options...)
    return Issue(result)
end

@api_default function issues(api::GitHubAPI, repo; options...)
    results, page_data = gh_get_paged_json(api, "/repos/$(name(repo))/issues"; options...)
    return map(Issue, results), page_data
end

@api_default function create_issue(api::GitHubAPI, repo; options...)
    result = gh_post_json(api, "/repos/$(name(repo))/issues"; options...)
    return Issue(result)
end

@api_default function edit_issue(api::GitHubAPI, repo, issue; options...)
    result = gh_patch_json(api, "/repos/$(name(repo))/issues/$(name(issue))"; options...)
    return Issue(result)
end


text:
"""
Modia module with electric component models (inspired from Modelica Standard Library).

* Developer: Hilding Elmqvist, Mogram AB
* Copyright (c) 2016-2018: Hilding Elmqvist, Toivo Henningsson, Martin Otter
* License: MIT (expat)

The building blocks for electric components are:

- `Pin` - The main connector representing an electrical node with variables `v` and `i`.
- `OnePort` - Base model for an electric device with two `Pin`s with variables
    `v`, `i`, `p`, and `n`.

The following functions define variables with appropriate units:

- `Voltage()` - Electric potential
- `Current()` - The main flow quantity
- `Resistance()`
- `Capacitance()`

"""
module Electric

using ..Blocks
using ..Synchronous: positive
using Unitful

using Modia

export Pin, Ground, OnePort, Resistor, Capacitor, Inductor,
  ConstantVoltage, StepVoltage, SignalVoltage, SineVoltage, IdealOpAmp3Pin, IdealDiode,
  Voltage, Current, Resistance, Capacitance

"Electric potential, volts"
Voltage(; args...) = Variable(;T=Unitful.V, size=(), start = 0.0, info = "Voltage", args...)
"Electric current, amperes"
Current(; args...) = Variable(;T=Unitful.A, size=(), start = 0.0, info = "Current", args...)
"Electric resistance, ohms"
Resistance(; args...) = Variable(;T=Unitful.Ω, size=(), args...)
"Electric capacitance, farads"
Capacitance(; args...) = Variable(;T=Unitful.F, size=(), args...)

@model Pin1 begin
    v = Float()
    i = Float(flow=true)
end

"An electric node for connections"
@model Pin begin
    v = Voltage(info = "Potential of the pin")
    i = Current(info = "Current into the pin", flow = true)
end

"Grounded `Pin` with zero voltage"
@model Ground begin
    p = Pin(info = "Grounded pin (zero volts)")
    @equations begin
        p.v = 0
    end
end

@model OnePort1 begin
    v = Float()
    i = Float()
    p = Pin()
    n = Pin()
    @equations begin
        v = p.v - n.v
        0 = p.i + n.i
        i = p.i
    end
end

"Base model for an electric device with two `Pin`s"
@model OnePort begin
    v = Voltage(info = "Voltage between `p` and `n`")
    i = Current(info = "Current from `p` to `n`")
    p = Pin(info = "Positive pin")
    n = Pin(info = "Negative pin")
    @equations begin
        v = p.v - n.v
        0 = p.i + n.i
        i = p.i
    end
end

"Ideal linear electrical resistor"
@model Resistor begin
    R = Parameter(1.0, start = 1.0, info = "Resistance", T = Unitful.Ω)
    @extends OnePort()
    @inherits i, v
    @equations begin
        R*i = v
    end
end

@model Resistor1 begin
    p = Pin()
    n = Pin()
    v = Float()
    i = Float()
    R = Parameter(description="Resistance")
    @equations begin
        v = p.v - n.v # Voltage drop
        0 = p.i + n.i # Kirchhoff's current law within component
        i = p.i
        R * i = v
    end
end

"Ideal linear electric capacitor"
@model Capacitor begin
    C = Parameter(1.0, start = 1.0, info = "Capacitance", T = Unitful.F)
    @extends OnePort(v = Voltage(start = 0.0))   # Setting state=false for v does not work with extends.
    @inherits i, v
    @equations begin
        C*der(v) = i
    end
end

@model Capacitor1 begin
    p = Pin()
    n = Pin()
    v = Float()
    i = Float()
    C = undefined
    @equations begin
        v = p.v - n.v # Voltage drop
        0 = p.i + n.i # Kirchhoff's current law within component
        i = p.i
        C * der(v) = i
    end
end

"Ideal linear electric inductor"
@model Inductor begin
    L = Parameter(1.0, start = 1.0, info = "Inductance", T = Unitful.H)
    @extends OnePort()
    @inherits i, v
    @equations begin
        L * der(i) = v
    end
end

"Constant voltage source"
@model ConstantVoltage begin
    V = Parameter(1.0, start = 1.0, info = "Voltage", T = Unitful.V)
    @extends OnePort()
    @inherits v
    @equations begin
        v = V
    end
end

"Step voltage source"
@model StepVoltage begin
    V = Parameter(1.0, start = 1.0, info = "Voltage", T = Unitful.V)
    startTime = Parameter(0.0, start = 0.0, info = "Start time", T = Unitful.s)
    t = Var(start = 0.0)
    @extends OnePort()
    @inherits v
    @equations begin
        v = if t < startTime; 0 else V end
        der(t) = 1
    end
end

"Generic voltage source using the input signal as source voltage"
@model SignalVoltage begin
    p = Pin()
    n = Pin()
    v = Float()
    i = Float()
    @equations begin
        v = p.v - n.v
        0 = p.i + n.i
        i = p.i
    end
end

@model VoltageSource begin
    @extends OnePort()
    @inherits v
    offset = 0 # Voltage offset
    startTime = 0 # Time offset
    signalSource = SignalSource(offset=offset, startTime=startTime)
    @equations begin
        v = signalSource.y
    end
end

@model SineVoltage1 begin
    # Sine voltage source
    V = Parameter() # Amplitude of sine wave
    phase = 0 # Phase of sine wave
    freqHz = Parameter() # Frequency of sine wave
    @extends VoltageSource(signalSource=Sine(amplitude=V, freqHz=freqHz, phase=phase))
end

"""
Sinusoidal voltage source
"""
@model SineVoltage begin
    # Sine voltage source
    V = Parameter() # Amplitude of sine wave
    phase = 0 # Phase of sine wave
    freqHz = Parameter() # Frequency of sine wave
    # @extends VoltageSource(signalSource = Sine())
    @extends VoltageSource(signalSource=Sine(amplitude=1, freqHz=1, offset=0, startTime=0))
    # @extends VoltageSource(signalSource = Sine(amplitude=V, freqHz=freqHz, phase=phase,
    # offset=offset, startTime=startTime))
    # @inherits offset, startTime
end

"""
Ideal operational amplifier (norator-nullator pair), but 3 pins
"""
@model IdealOpAmp3Pin begin
    in_p = Pin()
    in_n = Pin()
    out = Pin()
    @equations begin
        in_p.v = in_n.v
        in_p.i = 0
        in_n.i = 0
    end
end

"""
Ideal diode
"""
@model IdealDiode begin
    @extends OnePort()
    @inherits v, i
    s = Float(start=0.0, info = "Auxiliary variable") # For position on the diode characteristic
    #=
    s = 0: knee point
    s < 0: below knee point, diode conducting
    s > 0: above knee point, diode locking
    =#
    @equations begin
        v = if positive(s); 0 else s end
        i = if positive(s); s else 0 end
    end
end

@model IdealDiode2 begin # Ideal diode
    @extends OnePort()
    @inherits v, i
    Ron = 0 # 1.0E-5 # Forward state-on differential resistance (closed diode resistance)
    Goff = 0 # 1.0E-5 # Backward state-off conductance (opened diode conductance)
    Vknee = 0 # Forward threshold voltage
    # off = Variable(start=true) # Switching state
    s = Float(start=0.0) # Auxiliary variable for actual position on the ideal diode characteristic
    #=
    s = 0: knee point
    s < 0: below knee point, diode conducting
    s > 0: above knee point, diode locking
    =#
    @equations begin
        # off := s < 0
        v = s * if !positive(s); 1 else Ron end + Vknee
        i = s * if !positive(s); Goff else 1 end + Goff * Vknee
    end
end

@model IdealDiode1 begin # Ideal diode
    @extends OnePort()
    @inherits v, i
    Ron = 0 # 1.0E-5 # Forward state-on differential resistance (closed diode resistance)
    Goff = 0 # 1.0E-5 # Backward state-off conductance (opened diode conductance)
    Vknee = 0 # Forward threshold voltage
    off = Variable(start=true) # Switching state
    s = Float(start=0) # Auxiliary variable for actual position on the ideal diode characteristic
    #=
    s = 0: knee point
    s < 0: below knee point, diode conducting
    s > 0: above knee point, diode locking
    =#
    @equations begin
        off := positive(-s * 1.0)
        v = s * if off; 1 else Ron end + Vknee
        i = s * if off; Goff else 1 end + Goff * Vknee
    end
end

end


text:
export AbsDocOrNode, AbsDocument, AbsNode, AbsXml, AbsHtml, AbsAttribute, AbsNormal, AbsEmpty, AbsIgnore, AbsText, UN

UN{T}= Union{T, Nothing}

abstract type AbsDocOrNode end

abstract type AbsDocument <:AbsDocOrNode end
abstract type AbsNode <:AbsDocOrNode end

# Documents
abstract type AbsXml <: AbsDocument end
abstract type AbsHtml <: AbsDocument end

# Nodes
abstract type AbsNormal <: AbsNode end
abstract type AbsAttribute <: AbsNode end
abstract type AbsText <: AbsNode end

abstract type AbsEmpty <: AbsNormal end

# Ignore
abstract type AbsIgnore <: AbsDocOrNode end


text:
module LabViewXML

import LightXML

include("xmlparser.jl")

end # module


text:
using KDTrees
using StatsBase
using Plotly

function run_bench_build_tree(dim, knn, exps, rounds)
    println("Running build tree benchmark for 10^(", exps, ") points in ", dim, " dimensions.")
    n_points = [10^i for i in exps]

    times = zeros(length(n_points), rounds)

    timer = 0.0
    for (i, n_point) in enumerate(n_points)
        for (j , round) in enumerate(1:rounds)
            n_iters = 5
            println("Round ", j, " out of ", rounds, " for ", dim, "x", n_point, "...\r")
            while true
                timer = time_ns()
                for k in 1:n_iters
                  data = rand(dim, int(n_point))
                  tree = KDTree(data, 10, true)
                end
                timer = (time_ns() - float(timer)) / 10^9 # To seconds
                if timer < 1.0
                    n_iters *= 3
                    continue
                end
                break # Ends this round
            end
        times[i, j] = timer / n_iters
        end
        println("\n")
    end
    println("\nDone!")
    return mean_and_std(times, 2)
end

dim = 3
knn = 5
exps = 3:0.5:6
rounds = 3

times, stds = run_bench_build_tree(dim, knn, exps, rounds)
times = vec(times)
stds = vec(stds)

stderr = stds / sqrt(length(stds))
ymins = times - 1.96*stderr
ymaxs = times + 1.96*stderr
sizes = [10.0^i::Float64 for i in exps]

####################################################################
# Plotting
####################################################################
data = [
  [
    "x" => sizes,
    "y" => times,
    "type" => "scatter",
    "mode" => "lines+markers",
    "error_y" => [
        "type" => "data",
        "array" => stderr*1.96*2,
        "visible" => true
      ]
    ]
]

layout = [
  "title" => "Build time for KDTree (dim = 3)",
  "xaxis" => [
      "title" => "Number of points",
      "type" => "log"
  ],
  "yaxis" => [
      "title" => "Build time [s]",
      "type" => "log"
  ],
  "type" => "log",
  "autorange" => true
]

Plotly.signin("kcarlsson89", "lolololoololololo")

response = Plotly.plot(data, ["layout" => layout,
                              "filename" => "bench_build_x",
                              "fileopt" => "overwrite"])
plot_url = response["url"]

#=
using Gadfly
p = plot(x = sizes, y = times,  Geom.point, Geom.line, Geom.errorbar, Scale.x_log10, Scale.y_log10,
         ymin = ymins,
         ymax = ymaxs,
         Guide.xticks(ticks=collect(exps)),
         Guide.xlabel("Number of points"),
         Guide.ylabel("Build time [s]"),
         Guide.title("Build time for KDTree (dim = 3)"))
g = render(p)
img = PNG("build_plot.png", 5inch, 4inch)
draw(img, g)
=#


text:
using Weave

# weave("docs.jmd",out_path=:pwd,doctype="md2html")
weave(
    "docs.jmd",
    out_path=:pwd,
    doctype="md2pdf",
    template="julia_tex.tpl",
)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


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

mutable struct IoK8sApiRbacV1alpha1Role <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    rules::Any # spec type: Union{ Nothing, Vector{IoK8sApiRbacV1alpha1PolicyRule} } # spec name: rules

    function IoK8sApiRbacV1alpha1Role(;apiVersion=nothing, kind=nothing, metadata=nothing, rules=nothing)
        o = new()
        validate_property(IoK8sApiRbacV1alpha1Role, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiRbacV1alpha1Role, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiRbacV1alpha1Role, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiRbacV1alpha1Role, Symbol("rules"), rules)
        setfield!(o, Symbol("rules"), rules)
        o
    end
end # type IoK8sApiRbacV1alpha1Role

const _property_map_IoK8sApiRbacV1alpha1Role = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"), Symbol("rules")=>Symbol("rules"))
const _property_types_IoK8sApiRbacV1alpha1Role = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("rules")=>"Vector{IoK8sApiRbacV1alpha1PolicyRule}")
Base.propertynames(::Type{ IoK8sApiRbacV1alpha1Role }) = collect(keys(_property_map_IoK8sApiRbacV1alpha1Role))
Swagger.property_type(::Type{ IoK8sApiRbacV1alpha1Role }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiRbacV1alpha1Role[name]))}
Swagger.field_name(::Type{ IoK8sApiRbacV1alpha1Role }, property_name::Symbol) =  _property_map_IoK8sApiRbacV1alpha1Role[property_name]

function check_required(o::IoK8sApiRbacV1alpha1Role)
    true
end

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


text:
import Base.iterate
export FactSampler, trace

struct FactSampler{TF,T∇ϕ,Tc,Tu0,TG,Trng,Targs} <: PDMPSampler
    F::TF
    ∇ϕ::T∇ϕ
    c::Tc
    u0::Tu0
    G::TG
    rng::Trng

    args::Targs

    factor::Float64
    adapt::Bool
end

function FactSampler(∇ϕ, u0, c, G, F::Union{ZigZag,FactBoomerang}, args...;
    factor=1.8, adapt=false, seed=Seed())
    return FactSampler(F, ∇ϕ, c, u0, G, Rng(seed), args, factor, adapt)
end
function FactSampler(∇ϕ, u0, c, F::Union{ZigZag,FactBoomerang}, args...;
    factor=1.8, adapt=false)
    return FactSampler(F, ∇ϕ, c, u0, nothing, Rng(), args, factor, adapt)
end

function iterate(FS::FactSampler)
    t0, (x0, θ0) = FS.u0
    F = FS.F
    n = length(x0)
    t′ = t0
    t = fill(t′, size(θ0)...)
    t_old = copy(t)
    G1 = [i => rowvals(F.Γ)[nzrange(F.Γ, i)] for i in eachindex(θ0)]
    G = FS.G
    if G === nothing
        G = G1
    end
    @assert all(a.second ⊇ b.second for (a,b) in zip(G, G1))
    G2 = [i => setdiff(union((G1[j].second for j in G1[i].second)...), G[i].second) for i in eachindex(G1)]
    x, θ = copy(x0), copy(θ0)
    num = 0
    acc = zeros(Int, n)
    Q = SPriorityQueue{Int,Float64}()
    b = [ab(G1, i, x, θ, FS.c, F) for i in eachindex(θ)]
    for i in eachindex(θ)
        enqueue!(Q, i => poisson_time(b[i], rand(FS.rng)))
    end
    if hasrefresh(F)
        enqueue!(Q, (n + 1) => waiting_time_ref(F))
    end
    iterate(FS, ((t => (x, θ)), t_old, (acc, num), Q, b, G, G1, G2))
end

function iterate(FS::FactSampler, (u, t_old, (acc, num), Q, b, G, G1, G2))
    t, (x, θ) = u
    n = length(x)
    ev, t, x, θ, t′, (acc, num), _,  b, t_old = spdmp_inner!(FS.rng, G, G1, G2, FS.∇ϕ, t, x, θ, Q,
    FS.c, b, t_old, (acc, num), FS.F, FS.args...; factor=FS.factor, adapt=FS.adapt)
    u = t => (x, θ)
    return (t′ => ev), (u, t_old, (acc, num), Q, b, G, G1, G2)
end

function trace(FS::FactSampler, T)
    t0, (x0, θ0) = FS.u0
    Ξ = Trace(t0, x0, θ0, FS.F)
    ϕ = iterate(FS)
    ϕ === nothing && error("No events")
    ev, state = ϕ
    while true
        ϕ = iterate(FS, state)
        ϕ === nothing && return Ξ
        (t, ev), state = ϕ
        t > T && return Ξ
        push!(Ξ, ev)
    end
end


text:
using Test
using Measures

const imm = Length(:mm, 1)
const rmm = Length(:mm, 1//1)

const xx = Length(:xx, 1.0)

@testset "Bounding boxes" begin
    @testset "Valid bounding box constructors" begin
        @test BoundingBox((0mm, 1mm), 1mm, 2mm) ===
            BoundingBox((0mm, 1mm), 1mm, 2mm)
        @test BoundingBox((2mm, 1mm, 3mm), 1mm, 2mm, 4mm) ===
            BoundingBox((2mm, 1mm, 3mm), 1mm, 2mm, 4mm)
    end

    @testset "Invalid bounding box constructors" begin
        @test_throws MethodError BoundingBox((0mm,), 1mm, 2mm)
        @test_throws MethodError BoundingBox((0mm, 3mm, 8mm), 1mm, 4mm)
        @test_throws MethodError BoundingBox((0mm, 3mm), 1mm, 2mm, 4mm)
        @test_throws MethodError BoundingBox((0mm, 6mm, 3mm, 8mm), 1mm, 2mm, 4mm)
    end
end

@testset "Measure constructors" begin
    @testset "Length constructor" begin
        @test Length(:mm, 1) === Length{:mm, Int}(1)
        @test Length(:mm, 1.1) === Length{:mm, Float64}(1.1)
        @test 1mm |> typeof === Length{:mm, Float64}
        @test convert(Length{:mm, Float64}, 1imm) |> typeof === Length{:mm, Float64}
    end

end

function check_binop(∘) #\circ
    @test (1imm ∘ 2imm).value === 1 ∘ 2
    @test (1mm ∘ 2imm).value === 1.0 ∘ 2
    @test (2mm ∘ 1imm).value === 2.0 ∘ 1
    @test (typeof(2mm ∘ 3xx) <: Length) === false
    @test typeof(2mm ∘ 3xx) <: Measures.BinaryOp
end

function check_unaryop(>:)
    @test >:1imm.value === >:1
    @test >:1.0imm.value !== >:1
    @test >:1.0imm.value === >:1.0
end

@testset "operators" begin
    @testset "type promotion in +" begin
        @test typeof(1imm + 2imm) === Length{:mm, Int}
        @test typeof(1imm + 2mm) === Length{:mm, Float64}
        @test typeof(1mm + 3xx) === Measures.Add{Length{:mm, Float64}, Length{:xx, Float64}}
    end

    @testset "equality" begin
        @test 1imm == 1mm
        @test 1imm !== 1xx
        @test 1imm + 1xx !== 1mm + 1xx # perhaps make this true?
    end

    @testset "hash" begin
        @test hash(1imm) === hash(1mm)
        @test hash(1imm) !== hash(1xx)
        @test isequal(1imm, 1mm)
        @test isequal(1imm, 1xx) === false
    end

    @testset "+" begin
        check_unaryop(+)
        check_binop(+)
    end
    @testset "-" begin
        check_unaryop(-)
        check_binop(-)
    end
    @testset "max" begin
        check_binop(max)
    end
    @testset "min" begin
        check_binop(max)
    end
    @testset "broadcasting" begin
        a = [1mm, 2mm, 3mm]
        @test all((a .+ 1mm) .== [2mm, 3mm, 4mm])
    end
end


text:
export normlogpdf

using Distributions: DiscreteNonParametric, support, probs

# watch https://github.com/JuliaStats/Distributions.jl/issues/1183

"""
GPU automatic differentiable version for the logpdf function of normal distributions.
Adding an epsilon value to guarantee numeric stability if sigma is exactly zero
(e.g. if relu is used in output layer).
"""
function normlogpdf(μ, σ, x; ϵ = 1.0f-8)
    z = (x .- μ) ./ (σ .+ ϵ)
    -(z .^ 2 .+ log(2.0f0π)) / 2.0f0 .- log.(σ .+ ϵ)
end

# watch https://github.com/JuliaStats/Distributions.jl/pull/1184

Base.convert(::Type{DiscreteNonParametric{T,P,Ts,Ps}}, d::DiscreteNonParametric) where {T,P,Ts,Ps} = DiscreteNonParametric{T,P,Ts,Ps}(convert(Ts, support(d)), convert(Ps, probs(d)), check_args=false)


text:
### A Pluto.jl notebook ###
# v0.14.5

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 1e06178a-1fbf-11eb-32b3-61769a79b7c0
begin
	import Pkg
	Pkg.activate(mktempdir())
	Pkg.add([
			"Plots",
			"PlutoUI",
			"LaTeXStrings",
			"Distributions",
			"Random",
	])
	using LaTeXStrings
	using Plots
	using PlutoUI
	using Random, Distributions
end

# ╔═╡ 169727be-2433-11eb-07ae-ab7976b5be90
md"_Homework 10 Spring 2021, version 1_"

# ╔═╡ 21524c08-2433-11eb-0c55-47b1bdc9e459
md"""

# **Homework 10**: _Climate modeling I_
`18.S191`, spring 2021
"""

# ╔═╡ 23335418-2433-11eb-05e4-2b35dc6cca0e
# edit the code below to set your name and kerberos ID (i.e. email without @mit.edu)

student = (name = "Jazzy Doe", kerberos_id = "jazz")

# you might need to wait until all other cells in this notebook have completed running.
# scroll around the page to see what's up

# ╔═╡ 18be4f7c-2433-11eb-33cb-8d90ca6f124c
md"""

Submission by: **_$(student.name)_** ($(student.kerberos_id)@mit.edu)
"""

# ╔═╡ 253f4da0-2433-11eb-1e48-4906059607d3
md"_Let's create a package environment:_"

# ╔═╡ 87e68a4a-2433-11eb-3e9d-21675850ed71
html"""

"""

# ╔═╡ fe3304f8-2668-11eb-066d-fdacadce5a19
md"""
_Before working on the homework, make sure that you have watched the first lecture on climate modeling 👆. We have included the important functions from this lecture notebook in the next cell. Feel free to have a look!_
"""

# ╔═╡ 930d7154-1fbf-11eb-1c3a-b1970d291811
module Model

const S = 1368; # solar insolation [W/m^2]  (energy per unit time per unit area)
const α = 0.3; # albedo, or planetary reflectivity [unitless]
const B = -1.3; # climate feedback parameter [W/m^2/°C],
const T0 = 14.; # preindustrial temperature [°C]

absorbed_solar_radiation(; α=α, S=S) = S*(1 - α)/4; # [W/m^2]
outgoing_thermal_radiation(T; A=A, B=B) = A - B*T;

const A = S*(1. - α)/4 + B*T0; # [W/m^2].

greenhouse_effect(CO2; a=a, CO2_PI=CO2_PI) = a*log(CO2/CO2_PI);

const a = 5.0; # CO2 forcing coefficient [W/m^2]
const CO2_PI = 280.; # preindustrial CO2 concentration [parts per million; ppm];
CO2_const(t) = CO2_PI; # constant CO2 concentrations

const C = 51.; # atmosphere and upper-ocean heat capacity [J/m^2/°C]

function timestep!(ebm)
	append!(ebm.T, ebm.T[end] + ebm.Δt*tendency(ebm));
	append!(ebm.t, ebm.t[end] + ebm.Δt);
end;

tendency(ebm) = (1. /ebm.C) * (
	+ absorbed_solar_radiation(α=ebm.α, S=ebm.S)
	- outgoing_thermal_radiation(ebm.T[end], A=ebm.A, B=ebm.B)
	+ greenhouse_effect(ebm.CO2(ebm.t[end]), a=ebm.a, CO2_PI=ebm.CO2_PI)
);

begin
	mutable struct EBM
		T::Array{Float64, 1}

		t::Array{Float64, 1}
		Δt::Float64

		CO2::Function

		C::Float64
		a::Float64
		A::Float64
		B::Float64
		CO2_PI::Float64

		α::Float64
		S::Float64
	end;

	# Make constant parameters optional kwargs
	EBM(T::Array{Float64, 1}, t::Array{Float64, 1}, Δt::Real, CO2::Function;
		C=C, a=a, A=A, B=B, CO2_PI=CO2_PI, α=α, S=S) = (
		EBM(T, t, Δt, CO2, C, a, A, B, CO2_PI, α, S)
	);

	# Construct from float inputs for convenience
	EBM(T0::Real, t0::Real, Δt::Real, CO2::Function;
		C=C, a=a, A=A, B=B, CO2_PI=CO2_PI, α=α, S=S) = (
		EBM(Float64[T0], Float64[t0], Δt, CO2;
			C=C, a=a, A=A, B=B, CO2_PI=CO2_PI, α=α, S=S);
	);
end;

begin
	function run!(ebm::EBM, end_year::Real)
		while ebm.t[end] < end_year
			timestep!(ebm)
		end
	end;

	run!(ebm) = run!(ebm, 200.) # run for 200 years by default
end

CO2_hist(t) = CO2_PI * (1 .+ fractional_increase(t));
fractional_increase(t) = ((t .- 1850.)/220).^3;

begin
	CO2_RCP26(t) = CO2_PI * (1 .+ fractional_increase(t) .* min.(1., exp.(-((t .-1850.).-170)/100))) ;
	RCP26 = EBM(T0, 1850., 1., CO2_RCP26)
	run!(RCP26, 2100.)

	CO2_RCP85(t) = CO2_PI * (1 .+ fractional_increase(t) .* max.(1., exp.(((t .-1850.).-170)/100)));
	RCP85 = EBM(T0, 1850., 1., CO2_RCP85)
	run!(RCP85, 2100.)
end

end

# ╔═╡ 1312525c-1fc0-11eb-2756-5bc3101d2260
md"""## **Exercise 1** - _policy goals under uncertainty_
A recent ground-breaking [review paper](https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2019RG000678) produced the most comprehensive and up-to-date estimate of the *climate feedback parameter*, which they find to be

$B \approx \mathcal{N}(-1.3, 0.4),$

i.e. our knowledge of the real value is normally distributed with a mean value $\overline{B} = -1.3$ W/m²/K and a standard deviation $\sigma = 0.4$ W/m²/K. These values are not very intuitive, so let us convert them into more policy-relevant numbers.

**Definition:** *Equilibrium climate sensitivity (ECS)* is defined as the amount of warming $\Delta T$ caused by a doubling of CO₂ (e.g. from the pre-industrial value 280 ppm to 560 ppm), at equilibrium.

At equilibrium, the energy balance model equation is:

$0 = \frac{S(1 - α)}{4} - (A - BT_{eq}) + a \ln\left( \frac{2\;\text{CO}₂_{\text{PI}}}{\text{CO}₂_{\text{PI}}} \right)$

From this, we subtract the preindustrial energy balance, which is given by:

$0 = \frac{S(1-α)}{4} - (A - BT_{0}),$

The result of this subtraction, after rearranging, is our definition of $\text{ECS}$:

$\text{ECS} \equiv T_{eq} - T_{0} = -\frac{a\ln(2)}{B}$
"""

# ╔═╡ 7f961bc0-1fc5-11eb-1f18-612aeff0d8df
md"""The plot below provides an example of an "abrupt 2 × CO₂" experiment, a classic experimental treatment method in climate modelling which is used in practice to estimate ECS for a particular model. (Note: in complicated climate models the values of the parameters $a$ and $B$ are not specified *a priori*, but *emerge* as outputs of the simulation.)

The simulation begins at the preindustrial equilibrium, i.e. a temperature $T_{0} = 14$°C is in balance with the pre-industrial CO₂ concentration of 280 ppm until CO₂ is abruptly doubled from 280 ppm to 560 ppm. The climate responds by warming rapidly, and after a few hundred years approaches the equilibrium climate sensitivity value, by definition.
"""

# ╔═╡ fa7e6f7e-2434-11eb-1e61-1b1858bb0988
md"""
``B = `` $(@bind B_slider Slider(-2.5:.001:0; show_value=true, default=-1.3))
"""

# ╔═╡ 16348b6a-1fc2-11eb-0b9c-65df528db2a1
md"""
##### Exercise 1.1 - _Develop understanding for feedbacks and climate sensitivity_
"""

# ╔═╡ e296c6e8-259c-11eb-1385-53f757f4d585
md"""
👉 Change the value of $B$ using the slider above. What does it mean for a climate system to have a more negative value of $B$? Explain why we call $B$ the _climate feedback parameter_.
"""

# ╔═╡ a86f13de-259d-11eb-3f46-1f6fb40020ce
observations_from_changing_B = md"""
Hello world!
"""

# ╔═╡ 3d66bd30-259d-11eb-2694-471fb3a4a7be
md"""
👉 What happens when $B$ is greater than or equal to zero?
"""

# ╔═╡ 5f82dec8-259e-11eb-2f4f-4d661f44ef41
observations_from_nonnegative_B = md"""
Hello world!
"""

# ╔═╡ 56b68356-2601-11eb-39a9-5f4b8e580b87
md"Reveal answer: $(@bind reveal_nonnegative_B_answer CheckBox())"

# ╔═╡ 7d815988-1fc7-11eb-322a-4509e7128ce3
if reveal_nonnegative_B_answer
	md"""
This is known as the "runaway greenhouse effect", where warming self-amplifies so strongly through *positive feedbacks* that the warming continues forever (or until the oceans boil away and there is no longer a reservoir or water to support a *water vapor feedback*. This is thought to explain Venus' extremely hot and hostile climate, but as you can see is extremely unlikely to occur on present-day Earth.
"""
end

# ╔═╡ aed8f00e-266b-11eb-156d-8bb09de0dc2b
md"""
👉 Create a graph to visualize ECS as a function of B.
"""

# ╔═╡ b9f882d8-266b-11eb-2998-75d6539088c7

# ╔═╡ 269200ec-259f-11eb-353b-0b73523ef71a
md"""
#### Exercise 1.2 - _Doubling CO₂_

To compute ECS, we doubled the CO₂ in our atmosphere. This factor 2 is not entirely arbitrary: without substantial effort to reduce CO₂ emissions, we are expected to **at least** double the CO₂ in our atmosphere by 2100.

Right now, our CO₂ concentration is 415 ppm -- $(round(415 / 280, digits=3)) times the pre-industrial value of 280 ppm from 1850.

The CO₂ concentrations in the _future_ depend on human action. There are several models for future concentrations, which are formed by assuming different _policy scenarios_. A baseline model is RCP8.5 - a "worst-case" high-emissions scenario. In our notebook, this model is given as a function of ``t``.
"""

# ╔═╡ 2dfab366-25a1-11eb-15c9-b3dd9cd6b96c
md"""
👉 In what year are we expected to have doubled the CO₂ concentration, under policy scenario RCP8.5?
"""

# ╔═╡ 50ea30ba-25a1-11eb-05d8-b3d579f85652
expected_double_CO2_year = let

	missing
end

# ╔═╡ bade1372-25a1-11eb-35f4-4b43d4e8d156
md"""
#### Exercise 1.3 - _Uncertainty in B_

The climate feedback parameter ``B`` is not something that we can control– it is an emergent property of the global climate system. Unfortunately, ``B`` is also difficult to quantify empirically (the relevant processes are difficult or impossible to observe directly), so there remains uncertainty as to its exact value.

A value of ``B`` close to zero means that an increase in CO₂ concentrations will have a larger impact on global warming, and that more action is needed to stay below a maximum temperature. In answering such policy-related question, we need to take the uncertainty in ``B`` into account. In this exercise, we will do so using a Monte Carlo simulation: we generate a sample of values for ``B``, and use these values in our analysis.
"""

# ╔═╡ 02232964-2603-11eb-2c4c-c7b7e5fed7d1
B̅ = -1.3; σ = 0.4

# ╔═╡ c4398f9c-1fc4-11eb-0bbb-37f066c6027d
ECS(; B=B̅, a=Model.a) = -a*log(2.)./B;

# ╔═╡ 25f92dec-1fc4-11eb-055d-f34deea81d0e
let
	double_CO2(t) = if t >= 0
		2*Model.CO2_PI
	else
		Model.CO2_PI
	end

	# the definition of A depends on B, so we recalculate:
	A = Model.S*(1. - Model.α)/4 + B_slider*Model.T0
	# create the model
	ebm_ECS = Model.EBM(14., -100., 1., double_CO2, A=A, B=B_slider);
	Model.run!(ebm_ECS, 300)

	ecs = ECS(B=B_slider)

	p = plot(
		size=(500,250), legend=:bottomright,
		title="Transient response to instant doubling of CO₂",
		ylabel="temperature change [°C]", xlabel="years after doubling",
		ylim=(-.5, (isfinite(ecs) && ecs < 4) ? 4 : 10),
	)

	plot!(p, [ebm_ECS.t[1], ebm_ECS.t[end]], ecs .* [1,1],
		ls=:dash, color=:darkred, label="ECS")

	plot!(p, ebm_ECS.t, ebm_ECS.T .- ebm_ECS.T[1],
		label="ΔT(t) = T(t) - T₀")
end |> as_svg

# ╔═╡ 736ed1b6-1fc2-11eb-359e-a1be0a188670
B_samples = let
	B_distribution = Normal(B̅, σ)
	Nsamples = 5000

	samples = rand(B_distribution, Nsamples)
	# we only sample negative values of B
	filter(x -> x < 0, samples)
end

# ╔═╡ 49cb5174-1fc3-11eb-3670-c3868c9b0255
histogram(B_samples, size=(600, 250), label=nothing, xlabel="B [W/m²/K]", ylabel="samples")

# ╔═╡ f3abc83c-1fc7-11eb-1aa8-01ce67c8bdde
md"""
👉 Generate a probability distribution for the ECS based on the probability distribution function for $B$ above. Plot a histogram.
"""

# ╔═╡ 3d72ab3a-2689-11eb-360d-9b3d829b78a9
ECS_samples = missing

# ╔═╡ b6d7a362-1fc8-11eb-03bc-89464b55c6fc
md"**Answer:**"

# ╔═╡ 1f148d9a-1fc8-11eb-158e-9d784e390b24

# ╔═╡ cf8dca6c-1fc8-11eb-1f89-099e6ba53c22
md"It looks like the ECS distribution is **not normally distributed**, even though $B$ is.

👉 How does $\overline{\text{ECS}(B)}$ compare to $\text{ECS}(\overline{B})$? What is the probability that $\text{ECS}(B)$ lies above $\text{ECS}(\overline{B})$?
"

# ╔═╡ 02173c7a-2695-11eb-251c-65efb5b4a45f

# ╔═╡ 440271b6-25e8-11eb-26ce-1b80aa176aca
md"👉 Does accounting for uncertainty in feedbacks make our expectation of global warming better (less implied warming) or worse (more implied warming)?"

# ╔═╡ cf276892-25e7-11eb-38f0-03f75c90dd9e
observations_from_the_order_of_averaging = md"""
Hello world!
"""

# ╔═╡ 5b5f25f0-266c-11eb-25d4-17e411c850c9
md"""
#### Exercise 1.5 - _Running the model_

In the lecture notebook we introduced a _mutable struct_ `EBM` (_energy balance model_), which contains:
- the parameters of our climate simulation (`C`, `a`, `A`, `B`, `CO2_PI`, `α`, `S`, see details below)
- a function `CO2`, which maps a time `t` to the concentrations at that year. For example, we use the function `t -> 280` to simulate a model with concentrations fixed at 280 ppm.

`EBM` also contains the simulation results, in two arrays:
- `T` is the array of tempartures (°C, `Float64`).
- `t` is the array of timestamps (years, `Float64`), of the same size as `T`.
"""

# ╔═╡ 3f823490-266d-11eb-1ba4-d5a23975c335
html"""

Properties of an EBM obect:

Name	Description

A	Linearized outgoing thermal radiation: offset <soft>[W/m²]</soft>

B	Linearized outgoing thermal radiation: slope. or: climate feedback parameter <soft>[W/m²/°C]</soft>

α	Planet albedo, 0.0-1.0 <soft>[unitless]</soft>

S	Solar insulation <soft>[W/m²]</soft>

C	Atmosphere and upper-ocean heat capacity <soft>[J/m²/°C]</soft>

a	CO₂ forcing effect <soft>[W/m²]</soft>

CO2_PI	Pre-industrial CO₂ concentration <soft>[ppm]</soft>

"""

# ╔═╡ 971f401e-266c-11eb-3104-171ae299ef70
md"""

You can set up an instance of `EBM` like so:
"""

# ╔═╡ 746aa5bc-266c-11eb-14c9-63ccc313f5de
empty_ebm = Model.EBM(
	14.0, # initial temperature
	1850, # initial year
	1, # Δt
	t -> 280.0, # CO2 function
)

# ╔═╡ a919d584-2670-11eb-1cf9-2327c8135d6d
md"""
Have a look inside this object. We see that `T` and `t` are initialized to a 1-element array.

Let's run our model:
"""

# ╔═╡ bfb07a0a-2670-11eb-3938-772499c637b1
simulated_model = let
	ebm = Model.EBM(14.0, 1850, 1, t -> 280.0)
	Model.run!(ebm, 2020)
	ebm
end

# ╔═╡ 12cbbab0-2671-11eb-2b1f-038c206e84ce
md"""
Again, look inside `simulated_model` and notice that `T` and `t` have accumulated the simulation results.

In this simulation, we used `T0 = 14` and `CO2 = t -> 280`, which is why `T` is constant during our simulation. These parameters are the default, pre-industrial values, and our model is based on this equilibrium.

👉 Run a simulation with policy scenario RCP8.5, and plot the computed temperature graph. What is the global temperature in the year 2100?
"""

# ╔═╡ 9596c2dc-2671-11eb-36b9-c1af7e5f1089
simulated_rcp85_model = let

	missing
end

# ╔═╡ f94a1d56-2671-11eb-2cdc-810a9c7a8a5f

# ╔═╡ 4b091fac-2672-11eb-0db8-75457788d85e
md"""
Additional parameters can be set using keyword arguments. For example:

```julia
Model.EBM(14, 1850, 1, t -> 280.0; B=-2.0)
```
Creates the same model as before, but with `B = -2.0`.
"""

# ╔═╡ 9cdc5f84-2671-11eb-3c78-e3495bc64d33
md"""
👉 Write a function `temperature_response` that takes a function `CO2` and an optional value `B` as parameters, and returns the temperature at 2100 according to our model.
"""

# ╔═╡ f688f9f2-2671-11eb-1d71-a57c9817433f
function temperature_response(CO2::Function, B::Float64=-1.3)

	return missing
end

# ╔═╡ 049a866e-2672-11eb-29f7-bfea7ad8f572
temperature_response(t -> 280)

# ╔═╡ 09901de6-2672-11eb-3d50-05b176b729e7
temperature_response(Model.CO2_RCP85)

# ╔═╡ aea0d0b4-2672-11eb-231e-395c863827d3
temperature_response(Model.CO2_RCP85, -1.0)

# ╔═╡ 9c32db5c-1fc9-11eb-029a-d5d554de1067
md"""#### Exercise 1.6 - _Application to policy relevant questions_

We talked about two _emissions scenarios_: RCP2.6 (strong mitigation - controlled CO2 concentrations) and RCP8.5 (no mitigation - high CO2 concentrations). These are given by the following functions:
"""

# ╔═╡ ee1be5dc-252b-11eb-0865-291aa823b9e9
t = 1850:2100

# ╔═╡ e10a9b70-25a0-11eb-2aed-17ed8221c208
plot(t, Model.CO2_RCP85.(t),
	ylim=(0,1200), ylabel="CO2 concentration [ppm]")

# ╔═╡ 40f1e7d8-252d-11eb-0549-49ca4e806e16
@bind t_scenario_test Slider(t; show_value=true, default=1850)

# ╔═╡ 19957754-252d-11eb-1e0a-930b5208f5ac
Model.CO2_RCP26(t_scenario_test), Model.CO2_RCP85(t_scenario_test)

# ╔═╡ 06c5139e-252d-11eb-2645-8b324b24c405
md"""
We are interested in how the **uncertainty in our input** $B$ (the climate feedback paramter) *propagates* through our model to determine the **uncertainty in our output** $T(t)$, for a given emissions scenario. The goal of this exercise is to answer the following by using *Monte Carlo Simulation* for *uncertainty propagation*:

> 👉 What is the probability that we see more than 2°C of warming by 2100 under the low-emissions scenario RCP2.6? What about under the high-emissions scenario RCP8.5?

"""

# ╔═╡ f2e55166-25ff-11eb-0297-796e97c62b07

# ╔═╡ 1ea81214-1fca-11eb-2442-7b0b448b49d6
md"""
## **Exercise 2** - _How did Snowball Earth melt?_

In lecture 21 (see below), we discovered that increases in the brightness of the Sun are not sufficient to explain how Snowball Earth eventually melted.
"""

# ╔═╡ a0ef04b0-25e9-11eb-1110-cde93601f712
html"""

"""

# ╔═╡ 3e310cf8-25ec-11eb-07da-cb4a2c71ae34
md"""
We talked about a second theory -- a large increase in CO₂ (by volcanoes) could have caused a strong enough greenhouse effect to melt the Snowball. If we imagine that the CO₂ then decreased (e.g. by getting sequestered by the now liquid ocean), we might be able to explain how we transitioned from a hostile Snowball Earth to today's habitable "Waterball" Earth.

In this exercise, you will estimate how much CO₂ would be needed to melt the Snowball and visualize a possible trajectory for Earth's climate over the past 700 million years by making an interactive *bifurcation diagram*.

#### Exercise 2.1

In the [lecture notebook](https://github.com/hdrake/simplEarth/blob/master/2_ebm_multiple_equilibria.jl) (video above), we had a bifurcation diagram of $S$ (solar insolation) vs $T$ (temperature). We increased $S$, watched our point move right in the diagram until we found the tipping point. This time we will do the same, but we vary the CO₂ concentration, and keep $S$ fixed at its default (present day) value.
"""

# ╔═╡ d6d1b312-2543-11eb-1cb2-e5b801686ffb
md"""
Below we have an empty diagram, which is already set up with a CO₂ vs $T$ diagram, with a logarthmic horizontal axis. Now it's your turn! We have written some pointers below to help you, but feel free to do it your own way.
"""

# ╔═╡ 3cbc95ba-2685-11eb-3810-3bf38aa33231
md"""
We used two helper functions:
"""

# ╔═╡ 68b2a560-2536-11eb-0cc4-27793b4d6a70
function add_cold_hot_areas!(p)

	left, right = xlims(p)

	plot!(p,
		[left, right], [-60, -60],
		fillrange=[-10., -10.], fillalpha=0.3, c=:lightblue, label=nothing
	)
	annotate!(p,
		left+12, -19,
		text("completely\nfrozen", 10, :darkblue, :left)
	)

	plot!(p,
		[left, right], [10, 10],
		fillrange=[80., 80.], fillalpha=0.09, c=:red, lw=0., label=nothing
	)
	annotate!(p,
		left+12, 15,
		text("no ice", 10, :darkred, :left)
	)
end

# ╔═╡ 0e19f82e-2685-11eb-2e99-0d094c1aa520
function add_reference_points!(p)
	plot!(p,
		[Model.CO2_PI, Model.CO2_PI], [-55, 75],
		color=:grey, alpha=0.3, lw=8,
		label="Pre-industrial CO2"
	)
	plot!(p,
		[Model.CO2_PI], [Model.T0],
		shape=:circle, color=:orange, markersize=8,
		label="Our preindustrial climate"
	)
	plot!(p,
		[Model.CO2_PI], [-38.3],
		shape=:circle, color=:aqua, markersize=8,
		label="Alternate preindustrial climate"
	)
end

# ╔═╡ 1eabe908-268b-11eb-329b-b35160ec951e
md"""
👉 Create a slider for `CO2` between `CO2min` and `CO2max`. Just like the horizontal axis of our plot, we want the slider to be _logarithmic_.
"""

# ╔═╡ 1d388372-2695-11eb-3068-7b28a2ccb9ac

# ╔═╡ 4c9173ac-2685-11eb-2129-99071821ebeb
md"""
👉 Write a function `step_model!` that takes an existing `ebm` and `new_CO2`, which performs a step of our interactive process:
- Reset the model by setting the `ebm.t` and `ebm.T` arrays to a single element. _Which value?_
- Assign a new function to `ebm.CO2`. _What function?_
- Run the model.
"""

# ╔═╡ 736515ba-2685-11eb-38cb-65bfcf8d1b8d
function step_model!(ebm::Model.EBM, CO2::Real)

	# your code here

	return ebm
end

# ╔═╡ 8b06b944-268c-11eb-0bfc-8d4dd21e1f02
md"""
👉 Inside the plot cell, call the function `step_model!`.
"""

# ╔═╡ 09ce27ca-268c-11eb-0cdd-c9801db876f8
md"""
##### Parameters
"""

# ╔═╡ 298deff4-2676-11eb-2595-e7e22f613ea1
CO2min = 10

# ╔═╡ 2bbf5a70-2676-11eb-1085-7130d4a30443
CO2max = 1_000_000

# ╔═╡ de95efae-2675-11eb-0909-73afcd68fd42
Tneo = -48

# ╔═╡ 06d28052-2531-11eb-39e2-e9613ab0401c
ebm = Model.EBM(Tneo, 0., 5., Model.CO2_const)

# ╔═╡ 378aed18-252b-11eb-0b37-a3b511af2cb5
let
	p = plot(
		xlims=(CO2min, CO2max), ylims=(-55, 75),
		xaxis=:log,
		xlabel="CO2 concentration [ppm]",
		ylabel="Global temperature T [°C]",
		title="Earth's CO2 concentration bifurcation diagram",
		legend=:topleft
	)

	add_cold_hot_areas!(p)
	add_reference_points!(p)

	# your code here

	plot!(p,
		[ebm.CO2(ebm.t[end])], [ebm.T[end]],
		label=nothing,
		color=:black,
		shape=:circle,
	)

end |> as_svg

# ╔═╡ c78e02b4-268a-11eb-0af7-f7c7620fcc34
md"""
The albedo feedback is implemented by the methods below:
"""

# ╔═╡ d7801e88-2530-11eb-0b93-6f1c78d00eea
function α(T; α0=Model.α, αi=0.5, ΔT=10.)
	if T < -ΔT
		return αi
	elseif -ΔT <= T < ΔT
		return αi + (α0-αi)*(T+ΔT)/(2ΔT)
	elseif T >= ΔT
		return α0
	end
end

# ╔═╡ 607058ec-253c-11eb-0fb6-add8cfb73a4f
function Model.timestep!(ebm)
	ebm.α = α(ebm.T[end]) # Added this line
	append!(ebm.T, ebm.T[end] + ebm.Δt*Model.tendency(ebm));
	append!(ebm.t, ebm.t[end] + ebm.Δt);
end

# ╔═╡ 9c1f73e0-268a-11eb-2bf1-216a5d869568
md"""
If you like, make the visualization more informative! Like in the lecture notebook, you could add a trail behind the black dot, or you could plot the stable and unstable branches. It's up to you!
"""

# ╔═╡ 11096250-2544-11eb-057b-d7112f20b05c
md"""
#### Exercise 2.2

👉 Find the **lowest CO₂ concentration** necessary to melt the Snowball, programatically (i.e., using code).
"""

# ╔═╡ 9eb07a6e-2687-11eb-0de3-7bc6aa0eefb0
co2_to_melt_snowball = let

	missing
end

# ╔═╡ 36e2dfea-2433-11eb-1c90-bb93ab25b33c
if student.name == "Jazzy Doe" || student.kerberos_id == "jazz"
	md"""
	!!! danger "Before you submit"
	    Remember to fill in your **name** and **Kerberos ID** at the top of this notebook.
	"""
end

# ╔═╡ 36ea4410-2433-11eb-1d98-ab4016245d95
md"## Function library

Just some helper functions used in the notebook."

# ╔═╡ 36f8c1e8-2433-11eb-1f6e-69dc552a4a07
hint(text) = Markdown.MD(Markdown.Admonition("hint", "Hint", [text]))

# ╔═╡ 51e2e742-25a1-11eb-2511-ab3434eacc3e
hint(md"The function `findfirst` might be helpful.")

# ╔═╡ 53c2eaf6-268b-11eb-0899-b91c03713da4
hint(md"
```julia
@bind log_CO2 Slider(❓)
```

```julia
CO2 = 10^log_CO2
```

")

# ╔═╡ cb15cd88-25ed-11eb-2be4-f31500a726c8
hint(md"Use a condition on the albedo or temperature to check whether the Snowball has melted.")

# ╔═╡ 232b9bec-2544-11eb-0401-97a60bb172fc
hint(md"Start by writing a function `equilibrium_temperature(CO2)` which creates a new `EBM` at the Snowball Earth temperature T = $(Tneo) and returns the final temperature for a given CO2 level.")

# ╔═╡ 37061f1e-2433-11eb-3879-2d31dc70a771
almost(text) = Markdown.MD(Markdown.Admonition("warning", "Almost there!", [text]))

# ╔═╡ 371352ec-2433-11eb-153d-379afa8ed15e
still_missing(text=md"Replace `missing` with your answer.") = Markdown.MD(Markdown.Admonition("warning", "Here we go!", [text]))

# ╔═╡ 372002e4-2433-11eb-0b25-39ce1b1dd3d1
keep_working(text=md"The answer is not quite right.") = Markdown.MD(Markdown.Admonition("danger", "Keep working on it!", [text]))

# ╔═╡ 372c1480-2433-11eb-3c4e-95a37d51835f
yays = [md"Fantastic!", md"Splendid!", md"Great!", md"Yay ❤", md"Great! 🎉", md"Well done!", md"Keep it up!", md"Good job!", md"Awesome!", md"You got the right answer!", md"Let's move on to the next section."]

# ╔═╡ 3737be8e-2433-11eb-2049-2d6d8a5e4753
correct(text=rand(yays)) = Markdown.MD(Markdown.Admonition("correct", "Got it!", [text]))

# ╔═╡ 374522c4-2433-11eb-3da3-17419949defc
not_defined(variable_name) = Markdown.MD(Markdown.Admonition("danger", "Oopsie!", [md"Make sure that you define a variable called **$(Markdown.Code(string(variable_name)))**"]))

# ╔═╡ 37552044-2433-11eb-1984-d16e355a7c10
TODO = html"TODO"

# ╔═╡ Cell order:
# ╟─169727be-2433-11eb-07ae-ab7976b5be90
# ╟─18be4f7c-2433-11eb-33cb-8d90ca6f124c
# ╠═21524c08-2433-11eb-0c55-47b1bdc9e459
# ╠═23335418-2433-11eb-05e4-2b35dc6cca0e
# ╟─253f4da0-2433-11eb-1e48-4906059607d3
# ╠═1e06178a-1fbf-11eb-32b3-61769a79b7c0
# ╟─87e68a4a-2433-11eb-3e9d-21675850ed71
# ╟─fe3304f8-2668-11eb-066d-fdacadce5a19
# ╟─930d7154-1fbf-11eb-1c3a-b1970d291811
# ╟─1312525c-1fc0-11eb-2756-5bc3101d2260
# ╠═c4398f9c-1fc4-11eb-0bbb-37f066c6027d
# ╠═7f961bc0-1fc5-11eb-1f18-612aeff0d8df
# ╟─25f92dec-1fc4-11eb-055d-f34deea81d0e
# ╟─fa7e6f7e-2434-11eb-1e61-1b1858bb0988
# ╟─16348b6a-1fc2-11eb-0b9c-65df528db2a1
# ╟─e296c6e8-259c-11eb-1385-53f757f4d585
# ╠═a86f13de-259d-11eb-3f46-1f6fb40020ce
# ╟─3d66bd30-259d-11eb-2694-471fb3a4a7be
# ╠═5f82dec8-259e-11eb-2f4f-4d661f44ef41
# ╟─56b68356-2601-11eb-39a9-5f4b8e580b87
# ╟─7d815988-1fc7-11eb-322a-4509e7128ce3
# ╟─aed8f00e-266b-11eb-156d-8bb09de0dc2b
# ╠═b9f882d8-266b-11eb-2998-75d6539088c7
# ╟─269200ec-259f-11eb-353b-0b73523ef71a
# ╠═e10a9b70-25a0-11eb-2aed-17ed8221c208
# ╟─2dfab366-25a1-11eb-15c9-b3dd9cd6b96c
# ╠═50ea30ba-25a1-11eb-05d8-b3d579f85652
# ╟─51e2e742-25a1-11eb-2511-ab3434eacc3e
# ╟─bade1372-25a1-11eb-35f4-4b43d4e8d156
# ╠═02232964-2603-11eb-2c4c-c7b7e5fed7d1
# ╟─736ed1b6-1fc2-11eb-359e-a1be0a188670
# ╠═49cb5174-1fc3-11eb-3670-c3868c9b0255
# ╟─f3abc83c-1fc7-11eb-1aa8-01ce67c8bdde
# ╠═3d72ab3a-2689-11eb-360d-9b3d829b78a9
# ╟─b6d7a362-1fc8-11eb-03bc-89464b55c6fc
# ╠═1f148d9a-1fc8-11eb-158e-9d784e390b24
# ╟─cf8dca6c-1fc8-11eb-1f89-099e6ba53c22
# ╠═02173c7a-2695-11eb-251c-65efb5b4a45f
# ╟─440271b6-25e8-11eb-26ce-1b80aa176aca
# ╠═cf276892-25e7-11eb-38f0-03f75c90dd9e
# ╟─5b5f25f0-266c-11eb-25d4-17e411c850c9
# ╟─3f823490-266d-11eb-1ba4-d5a23975c335
# ╟─971f401e-266c-11eb-3104-171ae299ef70
# ╠═746aa5bc-266c-11eb-14c9-63ccc313f5de
# ╟─a919d584-2670-11eb-1cf9-2327c8135d6d
# ╠═bfb07a0a-2670-11eb-3938-772499c637b1
# ╟─12cbbab0-2671-11eb-2b1f-038c206e84ce
# ╠═9596c2dc-2671-11eb-36b9-c1af7e5f1089
# ╠═f94a1d56-2671-11eb-2cdc-810a9c7a8a5f
# ╟─4b091fac-2672-11eb-0db8-75457788d85e
# ╟─9cdc5f84-2671-11eb-3c78-e3495bc64d33
# ╠═f688f9f2-2671-11eb-1d71-a57c9817433f
# ╠═049a866e-2672-11eb-29f7-bfea7ad8f572
# ╠═09901de6-2672-11eb-3d50-05b176b729e7
# ╠═aea0d0b4-2672-11eb-231e-395c863827d3
# ╟─9c32db5c-1fc9-11eb-029a-d5d554de1067
# ╠═19957754-252d-11eb-1e0a-930b5208f5ac
# ╠═40f1e7d8-252d-11eb-0549-49ca4e806e16
# ╟─ee1be5dc-252b-11eb-0865-291aa823b9e9
# ╟─06c5139e-252d-11eb-2645-8b324b24c405
# ╠═f2e55166-25ff-11eb-0297-796e97c62b07
# ╟─1ea81214-1fca-11eb-2442-7b0b448b49d6
# ╟─a0ef04b0-25e9-11eb-1110-cde93601f712
# ╟─3e310cf8-25ec-11eb-07da-cb4a2c71ae34
# ╟─d6d1b312-2543-11eb-1cb2-e5b801686ffb
# ╠═378aed18-252b-11eb-0b37-a3b511af2cb5
# ╟─3cbc95ba-2685-11eb-3810-3bf38aa33231
# ╟─68b2a560-2536-11eb-0cc4-27793b4d6a70
# ╟─0e19f82e-2685-11eb-2e99-0d094c1aa520
# ╟─1eabe908-268b-11eb-329b-b35160ec951e
# ╠═1d388372-2695-11eb-3068-7b28a2ccb9ac
# ╟─53c2eaf6-268b-11eb-0899-b91c03713da4
# ╠═06d28052-2531-11eb-39e2-e9613ab0401c
# ╟─4c9173ac-2685-11eb-2129-99071821ebeb
# ╠═736515ba-2685-11eb-38cb-65bfcf8d1b8d
# ╟─8b06b944-268c-11eb-0bfc-8d4dd21e1f02
# ╟─09ce27ca-268c-11eb-0cdd-c9801db876f8
# ╟─298deff4-2676-11eb-2595-e7e22f613ea1
# ╟─2bbf5a70-2676-11eb-1085-7130d4a30443
# ╟─de95efae-2675-11eb-0909-73afcd68fd42
# ╟─c78e02b4-268a-11eb-0af7-f7c7620fcc34
# ╠═d7801e88-2530-11eb-0b93-6f1c78d00eea
# ╠═607058ec-253c-11eb-0fb6-add8cfb73a4f
# ╟─9c1f73e0-268a-11eb-2bf1-216a5d869568
# ╟─11096250-2544-11eb-057b-d7112f20b05c
# ╠═9eb07a6e-2687-11eb-0de3-7bc6aa0eefb0
# ╟─cb15cd88-25ed-11eb-2be4-f31500a726c8
# ╟─232b9bec-2544-11eb-0401-97a60bb172fc
# ╟─36e2dfea-2433-11eb-1c90-bb93ab25b33c
# ╟─36ea4410-2433-11eb-1d98-ab4016245d95
# ╟─36f8c1e8-2433-11eb-1f6e-69dc552a4a07
# ╟─37061f1e-2433-11eb-3879-2d31dc70a771
# ╟─371352ec-2433-11eb-153d-379afa8ed15e
# ╟─372002e4-2433-11eb-0b25-39ce1b1dd3d1
# ╟─372c1480-2433-11eb-3c4e-95a37d51835f
# ╟─3737be8e-2433-11eb-2049-2d6d8a5e4753
# ╟─374522c4-2433-11eb-3da3-17419949defc
# ╟─37552044-2433-11eb-1984-d16e355a7c10


text:
struct Options
    mode::String
    model_type::String
    x::Array{Data,1}
    n_total::Int64
    n_keep::Int64
    n_burn::Int64
    verbose::Bool
    use_hyperprior::Bool
    t_max::Int64
    # MFM options
    gamma::Float64
    log_pk::String
    # DPM options
    alpha_random::Bool
    alpha::Float64
    # Jain-Neal split-merge options
    use_splitmerge::Bool
    n_split::Int64
    n_merge::Int64
    # RJMCMC options
    k_max::Int64
    # Partition distribution values
    a::Float64
    b::Float64
    log_v::Array{Float64,1}
    # Other
    n::Int64
end

struct Result
    options::Options
    t::Array{Int8,1}
    N::Array{Int16,2}
    z::Array{Int8,2}
    theta::Array{Theta,2}
    keepers::Array{Int64,1}
    elapsed_time::Float64
    time_per_step::Float64
end


text:
module Groups

using AbstractAlgebra
import AbstractAlgebra: Group, GroupElem, Ring
import AbstractAlgebra: parent, parent_type, elem_type
import AbstractAlgebra: order, gens, matrix_repr

import Base: length, ==, hash, show, convert, eltype, iterate
import Base: inv, reduce, *, ^, power_by_squaring
import Base: findfirst, findnext
import Base: deepcopy_internal

export elements

using LinearAlgebra
using Markdown

###############################################################################
#
#   ParentType / ObjectType definition
#
###############################################################################

@doc doc"""
    ::GSymbol
> Abstract type which all group symbols of AbstractFPGroups should subtype. Each
> concrete subtype should implement fields:
> * `id` which is the `Symbol` representation/identification of a symbol
> * `pow` which is the (multiplicative) exponent of a symbol.

"""
abstract type GSymbol end

abstract type GWord{T<:GSymbol} <:GroupElem end

@doc doc"""
    W::GroupWord{T} <: GWord{T<:GSymbol} <:GroupElem
> Basic representation of element of a finitely presented group. `W.symbols`
> fieldname contains particular group symbols which multiplied constitute a
> group element, i.e. a word in generators.
> As reduction (inside group) of such word may be time consuming we provide
> `savedhash` and `modified` fields as well:
> hash (used e.g. in the `unique` function) is calculated by reducing the word,
> setting `modified` flag to `false` and computing the hash which is stored in
> `savedhash` field.
> whenever word `W` is changed `W.modified` is set to `false`;
> Future comparisons don't perform reduction (and use `savedhash`) as long as
> `modified` flag remains `false`.

"""
mutable struct GroupWord{T} <: GWord{T}
   symbols::Vector{T}
   savedhash::UInt
   modified::Bool
   parent::Group

   function GroupWord{T}(symbols::Vector{T}) where {T}
      return new{T}(symbols, hash(symbols), true)
   end
end

abstract type AbstractFPGroup <: Group end

###############################################################################
#
#   Includes
#
###############################################################################

include("FreeGroup.jl")
include("FPGroups.jl")
include("AutGroup.jl")

include("DirectPower.jl")
include("WreathProducts.jl")

###############################################################################
#
#   Type and parent object methods
#
###############################################################################

parent(w::GWord{T}) where {T<:GSymbol} = w.parent

###############################################################################
#
#   ParentType / ObjectType constructors
#
###############################################################################

GroupWord(s::T) where {T<:GSymbol} = GroupWord{T}(T[s])
GroupWord{T}(s::T) where {T<:GSymbol} = GroupWord{T}(T[s])
GroupWord(w::GroupWord{T}) where {T<:GSymbol} = w
convert(::Type{GroupWord{T}}, s::T) where {T<:GSymbol} = GroupWord{T}(T[s])

###############################################################################
#
#   Basic manipulation
#
###############################################################################

function hash(W::GWord, h::UInt)
    W.modified && reduce!(W)
    return xor(W.savedhash, h)
end

# WARNING: Due to specialised (constant) hash function of GWords this one is actually necessary!
function deepcopy_internal(W::T, dict::IdDict) where {T<:GWord}
    G = parent(W)
    return G(T(deepcopy(W.symbols)))
end

length(W::GWord) = sum([length(s) for s in W.symbols])

function deleteids!(W::GWord)
    to_delete = Int[]
    for i in 1:length(W.symbols)
        if W.symbols[i].pow == 0
           push!(to_delete, i)
        end
    end
    deleteat!(W.symbols, to_delete)
end

function freereduce!(W::GWord)
    reduced = true
    for i in 1:length(W.symbols) - 1
        if W.symbols[i].pow == 0
            continue
        elseif W.symbols[i].id == W.symbols[i+1].id
            reduced = false
            p1 = W.symbols[i].pow
            p2 = W.symbols[i+1].pow

            W.symbols[i+1] = change_pow(W.symbols[i], p1 + p2)
            W.symbols[i] = change_pow(W.symbols[i], 0)
        end
    end
    deleteids!(W)
    return reduced
end

function reduce!(W::GWord)
    if length(W) < 2
        deleteids!(W)
    else
        reduced = false
        while !reduced
            reduced = freereduce!(W)
        end
    end

    W.savedhash = hash(W.symbols, hash(typeof(W), hash(parent(W), zero(UInt))))
    W.modified = false

    return W
end

@doc doc"""
    reduce(W::GWord)
> performs reduction/simplification of a group element (word in generators).
> The default reduction is the free group reduction, i.e. consists of
> multiplying adjacent symbols with the same `id` identifier and deleting the
> identity elements from `W.symbols`.
> More specific procedures should be dispatched on `GWord`s type parameter.

"""
reduce(W::GWord) = reduce!(deepcopy(W))

@doc doc"""
    gens(G::AbstractFPGroups)
> returns vector of generators of `G`, as its elements.

"""
gens(G::AbstractFPGroup) = [G(g) for g in G.gens]

###############################################################################
#
#   String I/O
#
###############################################################################

@doc doc"""
    show(io::IO, W::GWord)
> The actual string produced by show depends on the eltype of `W.symbols`.

"""
function show(io::IO, W::GWord)
    if length(W) == 0
        print(io, "(id)")
    else
        join(io, [string(s) for s in W.symbols], "*")
    end
end

function show(io::IO, s::T) where {T<:GSymbol}
   if s.pow == 1
      print(io, string(s.id))
   else
      print(io, string((s.id))*"^$(s.pow)")
   end
end

###############################################################################
#
#   Comparison
#
###############################################################################

function (==)(W::GWord, Z::GWord)
    parent(W) == parent(Z) || return false

    W.modified && reduce!(W)
    Z.modified && reduce!(Z)

    if W.savedhash != Z.savedhash
        return false
    end

    return W.symbols == Z.symbols
end

function (==)(s::GSymbol, t::GSymbol)
   s.pow == t.pow || return false
   s.pow ==  0 && return true
   s.id == t.id || return false
   return true
end

###############################################################################
#
#   Binary operators
#
###############################################################################

function AbstractAlgebra.mul!(out::GWord, x::GWord, y::GWord; reduced::Bool=true)
    resize!(out.symbols, length(x.symbols)+length(y.symbols))
    for i in eachindex(x.symbols)
        out.symbols[i] = x.symbols[i]
    end
    for i in eachindex(y.symbols)
        out.symbols[length(x.symbols)+i] = y.symbols[i]
    end
    if reduced
        reduce!(out)
    end
    return out
end

function r_multiply!(W::GWord, x; reduced::Bool=true)
    if length(x) > 0
        append!(W.symbols, x)
    end
    if reduced
        reduce!(W)
    end
    return W
end

function l_multiply!(W::GWord, x; reduced::Bool=true)
    if length(x) > 0
        prepend!(W.symbols, x)
    end
    if reduced
        reduce!(W)
    end
    return W
end

r_multiply(W::GWord, x; reduced=true) =
    r_multiply!(deepcopy(W),x, reduced=reduced)
l_multiply(W::GWord, x; reduced=true) =
    l_multiply!(deepcopy(W),x, reduced=reduced)

(*)(W::GWord, Z::GWord) = r_multiply(W, Z.symbols)
(*)(W::GWord, s::GSymbol) = r_multiply(W, [s])
(*)(s::GSymbol, W::GWord) = l_multiply(W, [s])

function power_by_squaring(W::GWord, p::Integer)
    if p < 0
        return power_by_squaring(inv(W), -p)
    elseif p == 0
        return parent(W)()
    elseif p == 1
        return W
    elseif p == 2
        return W*W
    end
    W = deepcopy(W)
    t = trailing_zeros(p) + 1
    p >>= t
    while (t -= 1) > 0
        r_multiply!(W, W.symbols)
    end
    Z = deepcopy(W)
    while p > 0
        t = trailing_zeros(p) + 1
        p >>= t
        while (t -= 1) >= 0
            r_multiply!(W, W.symbols)
        end
        r_multiply!(Z, W.symbols)
    end
    return Z
end

(^)(x::GWord, n::Integer) = power_by_squaring(x,n)

###############################################################################
#
#   Inversion
#
###############################################################################

function inv(W::T) where {T<:GWord}
    if length(W) == 0
        return W
    else
        G = parent(W)
        w = T(reverse([inv(s) for s in W.symbols]))
        w.modified = true
        return G(w)
    end
end

###############################################################################
#
#   Replacement of symbols / sub-words
#
###############################################################################

issubsymbol(s::GSymbol, t::GSymbol) =
    s.id == t.id && (0 ≤ s.pow ≤ t.pow || 0 ≥ s.pow ≥ t.pow)

"""doc
Find the first linear index k>=i such that Z < W.symbols[k:k+length(Z)-1]
"""
function findnext(W::GWord, Z::GWord, i::Int)
    n = length(Z.symbols)
   if n == 0
      return 0
   elseif n == 1
      for idx in i:lastindex(W.symbols)
         if issubsymbol(Z.symbols[1], W.symbols[idx])
            return idx
         end
      end
      return 0
   else
      for idx in i:lastindex(W.symbols) - n + 1
         foundfirst = issubsymbol(Z.symbols[1], W.symbols[idx])
         lastmatch = issubsymbol(Z.symbols[end], W.symbols[idx+n-1])
         if foundfirst && lastmatch
            # middles match:
            if view(Z.symbols, 2:n-1) == view(W.symbols, idx+1:idx+n-2)
               return idx
            end
         end
      end
   end
   return 0
end

findfirst(W::GWord, Z::GWord) = findnext(W, Z, 1)

function replace!(W::GWord, index, toreplace::GWord, replacement::GWord; check=true)
   n = length(toreplace.symbols)
   if n == 0
      return reduce!(W)

   elseif n == 1
      if check
         @assert issubsymbol(toreplace.symbols[1], W.symbols[index])
      end

      first = change_pow(W.symbols[index],
         W.symbols[index].pow - toreplace.symbols[1].pow)
      last = change_pow(W.symbols[index], 0)

   else
      if check
         @assert issubsymbol(toreplace.symbols[1], W.symbols[index])
         @assert W.symbols[index+1:index+n-2] == toreplace.symbols[2:end-1]
         @assert issubsymbol(toreplace.symbols[end], W.symbols[index+n-1])
      end

      first = change_pow(W.symbols[index],
      W.symbols[index].pow - toreplace.symbols[1].pow)
      last = change_pow(W.symbols[index+n-1],
      W.symbols[index+n-1].pow - toreplace.symbols[end].pow)
   end

   replacement = first * replacement * last
   splice!(W.symbols, index:index+n-1, replacement.symbols)
   return reduce!(W)
end

function replace(W::GWord, index, toreplace::GWord, replacement::GWord)
    replace!(deepcopy(W), index, toreplace, replacement)
end

function replace_all!(W::T,subst_dict::Dict{T,T}) where {T<:GWord}
    modified = false
    for toreplace in reverse!(sort!(collect(keys(subst_dict)), by=length))
        replacement = subst_dict[toreplace]
        i = findfirst(W, toreplace)
        while i ≠ 0
            modified = true
            replace!(W,i,toreplace, replacement)
            i = findnext(W, toreplace, i)
        end
    end
    return modified
end

function replace_all(W::T, subst_dict::Dict{T,T}) where {T<:GWord}
   W = deepcopy(W)
   replace_all!(W, subst_dict)
   return W
end

###############################################################################
#
#   Misc
#
###############################################################################

function generate_balls(S::AbstractVector{T}, Id::T=parent(first(S))();
        radius=2, op=*) where T<:GroupElem
    sizes = Int[]
    B = [Id]
    for i in 1:radius
        BB = [op(i,j) for (i,j) in Base.product(B,S)]
        B = unique([B; vec(BB)])
        push!(sizes, length(B))
    end
    return B, sizes
end

function generate_balls(S::AbstractVector{T}, Id::T=one(parent(first(S)));
        radius=2, op=*) where {T<:NCRingElem}
    sizes = Int[]
    B = [Id]
    for i in 1:radius
        BB = [op(i,j) for (i,j) in Base.product(B,S)]
        B = unique([B; vec(BB)])
        push!(sizes, length(B))
    end
    return B, sizes
end

########### iteration for GFField

length(F::AbstractAlgebra.GFField) = order(F)

function iterate(F::AbstractAlgebra.GFField, s=0)
   if s >= order(F)
      return nothing
   else
      return F(s), s+1
   end
end

eltype(::Type{AbstractAlgebra.GFField{I}}) where I = AbstractAlgebra.gfelem{I}

end # of module Groups


text:
function paint!(r::Robot)
    move_till_border!(r)
    putmarkers_border!(r)
    show!(r)
end

function move_till_border!(r::Robot)
    while isborder(r, Nord) == false
        move!(r, Nord)
    end
    while isborder(r, West) == false
        move!(r, West)
    end
end

function putmarkers_border!(r::Robot)
    while ismarker(r) == false
        while isborder(r, Sud) == false
            putmarker!(r)
            move!(r, Sud)
        end
        if isborder(r, Ost) == false
            putmarker!(r)
            move!(r, Ost)
        end
        while isborder(r, Nord) == false
            putmarker!(r)
            move!(r, Nord)
        end
        if isborder(r, Ost) == false
            putmarker!(r)
            move!(r, Ost)
        end
    end
end


text:
"""
Provide a set of fast math functions
"""
module FastMath
using LinearAlgebra
# using Cuba
export dot, squaredNorm, norm

# include("Yeppp.jl")
# using .Yeppp

# @inline dot(x, y) = Yeppp.dot(x, y)
# @inline squaredNorm(x) = Yeppp.dot(x, x)
# @inline norm(x) = sqrt(Yeppp.dot(x, x))

@inline dot(x, y) = dot(x, y)
@inline squaredNorm(x) = dot(x, x)
@inline norm(x) = sqrt(dot(x, x))

"""
    invsqrt(x)

The Legendary Fast Inverse Square Root
See the following links: [wikipedia](https://en.wikipedia.org/wiki/Fast_inverse_square_root) and [thesis](https://cs.uwaterloo.ca/~m32rober/rsqrt.pdf)
"""
@inline function invsqrt(x::Float64)
    #   y = x
    x2::Float64 = x * 0.5
    i::Int64 = reinterpret(Int64, x)
    # The magic number is for doubles is from https://cs.uwaterloo.ca/~m32rober/rsqrt.pdf
    i = 0x5fe6eb50c7b537a9 - (i >> 1)
    y = reinterpret(Float64, i)
    y = y * (1.5 - (x2 * y * y)) # 1st iteration
    y = y * (1.5 - (x2 * y * y))  # 2nd iteration, this can be removed
    return y
end

@inline function invsqrt(x::Float32)
    x2 = x * 0.5f0
    i::Int32 = reinterpret(Int32, x)  # evil floating point bit level hacking
    i = 0x5f3759df - (i >> 1)
    y = reinterpret(Float32, i)
    y = y * (1.5f0 - (x2 * y * y)) # 1st iteration
    y = y * (1.5f0 - (x2 * y * y))  # 2nd iteration, this can be removed
    return y
end

# function integrator1D(type, f, xmin, xmax, rtol)
#     function integrand(x, g)
#         if xmin != -Inf && xmax != Inf
#             a, b = xmin, xmax
#             g[1] = f(a + (b - a) * x[1]) * (b - a)
#         elseif xmin != -Inf && xmax == Inf
#             a = xmin
#             g[1] = f(a + x[1] / (1 - x[1])) / (1 - x[1])^2
#         elseif xmin == -Inf && xmax != Inf
#             b = xmax
#             g[1] = - f(b + 1 - 1 / x[1]) / (x[1])^2
#         else
#             y = x[1]
#             denorm = (1 - y) * y
#             g[1] = - f((2y - 1) / denorm) * (2y^2 - 2y + 1) / denorm^2
#         end
#     end

#     if type == :vegas
#         result, err = Cuba.vegas(integrand, rtol=rtol)
#         return result[1], err[1]
#     elseif type == :cuhre
#         result, err = Cuba.cuhre(integrand, rtol=rtol)
#         return result[1], err[1]
#     else
#         @error "Not implemented!"
#     end
# end

end


text:
function makeconvergencegraphs(E::Estimation;
  convergencewindows::Vector{Int} = CONVERGENCE_WINDOWS,
  prefix::String = "",  outgraphpath::String=OUT_GRAPH_PATH, display::Bool = false)

  local plotnames::Vector{String} = Vector{String}()
  local plots::Vector{PlotContainer} = Vector{PlotContainer}()
  local subplots::Vector{PlotContainer} = Vector{PlotContainer}()

  df::DataFrame = paramdf(E::Estimation)
  df[:zetaGP_abs] = (abs).(df[:zetaGP])
  df[:SGP_abs] = (abs).(df[:SGP])
  paramnames::Vector{Symbol} = [(Symbol).(names(E.Θ)); :zetaGP_abs; :SGP_abs]

  sort!(df, :pass)
  for f::Symbol ∈ [:sigma2G, :SG, :zeta2G, :zeta2P, :zetaGP_abs]
    for convergencewindow ∈ convergencewindows
      dflong::DataFrame = melt(df[df[:pass] .≤ convergencewindow,:], :pass, paramnames)
      df2Gplot::SubDataFrame = view(dflong, (s->s∈[:pass, f]).(dflong[:variable]),:)
      push!(subplots,
        plot(df2Gplot,
        x=:pass, color=:variable,y=:value,
        Scale.y_log10,
        #Guide.title(plotnames[end]),
        Guide.ylabel("Value"),Guide.xlabel(nothing),
        Guide.Theme(key_position=:none, line_width=0.2mm),
        Geom.line, Coord.cartesian(ymin=-9, ymax = -3, xmin=1, xmax=convergencewindow)))
      end

      push!(plotnames, "$(f)_convergence")
      push!(plots, vstack([subplots[1], subplots[2], subplots[3], subplots[4]]))

      for i ∈ 1:length(plots) #write the graphs
        draw(PNG("$outgraphpath\\$(prefix)_$(plotnames[i]).png", 9inch, 7inch), plots[i])
      end

      plots = Vector{PlotContainer}()
      subplots = Vector{PlotContainer}()
      plotnames = Vector{String}()
    end

# NOTE TODO: Uncomment THIS!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  for convergencewindow ∈ convergencewindows
    subdf::SubDataFrame  = view(df, df[:pass] .≤ convergencewindow, :)
    push!(subplots,
      plot(subdf,
      x=:sigma2G, y=:SG,
      Scale.y_log10, Scale.x_log10,
      Guide.title("Minvar sample and pop var i=$convergencewindow"),
      Guide.xlabel("Pop var"), Guide.ylabel("Minvar sample var"),
      Gadfly.Theme(line_width=0.2mm),
      Coord.cartesian(ymin=-9, ymax = -3, xmin=-9, xmax = -3),
      Geom.path))
    push!(subplots,
      plot(subdf,
      x=:zeta2G, y=:zeta2P,
      Scale.y_log10, Scale.x_log10,
      Guide.title("Var of Var i=$convergencewindow"),
      Guide.xlabel("Minvar Var of Sample Var"), Guide.ylabel("Test Port. Var of Sample Var"),
      Gadfly.Theme(line_width=0.2mm),
      Coord.cartesian(ymin=-8, ymax = -2, xmin=-8, xmax = -2),
      Geom.path))
  end

  push!(plotnames, "sample_pop_var_convergence")
  push!(plots,
    gridstack([subplots[1] subplots[3]; subplots[5] subplots[7]]))

  push!(plotnames, "var_of_var_convergence")
  push!(plots,
    gridstack([subplots[2] subplots[4]; subplots[6] subplots[8]]))

  for i ∈ 1:length(plots) #write the graphs
    draw(PNG("$outgraphpath\\$(prefix)_$(plotnames[i]).png", 9inch, 7inch), plots[i])
  end

  println("Convergence graphs drawn")
  return nothing
end

function makeconvergencegraphs(;outpath::String = OUT_PATH,
  convergencewindows::Vector{Int} = CONVERGENCE_WINDOWS,
  outname::String = OUT_ESTIMATE_NAME,
  prefix="")::Nothing

  makeconvergencegraphs(deserialize("$outpath\\$(prefix)_$(outname).jls"),
    prefix=prefix, convergencewindows=convergencewindows)

  return
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
export horizontalbezier, perpendicularbezier
export bracket, arrow, arcarrow
export reversed, concat
export rotate
export centeredin

Base.:*(b::Bezier, r::Real) = Bezier(r * b.from, r * b.c1, r * b.c2, r * b.to)
Base.:+(b::Bezier, p::Point) = Bezier(b.from + p, b.c1 + p, b.c2 + p, b.to + p)
Base.:-(b::Bezier, p::Point) = Bezier(b.from - p, b.c1 - p, b.c2 - p, b.to - p)
rotate(b::Bezier, ang::Angle) = Bezier(rotate.((b.from, b.c1, b.c2, b.to), ang)...)

function horizontalbezier(p1::Point, p2::Point, strength=1)
    diff = p1 → p2
    Bezier(
        p1,
        Point(p1.x + 0.5 * strength * diff.x, p1.y),
        Point(p2.x - 0.5 * strength * diff.x, p2.y),
        p2
    )
end

function perpendicularbezier(l1::Line, l2::Line, fraction1=0.5, fraction2=0.5; strength=1, reverse1=false, reverse2=false)
    fpoint1 = fraction(l1, fraction1)
    fpoint2 = fraction(l2, fraction2)
    perpendicular1 = perpendicular(l1, reverse1)
    perpendicular2 = perpendicular(l2, reverse2)
    dist = magnitude(fpoint1 → fpoint2)
    c1 = fpoint1 + perpendicular1 * 0.5dist * strength
    c2 = fpoint2 + perpendicular2 * 0.5dist * strength
    Bezier(fpoint1, c1, c2, fpoint2)
end

needed_attributes(::Type{Bezier}) = (Linewidth, Stroke, Linestyle, Fill, Visible)

# function Path(closed::Bool, segments::Vararg{<: BezierSegment, N}) where N
#     Path([segments...], closed)
# end

Base.:+(path::Path, p::Point) = Path(path.commands .+ p)
Base.:-(path::Path, p::Point) = Path(path.commands .- p)
Base.:*(path::Path, r::Real) = Path(path.commands .* r)
Base.:*(r::Real, path::Path) = Path(path.commands .* r)

function bracket(p1::Point, p2::Point, widthscale::Real = 0.1, innerstrength=1, outerstrength=1; flip=false)
    l = Line(p1, p2)
    perp1 = perpendicular(l, flip)
    tipdist = widthscale * length(l)
    tipvec = perp1 * tipdist
    tip = fraction(l, 0.5) + tipvec
    bez1 = Bezier(p1, p1 + 0.5tipvec * outerstrength, tip - 0.5tipvec * innerstrength, tip)
    bez2 = Bezier(tip, tip - 0.5tipvec * innerstrength, p2 + 0.5tipvec * outerstrength, p2)
    Path([bez1, bez2], false)
end

needed_attributes(::Type{Path}) = (Visible, Linewidth, Stroke, Linestyle, Fill)

function arrow(from::Point, to::Point, tiplength, tipwidth, shaftwidthback, shaftwidthfront, tipretraction = 0)
    vector = from → to
    tipconnection = to - normalize(vector) * tiplength
    tipconnection_retracted = to - normalize(vector) * tiplength * (1-tipretraction)
    ortholeft = normalize(rotate(vector, deg(90)))
    tipleft = tipconnection + 0.5tipwidth * ortholeft
    tipright = tipconnection - 0.5tipwidth * ortholeft
    endleft = from + 0.5shaftwidthback * ortholeft
    endright = from - 0.5shaftwidthback * ortholeft
    tipconnleft = tipconnection_retracted + 0.5shaftwidthfront * ortholeft
    tipconnright = tipconnection_retracted - 0.5shaftwidthfront * ortholeft
    Path([
        Move(endleft),
        Lineto(tipconnleft),
        Lineto(tipleft),
        Lineto(to),
        Lineto(tipright),
        Lineto(tipconnright),
        Lineto(endright),
        Close()
    ])
end

arrow(from::Point, to::Point, tipsize, shaftsize, tipretraction = 0) =
    arrow(from, to, tipsize, tipsize, shaftsize, shaftsize, tipretraction)

arrow(l::Line, args...) = arrow(l.from, l.to, args...)

reversed(b::Bezier) = Bezier(b.to, b.c2, b.c1, b.from)
reversed(b::Path) = Path(reverse!(reversed.(b.segments)), b.closed)

function arcarrow(from::Point, to::Point, radiusfraction::Real, tiplength::Real, tipwidth::Real, tipretraction::Real=0)
    arc = Arc(from, to, radiusfraction)
    alength = arclength(arc)

    tipconnection = fraction(arc, 1 - tiplength / alength)
    tipconnection_retracted = fraction(arc, 1 - tiplength * (1 - tipretraction) / alength)

    arc_retracted = lengthen(arc, (arc.end_angle - arc.start_angle) * -tiplength * (1 - tipretraction) / alength)
    arcbezier = Base.convert(Path, arc_retracted)

    ortholeft = normalize(rotate(tipconnection → to, deg(90)))
    tipleft = tipconnection + 0.5tipwidth * ortholeft
    tipright = tipconnection - 0.5tipwidth * ortholeft

    segments = BezierSegment[
        arcbezier.segments...,
        Line(tipconnection_retracted, tipright),
        Line(tipright, to),
        Line(to, tipleft),
        Line(tipleft, tipconnection_retracted),
        reversed(arcbezier).segments...
    ]
    Path(segments, false)
end

function concat(bp::Path, paths...)
    segments = vcat((b.segments for b in (bp, paths...))...)
    Path(
        segments,
        false
    )
end

function rotate(p::Path, ang::Angle)
    Path(rotate.(p.commands, ang))
end

continuecommands(l::Line) = PathCommand[Lineto(stop(l))]
continuecommand(b::Bezier) = CurveTo(b.c1, b.c2, b.to)
continuecommands(b::Bezier) = PathCommand[continuecommand(b)]
function continuecommands(a::Arc)
    # https://stackoverflow.com/a/44829356/2279303

    function uptoquarter(center, p1, p2)
        x1, y1 = p1.xy
        x4, y4 = p2.xy
        xc, yc = center.xy

        ax = x1 - xc
        ay = y1 - yc
        bx = x4 - xc
        by = y4 - yc
        q1 = ax * ax + ay * ay
        q2 = q1 + ax * bx + ay * by
        k2 = 4/3 * (√(2 * q1 * q2) - q2) / (ax * by - ay * bx)

        x2 = xc + ax - k2 * ay
        y2 = yc + ay + k2 * ax
        x3 = xc + bx + k2 * by
        y3 = yc + by - k2 * bx
        Bezier(p1, P(x2, y2), P(x3, y3), p2)
    end

    quarterfractions = (a.end_angle - a.start_angle) / deg(90)
    nsegments=ceil(abs(quarterfractions))

    points = [fraction(a, i / nsegments) for i in 0:nsegments]
    segments = [uptoquarter(a.center, p1, p2) for (p1, p2) in zip(points[1:end-1], points[2:end])]
    PathCommand[continuecommand.(segments)...]
end

function Path(segments::Vector{T}) where T <:Union{Line, Bezier, Arc}

    laststop = nothing
    commands = PathCommand[]

    i = 1
    while i <= length(segments)
        seg = segments[i]

        if i == 1 || !(start(seg) ≈ laststop)
            push!(commands, Move(start(seg)))
        end

        append!(commands, continuecommands(seg))
        laststop = stop(seg)

        i += 1

    end

    if laststop ≈ start(segments[1])
        # replace first value so start and stop are exactly the same
        commands[1] = Move(laststop)
        push!(commands, Close())
    end

    Path(commands)

end

function Path(svg::AbstractString)

    # args = split(svg, r"[\s,]+", keepempty=false)
    # args = split(svg, r"((?<=[a-zA-Z])(?=\d)|(?<=\d)(?=[a-zA-Z])|([\s,]+)|((?<=\d)(?=\-))|((?<=[a-zA-Z])(?=\-)))|((?<=\.\d+))", keepempty=false)
    args = [e.match for e in eachmatch(r"([a-zA-Z])|(\-?\d*\.?\d+)", svg)]

    i = 1

    commands = PathCommand[]
    lastcomm = nothing
    lastp() = isnothing(lastcomm) ? O : commands[end].p

    while i <= length(args)

        comm = args[i]

        # command letter is omitted, use last command
        if isnothing(match(r"[a-zA-Z]", comm))
            comm = lastcomm
            i -= 1
        end

        if comm == "M"
            x, y = parse.(Float64, args[i+1:i+2])
            push!(commands, Move(P(x, y)))
            i += 3
        elseif comm == "m"
            x, y = parse.(Float64, args[i+1:i+2])
            push!(commands, Move(P(x, y) + lastp()))
            i += 3
        elseif comm == "L"
            x, y = parse.(Float64, args[i+1:i+2])
            push!(commands, Lineto(P(x, y)))
            i += 3
        elseif comm == "l"
            x, y = parse.(Float64, args[i+1:i+2])
            push!(commands, Lineto(P(x, y) + lastp()))
            i += 3
        elseif comm == "H"
            x = parse(Float64, args[i+1])
            push!(commands, Lineto(P(x, lastp().y)))
            i += 2
        elseif comm == "h"
            x = parse(Float64, args[i+1])
            push!(commands, Lineto(X(x) + lastp()))
            i += 2
        elseif comm == "Z"
            push!(commands, Close())
            i += 1
        elseif comm == "z"
            push!(commands, Close())
            i += 1
        elseif comm == "C"
            x1, y1, x2, y2, x3, y3 = parse.(Float64, args[i+1:i+6])
            push!(commands, CurveTo(P(x1, y1), P(x2, y2), P(x3, y3)))
            i += 7
        elseif comm == "c"
            x1, y1, x2, y2, x3, y3 = parse.(Float64, args[i+1:i+6])
            l = lastp()
            push!(commands, CurveTo(P(x1, y1) + l, P(x2, y2) + l, P(x3, y3) + l))
            i += 7
        elseif comm == "S"
            x1, y1, x2, y2 = parse.(Float64, args[i+1:i+4])
            prev_control = commands[end-1]
            prev = commands[end]
            reflected = prev.p + (prev_control → prev.p)
            push!(commands, CurveTo(reflected, P(x1, y1), P(x2, y2)))
            i += 5
        elseif comm == "s"
            x1, y1, x2, y2 = parse.(Float64, args[i+1:i+4])
            prev_control = commands[end-1]
            prev = commands[end]
            reflected = prev.p + (prev_control.p → prev.p)
            l = lastp()
            push!(commands, CurveTo(reflected, P(x1, y1) + l, P(x2, y2) + l))
            i += 5
        else
            for c in commands
                println(c)
            end
            error("Parsing $comm not implemented.")
        end

        lastcomm = comm

    end

    Path(commands)

end

function Path(ps::AbstractArray{Point}, closed::Bool=false)
    segs = [Move(ps[1]), Lineto.(ps[2:end])...]
    closed && push!(segs, Close())
    Path(segs)
end

function Path(xs::AbstractArray, ys::AbstractArray, closed::Bool=false)
    Path(P.(xs, ys), closed)
end

function bbox(b::Bezier)

    p0 = [b.from.xy...]
    p1 = [b.c1.xy...]
    p2 = [b.c2.xy...]
    p3 = [b.to.xy...]

    mi = min.(p0, p3)
    ma = max.(p0, p3)

    c = -p0 + p1
    b =  p0 - 2p1 + p2
    a = -p0 + 3p1 - 3p2 + 1p3

    h = b.*b - a.*c

    if h[1] > 0
        h[1] = sqrt(h[1])
        t = (-b[1] - h[1]) / a[1]
        if t > 0 && t < 1
            s = 1.0-t
            q = s*s*s*p0[1] + 3.0*s*s*t*p1[1] + 3.0*s*t*t*p2[1] + t*t*t*p3[1]
            mi[1] = min(mi[1],q)
            ma[1] = max(ma[1],q)
        end
        t = (-b[1] + h[1])/a[1]
        if t>0 && t<1
            s = 1.0-t
            q = s*s*s*p0[1] + 3.0*s*s*t*p1[1] + 3.0*s*t*t*p2[1] + t*t*t*p3[1]
            mi[1] = min(mi[1],q)
            ma[1] = max(ma[1],q)
        end
    end

    if h[2]>0.0
        h[2] = sqrt(h[2])
        t = (-b[2] - h[2])/a[2]
        if t>0.0 && t<1.0
            s = 1.0-t
            q = s*s*s*p0[2] + 3.0*s*s*t*p1[2] + 3.0*s*t*t*p2[2] + t*t*t*p3[2]
            mi[2] = min(mi[2],q)
            ma[2] = max(ma[2],q)
        end
        t = (-b[2] + h[2])/a[2]
        if t>0.0 && t<1.0
            s = 1.0-t
            q = s*s*s*p0[2] + 3.0*s*s*t*p1[2] + 3.0*s*t*t*p2[2] + t*t*t*p3[2]
            mi[2] = min(mi[2],q)
            ma[2] = max(ma[2],q)
        end
    end

    BBox(P(mi...), P(ma...))
end

function bbox(bboxes::Array{BBox})
    mi = [bboxes[1].from.xy...]
    ma = [bboxes[1].to.xy...]

    for i in 2:length(bboxes)
        mi .= min.(mi, bboxes[i].from.xy)
        ma .= max.(ma, bboxes[i].to.xy)
    end
    BBox(P(mi...), P(ma...))
end

function segments(p::Path)
    segs = []
    last = nothing
    for c in p.commands
        if typeof(c) <: Move
            last = c.p
        elseif typeof(c) <: Lineto
            push!(segs, Line(last, c.p))
            last = c.p
        elseif typeof(c) <: CurveTo
            push!(segs, Bezier(last, c.c1, c.c2, c.p))
            last = c.p
        elseif typeof(c) <: Close
            # what doing
        end
    end
    segs
end

bbox(p::Path) = bbox(bbox.(segments(p)))

function centeredin(p::Path, squarelength::Real)
    bb = bbox(p)
    r = Rect(bb)
    factor = squarelength / max(r.width, r.height)
    (p - r.center) * factor
end

Base.:+(m::Move, p::Point) = Move(m.p + p)
Base.:+(l::Lineto, p::Point) = Lineto(l.p + p)
Base.:+(c::CurveTo, p::Point) = CurveTo(c.c1 + p, c.c2 + p, c.p + p)
Base.:+(c::Close, p::Point) = c

Base.:-(m::Move, p::Point) = Move(m.p - p)
Base.:-(l::Lineto, p::Point) = Lineto(l.p - p)
Base.:-(c::CurveTo, p::Point) = CurveTo(c.c1 - p, c.c2 - p, c.p - p)
Base.:-(c::Close, p::Point) = c

Base.:*(m::Move, s::Real) = Move(m.p * s)
Base.:*(l::Lineto, s::Real) = Lineto(l.p * s)
Base.:*(c::CurveTo, s::Real) = CurveTo(c.c1 * s, c.c2 * s, c.p * s)
Base.:*(c::Close, s::Real) = c

Base.:*(s::Real, m::Move) = m * s
Base.:*(s::Real, l::Lineto) = l * s
Base.:*(s::Real, c::CurveTo) = c * s
Base.:*(s::Real, c::Close) = c * s

rotate(m::Move, a::Angle) = Move(rotate(m.p, a))
rotate(l::Lineto, a::Angle) = Lineto(rotate(l.p, a))
rotate(c::CurveTo, a::Angle) = CurveTo(rotate.([c.c1, c.c2, c.p], a)...)
rotate(c::Close, a::Angle) = c


text:
#include("PrivateMultiplicativeWeights.jl")

#using PrivateMultiplicativeWeights

using LinearAlgebra

"""
return the lower bound for the "Query andWorkload Error"
    according to the article:
    "Ecient Batch Query Answering Under Dierential Privacy"
"""
function getLowerBound(queries::HistogramQueries, eps, delta)
    m,n = size(queries.queries)
    P_eps_delta = (2*log2(2/delta))/(eps^2)
    singular_values = svd(queries.queries).S
    ss_svd = (1/n)*(sum(singular_values)^2)

    ss_svd*P_eps_delta
end

"""
return the lower bound for the "Query andWorkload Error"
    according to the article:
    "Ecient Batch Query Answering Under Dierential Privacy"
"""
function getLowerBound(queries::RangeQueries, eps, delta)

    query_matrix = zeros(Int64, length(queries.intervals), queries.domain)

    for (index, interval) in enumerate(queries.intervals)
        query_matrix[index,:] = get_query_vector_from_interval(interval, queries.domain)
    end

    m,n = size(query_matrix)
    P_eps_delta = (2*log2(2/delta))/(eps^2)
    singular_values = svd(query_matrix).S
    ss_svd = (1/n)*(sum(singular_values)^2)

    ss_svd*P_eps_delta
end

max_stars_repo_name: giladroyz/MWEM-project

max_issues_repo_name: giladroyz/MWEM-project

max_forks_repo_name: giladroyz/MWEM-project