text:
const __bodyfunction__ = Dict{Method,Any}()

# Find keyword "body functions" (the function that contains the body
# as written by the developer, called after all missing keyword-arguments
# have been assigned values), in a manner that doesn't depend on
# gensymmed names.
# `mnokw` is the method that gets called when you invoke it without
# supplying any keywords.
function __lookup_kwbody__(mnokw::Method)
    function getsym(arg)
        isa(arg, Symbol) && return arg
        @assert isa(arg, GlobalRef)
        return arg.name
    end

    f = get(__bodyfunction__, mnokw, nothing)
    if f === nothing
        fmod = mnokw.module
        # The lowered code for `mnokw` should look like
        #   %1 = mkw(kwvalues..., #self#, args...)
        #        return %1
        # where `mkw` is the name of the "active" keyword body-function.
        ast = Base.uncompressed_ast(mnokw)
        if isa(ast, Core.CodeInfo) && length(ast.code) >= 2
            callexpr = ast.code[end-1]
            if isa(callexpr, Expr) && callexpr.head == :call
                fsym = callexpr.args[1]
                if isa(fsym, Symbol)
                    f = getfield(fmod, fsym)
                elseif isa(fsym, GlobalRef)
                    if fsym.mod === Core && fsym.name === :_apply
                        f = getfield(mnokw.module, getsym(callexpr.args[2]))
                    elseif fsym.mod === Core && fsym.name === :_apply_iterate
                        f = getfield(mnokw.module, getsym(callexpr.args[3]))
                    else
                        f = getfield(fsym.mod, fsym.name)
                    end
                else
                    f = missing
                end
            else
                f = missing
            end
        else
            f = missing
        end
        __bodyfunction__[mnokw] = f
    end
    return f
end

function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.eyeM)),NamedTuple{(:like,),Tuple{Array{Int8,2}}},typeof(eyeM),Int64,Int64})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.freqspaceM)),NamedTuple{(:dim,),Tuple{Int64}},typeof(freqspaceM),Int64})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.ndgridM)),NamedTuple{(:dim,),Tuple{Int64}},typeof(ndgridM),StepRange{Int64,Int64}})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.onesM)),NamedTuple{(:like,),Tuple{Array{Float64,2}}},typeof(onesM),Array{Int64,1}})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.onesM)),NamedTuple{(:like,),Tuple{Array{Int32,2}}},typeof(onesM),Symbol,Int64})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.onesM)),NamedTuple{(:like,),Tuple{Array{Int64,2}}},typeof(onesM),Int64})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.randM)),NamedTuple{(:like,),Tuple{Array{Float64,2}}},typeof(randM),Array{Int64,1}})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.randM)),NamedTuple{(:like,),Tuple{Array{Int32,2}}},typeof(randM),Symbol,Int64})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.randM)),NamedTuple{(:like,),Tuple{Array{Int64,2}}},typeof(randM),Int64})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.sortM)),NamedTuple{(:I,),Tuple{Bool}},typeof(sortM),Array{Int64,1}})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.sortM)),NamedTuple{(:I,),Tuple{Bool}},typeof(sortM),Array{Int64,2},Int64,Symbol})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.sortperm)),NamedTuple{(:dims,),Tuple{Int64}},typeof(sortperm),Array{Int64,2}})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.zerosM)),NamedTuple{(:like,),Tuple{Array{Float64,2}}},typeof(zerosM),Array{Int64,1}})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.zerosM)),NamedTuple{(:like,),Tuple{Array{Int32,2}}},typeof(zerosM),Symbol,Int64})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.zerosM)),NamedTuple{(:like,),Tuple{Array{Int64,2}}},typeof(zerosM),Int64})
    Base.precompile(Tuple{typeof(blkdiagM),Array{Float64,2},Array{Float64,2},Vararg{Array{Float64,2},N} where N})
    Base.precompile(Tuple{typeof(blkdiagM),Symbol,Array{Float64,2},Vararg{Array{Float64,2},N} where N})
    Base.precompile(Tuple{typeof(catM),Int64,Array{Float64,2},Vararg{Array{Float64,2},N} where N})
    Base.precompile(Tuple{typeof(catM),Int64,Array{Float64,3},Vararg{Array{Float64,3},N} where N})
    Base.precompile(Tuple{typeof(clcM)})
    Base.precompile(Tuple{typeof(diagM),Array{Int64,1}})
    Base.precompile(Tuple{typeof(diagM),Array{Int64,2}})
    Base.precompile(Tuple{typeof(doubleM),Float64})
    Base.precompile(Tuple{typeof(eyeM),Array{Int64,1}})
    Base.precompile(Tuple{typeof(eyeM),Int64,Int64})
    Base.precompile(Tuple{typeof(eyeM),Int64})
    Base.precompile(Tuple{typeof(eyeM),Tuple{Int64,Int64}})
    Base.precompile(Tuple{typeof(eyeM),Type{T} where T,Array{Int64,1}})
    Base.precompile(Tuple{typeof(eyeM),Type{T} where T,Int64,Int64})
    Base.precompile(Tuple{typeof(eyeM),Type{T} where T,Tuple{Int64,Int64}})
    Base.precompile(Tuple{typeof(falseM),Array{Int64,1}})
    Base.precompile(Tuple{typeof(falseM),Int64})
    Base.precompile(Tuple{typeof(falseM),Symbol,Int64})
    Base.precompile(Tuple{typeof(falseM),Tuple{Int64,Int64}})
    Base.precompile(Tuple{typeof(flipM),Array{Int64,1}})
    Base.precompile(Tuple{typeof(flipM),Array{Int64,2}})
    Base.precompile(Tuple{typeof(flipM),String})
    Base.precompile(Tuple{typeof(int16M),Int64})
    Base.precompile(Tuple{typeof(int32M),Int64})
    Base.precompile(Tuple{typeof(int64M),Int64})
    Base.precompile(Tuple{typeof(int8M),Int64})
    Base.precompile(Tuple{typeof(iscolumnM),Array{Float64,1}})
    Base.precompile(Tuple{typeof(iscolumnM),Array{Float64,2}})
    Base.precompile(Tuple{typeof(iscolumnM),Array{Int64,1}})
    Base.precompile(Tuple{typeof(iscolumnM),Array{Int64,2}})
    Base.precompile(Tuple{typeof(iscolumnM),Array{String,1}})
    Base.precompile(Tuple{typeof(iscolumnM),Array{String,2}})
    Base.precompile(Tuple{typeof(iscolumnM),Int64})
    Base.precompile(Tuple{typeof(iscolumnM),String})
    Base.precompile(Tuple{typeof(iscolumnM),Symbol,Array{Float64,1}})
    Base.precompile(Tuple{typeof(iscolumnM),Symbol,Array{Float64,2}})
    Base.precompile(Tuple{typeof(iscolumnM),Symbol,Array{Int64,1}})
    Base.precompile(Tuple{typeof(iscolumnM),Symbol,Array{Int64,2}})
    Base.precompile(Tuple{typeof(iscolumnM),Symbol,Array{String,1}})
    Base.precompile(Tuple{typeof(iscolumnM),Symbol,Int64})
    Base.precompile(Tuple{typeof(iscolumnM),Symbol,String})
    Base.precompile(Tuple{typeof(ismatrixM),Array{Float64,1}})
    Base.precompile(Tuple{typeof(ismatrixM),Array{Float64,2}})
    Base.precompile(Tuple{typeof(ismatrixM),Array{Int64,1}})
    Base.precompile(Tuple{typeof(ismatrixM),Array{Int64,2}})
    Base.precompile(Tuple{typeof(ismatrixM),Array{String,1}})
    Base.precompile(Tuple{typeof(ismatrixM),Array{String,2}})
    Base.precompile(Tuple{typeof(ismatrixM),String})
    Base.precompile(Tuple{typeof(ismatrixM),Symbol,Array{Float64,1}})
    Base.precompile(Tuple{typeof(ismatrixM),Symbol,Array{Float64,2}})
    Base.precompile(Tuple{typeof(ismatrixM),Symbol,Array{Int64,1}})
    Base.precompile(Tuple{typeof(ismatrixM),Symbol,Array{Int64,2}})
    Base.precompile(Tuple{typeof(ismatrixM),Symbol,Array{String,1}})
    Base.precompile(Tuple{typeof(ismatrixM),Symbol,String})
    Base.precompile(Tuple{typeof(isrowM),Array{Float64,1}})
    Base.precompile(Tuple{typeof(isrowM),Array{Float64,2}})
    Base.precompile(Tuple{typeof(isrowM),Array{Int64,1}})
    Base.precompile(Tuple{typeof(isrowM),Array{Int64,2}})
    Base.precompile(Tuple{typeof(isrowM),Array{String,1}})
    Base.precompile(Tuple{typeof(isrowM),Array{String,2}})
    Base.precompile(Tuple{typeof(isrowM),Int64})
    Base.precompile(Tuple{typeof(isrowM),String})
    Base.precompile(Tuple{typeof(isrowM),Symbol,Array{Float64,1}})
    Base.precompile(Tuple{typeof(isrowM),Symbol,Array{Float64,2}})
    Base.precompile(Tuple{typeof(isrowM),Symbol,Array{Int64,1}})
    Base.precompile(Tuple{typeof(isrowM),Symbol,Array{Int64,2}})
    Base.precompile(Tuple{typeof(isrowM),Symbol,Array{String,1}})
    Base.precompile(Tuple{typeof(isrowM),Symbol,Int64})
    Base.precompile(Tuple{typeof(isrowM),Symbol,String})
    Base.precompile(Tuple{typeof(isscalarM),Array{Float64,3}})
    Base.precompile(Tuple{typeof(isscalarM),Array{String,1}})
    Base.precompile(Tuple{typeof(isscalarM),Array{String,2}})
    Base.precompile(Tuple{typeof(isscalarM),Bool})
    Base.precompile(Tuple{typeof(isscalarM),Char})
    Base.precompile(Tuple{typeof(isscalarM),Int64})
    Base.precompile(Tuple{typeof(isscalarM),Symbol,Array{Float64,3}})
    Base.precompile(Tuple{typeof(isscalarM),Symbol,Array{Int64,1}})
    Base.precompile(Tuple{typeof(isscalarM),Symbol,Array{String,1}})
    Base.precompile(Tuple{typeof(isscalarM),Symbol,Array{String,2}})
    Base.precompile(Tuple{typeof(isscalarM),Symbol,Bool})
    Base.precompile(Tuple{typeof(isscalarM),Symbol,Char})
    Base.precompile(Tuple{typeof(isscalarM),Symbol,Int64})
    Base.precompile(Tuple{typeof(isscalarM),Symbol,String})
    Base.precompile(Tuple{typeof(isvectorM),Array{Float64,1}})
    Base.precompile(Tuple{typeof(isvectorM),Array{Float64,2}})
    Base.precompile(Tuple{typeof(isvectorM),Array{Int64,1}})
    Base.precompile(Tuple{typeof(isvectorM),Array{Int64,2}})
    Base.precompile(Tuple{typeof(isvectorM),Array{String,1}})
    Base.precompile(Tuple{typeof(isvectorM),Array{String,2}})
    Base.precompile(Tuple{typeof(isvectorM),Int64})
    Base.precompile(Tuple{typeof(isvectorM),String})
    Base.precompile(Tuple{typeof(isvectorM),Symbol,Array{Float64,1}})
    Base.precompile(Tuple{typeof(isvectorM),Symbol,Array{Float64,2}})
    Base.precompile(Tuple{typeof(isvectorM),Symbol,Array{Int64,1}})
    Base.precompile(Tuple{typeof(isvectorM),Symbol,Array{Int64,2}})
    Base.precompile(Tuple{typeof(isvectorM),Symbol,Array{String,1}})
    Base.precompile(Tuple{typeof(isvectorM),Symbol,Int64})
    Base.precompile(Tuple{typeof(isvectorM),Symbol,String})
    Base.precompile(Tuple{typeof(lengthM),Array{Any,2}})
    Base.precompile(Tuple{typeof(lengthM),Array{Int64,2}})
    Base.precompile(Tuple{typeof(linspaceM),Int64,Int64,Int64})
    Base.precompile(Tuple{typeof(linspaceM),Int64,Int64})
    Base.precompile(Tuple{typeof(linspaceM),Symbol,Int64,Int64,Int64})
    Base.precompile(Tuple{typeof(linspaceM),Symbol,Int64,Int64})
    Base.precompile(Tuple{typeof(logspaceM),Int64,Int64,Int64})
    Base.precompile(Tuple{typeof(logspaceM),Int64,Int64,Symbol})
    Base.precompile(Tuple{typeof(logspaceM),Int64,Int64})
    Base.precompile(Tuple{typeof(meshgridM),StepRange{Int64,Int64}})
    Base.precompile(Tuple{typeof(meshgridM),UnitRange{Int64},UnitRange{Int64},UnitRange{Int64}})
    Base.precompile(Tuple{typeof(ndgridM),UnitRange{Int64},UnitRange{Int64},UnitRange{Int64}})
    Base.precompile(Tuple{typeof(onesM),Array{Int64,1}})
    Base.precompile(Tuple{typeof(onesM),Int64})
    Base.precompile(Tuple{typeof(onesM),Symbol,Int64})
    Base.precompile(Tuple{typeof(onesM),Tuple{Int64,Int64}})
    Base.precompile(Tuple{typeof(onesM),Type{Int32},Array{Int64,1}})
    Base.precompile(Tuple{typeof(onesM),Type{T} where T,Int64,Int64})
    Base.precompile(Tuple{typeof(onesM),Type{T} where T,Tuple{Int64,Int64}})
    Base.precompile(Tuple{typeof(randM),Array{Int64,1}})
    Base.precompile(Tuple{typeof(randM),Int64})
    Base.precompile(Tuple{typeof(randM),Symbol,Int64})
    Base.precompile(Tuple{typeof(randM),Tuple{Int64,Int64}})
    Base.precompile(Tuple{typeof(randM),Type{Int32},Array{Int64,1}})
    Base.precompile(Tuple{typeof(randM),Type{T} where T,Int64,Int64})
    Base.precompile(Tuple{typeof(randM),Type{T} where T,Tuple{Int64,Int64}})
    Base.precompile(Tuple{typeof(repelemM),Array{Int64,2},Int64})
    Base.precompile(Tuple{typeof(repmatM),Array{Int64,1},Array{Int64,1}})
    Base.precompile(Tuple{typeof(repmatM),Array{Int64,1},Int64,Int64})
    Base.precompile(Tuple{typeof(repmatM),Array{Int64,1},Int64})
    Base.precompile(Tuple{typeof(repmatM),Array{Int64,1},Tuple{Int64,Int64}})
    Base.precompile(Tuple{typeof(repmatM),Int64,Int64,Int64})
    Base.precompile(Tuple{typeof(singleM),Float64})
    Base.precompile(Tuple{typeof(sizeM),Array{Float64,3},Int64,Int64})
    Base.precompile(Tuple{typeof(sizeM),Array{Int64,2},Int64})
    Base.precompile(Tuple{typeof(sizeM),Symbol,Array{Float64,3},Int64,Int64})
    Base.precompile(Tuple{typeof(sizeM),Symbol,Array{Int64,2},Int64})
    Base.precompile(Tuple{typeof(sizeM),Symbol,Array{Int64,2}})
    Base.precompile(Tuple{typeof(sortM),Array{Int64,1}})
    Base.precompile(Tuple{typeof(sortM),Array{Int64,2}})
    Base.precompile(Tuple{typeof(sortM),Array{Integer,3},Int64})
    Base.precompile(Tuple{typeof(squeezeM),Array{Float64,3}})
    Base.precompile(Tuple{typeof(transposeM),Array{Int64,2}})
    Base.precompile(Tuple{typeof(transposeM),Symbol,Array{Int64,2}})
    Base.precompile(Tuple{typeof(trueM),Array{Int64,1}})
    Base.precompile(Tuple{typeof(trueM),Int64})
    Base.precompile(Tuple{typeof(trueM),Symbol,Int64})
    Base.precompile(Tuple{typeof(trueM),Tuple{Int64,Int64}})
    Base.precompile(Tuple{typeof(uint16M),Int64})
    Base.precompile(Tuple{typeof(uint32M),Int64})
    Base.precompile(Tuple{typeof(uint64M),Int64})
    Base.precompile(Tuple{typeof(uint8M),Int64})
    Base.precompile(Tuple{typeof(zerosM),Array{Int64,1}})
    Base.precompile(Tuple{typeof(zerosM),Int64})
    Base.precompile(Tuple{typeof(zerosM),Symbol,Int64})
    Base.precompile(Tuple{typeof(zerosM),Tuple{Int64,Int64}})
    Base.precompile(Tuple{typeof(zerosM),Type{Int32},Array{Int64,1}})
    Base.precompile(Tuple{typeof(zerosM),Type{T} where T,Int64,Int64,Int64})
    Base.precompile(Tuple{typeof(zerosM),Type{T} where T,Int64,Int64})
    Base.precompile(Tuple{typeof(zerosM),Type{T} where T,Tuple{Int64,Int64}})
    isdefined(MatLang, Symbol("#32#33")) && Base.precompile(Tuple{getfield(MatLang, Symbol("#32#33")),Array{Int64,1}})
    let fbody = try __lookup_kwbody__(which(sortperm, (Array{Int64,2},))) catch missing end
        if !ismissing(fbody)
            precompile(fbody, (Int64,Bool,typeof(sortperm),Array{Int64,2},))
        end
    end
end

max_stars_repo_name: juliamatlab/MatLang

max_issues_repo_name: juliamatlab/MatLang

max_forks_repo_name: juliamatlab/MatLang


text:
"""
    Wirtinger(primal::Union{Number,AbstractDifferential},
              conjugate::Union{Number,AbstractDifferential})

Returns a `Wirtinger` instance representing the complex differential:

```
df = ∂f/∂z * dz + ∂f/∂z̄ * dz̄
```

where `primal` corresponds to `∂f/∂z * dz` and `conjugate` corresponds to `∂f/∂z̄ * dz̄`.

The two fields of the returned instance can be accessed generically via the
[`wirtinger_primal`](@ref) and [`wirtinger_conjugate`](@ref) methods.
"""
struct Wirtinger{P,C} <: AbstractDifferential
    primal::P
    conjugate::C
    function Wirtinger(primal::Union{Number,AbstractDifferential},
                       conjugate::Union{Number,AbstractDifferential})
        return new{typeof(primal),typeof(conjugate)}(primal, conjugate)
    end
end

wirtinger_primal(x::Wirtinger) = x.primal
wirtinger_primal(x) = x

wirtinger_conjugate(x::Wirtinger) = x.conjugate
wirtinger_conjugate(::Any) = Zero()

extern(x::Wirtinger) = throw(ArgumentError("`Wirtinger` cannot be converted to an external type."))

Base.Broadcast.broadcastable(w::Wirtinger) = Wirtinger(broadcastable(w.primal),
                                                       broadcastable(w.conjugate))

Base.iterate(x::Wirtinger) = (x, nothing)
Base.iterate(::Wirtinger, ::Any) = nothing

# TODO: define `conj` for` `Wirtinger`
Base.conj(x::Wirtinger) = throw(MethodError(conj, x))

function refine_differential(::Type{<:Union{<:Real, AbstractArray{<:Real}}}, w::Wirtinger)
    return wirtinger_primal(w) + wirtinger_conjugate(w)
end

ext: jl


text:
_opmap = Dict{Tuple{UInt32, Symbol}, Tuple{SOperator, UInt32}}()

clearops() = empty!(_opmap)

"""
    @simulatable

Make any operator simulatable by providing:
- A composite (struct) type corresponding to the hardware evaluation
- An function that computes the floating-point value

# Examples

The following makes a simulatable stochastic add operation.
```julia
@simulatable(SSignedAdder, +(x::SBit, y::SBit) = x.value + y.value)
```

The following makes a simulatable stochastic divide operation.
```julia
@simulatable(SSignedDivider,
function /(x::SBit, y::SBit)
    if y.value <= 0
        error("SBit only supports divisors > 0.")
    end

    x.value / y.value
end)
```

Some operators require arguments to the constructor. You can optionally specify these
to `@simulatable` as a third argument.
```julia
@simulatable(SSignedMatMultiplier,
    *(x::VecOrMat{SBit}, y::VecOrMat{SBit}) = float.(x) * float.(y), (size(x, 1), size(y, 2)))
```
"""
macro simulatable(optype, fdef, args = ())
    def = splitdef(fdef)
    opsym = "$(def[:name]) -> $optype"
    idstrs = []
    opargs = []
    for arg in splitarg.(def[:args])
        push!(idstrs, :(_getidstr($(arg[1]))))
        push!(opargs, :($(arg[1])))
    end
    hashkey = :(vcat($(idstrs...)))

    newbody = quote
        key = (hashn!($hashkey, 4), Symbol($opsym))
        if haskey(_opmap, key)
            op, id = _opmap[key]
            value = $(def[:body])

            SBit(op($(opargs...)), value, id)
        else
            op = $(optype)($args...)
            id = _genid()
            value = $(def[:body])
            outbit = SBit(op($(opargs...)), value, id)
            _opmap[key] = (op, id)

            return outbit
        end
    end
    def[:body] = newbody

    return MacroTools.combinedef(def)
end


text:
function parse_ages(io)
    lines = readlines(io)
    return parse_ints_from_row(lines[1]; sep=",")
end

function ages_to_counters(ages)
    return [sum(ages .== i) for i in 0:8]
end

function step_ages(ages)
    counters = ages_to_counters(ages)
    return step_counters(counters)
end

function step_counters(counters)
    n_new = counters[1]
    for i in 1:8
        counters[i] = counters[i+1]
    end
    counters[9] = 0
    if n_new > 0
        counters[7] += n_new
        counters[9] = n_new
    end
    return counters
end

function step_n_ages(days, ages)
    counters = ages_to_counters(ages)
    return step_n_counters(days, counters)
end

function step_n_counters(days, counters)
    for i in 1:days
        counters = step_counters(counters)
        # println(counters)
    end
    return counters
end

function day06a(counters)
    counters = step_n_counters(80, counters)
    return sum(counters)
end

function day06b(counters)
    counters = step_n_counters(256, counters)
    return sum(counters)
end

function day06(io)
    ages = parse_ages(io)
    counters = ages_to_counters(ages)
    return (day06a(copy(counters)), day06b(copy(counters)))
end


text:
## Basic understanding of logit regression.

using DataFrames
using CSV
using GLM
using FreqTables

obs = CSV.File("./logit_regression/logit_regression_data.csv", header=true) |> DataFrame

## summary of methods.
## want to model variable p in range (0, 1). This is really hard to do.
## probability is linked with odds ratio = prob of success / prob of failure
## taking the log of this is the log odds ratio.
## how do we transform p ∈ (0, 1) to odds ∈ (0, ∞)
## then take the log of odds to map to (-∞, ∞)
## This transformation is called logit transformation. The other common choice is the probit transformation,

## so regression of p is
## log(p/(1-p)) =  linear predictors

## now lets see what this means in the sample data.
# for the purpose of illustration.
# The data set has 200 observations and the outcome variable used will be hon, indicating if a student is in
# an honors class or not. So our p = prob(hon=1).

# Model 1, no predictors
ols = glm(@formula(hon ~ 1), obs, Binomial(), LogitLink())
coef(ols) # so the intercept here is the log odds -1.12
# taking the exponential gives us the odds ratio. I.e. what are the odds of getting into a honors class
exp.(coef(ols))
exp.(coef(ols)) ./ (1 .+ exp.(coef(ols))) # we can get p back as well

# we can verify all of this by looking at a frequnecy table
# lets look at the FreqTables
#freqtable(obs, :hon)

# Model 2, a single predictor: female
ols = glm(@formula(hon ~ female), obs, Binomial(), LogitLink())
coef(ols)   ## gets the log odds
exp.(coef(ols)) ## gets the odds ratio
## and can even calculate the prob given whether person is male or female

# here theintercept of -1.471 is the log odds for males since male is the reference group (female = 0).
# so this means the odds ratio is 0.23

# Using the frequency for males, we can confirm this: log(.23) = -1.47.

# Model 3, using math
ols = glm(@formula(hon ~ math), obs, Binomial(), LogitLink())
coef(ols)   ## gets the log odds
exp.(coef(ols)) ## gets the odds ratio
## and can even calculate the prob given math score.

# Model 3, using combination
ols = glm(@formula(hon ~ math + female + read), obs, Binomial(), LogitLink())
coef(ols)   ## gets the log odds
exp.(coef(ols)) ## gets the odds ratio
## and can even calculate the prob given a different combination of predictors

# This fitted model says that, holding math and reading at a fixed value, the odds of getting into an
# honors class for females (female = 1)over the odds of getting into an honors class for males (female = 0)
# is exp(.979948) = 2.66

# Model 4, interaction term
ols = glm(@formula(hon ~ math + female + female*math), obs, Binomial(), LogitLink())
coef(ols)   ## gets the log odds
exp.(coef(ols)) ## gets the odds ratio

max_stars_repo_name: affans/notebooks

max_issues_repo_name: affans/notebooks

max_forks_repo_name: affans/notebooks


text:
using GNSSTools
using JLD

directory = "/home/sjbilardi/Dropbox/Apps/Overleaf/MastersThesis/"
doppler_file_name = string(directory, "ch1_doppler_curves.jld")
doppler_raw_gps, doppler_rate_raw_gps, sig_freq = load(doppler_file_name,
                                                       "doppler_raw_gps",
                                                       "doppler_rate_raw_gps",
                                                       "freq")
doppler_raw_iridium, doppler_rate_raw_iridium = load(doppler_file_name,
                                                     "doppler_raw_iridium",
                                                     "doppler_rate_raw_iridium")
doppler_raw_starlink, doppler_rate_raw_starlink = load(doppler_file_name,
                                                       "doppler_raw_starlink",
                                                       "doppler_rate_raw_starlink")
doppler_raw_oneweb, doppler_rate_raw_oneweb = load(doppler_file_name,
                                                   "doppler_raw_oneweb",
                                                   "doppler_rate_raw_oneweb")

Ts = [1e-3 2e-3 5e-3 10e-3 20e-3 50e-3 100e-3]
constellation_dopplers = [doppler_raw_gps,
                          doppler_raw_iridium,
                          doppler_raw_starlink,
                          doppler_raw_oneweb]
constellation_names = ["GPS            ",
                       "Iridium        ",
                       "Starlink       ",
                       "OneWeb         "]
header = "Constellation"
for i in 1:length(Ts)
    global header = string(header, "\tT=$(floor(Int, Ts[i]*1000))")
end
println(header)
println("----------------------------------------------------------------------")
for i in 1:length(constellation_names)
    row = constellation_names[i]
    dopplers = constellation_dopplers[i]
    fd_range = max(maximum(dopplers), abs(minimum(dopplers)))
    for j in 1:length(Ts)
        Δf = 1/Ts[j]
        M = 2*ceil(Int, fd_range/Δf) + 1
        row = string(row, "\t$(M)")
    end
    println(row)
end
println("----------------------------------------------------------------------")

println("\n\n\nConstellation\tMax f_d (kHz)")
println("---------------------------------")
for i in 1:length(constellation_names)
    row = constellation_names[i]
    dopplers = constellation_dopplers[i]
    fd_range = max(maximum(dopplers), abs(minimum(dopplers)))
    fd_range = round(fd_range/1e3, digits=1)
    row = string(row, "\t$(fd_range)kHz")
    println(row)
end
println("------------------------------")


text:
#####################
# Probabilistic circuits
#####################
abstract type ProbΔNode{O} <: DecoratorΔNode{O} end
abstract type ProbLeafNode{O} <: ProbΔNode{O} end
abstract type ProbInnerNode{O} <: ProbΔNode{O} end

mutable struct ProbLiteral{O} <: ProbLeafNode{O}
    origin::O
    data
    bit::Bool
    ProbLiteral(n) = new{node_type(n)}(n, nothing, false)
end

mutable struct Prob⋀{O} <: ProbInnerNode{O}
    origin::O
    children::Vector{<:ProbΔNode{<:O}}
    data
    bit::Bool
    Prob⋀(n, children) = begin
        new{node_type(n)}(n, convert(Vector{ProbΔNode{node_type(n)}},children), nothing, false)
    end
end

mutable struct Prob⋁{O} <: ProbInnerNode{O}
    origin::O
    children::Vector{<:ProbΔNode{<:O}}
    log_thetas::Vector{Float64}
    data
    bit::Bool
    Prob⋁(n, children) = new{node_type(n)}(n, convert(Vector{ProbΔNode{node_type(n)}},children), some_vector(Float64, length(children)), nothing, false)
end

const ProbΔ{O} = AbstractVector{<:ProbΔNode{<:O}}

Base.eltype(::Type{ProbΔ{O}}) where {O} = ProbΔNode{<:O}

#####################
# traits
#####################

import LogicCircuits.GateType # make available for extension
import LogicCircuits.node_type

@inline GateType(::Type{<:ProbLiteral}) = LiteralGate()
@inline GateType(::Type{<:Prob⋀}) = ⋀Gate()
@inline GateType(::Type{<:Prob⋁}) = ⋁Gate()

@inline node_type(::ProbΔNode) = ProbΔNode

#####################
# constructors and conversions
#####################

const ProbCache = Dict{ΔNode, ProbΔNode}

function ProbΔ2(circuit::Δ)::ProbΔ
    node2dag(ProbΔ2(circuit[end]))
end

function ProbΔ2(circuit::ΔNode)::ProbΔNode
    f_con(n) = error("Cannot construct a probabilistic circuit from constant leafs: first smooth and remove unsatisfiable branches.")
    f_lit(n) = ProbLiteral(n)
    f_a(n, cn) = Prob⋀(n, cn)
    f_o(n, cn) = Prob⋁(n, cn)
    foldup_aggregate(circuit, f_con, f_lit, f_a, f_o, ProbΔNode{node_type(circuit)})
end

function ProbΔ(circuit::Δ, cache::ProbCache = ProbCache())

    sizehint!(cache, length(circuit)*4÷3)

    pc_node(::LiteralGate, n::ΔNode) = ProbLiteral(n)
    pc_node(::ConstantGate, n::ΔNode) = error("Cannot construct a probabilistic circuit from constant leafs: first smooth and remove unsatisfiable branches.")

    pc_node(::⋀Gate, n::ΔNode) = begin
        children = map(c -> cache[c], n.children)
        Prob⋀(n, children)
    end

    pc_node(::⋁Gate, n::ΔNode) = begin
        children = map(c -> cache[c], n.children)
        Prob⋁(n, children)
    end

    map(circuit) do node
        pcn = pc_node(GateType(node), node)
        cache[node] = pcn
        pcn
    end
end

#####################
# methods
#####################

import LogicCircuits: literal, children # make available for extension

@inline literal(n::ProbLiteral)::Lit  = literal(n.origin)
@inline children(n::ProbInnerNode) = n.children

num_parameters(n::Prob⋁) = num_children(n)
num_parameters(c::ProbΔ) = sum(n -> num_parameters(n), ⋁_nodes(c))

"Return the first origin that is a probabilistic circuit node"
prob_origin(n::DecoratorΔNode)::ProbΔNode = origin(n, ProbΔNode)

"Return the first origin that is a probabilistic circuit"
prob_origin(c::DecoratorΔ)::ProbΔ = origin(c, ProbΔNode)

function estimate_parameters2(pc::ProbΔ, data::XData{Bool}; pseudocount::Float64)
    Logical.pass_up_down2(pc, data)
    w = (data isa PlainXData) ? nothing : weights(data)
    estimate_parameters_cached2(pc, w; pseudocount=pseudocount)
end

function estimate_parameters_cached2(pc::ProbΔ, w; pseudocount::Float64)
    flow(n) = Float64(sum(sum(n.data)))
    children_flows(n) = sum.(map(c -> c.data[1] .& n.data[1], children(n)))

    if issomething(w)
        flow_w(n) = sum(Float64.(n.data[1]) .* w)
        children_flows_w(n) = sum.(map(c -> Float64.(c.data[1] .& n.data[1]) .* w, children(n)))
        flow = flow_w
        children_flows = children_flows_w
    end

    estimate_parameters_node2(n::ProbΔNode) = ()
    function estimate_parameters_node2(n::Prob⋁)
        if num_children(n) == 1
            n.log_thetas .= 0.0
        else
            smoothed_flow = flow(n) + pseudocount
            uniform_pseudocount = pseudocount / num_children(n)
            n.log_thetas .= log.((children_flows(n) .+ uniform_pseudocount) ./ smoothed_flow)
            @assert isapprox(sum(exp.(n.log_thetas)), 1.0, atol=1e-6) "Parameters do not sum to one locally"
            # normalize away any leftover error
            n.log_thetas .- logsumexp(n.log_thetas)
        end
    end

    foreach(estimate_parameters_node2, pc)
end

function log_likelihood_per_instance2(pc::ProbΔ, data::XData{Bool})
    Logical.pass_up_down2(pc, data)
    log_likelihood_per_instance_cached(pc, data)
end

function log_likelihood_per_instance_cached(pc::ProbΔ, data::XData{Bool})
    log_likelihoods = zeros(num_examples(data))
    indices = some_vector(Bool, num_examples(data))::BitVector
    for n in pc
         if n isa Prob⋁ && num_children(n) != 1 # other nodes have no effect on likelihood
            foreach(n.children, n.log_thetas) do c, log_theta
                indices = n.data[1] .& c.data[1]
                view(log_likelihoods, indices::BitVector) .+=  log_theta # see MixedProductKernelBenchmark.jl
            end
         end
    end
    log_likelihoods
end

import LogicCircuits: conjoin_like, disjoin_like, literal_like, copy_node, normalize, replace_node # make available for extension

"Conjoin nodes in the same way as the example"
@inline function conjoin_like(example::ProbΔNode, arguments::Vector)
    if isempty(arguments)
        # @assert false "Probabilistic circuit does not have anonymous true node"
        nothing
    elseif example isa Prob⋀ && children(example) == arguments
        example
    else
        n = conjoin_like(origin(example), origin.(arguments))
        Prob⋀(n, arguments)
    end
end

"Disjoin nodes in the same way as the example"
@inline function disjoin_like(example::ProbΔNode, arguments::Vector)
    if isempty(arguments)
        # @assert false "Probabilistic circuit does not have false node"
        nothing
    elseif example isa Prob⋁ && children(example) == arguments
        example
    else
        n = disjoin_like(origin(example), origin.(arguments))
        # normalize parameters
        thetas = zeros(Float64, length(arguments))
        flag = falses(length(arguments))
        for (i, c) in enumerate(arguments)
            ind = findfirst(x -> x == c, children(example))
            if issomething(ind)
                thetas[i] = exp(example.log_thetas[ind])
                flag[i] = true
            end
        end
        if all(flag)
            thetas = thetas / sum(thetas)
        end
        p = Prob⋁(n, arguments)
        p.log_thetas .= log.(thetas)
        p
    end
end

"Construct a new literal node like the given node's type"
@inline literal_like(::ProbΔNode, lit::Lit) = ProbLiteral(lit)

@inline copy_node(n::Prob⋁, cns) = begin
    orig = copy_node(origin(n), origin.(cns))
    p = Prob⋁(orig, cns)
    p.log_thetas .= copy(n.log_thetas)
    p
end

@inline copy_node(n::Prob⋀, cns) = begin
    orig = copy_node(origin(n), origin.(cns))
    Prob⋀(orig, cns)
end

import LogicCircuits.normalize

@inline normalize(n::Prob⋁, old_n::Prob⋁, kept::Union{Vector{Bool}, BitArray}) = begin
     thetas = exp.(old_n.log_thetas[kept])
     n.log_thetas .= log.(thetas / sum(thetas))
end

function estimate_parameters(pc::ProbΔ, data::XBatches{Bool}; pseudocount::Float64)
    estimate_parameters(AggregateFlowΔ(pc, aggr_weight_type(data)), data; pseudocount=pseudocount)
end

function estimate_parameters(afc::AggregateFlowΔ, data::XBatches{Bool}; pseudocount::Float64)
    @assert feature_type(data) == Bool "Can only learn probabilistic circuits on Bool data"
    @assert (afc[end].origin isa ProbΔNode) "AggregateFlowΔ must originate in a ProbΔ"
    collect_aggr_flows(afc, data)
    estimate_parameters_cached(afc; pseudocount=pseudocount)
    afc
end

function estimate_parameters(fc::FlowΔ, data::XBatches{Bool}; pseudocount::Float64)
    @assert feature_type(data) == Bool "Can only learn probabilistic circuits on Bool data"
    @assert (prob_origin(afc[end]) isa ProbΔNode) "FlowΔ must originate in a ProbΔ"
    collect_aggr_flows(fc, data)
    estimate_parameters_cached(origin(fc); pseudocount=pseudocount)
end

 # turns aggregate statistics into theta parameters
function estimate_parameters_cached(afc::AggregateFlowΔ; pseudocount::Float64)
    foreach(n -> estimate_parameters_node(n; pseudocount=pseudocount), afc)
end

estimate_parameters_node(::AggregateFlowΔNode; pseudocount::Float64) = () # do nothing
function estimate_parameters_node(n::AggregateFlow⋁; pseudocount)
    origin = n.origin::Prob⋁
    if num_children(n) == 1
        origin.log_thetas .= 0.0
    else
        smoothed_aggr_flow = (n.aggr_flow + pseudocount)
        uniform_pseudocount = pseudocount / num_children(n)
        origin.log_thetas .= log.( (n.aggr_flow_children .+ uniform_pseudocount) ./ smoothed_aggr_flow )
        @assert isapprox(sum(exp.(origin.log_thetas)), 1.0, atol=1e-6) "Parameters do not sum to one locally: $(exp.(origin.log_thetas)), estimated from $(n.aggr_flow) and $(n.aggr_flow_children). Did you actually compute the aggregate flows?"
        #normalize away any leftover error
        origin.log_thetas .- logsumexp(origin.log_thetas)
    end
end

# compute log likelihood
function compute_log_likelihood(pc::ProbΔ, data::XBatches{Bool})
    compute_log_likelihood(AggregateFlowΔ(pc, aggr_weight_type(data)))
end

# compute log likelihood, reusing AggregateFlowΔ but ignoring its current aggregate values
function compute_log_likelihood(afc::AggregateFlowΔ, data::XBatches{Bool})
    @assert feature_type(data) == Bool "Can only test probabilistic circuits on Bool data"
    collect_aggr_flows(afc, data)
    ll = log_likelihood(afc)
    (afc, ll)
end

# return likelihoods given current aggregate flows.
function log_likelihood(afc::AggregateFlowΔ)
    sum(n -> log_likelihood(n), afc)
end

log_likelihood(::AggregateFlowΔNode) = 0.0
log_likelihood(n::AggregateFlow⋁) = sum(n.origin.log_thetas .* n.aggr_flow_children)

"""
Calculates log likelihood for a batch of fully observed samples.
(Also retures the generated FlowΔ)
"""
function log_likelihood_per_instance(pc::ProbΔ, batch::PlainXData{Bool})
    fc = FlowΔ(pc, num_examples(batch), Bool)
    (fc, log_likelihood_per_instance(fc, batch))
end

function log_proba(pc::ProbΔ, batch::PlainXData{Bool})
    log_likelihood_per_instance(pc, batch)[2]
end

function log_proba(pc::ProbΔ, batch::PlainXData{Int8})
    marginal_log_likelihood_per_instance(pc, batch)[2]
end
"""
Calculate class dominance considering missing values
"""

function dominance_missing(pc::ProbΔ, batch::PlainXData{Int8}, y::Integer)
    dominance_per_instance(pc, batch, y)[2]
end

"""
Calculate log likelihood per instance for batches of samples.
"""
function log_likelihood_per_instance(pc::ProbΔ, batches::XBatches{Bool})::Vector{Float64}
    mapreduce(b -> log_likelihood_per_instance(pc, b)[2], vcat, batches)
end

"""
Calculate log likelihood for a batch of fully observed samples.
(This is for when you already have a FlowΔ)
"""
function log_likelihood_per_instance(fc::FlowΔ, batch::PlainXData{Bool})
    @assert (prob_origin(fc[end]) isa ProbΔNode) "FlowΔ must originate in a ProbΔ"
    pass_up_down(fc, batch)
    log_likelihoods = zeros(num_examples(batch))
    indices = some_vector(Bool, flow_length(fc))::BitVector
    for n in fc
         if n isa DownFlow⋁ && num_children(n) != 1 # other nodes have no effect on likelihood
            origin = prob_origin(n)::Prob⋁
            foreach(n.children, origin.log_thetas) do c, log_theta
                #  be careful here to allow for the Boolean multiplication to be done using & before switching to float arithmetic, or risk losing a lot of runtime!
                # log_likelihoods .+= prod_fast(downflow(n), pr_factors(c)) .* log_theta
                assign_prod(indices, downflow(n), pr_factors(c))
                view(log_likelihoods, indices::BitVector) .+=  log_theta # see MixedProductKernelBenchmark.jl
                # TODO put the lines above in Utils in order to ensure we have specialized types
            end
         end
    end
    log_likelihoods
end

"""
Calculate log likelihood for a batch of samples with partial evidence P(e).
(Also returns the generated FlowΔ)

To indicate a variable is not observed, pass -1 for that variable.
"""
function marginal_log_likelihood_per_instance(pc::ProbΔ, batch::PlainXData{Int8})
    opts = (flow_opts★..., el_type=Float64, compact⋀=false, compact⋁=false)
    fc = UpFlowΔ(pc, num_examples(batch), Float64, opts)
    (fc, marginal_log_likelihood_per_instance(fc, batch))
end

"""
Calculate class dominance for a batch of samples with missing values.

To indicate a variable is not observed, pass -1 for that variable.
"""
function dominance_per_instance(pc::ProbΔ, batch::PlainXData{Int8}, y::Integer)
    opts = (compact⋀=false, compact⋁=false)
    fc = UpFlowΔ(pc, 4, Tuple{Float64,Float64,Float64},opts);
    # opts = (flow_opts★..., el_type=Float64, compact⋀=false, compact⋁=false)
    # fc = UpFlowΔ(pc, num_examples(batch), Float64, opts)
    (fc, dominance_per_instance(fc, batch, y))
end

"""
Calculate log likelihood for a batch of samples with partial evidence P(e).
(If you already have a FlowΔ)

To indicate a variable is not observed, pass -1 for that variable.
"""
function marginal_log_likelihood_per_instance(fc::UpFlowΔ, batch::PlainXData{Int8})
    @assert (prob_origin(fc[end]) isa ProbΔNode) "FlowΔ must originate in a ProbΔ"
    marginal_pass_up(fc, batch)
    pr(fc[end])
end
"""
Calculate class dominance for a batch of samples with missing values.

To indicate a variable is not observed, pass -1 for that variable.
"""
function dominance_per_instance(fc::UpFlowΔ, batch::PlainXData{Int8}, y::Integer)
    @assert (prob_origin(fc[end]) isa ProbΔNode) "FlowΔ must originate in a ProbΔ"
    dominance_pass_up(fc, batch,y)
    # print(pr(fc[end]))
    pr(fc[end])
end

function check_parameter_integrity(circuit::ProbΔ)
    for node in filter(n -> GateType(n) isa Prob⋁, circuit)
        @assert all(θ -> !isnan(θ), node.log_thetas) "There is a NaN in one of the log_thetas"
    end
    true
end

##################
# Sampling from a psdd
##################

"""
Sample from a PSDD without any evidence
"""
function sample(circuit::ProbΔ)::AbstractVector{Bool}
    inst = Dict{Var,Int64}()
    simulate(circuit[end], inst)
    len = length(keys(inst))
    ans = Vector{Bool}()
    for i = 1:len
        push!(ans, inst[i])
    end
    ans
end

# Uniformly sample based on the probability of the items
# and return the selected index
function sample(probs::AbstractVector{<:Number})::Int32
    z = sum(probs)
    q = rand() * z
    cur = 0.0
    for i = 1:length(probs)
        cur += probs[i]
        if q <= cur
            return i
        end
    end
    return length(probs)
end

function simulate(node::ProbLiteral, inst::Dict{Var,Int64})
    if positive(node)
        inst[variable(node.origin)] = 1
    else
        inst[variable(node.origin)] = 0
    end
end

function simulate(node::Prob⋁, inst::Dict{Var,Int64})
    idx = sample(exp.(node.log_thetas))
    simulate(node.children[idx], inst)
end
function simulate(node::Prob⋀, inst::Dict{Var,Int64})
    for child in node.children
        simulate(child, inst)
    end
end

"""
Sampling with Evidence from a psdd.
Internally would call marginal pass up on a newly generated flow circuit.
"""
function sample(circuit::ProbΔ, evidence::PlainXData{Int8})::AbstractVector{Bool}
    opts= (compact⋀=false, compact⋁=false)
    flow_circuit = UpFlowΔ(circuit, 1, Float64, opts)
    marginal_pass_up(flow_circuit, evidence)
    sample(flow_circuit)
end

"""
Sampling with Evidence from a psdd.
Assuming already marginal pass up has been done on the flow circuit.
"""
function sample(circuit::UpFlowΔ)::AbstractVector{Bool}
    inst = Dict{Var,Int64}()
    simulate2(circuit[end], inst)
    len = length(keys(inst))
    ans = Vector{Bool}()
    for i = 1:len
        push!(ans, inst[i])
    end
    ans
end

function simulate2(node::UpFlowLiteral, inst::Dict{Var,Int64})
    if positive(node)
        #TODO I don't think we need these 'grand_origin' parts below
        inst[variable(grand_origin(node))] = 1
    else
        inst[variable(grand_origin(node))] = 0
    end
end

function simulate2(node::UpFlow⋁, inst::Dict{Var,Int64})
    prs = [ pr(ch)[1] for ch in children(node) ]
    idx = sample(exp.(node.origin.log_thetas .+ prs))
    simulate2(children(node)[idx], inst)
end

function simulate2(node::UpFlow⋀, inst::Dict{Var,Int64})
    for child in children(node)
        simulate2(child, inst)
    end
end

##################
# Most Probable Explanation MPE of a psdd
#   aka MAP
##################

@inline function MAP(circuit::ProbΔ, evidence::PlainXData{Int8})::Matrix{Bool}
    MPE(circuit, evidence)
end

function MPE(circuit::ProbΔ, evidence::PlainXData{Int8})::Matrix{Bool}
    # Computing Marginal Likelihood for each node
    fc, lls = marginal_log_likelihood_per_instance(circuit, evidence)

    ans = Matrix{Bool}(zeros(size(evidence.x)))
    active_samples = Array{Bool}(ones( num_examples(evidence) ))

    mpe_simulate(fc[end], active_samples, ans)
    ans
end

"""
active_samples: bool vector indicating which samples are active for this node during mpe
result: Matrix (num_samples, num_variables) indicating the final result of mpe
"""
function mpe_simulate(node::UpFlowLiteral, active_samples::Vector{Bool}, result::Matrix{Bool})
    if positive(node)
        result[active_samples, variable(node)] .= 1
    else
        result[active_samples, variable(node)] .= 0
    end
end
function mpe_simulate(node::UpFlow⋁, active_samples::Vector{Bool}, result::Matrix{Bool})
    prs = zeros( length(node.children), size(active_samples)[1] )
    @simd  for i=1:length(node.children)
        prs[i,:] .= pr(node.children[i]) .+ (node.origin.log_thetas[i])
    end

    max_child_ids = [a[1] for a in argmax(prs, dims = 1) ]
    @simd for i=1:length(node.children)
        ids = Vector{Bool}( active_samples .* (max_child_ids .== i)[1,:] )  # Only active for this child if it was the max for that sample
        mpe_simulate(node.children[i], ids, result)
    end
end
function mpe_simulate(node::UpFlow⋀, active_samples::Vector{Bool}, result::Matrix{Bool})
    for child in node.children
        mpe_simulate(child, active_samples, result)
    end
end


text:
using Plots,Statistics,LaTeXStrings,JLD,LsqFit
function NeighborReturner(network_, i, j)
    neighbors=[]
    if j!=1 && network_[i,j-1]!=0
        push!(neighbors, network_[i,j-1]) # the left neighbor--> neighbor[1]
    end
    if i!=1 && network_[i-1,j]!=0 #if true --> the neighbors[2] is the upper neighbor
        push!(neighbors,network_[i-1,j])
    end
    return neighbors
end
function LabelFinder(a, L)
    while L[a] != L[L[a]]
        a = L[a]
    end
    return L[a]
end
function percolation_check(network_,L,dim)
    first_col=[]
    last_col=[]
    for i in 1:dim
        if network_[i,1]!=0
            push!(first_col,LabelFinder(network_[i,1],L))
        end
    end
    for i in 1:dim
        if network_[i,dim]!=0
            push!(last_col,LabelFinder(network_[i,dim],L))
        end
    end
    return intersect(first_col,last_col)
end
function percolation(dim,p)
    network_=zeros(Int, dim, dim)
    L=[]
    S=[]
    counter=1
    for col in 1:dim
        for row in 1:dim
            if p>rand()
                neighbors=NeighborReturner(network_, row, col)
                if length(neighbors)==0
                    network_[row,col]=counter
                    push!(L, counter)
                    push!(S, 1)
                    #println(S)
                    counter+=1
                elseif length(neighbors)==1
                    network_[row,col]=LabelFinder(neighbors[1], L)
                    S[LabelFinder(neighbors[1], L)]+=1
                elseif length(neighbors)==2 && LabelFinder(neighbors[1],L)==LabelFinder(neighbors[2],L)
                    network_[row,col]=LabelFinder(neighbors[1], L)
                    S[LabelFinder(neighbors[1], L)]+=1
                else
                    S[LabelFinder(neighbors[1], L)]+=1+S[LabelFinder(neighbors[2], L)]
                    S[LabelFinder(neighbors[2], L)]=0
                    network_[row,col]=LabelFinder(neighbors[1], L)
                    L[neighbors[2]]=LabelFinder(neighbors[1], L)
                end
            end
        end
    end
    return network_,L,S
end
function RadiusOfGyration(network_,L,S,dim)
    if length(Set(S)) < 3
        return 0.0
    end
    if findall(x->x==maximum(S),S)[1] ∈ percolation_check(network_,L,dim)
        S[findall(x->x==maximum(S),S)[1]]=0 #the infinit cluster is omitted from S in this way.
    end
    cluster_size=maximum(S)
    S_max=findall(x->x==maximum(S),S)[1]
    numerator=0
    i_nums=[]
    j_nums=[]
    for coordinate_num in findall(x->x==S_max,network_)
        push!(i_nums, coordinate_num[1])
        push!(j_nums, coordinate_num[2])
    end
    i_com=mean(i_nums)
    j_com=mean(j_nums)
    TheFraction_list=[]
    for i in 1:dim
        for j in 1:dim
             if network_[i, j]==S_max
                push!(TheFraction_list,((i-i_com)^2 + (j-j_com)^2))
            end
        end
    end
    RadiusOfGyration=sqrt(mean(TheFraction_list))
    return RadiusOfGyration
end
dim_list=[10,20,40,80,160]
SavedData_std=[]
SavedData_mean=[]
for n in 1:5
    dim=dim_list[n]
    Meanlist = []
    STDlist = []
    probability=[hcat(0:0.03:0.25)...,hcat(0.25:0.02:0.75)...,hcat(0.75:0.03:1)...]
    for p in probability
        xi=[]
        for run_num in 1:1000
            network_,L,S=percolation(dim,p)
            push!(xi,RadiusOfGyration(network_,L,S,dim))
        end
        push!(STDlist, std(xi))
        push!(Meanlist, mean(xi))
    end
    push!(SavedData_std, STDlist)
    push!(SavedData_mean, Meanlist)
end
save("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter4\\4.5\\AllData_HK.jld",
"std_data", SavedData_std,
"mean_data", SavedData_mean)

#creating the figures:
load("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter4\\4.5\\AllData_HK.jld")
#plot withough error bars:
plot(dpi=400)
for i in 1:5
    dim=dim_list[i]
    scatter!(probability, SavedData_mean[i][1:43],label=nothing,markersize=3,c=:black,alpha=0.3)
    plot!(probability, SavedData_mean[i][1:43],label=L"L=%$dim")

end
scatter!(xlabel="P", ylabel=L"\xi", title=L"\xi\_ P")
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter4\\Fig\\4.5_plot.png")
#scatter plot including error bars:
plot(dpi=400)
for i in 1:5
    dim=dim_list[i]
    scatter!(probability, SavedData_mean[i][1:43],yerr=SavedData_std[i][1:43],markersize=3,
    alpha=0.7,label=L"L=%$dim")

end
scatter!(xlabel="P", ylabel=L"\xi", title=L"\xi\_ P")
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter4\\Fig\\4.5_scatter.png")
#getting Pc(∞) and ν:
xdata=[]
ydata=dim_list
for i in 1:5
    Pc_L= probability[findall(x->x==maximum(SavedData_mean[i]),SavedData_mean[i])[1]]
    push!(xdata, Pc_L)

end
println(xdata)
println(ydata)
@.model(x,p)=abs(x-p[1])^(-p[2])
p0=[1.3, 0.59]
fit=curve_fit(model, xdata, ydata, p0)
#result: ν-->0.61 , Pc(∞)-->1.02


text:
module ParticleFilterCS373
using Distributions

world_size = 100
landmarks  = [(20.0, 20.0), (80.0, 80.0), (20.0, 80.0), (80.0, 20.0)]

export Robot, set!, set_noise!, evaluate, sense, move, Gaussian, measurement_prob

type Robot
  x::Float64
  y::Float64
  orientation::Float64
  forward_noise::Float64
  turn_noise::Float64
  sense_noise::Float64
end

function Robot()
  return Robot(rand()*world_size, rand()*world_size, rand()*2*π, 0., 0., 0.)
end

function set!(rob::Robot, new_x, new_y, new_orientation)
  new_x = abs(new_x % world_size)
  new_y = abs(new_y % world_size)

  if(new_x < 0 || new_x > world_size)
    throw(ArgumentError(@sprintf("X coordinate out of bounds: %s", new_x)))
  end
  if(new_y < 0 || new_y > world_size)
    throw(ArgumentError(@sprintf("Y coordinate out of bounds: %s", new_y)))
  end

  rob.x = new_x
  rob.y = new_y
  rob.orientation = new_orientation

  rob
end

function set_noise!(rob::Robot, new_f_noise, new_t_noise, new_s_noise)
  rob.forward_noise = new_f_noise
  rob.turn_noise = new_t_noise
  rob.sense_noise = new_s_noise

  rob
end

function evaluate(r::Robot, p)
  sum = 0.0
  for i = 1:length(p) # calculate mean error
    dx = (p[i].x - r.x + (world_size/2.0)) % world_size - (world_size/2.0)
    dy = (p[i].y - r.y + (world_size/2.0)) % world_size - (world_size/2.0)
    err = sqrt(dx * dx + dy * dy)
    sum += err
  end
  return sum / float(len(p))
end

function sense(rob::Robot)
  Z = Array(Float64, length(landmarks))
  for ((x, y), i) in zip(landmarks, [1:length(landmarks)])
    println(i)
    dist = sqrt((rob.x - x)^2 + (rob.y - y)^2)
    dist += rob.sense_noise == 0 ? 0 : rand(Normal(0, rob.sense_noise), 1)[1]
    Z[i] = dist
  end
  return Z
end

function move(rob::Robot, turn, forward)
  if(forward < 0)
    throw(ArgumentError("Robot cannot move backwards"))
  end

  # turn, and add randomness to the turning command
  turn_noise   = rob.turn_noise == 0 ? 0 : rand(Normal(0, rob.turn_noise), 1)[1]
  orientation  = rob.orientation + turn + turn_noise
  orientation %= 2π

  # move, and add randomness to the motion command
  forward_noise = rob.forward_noise == 0 ? 0 : rand(Normal(0, rob.forward_noise), 1)[1]
  dist = forward + forward_noise
  x = rob.x + cos(orientation) * dist
  y = rob.y + sin(orientation) * dist
  x %= world_size
  y %= world_size

  # set particle
  res = Robot()
  set!(res, x, y, orientation)
  set_noise!(res, rob.forward_noise, rob.turn_noise, rob.sense_noise)

  return res
end

function Gaussian(μ, σ, x)
  return exp(- ((μ - x)^2) / σ^2 / 2 ) / sqrt(2π*σ^2)
end

function measurement_prob(rob::Robot, measurements)
  prob = 1.0
  for ((x, y), m) in zip(landmarks, measurements)
    dist = sqrt((rob.x - x)^2 + (rob.y - y)^2)
    prob *= Gaussian(dist, rob.sense_noise, m)
  end
  return prob
end
end


text:
function solve()
    a, b, c, k = [parse(Int, x) for x in split(readline())]
    if k <= a
        k
    elseif k <= a+b
        a
    else
        a - (k-a-b)
    end
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
# The nth term of the sequence of triangle numbers is given by, tn = ½n(n+1);
# so the first ten triangle numbers are:
#
# 1, 3, 6, 10, 15, 21, 28, 36, 45, 55, ...
#
# By converting each letter in a word to a number corresponding to its
# alphabetical position and adding these values we form a word value. For
# example, the word value for SKY is 19 + 11 + 25 = 55 = t10. If the word value
# is a triangle number then we shall call the word a triangle word.
#
# Using words.txt (right click and 'Save Link/Target As...'), a 16K text file
# containing nearly two-thousand common English words, how many are triangle
# words?

using ProjectEulerSolutions

function p042solution()::Integer

    path = "data/p042_words.txt"

    # Store reasonable set of triangular numbers
    max_triangular_number = 26*20

    trinumbers = Set{Integer}()
    i = 1
    trinumber = 1
    while trinumber < max_triangular_number
        trinumber = i * (i + 1) / 2
        push!(trinumbers, trinumber)
        i += 1
    end

    allwords = Vector{String}()
    if isfile(path)
        open(path) do file
            for row in eachline(file)
                rownames = split(row, ",")
                append!(allwords, [n[2:end-1] for n in rownames]) # Remove quotes
            end
        end
    end

    triangular_words = 0
    for word in allwords
        if sum(collect(word) .- 'A' .+ 1) in trinumbers
            triangular_words += 1
        end
    end

    return triangular_words
end

p042 = Problems.Problem(p042solution)

Problems.benchmark(p042)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
module HerbSWIPL

using Base: UInt16, String
using Julog

include("structs.jl")
include("swipl_bindings.jl")
include("swipl_interface.jl")

export Swipl
export start, stop, cleanup, is_initialised, asserta, assertz, retract, resolve, swipl_resolve, register_foreign

end # module


text:
function createDataFrame()
    DataFrame(E=[0], c1=[0], c2=[0], c3=[0], c4=[0], c5=[0], c6=[0], c7=[0], c8=[0], c9=[0], c10=[0])
end

function exportData(dir, data, network, nodes)
    #Save log
    CSV.write("Data/$dir/data.csv", data)

    #Save nodes table
    nodes_df = DataFrame(ID = Int[], Vote = Int[], σ_1 = Int[], σ_2 = Int[], σ_3 = Int[],
        σ_4 = Int[], σ_5 = Int[], E = Int[])

    for i in 1:length(nodes)
        push!(nodes_df, (nodes[i].id, nodes[i].vote[], nodes[i].values[1], nodes[i].values[2],
            nodes[i].values[3], nodes[i].values[4], nodes[i].values[5], nodes[i].energy[]))
    end

    CSV.write("Data/$dir/nodes.csv", nodes_df)

    #Exporting the graph
    savegraph("Data/$dir/graph.net", network, "Network", GraphIO.NET.NETFormat())
end

function exportNetwork(dir, network, nodes, step)
    #Save nodes table
    nodes_df = DataFrame(ID = Int[], Vote = Int[], σ_1 = Int[], σ_2 = Int[], σ_3 = Int[],
        σ_4 = Int[], σ_5 = Int[], E = Int[])

    for i in 1:length(nodes)
        push!(nodes_df, (nodes[i].id, nodes[i].vote[], nodes[i].values[1], nodes[i].values[2],
            nodes[i].values[3], nodes[i].values[4], nodes[i].values[5], nodes[i].energy[]))
    end

    CSV.write("Data/$dir/Network/nodes_$step.csv", nodes_df)

    #Exporting the graph
    savegraph("Data/$dir/Network/graph_$step.net", network, "Network", GraphIO.NET.NETFormat())
end

function exportEquilibriumLog(T_step, equilibriumTime, dir)
    eq_df = DataFrame(T = Int[], eqTime = Int[])

    for i in 1:length(equilibriumTime)
        push!(eq_df, (equilibriumTime[i], i*T_step))
    end

    CSV.write("Data/$dir/Equilibrium_log.csv", eq_df)
end

max_stars_repo_name: luca-mm/Bachelor-Thesis

max_issues_repo_name: luca-mm/Bachelor-Thesis

max_forks_repo_name: luca-mm/Bachelor-Thesis


text:
# ---
# title: ElasticDistortion
# cover: elasticdistortion.gif
# id: op_elastic
# ---

# Smoothed random distortion

using Augmentor
using ImageShow, ImageCore
using Random

img_in = testpattern(RGB, ratio=0.5)

mosaicview(
    img_in,
    augment(img_in, ElasticDistortion(15,15,0.1)),
    augment(img_in, ElasticDistortion(10,10,0.2,4,3,true));
    fillvalue=colorant"white", nrow=1, npad=10
)

# ## Reference

#md # ```@docs
#md # ElasticDistortion
#md # ```

## save covers #src
using ImageMagick #src
include(joinpath("..", "assets", "utilities.jl")) #src
cover = make_gif(testpattern(RGB, ratio=0.5), ElasticDistortion(15,15,0.1), 10) #src
ImageMagick.save("elasticdistortion.gif", cover; fps=2) #src


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

"""
    renumbering(neworder::Array{Int})

Renumber the RCM ordered nodes. A part of the Reverse Cuthill-McKee Algorithm.

# Examples

```jldoctest
julia> neworder = Int[6, 7, 2, 5, 3, 8, 1, 4, 9];

julia> renumbering(neworder)
Dict{Int64,Int64} with 9 entries:
  7 => 2
  9 => 9
  4 => 8
  2 => 3
  3 => 5
  5 => 4
  8 => 6
  6 => 1
  1 => 7

```

# References

* Wikipedia contributors. "Cuthill–McKee algorithm". Wikipedia, The Free Encyclopedia. Wikipedia, The Free Encyclopedia, 11 Jul. 2017. Web. 20 Jul. 2017. https://en.wikipedia.org/wiki/Cuthill%E2%80%93McKee_algorithm
"""
function renumbering(neworder::Array{Int})
    finalnodes = Dict{Int, Int}()
    n=length(neworder)
    for i=1:n
      neworder[i]
      finalnodes[i] = neworder[i]
    end
Final = Dict((v => k for (k,v) in finalnodes))
return Final
end


text:
using LinearAlgebra, Statistics, StatsBase, Plots; pyplot()

n, N = 5, 10^6
P = diagm(-1 => fill(1/3,n-1),
           0 => fill(1/3,n),
           1 => fill(1/3,n-1))
P[1,n], P[n,1] = 1/3, 1/3

A = UpperTriangular(ones(n,n))
C = P*A

function f1(x,u)
    for xNew in 1:n
        if u <= C[x+1,xNew]
            return xNew-1
        end
    end
end

f2(x,xi) = mod(x + xi , n)

function countTau(f,rnd)
    t = 0
    visits = fill(false,n)
    state = 0
    while sum(visits) < n
        state = f(state,rnd())
        visits[state+1] |= true
        t += 1
    end
    return t-1
end

data1 = [countTau(f1,rand) for _ in 1:N]
data2 = [countTau(f2,()->rand([-1,0,1]) ) for _ in 1:N]
est1, est2 = mean(data1), mean(data2)
c1, c2 = counts(data1)/N,counts(data2)/N
println("Estimated mean value of tau using f1: ",est1)
println("Estimated mean value of tau using f2: ",est2)
println("\nThe matrix P:", P)
scatter(4:33,c1[1:30],
	c=:blue, ms=5, msw=0,
	label="Transition probability matrix")
scatter!(4:33,c2[1:30],
	c=:red, ms=5, msw=0, shape=:cross,
	label="Stochastic recursive formula", xlabel="Time", ylabel="Probability")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using Test
using TwentyFour

@test solve(5,5,3,1) == "No solution"
x = solve(5,5,2,1)
@test eval(Meta.parse(x)) == 24


text:
# sexp to Array{Union{T, Missing}}
for S in (:IntSxp, :RealSxp, :CplxSxp, :LglSxp, :StrSxp)
    @eval begin
        function rcopy(::Type{Array{Union{T, Missing}}}, s::Ptr{$S}) where T
            protect(s)
            v = try
                Union{T, Missing}[isNA(e) ? missing : rcopy(T,e) for e in s]
            finally
                unprotect(1)
            end
            reshape(v,size(s))
        end
        function rcopy(::Type{Vector{Union{T, Missing}}}, s::Ptr{$S}) where T
            protect(s)
            v = try
                Union{T, Missing}[isNA(e) ? missing : rcopy(T,e) for e in s]
            finally
                unprotect(1)
            end
            v
        end
    end
end

## Array{Union{T, Missing}} to sexp conversion.

for (J, S, C) in ((:Integer, :IntSxp, :integer),
                 (:AbstractFloat, :RealSxp, :numeric),
                 (:Complex, :CplxSxp, :complex),
                 (:Bool, :LglSxp, :logical),
                 (:AbstractString, :StrSxp, :character))
    @eval begin
        function sexp(::Type{RClass{$(QuoteNode(C))}}, v::AbstractArray{Union{T, Missing}}) where T <: $J
            rv = protect(allocArray($S, size(v)...))
            try
                for (i, x) in enumerate(v)
                    if ismissing(x)
                        rv[i] = naeltype($S)
                    else
                        rv[i] = x
                    end
                end
            finally
                unprotect(1)
            end
            rv
        end
    end
end


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

name = "NetCDF"
version = v"4.7.4"

# Collection of sources required to build NetCDF
sources = [
    ArchiveSource("https://github.com/Unidata/netcdf-c/archive/v$(version).zip",
                  "170c9c9020f8909811b06e1034d5ea9288b3d5bd90793e3dd27490191faa7566")
]

# HDF5.h in /workspace/artifacts/805ccba77cd286c1afc127d1e45aae324b507973/include
# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/netcdf-c-*

export CPPFLAGS="-I$prefix/include"
export LDFLAGS="-L${libdir}"

if [[ ${target} == *-mingw* ]]; then
   export LDFLAGS="-L${libdir} -lhdf5-0 -lhdf5_hl-0 -lcurl-4"

   if [ ! -f ${WORKSPACE}/destdir/bin/libzlib1.dll ]; then
       ln -s ${WORKSPACE}/destdir/bin/libz.dll ${WORKSPACE}/destdir/bin/libzlib1.dll
   fi

   ./configure --prefix=$prefix --build=${MACHTYPE} --host=${target}  --disable-utilities --enable-shared --disable-static

   # linking fails with: "libtool:   error: can't build x86_64-w64-mingw32 shared library unless -no-undefined is specified"
   # unless -no-undefined is added to LDFLAGS
   make LDFLAGS="-no-undefined -L$prefix/bin -lhdf5-0 -lhdf5_hl-0 -lcurl-4" -j${nproc}

elif [[ "${target}" == *-apple-* ]]; then
    # this file is referenced by hdf.h by not installed
    touch /workspace/destdir/include/features.h

    if [ ! -f ${libdir}/libhdf5.dylib ]; then
        ln -s ${libdir}/libhdf5.*.dylib ${libdir}/libhdf5.dylib
    fi

    if [ ! -f ${libdir}/libhdf5_hl.dylib ]; then
        ln -s ${libdir}/libhdf5_hl.*.dylib ${libdir}/libhdf5_hl.dylib
    fi

   ./configure --prefix=$prefix --build=${MACHTYPE} --host=${target}  --disable-utilities --enable-shared --disable-static
   make -j${nproc}
else
    # do not exist on Platform("x86_64", "linux")
    if [ ! -f ${libdir}/libhdf5.so ]; then
        ln -s ${libdir}/libhdf5.so.* ${libdir}/libhdf5.so
    fi

    if [ ! -f ${libdir}/libhdf5_hl.so ]; then
        ln -s ${libdir}/libhdf5_hl.so.* ${libdir}/libhdf5_hl.so
    fi

   ./configure --prefix=$prefix --build=${MACHTYPE} --host=${target}  --disable-utilities --enable-shared --disable-static
   make -j${nproc}
fi

make install
nc-config --all
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
# Set equal to the supported platforms in HDF5
platforms = [
    Platform("x86_64", "linux"),
    # HDF5_jll on armv7l should use the same glibc as the root filesystem
    # before it can be used
    # https://github.com/JuliaPackaging/Yggdrasil/pull/1090#discussion_r432683488
    # Platform("armv7l", "linux"; libc="glibc"),
    Platform("aarch64", "linux"; libc="glibc"),
    Platform("x86_64", "macos"),
    Platform("x86_64", "windows"),
    Platform("i686", "windows"),
]

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

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="HDF5_jll", version="1.12.0")),
    Dependency("Zlib_jll"),
    Dependency("LibCURL_jll", v"7.71.1"),
    # The following libraries are dependencies of LibCURL_jll which is now a
    # stdlib, but the stdlib doesn't explicitly list its dependencies
    Dependency("LibSSH2_jll"),
    Dependency("MbedTLS_jll"),
    Dependency("nghttp2_jll"),
]

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


text:
# uniform random walk
function proposal1(current, tuning, lb, ub)
    trial = copy(current)
    if rand() > 0.1
        i = rand(1:size(trial,1))
        trial[i] += tuning[i].*randn()
    else
        trial = lb + (ub - lb).*rand(size(lb,1))
    end
    return trial
end

max_stars_repo_name: mcreel/SV

max_issues_repo_name: mcreel/SV

max_forks_repo_name: mcreel/SV


text:
@testset "Hermite 5D" begin

    f(x, y, z, v, w) = cos(x^2 + y^2 + z^2 + v^2 + w^2)

    xmin, xmax, nx = -1, 1, 32
    ymin, ymax, ny = -1, 1, 32
    zmin, zmax, nz = -1, 1, 32
    vmin, vmax, nv = -1, 1, 32
    wmin, wmax, nw = -1, 1, 32

    ϵ = 0.1
    γ = 3.0

    hermite_x = Hermite(Chebyshev( xmin, xmax, nx ), ϵ, γ)
    hermite_y = Hermite(Chebyshev( ymin, ymax, ny ), ϵ, γ)
    hermite_z = Hermite(Chebyshev( zmin, zmax, nz ), ϵ, γ)
    hermite_v = Hermite(Chebyshev( vmin, vmax, nv ), ϵ, γ)
    hermite_w = Hermite(Chebyshev( wmin, wmax, nw ), ϵ, γ)

    xk = hermite_x.nodes
    yk = hermite_y.nodes
    zk = hermite_z.nodes
    vk = hermite_v.nodes
    wk = hermite_w.nodes

    nxe, nye, nze, nve, nwe = 32, 32, 32, 32, 32

    xe = collect(range(xmin, stop=xmax, length=nxe))
    ye = collect(range(ymin, stop=ymax, length=nye))
    ze = collect(range(zmin, stop=zmax, length=nze))
    ve = collect(range(vmin, stop=vmax, length=nve))
    we = collect(range(wmin, stop=wmax, length=nwe))
    dx = abs(xe[2] - xe[1])
    dy = abs(ye[2] - ye[1])
    dz = abs(ze[2] - ze[1])
    dv = abs(ve[2] - ve[1])
    dw = abs(we[2] - we[1])

    fk = [f(x, y, z, v, w) for x in xk, y in yk, z in zk, v in vk, w in wk]
    fe = [f(x, y, z, v, w) for x in xe, y in ye, z in ze, v in ve, w in we]

    @time s  = interpolate(hermite_x, hermite_y, hermite_z, hermite_v, hermite_w, fk, xe, ye, ze, ve, we)

    max_error = maximum(abs.(s .- fe))

    l2_error  = sqrt(trapz((s - fe).^2, dx, dy, dz, dv, dw))

    println(" Hermite 5d max error : $max_error ")
    println(" Hermite 5d L2  error : $l2_error  ")

    @test true

end


text:
"""
The enhanced GJK algorithm relies on a pre-computed set of neighbors for each
vertex in the mesh. In order to use those neighbors, we have to know from which
vertex to start. Specifically, we need to know the index of the vertex
corresponding to the point in the GJK simplex which we are trying to improve.
To do that, we introduce the notion of a `Tagged` point. A tagged point is
just a point and some arbitrary additional data field. All of the
`any_inside` and `support_vector_max` functions in this package return tagged
points. For most geometries, that tag is empty (`nothing`). But for our
NeighborMesh type, the tag is the linear index into the vertices of the mesh,
which lets us look up that mesh's neighbors faster later on.
"""
struct Tagged{P, T}
    point::P
    tag::T
end

Tagged(point::P) where {P} = Tagged(point, nothing)

Base.:*(n::Number, t::Tagged) = n * value(t)

value(t::Tagged) = t.point

function any_inside(geometry)
    point = gt.any_inside(geometry)
    Tagged(SVector(point))
end

function any_inside(mesh::gt.AbstractMesh{gt.Point{N, T}}) where {N,T}
    Tagged(SVector(first(gt.vertices(mesh))))
end

function any_inside(point::gt.Point)
    Tagged(SVector(point))
end


text:
"""
f(x,y) = x+y
"""
function model_1(n_train = 100, n_test = 30)
    n_train = 100
    xtrain=rand(2,n_train)
    ytrain=sum(xtrain, dims=1)

    xtest=10rand(2,n_test)
    ytest=sum(xtest, dims=1)

    return (xtrain=xtrain, ytrain=ytrain, xtest=xtest, ytest=ytest)
end


text:
module Runner

export execute

function execute(args::Dict{Symbol,Any})
  println("Running...")
end

end # module


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

import Statistics

const _abs_err_orientation_2d = 2*eps(Float64)

"""
    ConvexPolygon{N, T<:Real} <: PlanarShape{T}

General Convex Polygon surface, not a valid CSG object.
The rotation of the polygon around its normal is defined by `rotationvec`.
`rotationvec×surfacenormal` is taken as the vector along the u axis.

```julia
ConvexPolygon(local_frame::Transform{T}, local_polygon_points::Vector{SVector{2, T}}, interface::NullOrFresnel{T} = nullinterface(T))
```

The local frame defines the plane (spans by the right and up vectors) with the plane normal given by the forward vector.
the local_polygon_points are given with respect to the local frame and are 2D points.
NOTE: This class uses static vectors to hold the points which will lead to more efficient performance, but should not be used with polygons with more than 20-30 points.
"""
struct ConvexPolygon{N,T<:Real}  <: PlanarShape{T}
    plane::Plane{T,3}
    local_frame::Transform{T}
    # local_points::Vector{SVector{2, T}}
    local_points::SMatrix{2,N,T}
    # for efficency
    _local_frame_inv::Transform{T}                                  # cache the inverse matrix to avoid computing it for every intersection test
    _local_lines::Vector{SVector{3, SVector{2, T}}}                 # defines the edge points + a third point representing the slopes in order to save some calculationsduring ray checking
    _length::Int64                                                  # cache the length of the polygon

    function ConvexPolygon(
            local_frame::Transform{T},
            local_polygon_points::Vector{SVector{2, T}},
            interface::NullOrFresnel{T} = NullInterface(T)
        ) where {T<:Real}

        # need at least 3 points to define apolygon
        @assert length(local_polygon_points) > 3

        local_center = Statistics.mean(local_polygon_points)
        world_center = local2world(local_frame) * Vec3(local_center[1], local_center[2], zero(T))

        # poly_points = [GeometricalPredicates.Point2D(p[1], p[2]) for p in local_polygon_points]
        # poly = GeometricalPredicates.Polygon2D(poly_points...)
        pts = local_polygon_points
        local_lines = SVector(
            [SVector(
                pts[i],                                             # A
                pts[mod(i, length(pts))+1],                         # B
                                                                    # bx = Bx - Ax,  by = By - Ay
                SVector(pts[mod(i, length(pts))+1][1] - pts[i][1], pts[mod(i, length(pts))+1][2] - pts[i][2]))
                for i in 1:length(pts)]...
        )

        plane = Plane(forward(local_frame), world_center, interface = interface)
        N = length(local_polygon_points)
        temp = MMatrix{2,N,T}(undef)
        for (i,pt) in pairs(local_polygon_points)
            temp[:,i] = pt
        end

        N2 = 2*N
        new{N,T}(plane, local_frame, SMatrix{2,N,T,N2}(temp), inv(local_frame), local_lines, length(local_lines))
    end
end
export ConvexPolygon

centroid(poly::ConvexPolygon) = poly.plane.pointonplane

#function barrier to make vertices allocate less and be faster.
function to3d(pts::SMatrix{2,N,T,L}) where{N,L,T}
    temp = MMatrix{3,N,T}(undef)
    for row in 1:2
        for col in 1:N
            temp[row,col] = pts[row,col]
        end
    end

    for col in 1:N
        temp[3,col] = T(0)
    end

    return SMatrix{3,N,T}(temp)
return temp
end

#this function allocates. Don't know why, it shouldn't but it does.
function vertices(poly::ConvexPolygon{N,T}) where{N,T<:Real}
   return poly.local_frame * to3d(poly.local_points)
end

function surfaceintersection(poly::ConvexPolygon{N,T}, r::AbstractRay{T,3}) where {N,T<:Real}
    interval = surfaceintersection(poly.plane, r)
    if interval isa EmptyInterval{T} || isinfiniteinterval(interval)
        return EmptyInterval(T) # no ray plane intersection or inside plane but no hit
    else
        intersect = halfspaceintersection(interval)
        p = point(intersect)

        local_p = poly._local_frame_inv * p

        @inline function orientation(l::SVector{3, SVector{2, T}})::Int8 where {T<:Real}
            cx = local_p[1] - l[1][1]                         # getx(p) - getx(geta(l))
            cy = local_p[2] - l[1][2]                         # gety(p) - gety(geta(l))
            _pr2 = -(l[3][1]*cy) + l[3][2]*cx           # _pr2 = -l._bx*cy + l._by*cx
            if _pr2 < -_abs_err_orientation_2d
                1
            elseif _pr2 > _abs_err_orientation_2d
                -1
            else
                0   #  can implement something more accurate in the future to minimize numerical errors
            end
        end

        @inline function inpolygon()
            side = orientation(poly._local_lines[1])
            for i = 2:poly._length
                orientation(poly._local_lines[i]) == side || return false
            end
            return true
        end

        if !inpolygon()
            return EmptyInterval(T)
        else
            if dot(normal(poly), direction(r)) < zero(T)
                return positivehalfspace(intersect)
            else
                return rayorigininterval(intersect)
            end
        end
    end
end

"""
    makemesh(poly::ConvexPolygon{N, T}, ::Int = 0) where {N, T<:Real} -> TriangleMesh

Create a triangle mesh that can be rendered by iterating on the polygon's edges and for each edge use the centroid as the third vertex of the triangle.
"""
function makemesh(poly::ConvexPolygon{N,T}, ::Int = 0) where {N,T<:Real}
    c = centroid(poly)

    l2w = local2world(poly.local_frame)
    len = Size(poly.local_points)[2]

    triangles = []
    for i in 1:len
        p1 = poly.local_points[:,i]
        p2 = poly.local_points[:,mod(i,len) + 1]

        tri = Triangle(
            Vector(l2w * Vec3(p2[1], p2[2], zero(T))),
            Vector(c),
            Vector(l2w * Vec3(p1[1], p1[2], zero(T))))
        push!(triangles, tri)
    end

    triangles = Vector{Triangle{T}}(triangles)

    return TriangleMesh(triangles)
end


text:
using Plots, JLD2

cd(@__DIR__)
begin
    u0 = []
    @load "exp.jld2" x sol
    push!(u0, sol)
    @load "l2.jld2" x sol
    push!(u0, sol)
    @load "lasso.jld2" x sol
    push!(u0, sol)
    @load "nofilter.jld2" x sol
    push!(u0, sol)

    u = []
    @load "l2_apt.jld2" x sol
    push!(u, sol)
end

plot(x, u0[end][:, 2], label="No filter", lw=1.5, color=:gray27, xlabel="x", ylabel="standard deviation")
plot!(x, u0[1][:, 2], label="Exp", lw=1.5, color=1)
plot!(x, u0[3][:, 2], label="Lasso", lw=1.5, color=2)
plot!(x, u0[2][:, 2], label="L²", lw=1.5, color=3)
plot!(x, u[1][:, 2], label="Adaptive L²", lw=1.5, color=4)
savefig("burgers_std.pdf")

plot(x, u0[end][:, 1], label="No filter", lw=1.5, color=:gray27, xlabel="x", ylabel="expected value")
plot!(x, u0[1][:, 1], label="Exp", lw=1.5, color=1)
plot!(x, u0[3][:, 1], label="Lasso", lw=1.5, color=2)
plot!(x, u0[2][:, 1], label="L²", lw=1.5, color=3)
plot!(x, u[1][:, 1], label="Adaptive L²", lw=1.5, color=4)
savefig("burgers_mean.pdf")

max_stars_repo_path: src/burgers/6th order/plot.jl

max_stars_repo_name: CSMMLab/FRSG

max_issues_repo_path: src/burgers/6th order/plot.jl

max_issues_repo_name: CSMMLab/FRSG

max_forks_repo_path: src/burgers/6th order/plot.jl

max_forks_repo_name: CSMMLab/FRSG


text:
using DataStructures

"""
    stack_to_rational(st::Stack)

Take a stack consisting of the convergents of the continued
fraction expansion of a number and return the rational approximation.

Note: This function makes a deep copy of the given stack.

Ex.
    [3, 7, 15] ---> 333/106

julia> st = Stack(Int)
DataStructures.Stack{Int64}(Deque [Int64[]])

julia> push!(st, 3)
DataStructures.Stack{Int64}(Deque [[3]])

julia> push!(st, 7)
DataStructures.Stack{Int64}(Deque [[3,7]])

julia> push!(st, 15)
DataStructures.Stack{Int64}(Deque [[3,7,15]])

julia> stack_to_rational(st)
333//106
"""
function stack_to_rational(st)
    tmp = deepcopy(st)
    x = pop!(tmp)

    while (length(tmp) > 0)
        x = pop!(tmp) + 1 // x
    end

    return x
end

"""
    diophantine(x [, ε = 10.0-5])

Return a rational approximation of x with an error tolerance of ε (default value of 10^-5).

Using notation from Bremner (pg 146).

julia> diophantine(float(π))
355//113
"""
function diophantine(x, ε = 10.0^-5)
    st = Stack(Int)

    # step i=2
    γ = x
    qi = floor(Int, γ)
    push!(st, qi)

    # step i=3,4,...
    while (abs(stack_to_rational(st) - x) > ε)
        γ = 1 / (γ - qi)
        qi = floor(Int, γ)
        push!(st, qi)
    end

    return stack_to_rational(st)
end


text:
import Nemo: nf_elem, AnticNumberField, FmpzPolyRing, degree,
             lg, prime_decomposition, parent, length,
             norm, real, imag, inv, rows, getindex!, lll, hnf, cols, basis,
             trace, factor, mod, zero,
             hash, PolynomialRing, coeff,
             var, abs, min, iszero, one, sqrt, isone, deepcopy, rank, in,
             discriminant, log, sub, lift, FlintQQ, FlintZZ, elem_type,
             elem_from_mat_row, elem_to_mat_row!, order, signature,
             base_ring, compose, root, arf_struct, fmpq, valuation, remove,
             Ring, prec, conj, mul!, gen, divexact, derivative, zero!, divrem,
             resultant, evaluate, setcoeff!, div, isodd, iseven, max, floor,
             ceil, //, setindex!, transpose, colon, nf_elem, isreal,
             MatrixSpace, contains, overlaps, solve, unique_integer, gcd,
             minpoly, charpoly, det, howell_form, needs_parentheses,
             isnegative, characteristic, order intersect, lcm, fmpq


text:
"""
    AbstractIndexedGraph{T} <: AbstractGraph{T}
An abstract type representing an indexed graph.
`AbstractIndexedGraph`s must have the following elements:
  - `A::SparseMatrixCSC` adjacency matrix
"""
abstract type AbstractIndexedGraph{T} <: AbstractGraph{T} end

function Base.show(io::IO, g::AbstractIndexedGraph)
    s = is_directed(g) ? "directed" : "undirected"
    println(io, "{$(nv(g)), $(ne(g))} ", s, " sparse ", eltype(g), " graph")
    Base.print_array(io, adjacency_matrix(g))
end

Base.eltype(g::AbstractIndexedGraph{T}) where T = T

edgetype(g::AbstractIndexedGraph{T}) where T = IndexedEdge{T}

has_vertex(g::AbstractIndexedGraph, i::Integer) = i ≤ size(g.A, 2)

has_edge(g::AbstractIndexedGraph, i::Integer, j::Integer) = g.A[i,j] != 0

nv(g::AbstractIndexedGraph) = size(g.A, 2)

vertices(g::AbstractIndexedGraph) = 1:size(g.A, 2)

outneighbors(g::AbstractIndexedGraph, i::Integer) = @inbounds @view g.A.rowval[nzrange(g.A,i)]

# Returns sparse adj matrix. Elements default to Int (to match Graphs)
function adjacency_matrix(g::AbstractIndexedGraph, T::DataType=Int)
    SparseMatrixCSC(g.A.m, g.A.n, g.A.colptr, g.A.rowval, ones(T, nnz(g.A)))
end


text:
using ModifiedArrays
using Test

include("offset.jl")
include("adjtrans.jl")
include("eltype.jl")


text:
"""
Abstract types and methods for codes, error models and decoders.
"""
module Model

# imports
using ..Qecsim:QecsimError
using ..PauliTools:bsp
using LinearAlgebra:I

# exports
# - abstract model
export AbstractModel
export label
# - stabilizer code
export StabilizerCode
export stabilizers, logical_xs, logical_zs, logicals, nkd, validate
# - error model
export ErrorModel
export generate, probability_distribution
# - decoder
export Decoder
export DecodeResult
export decode

# AbstractModel

"""
Abstract supertype for models.
"""
abstract type AbstractModel end

"""
    label(model) -> String

Return a label suitable for use in plots and for grouping results.

!!! note "Abstract method"

    This method should be implemented for concrete subtypes or duck-typed implementations of
    [`AbstractModel`](@ref).
"""
function label end

# StabilizerCode

"""
    StabilizerCode <: AbstractModel

Abstract supertype for stabilizer codes.
"""
abstract type StabilizerCode <: AbstractModel end

"""
    stabilizers(code) -> BitMatrix

Return the stabilizers in binary symplectic form.

Each row is a stabilizer generator. An overcomplete set of generators can be included to
simplify decoding.

!!! note "Abstract method"

    This method should be implemented for concrete subtypes or duck-typed implementations of
    [`StabilizerCode`](@ref).
"""
function stabilizers end

"""
    logical_xs(code) -> BitMatrix

Return the logical X operators in binary symplectic form.

Each row is an operator. The order should match that of [`logical_zs`](@ref).

!!! note "Abstract method"

    This method should be implemented for concrete subtypes or duck-typed implementations of
    [`StabilizerCode`](@ref).
"""
function logical_xs end

"""
    logical_zs(code) -> BitMatrix

Return the logical Z operators in binary symplectic form.

Each row is an operator. The order should match that of [`logical_xs`](@ref).

!!! note "Abstract method"

    This method should be implemented for concrete subtypes or duck-typed implementations of
    [`StabilizerCode`](@ref).
"""
function logical_zs end

"""
    logicals(code) -> BitMatrix

Return the logical operators in binary symplectic form.

Each row is an operator. X operators are stacked above Z operators in the order given by
[`logical_xs`](@ref) and [`logical_zs`](@ref).
"""
function logicals(code)
    return vcat(logical_xs(code), logical_zs(code))
end

"""
    nkd(code) -> Tuple{Int,Int,Union{Int,Missing}}

Return a descriptor in the format `(n, k, d)`, where `n` is the number of physical qubits,
`k` is the number of logical qubits, and `d` is the distance of the code (or `missing` if
unknown).

!!! note "Abstract method"

    This method should be implemented for concrete subtypes or duck-typed implementations of
    [`StabilizerCode`](@ref).
"""
function nkd end

@doc raw"""
    validate(code)

Perform sanity checks.

If any of the following fail then a [`QecsimError`](@ref) is thrown:

* ``S ⊙ S^T = 0``
* ``S ⊙ L^T = 0``
* ``L ⊙ L^T = Λ``

where ``S`` and ``L`` are the code [`stabilizers`](@ref) and [`logicals`](@ref),
respectively, and ``⊙`` and ``Λ`` are defined in [`PauliTools.bsp`](@ref bsp).
"""
function validate(code)
    s, l = stabilizers(code), logicals(code)
    !any(bsp(s, transpose(s))) || throw(QecsimError("stabilizers do not mutually commute"))
    !any(bsp(s, transpose(l))) || throw(QecsimError(
        "stabilizers do not commute with logicals"))
    # twisted identity with same size as logicals
    nlogicals = size(l, 1)
    twistedI = circshift(Matrix(I, nlogicals, nlogicals), nlogicals / 2)
    (bsp(l, transpose(l)) == twistedI) || throw(QecsimError(
        "logicals do not mutually twist commute"))
    return nothing
end

# ErrorModel

"""
    ErrorModel <: AbstractModel

Abstract supertype for error models.
"""
abstract type ErrorModel <: AbstractModel end

"""
    generate(error_model, code, p::Real, [rng::AbstractRNG=GLOBAL_RNG]) -> BitVector

Generate a new error in binary symplectic form according to the `error_model` and `code`,
where `p` is typically the probability of an error on a single qubit.

!!! note "Abstract method"

    This method should be implemented for concrete subtypes or duck-typed implementations of
    [`ErrorModel`](@ref).
"""
function generate end

"""
    probability_distribution(error_model, p::Real) -> NTuple{4,Real}

Return the single-qubit probability distribution amongst Pauli I, X, Y and Z, where `p` is
the overall probability of an error on a single qubit.

!!! note "Abstract method [optional]"

    This method is **not** invoked by any core modules. Since it is often useful for
    decoders, it is provided as a template and concrete subtypes or duck-typed
    implementations of [`ErrorModel`](@ref) are encouraged to implement it when appropriate,
    particularly for IID error models.
"""
function probability_distribution end

# Decoder
"""
    Decoder <: AbstractModel

Abstract supertype for decoders.
"""
abstract type Decoder <: AbstractModel end

"""
    decode(decoder, code, syndrome::AbstractVector{Bool}; kwargs...) -> DecodeResult

Resolve a recovery operation for the given `code` and `syndrome`, or evaluate the success of
decoding, as encapsulated in the decode result.

The syndrome has length equal to the number of code stabilizers, and element values of 0
or 1 (false or true) indicate whether the corresponding stabilizer does or does not commute
with the error, respectively.

Keyword parameters `kwargs` may be provided by the client, e.g. [`App`](@ref App), with
context values such as `error_model`, `error_probability` and `error`. Most implementations
will ignore such parameters; however, if they are used, implementations should declare them
explicitly and treat them as optional.

See also [`DecodeResult`](@ref).

!!! note "Abstract method"

    This method should be implemented for concrete subtypes or duck-typed implementations of
    [`Decoder`](@ref).
"""
function decode end

"""
    DecodeResult(
        success::Union{Nothing,Bool},
        recovery::Union{Nothing,AbstractVector{Bool}},
        logical_commutations::Union{Nothing,AbstractVector{Bool}}
        custom_values::Union{Nothing,AbstractVector}
    )
    DecodeResult(;
        success::Union{Nothing,Bool}=nothing,
        recovery::Union{Nothing,AbstractVector{Bool}}=nothing,
        logical_commutations::Union{Nothing,AbstractVector{Bool}}=nothing
        custom_values::Union{Nothing,AbstractVector}=nothing
    )

Construct a decoding result as returned by [`decode`](@ref).

Typically decoders will provide a `recovery` operation and delegate the evaluation of
`success` and `logical_commutations` to the client, e.g. [`App`](@ref App). Optionally,
`success` and/or `logical_commutations` may be provided as overrides. At least one of
`recovery` or `success` must be specified to allow a success value to be resolved.
Additionally `custom_values` may be specified. If `logical_commutations` and/or
`custom_values` are provided then they should be of consistent type and size over
identically parameterized simulation runs. For example, [`App`](@ref App) will sum
`logical_commutations` across runs and, similarly, if `custom_values` are numbers they will
be summed across runs, and if they are arrays they will be concatenated using `vcat`.

See also [`decode`](@ref).

# Examples
```jldoctest
julia> DecodeResult(; recovery=BitVector([1, 0, 1, 0, 1, 1]))  # typical use-case
DecodeResult{Nothing}(nothing, Bool[1, 0, 1, 0, 1, 1], nothing, nothing)

julia> DecodeResult(; success=true)  # override success
DecodeResult{Nothing}(true, nothing, nothing, nothing)

julia> DecodeResult(; success=false, logical_commutations=BitVector([1, 0]))  # override all
DecodeResult{Nothing}(false, nothing, Bool[1, 0], nothing)

julia> DecodeResult(; success=true, custom_values=[2.3, 4.1])  # custom values (numbers)
DecodeResult{Vector{Float64}}(true, nothing, nothing, [2.3, 4.1])

julia> DecodeResult(; success=true, custom_values=[[2.3], [4.1]])  # custom values (arrays)
DecodeResult{Vector{Vector{Float64}}}(true, nothing, nothing, [[2.3], [4.1]])

julia> DecodeResult(true, nothing, nothing, [2.3, 4.1])  # positional parameters
DecodeResult{Vector{Float64}}(true, nothing, nothing, [2.3, 4.1])

julia> DecodeResult(; success=nothing, recovery=nothing)  # too few specified parameters
ERROR: QecsimError: at least one of 'success' or 'recovery' must be specified
Stacktrace:
[...]
```
"""
struct DecodeResult{T<:Union{Nothing,Vector}}
    # Parametric type to avoid abstract types in struct, see:
    # https://docs.julialang.org/en/v1/manual/performance-tips/#Type-declarations
    # https://discourse.julialang.org/t/union-of-nothing-and-parametric-types/11986/5
    success::Union{Nothing,Bool}
    recovery::Union{Nothing,BitVector}
    logical_commutations::Union{Nothing,BitVector}
    custom_values::T
    # lenient constructor that infers types
    function DecodeResult(success, recovery, logical_commutations,
                          custom_values::Union{Nothing,AbstractVector{S}}) where S
        if isnothing(success) && isnothing(recovery)
            throw(QecsimError("at least one of 'success' or 'recovery' must be specified"))
        end
        if isnothing(custom_values)
            new{Nothing}(success, recovery, logical_commutations, custom_values)
        else
            new{Vector{S}}(success, recovery, logical_commutations, collect(custom_values))
        end
    end
end
# keyword convenience constructor
function DecodeResult(; success=nothing, recovery=nothing, logical_commutations=nothing,
                      custom_values=nothing)
    DecodeResult(success, recovery, logical_commutations, custom_values)
end

end


text:
using Clang.Generators
using assimp_jll

cd(@__DIR__)

const ASSIMP_INCLUDE = joinpath(assimp_jll.artifact_dir, "include") |> normpath
const C_INTERFACE_HEADERS = ["cfileio.h",  "cimport.h", "cexport.h", "scene.h", "postprocess.h", "pbrmaterial.h", "version.h"]
const ASSIMP_HEADERS = [joinpath(ASSIMP_INCLUDE, "assimp", h) for h in C_INTERFACE_HEADERS]

options = load_options(joinpath(@__DIR__, "generator.toml"))

args = get_default_args()
push!(args, "-I$ASSIMP_INCLUDE")

ctx = create_context(ASSIMP_HEADERS, args, options)

build!(ctx)


text:
import ASTInterpreter: Interpreter, enter_call_expr, determine_line_and_file, next_line!,
  evaluated!, finish!, step_expr, idx_stack

import ..Atom: fullpath, handle, @msg, wsitem, Inline
using Media

function fileline(i::Interpreter)
  file, line = determine_line_and_file(i, idx_stack(i))[end]
  Atom.fullpath(string(file)), line
end

debugmode(on) = @msg debugmode(on)
stepto(file, line, text) = @msg stepto(file, line, text)
stepto(i::Interpreter) = stepto(fileline(i)..., stepview(i.next_expr[2]))
stepto(::Void) = debugmode(false)

function stepview(ex)
  @capture(ex, f_(as__)) || return render(Inline(), Text(string(ex)))
  render(Inline(), span(c(render(Inline(), f),
                          "(",
                          interpose([render(Inline(), a) for a in as], ", ")...,
                          ")")))
end

interp = nothing

const cond = Condition()

isdebugging() = interp ≠ nothing

validcall(x) =
  @capture(x, f_(args__)) &&
  !isa(f, Core.IntrinsicFunction) &&
  f ∉ [tuple, getfield]

function tocall!(interp)
  while !validcall(interp.next_expr[2])
    step_expr(interp) || return false
  end
  return true
end

function err(e)
  debugmode(false)
  interp = nothing
  notify(cond, e, error = true)
end

macro errs(ex)
  :(try
      $(esc(ex))
    catch e
      err(e)
    end)
end

function done(interp)
  stack, val = interp.stack, interp.retval
  stack = filter(x -> isa(x, Interpreter), stack)
  if stack[1] == interp
    debugmode(false)
    interp = nothing
    notify(cond)
  else
    i = findfirst(stack, interp)
    resize!(stack, i-1)
    interp = stack[end]
    evaluated!(interp, val)
    tocall!(interp)
  end
  return interp
end

handle("nextline") do
  @errs begin
    global interp = next_line!(interp) && tocall!(interp) ? interp : done(interp)
    stepto(interp)
  end
end

handle("stepin") do
  global interp
  isexpr(interp.next_expr[2], :call) || return
  new = enter_call_expr(interp, interp.next_expr[2])
  if new ≠ nothing
    interp = new
    tocall!(interp)
    stepto(interp)
  end
end

handle("finish") do
  global interp
  @errs begin
    finish!(interp)
    interp = done(interp)
    stepto(interp)
  end
end

handle("stepexpr") do
  @errs begin
    step_expr(interp)
    stepto(interp)
  end
end

contexts(i::Interpreter = interp) =
  reverse!([d(:context => i.linfo.def.name, :items => context(i)) for i in i.stack])

import Gallium: JuliaStackFrame

type Undefined end

@render Inline u::Undefined span(".fade", "<undefined>")

Atom.wsicon(::Undefined) = "icon-circle-slash"

function context(i::Union{Interpreter,JuliaStackFrame})
  items = []
  for (k, v) in zip(i.linfo.sparam_syms, i.env.sparams)
    push!(items, wsitem(k, v))
  end
  isdefined(i.linfo, :slotnames) || return items
  for (k, v) in zip(i.linfo.slotnames, i.env.locals)
    k in (Symbol("#self#"), Symbol("#unused#")) && continue
    push!(items, wsitem(k, isnull(v) ? Undefined() : get(v)))
  end
  return items
end

context(i) = []

function interpret(code::AbstractString, i::Interpreter = interp)
  code = parse(code)
  ok, result = ASTInterpreter.eval_in_interp(i, code)
  return ok ? result : Atom.EvalError(result)
end


text:
using NBodySimulator

function generate_bodies_in_line(n::Int, m::Real, v_dev::Real, L::Real)
    dL = L / (ceil(n^(1 / 3)))
    n_line = floor(Int, L / dL)
    rng = MersenneTwister(n);
    velocities = v_dev * randn(rng, Float64, (3, n_line))
    bodies = MassBody[]
    x = y = L / 2
    for i = [3] #1:n_line-1
        r = SVector(x, y, i * dL)
        v = SVector{3}(velocities[:,i])
        body = MassBody(r, v, m)
        push!(bodies, body)
    end
    return bodies
end

const qe = 1.6e-19
const Na = 6.022e23
const T = 298.16 # °K
const T0 = 275 # °K
const kb = 1.38e-23 # J/K
const ϵOO = 0.1554253*4184/Na
const σOO = 3.165492e-10 # m
const ρ = 216 # kg/m^3
const mO = 15.999 * 1.6747 * 1e-27 # kg
const mH = 1.00794 * 1.6747 * 1e-27#
const mH2O = mO+2*mH
const N = 216#floor(Int, ρ * L^3 / m)
const L = (mH2O*N/ρ)^(1/3)#10.229σ
const R = 9e-10 # 3*σOO
const Rel = 0.49*L
const v_dev = sqrt(kb * T /mH2O)
const τ = 0.5e-15 # σ/v
const t1 = 0τ
const t2 = 6000τ
const k_bond = 1059.162*4184*1e20/Na # J/m^2
const k_angle = 75.90*4184/Na # J/rad^2
const rOH = 1.012e-10 # m
const ∠HOH = 113.24*pi/180 # rad
const qH = 0.41*qe
const qO = -0.82*qe
const k = 9e9 #
#bodies = generate_bodies_in_cell_nodes(N, mH2O, v_dev, L)
bodies = load_water_molecules_from_pdb("C:/Users/Michael/Desktop/GSoC/pdbs/output_4.pdb")
jl_parameters = LennardJonesParameters(ϵOO, σOO, R)
e_parameters = ElectrostaticParameters(k, Rel)
spc_paramters = SPCFwParameters(rOH, ∠HOH, k_bond, k_angle)
pbc = CubicPeriodicBoundaryConditions(L)
water = WaterSPCFw(bodies, mH, mO, qH, qO,  jl_parameters, e_parameters, spc_paramters);
#thermostat = BerendsenThermostat(T0, 200τ)
simulation = NBodySimulation(water, (t1, t2), pbc, kb);
#result = run_simulation(simulation, Tsit5())
result = @time run_simulation(simulation, VelocityVerlet(), dt=τ)

#using JLD
#save("D:/water $Nactual molecules $timesteps steps $time_now.jld", "rs", rs, "grf", grf, "ts", ts, "dr2", dr2)
#save("D:/water $Nactual molecules $timesteps steps $time_now.jld", "rs", rs, "grf", grf, "ts", ts, "dr2", dr2, "e_tot", e_tot, "e_kin", e_kin, "e_pot", e_pot)

#using Plots
#import GR
#plot(rs/1e-10, grf, xlim=[0, 0.4999L/1e-10], label=["Radial distribution function"],ylabel="g(r)", xlabel="r, Å")

#@time animate(result, "D:/$Nactual H20 particles with $timesteps timesteps $time_now.gif")

#=
(rs, grf) = rdf(result)
(ts, dr2) = msd(result)

t = t1:τ:result.solution.t[end-1]
temper = temperature.(result, t)

time_now = Dates.format(now(), "yyyy_mm_dd_HH_MM_SS")
Nactual = length(bodies)
timesteps = round(length(result.solution.t))

using MAT
file = matopen("d:/water real units without thermostat and $timesteps timesteps_$time_now.mat", "w")
write(file, "t", collect(t))
write(file, "temper", temper)
write(file, "rs", rs)
write(file, "ts", ts)
write(file, "grf", grf)
write(file, "dr2", dr2)
close(file)

@time NBodySimulator.save_to_pdb(result, "d:/water in real units from PDB $timesteps timesteps $time_now.pdb" )
=#


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

#=

function _split_dims(dims::NTuple{N,Integer}) where {N}
    int_dims = Int.(dims)
    Base.front(int_dims), int_dims[end]
end

=#

@inline _tail_impl(x, ys...) = (ys...,)
@inline _tail(x) = _tail_impl(x...)

Base.@pure _ncolons(::Val{N}) where N = ntuple(_ -> Colon(), Val{N}())

Base.@propagate_inbounds front_tuple(x::NTuple{N,Any}, ::Val{M}) where {N,M} =
    Base.ntuple(i -> x[i], Val{M}())

Base.@propagate_inbounds back_tuple(x::NTuple{N,Any}, ::Val{M}) where {N,M} =
    Base.ntuple(i -> x[i + N - M], Val{M}())

Base.@propagate_inbounds split_tuple(x::NTuple{N,Any}, ::Val{M}) where {N,M} =
    (front_tuple(x, Val{M}()), back_tuple(x, Val{N - M}()))

Base.@propagate_inbounds swap_front_back_tuple(x::NTuple{N,Any}, ::Val{M}) where {N,M} =
    (back_tuple(x, Val{N - M}())..., front_tuple(x, Val{M}())...)

_convert_elype(::Type{T}, A::AbstractArray{T}) where {T} = A

_convert_elype(::Type{T}, A::AbstractArray{U}) where {T,U} = broadcast(x -> convert(T, x), A)

Base.@pure _add_vals(Val_M::Val{M}, Val_N::Val{N}) where {M,N} =
    Val{length((ntuple(identity, Val_M)..., ntuple(identity, Val_N)...))}()

Base.@pure require_ndims(A::AbstractArray{T,N}, Val_N::Val{N}) where {T,N} =
    nothing

Base.@pure require_ndims(A::AbstractArray{T,M}, Val_N::Val{N}) where {T,M,N} =
    throw(ArgumentError("Require an array with $N dimensions"))