text:
struct TrackedArray{T,N,A<:AbstractArray{T,N}} <: AbstractArray{T,N}
  tracker::Tracked{A}
  data::A
  grad::A
  TrackedArray{T,N,A}(t::Tracked{A}, data::A) where {T,N,A} = new(t, data)
  TrackedArray{T,N,A}(t::Tracked{A}, data::A, grad::A) where {T,N,A} = new(t, data, grad)
end

data(x::TrackedArray) = x.data
tracker(x::TrackedArray) = x.tracker

TrackedVector{T,A} = TrackedArray{T,1,A}
TrackedMatrix{T,A} = TrackedArray{T,2,A}
TrackedVecOrMat{T,A} = Union{TrackedVector{T,A},TrackedMatrix{T,A}}

track(c::Call, x::AbstractArray) = TrackedArray(c, x)

TrackedArray(c::Call, x::A) where A <: AbstractArray =
  TrackedArray{eltype(A),ndims(A),A}(Tracked{A}(c), x)

TrackedArray(c::Call, x::A, Δ::A) where A <: AbstractArray =
  TrackedArray{eltype(A),ndims(A),A}(Tracked{A}(c, Δ), x, Δ)

TrackedArray(x::AbstractArray) = TrackedArray(Call(), x, zeros(x))

Base.eltype(x::Type{<:TrackedArray{T}}) where T <: Real = TrackedReal{T}

Base.show(io::IO, ::Type{TrackedArray{T,N,A}}) where {T,N,A<:AbstractArray{T,N}} =
  print(io, "TrackedArray{…,$A}")

function Base.showarray(io::IO, X::TrackedArray, repr::Bool = true; header = true)
  if repr
    print(io, "param(")
    Base.showarray(io, data(X), true)
    print(io, ")")
  else
    header && print(io, "Tracked ")
    Base.showarray(io, data(X), false, header = header)
  end
end

Base.setindex!(xs::TrackedArray, v, i...) =
  error("Can't differentiate `setindex!`")

back!(::TrackedArray) = error("Value is not scalar; use `back!(sum(x))` or `back!(x, Δ)`")

# Fallthrough methods

for f in :[Base.size, Base.ndims].args
  @eval @inline $f(x::TrackedArray, a...) = $f(data(x), a...)
end

Base.size(x::TrackedArray, i::Integer, j::Integer, is::Integer...) =
  size(data(x), i, j, is...)

Base.similar(x::TrackedArray, dims::Union{AbstractUnitRange,Integer}...) =
  similar(data(x), dims...)

Base.similar(x::TrackedArray, T::Type) = similar(data(x), T)

Base.:(==)(x::TrackedArray, y) = data(x) == y
Base.:(==)(y, x::TrackedArray) = y == data(x)
Base.:(==)(x::TrackedArray, y::TrackedArray) = data(x) == data(y)

# Array Stdlib

Base.getindex(xs::TrackedArray, i...) = track(getindex, xs, i...)

@grad function getindex(xs::AbstractArray, i...)
  data(xs)[i...], function (Δ)
    Δ′ = zero(xs)
    Δ′[i...] = data(Δ)
    (nobacksies(:getindex, Δ′), map(_->nothing, i)...)
  end
end

Base.:-(xs::TrackedArray) = track(-, xs)

@grad -(xs) = -data(xs), Δ -> (-Δ,)

Base.transpose(xs::TrackedArray) = track(transpose, xs)
Base.ctranspose(xs::TrackedArray) = track(ctranspose, xs)

@grad transpose(xs) = data(xs).', Δ -> (reshape(Δ.', size(xs)),)
@grad ctranspose(xs) = data(xs)', Δ -> (reshape(Δ', size(xs)),)

Base.repmat(x::TrackedVecOrMat, a::Integer...) = track(repmat, x, a...)
Base.repmat(x::TrackedVecOrMat, a::Int64...) = track(repmat, x, a...)

@grad function repmat(xs, m, n = 1)
  repmat(data(xs), m, n), function (Δ)
    Δ′ = similar(xs)
    S = size(xs)
    for (i,v) in enumerate(data(Δ))
        d1 = divrem(i-1, S[1]*m)
        x = d1[2] % S[1]+1
        y = d1[1] % S[2]+1
        Δ′[x, y] += v
    end
    return (nobacksies(:repmat, Δ′), nothing, nothing)
  end
end

Base.repeat(A::TrackedArray; kw...) = track_kw(repeat, A; kw...)

@grad function repeat(xs; inner=ntuple(x->1, ndims(A)), outer=ntuple(x->1, ndims(A)))
  repeat(data(xs), inner = inner, outer = outer), function (Δ)
    Δ′ = zero(xs)
    S = size(xs)

    # Loop through each element of Δ, calculate source dimensions, accumulate into Δ′
    for (dest_idx, val) in enumerate(IndexCartesian(), data(Δ))
        # First, round dest_idx[dim] to nearest gridpoint defined by inner[dim], then
        # wrap around based on original size S.
        src_idx = [mod1(div(dest_idx[dim] - 1, inner[dim]) + 1, S[dim]) for dim in 1:length(S)]
        Δ′[src_idx...] += val
    end
    (nobacksies(:repeat, Δ′),)
  end
end

for f in [:vcat, :hcat]
  @eval begin
    # This section is a bit of a hack since julia doesn't have a standardised
    # promotion mechanism for concatenation yet
    # https://github.com/JuliaLang/julia/pull/20815

    # It should support tracked concatenation with rank ∈ (1,2) with a
    # TrackedArray anywhere among the arguments This works as long as base has
    # other functions that captures `(::Union{Vector,RowVector,Matrix}...)`.
    Base.$f(a::Union{TrackedArray,Vector,RowVector,Matrix}...) = track($f, a...)

    # It should support tracked concatenation with rank>2 if the TrackedArray is
    # first
    Base.$f(a::TrackedArray, b::AbstractArray...) = track($f, a, b...)
    Base.$f(a::TrackedArray, b::Union{TrackedArray,Vector,RowVector,Matrix}...) = track($f, a, b...) # resolves ambiguity introduced by previous row

    # It should support tracked concatenation with rank>2 if the TrackedArray is
    # second
    Base.$f(a::Array, b::TrackedArray, c::AbstractArray...) = track($f, a, b, c...)
    Base.$f(a::Union{Vector,RowVector,Matrix}, b::TrackedArray,
            c::Union{TrackedArray,Vector,RowVector,Matrix}...) =
      track($f, a, b, c...) # resolves ambiguity introduced by previous row
  end
end

@grad function vcat(xs...)
  vcat(data.(xs)...), function (Δ)
    start = 0
    Δs = [begin
      i = map(_ -> :, size(xsi)) |> Base.tail
      d = Δ[start+1:start+size(xsi,1), i...]
      start += size(xsi, 1)
      d
    end for xsi in xs]
    return (Δs...,)
  end
end

@grad function hcat(xs...)
  hcat(data.(xs)...), function (Δ)
    start = 0
    Δs = [begin
      d = if ndims(xsi) == 1
        Δ[:, start+1]
      else
        i = map(_ -> :, size(xsi)) |> Base.tail |> Base.tail
        Δ[:, start+1:start+size(xsi,2), i...]
      end
      start += size(xsi, 2)
      d
    end for xsi in xs]
    return (Δs...,)
  end
end

Base.cat(dims, a::TrackedArray, b::AbstractArray...) = track(cat, dims, a, b...)
Base.cat(dims, a::Union{RowVector,Array}, b::TrackedArray, c::AbstractArray...) = track(cat, dims, a, b, c...)

@grad function cat(dims, Xs...)
  cat(dims, data.(Xs)...), function (Δ)
    start = ntuple(i -> 0, Val{ndims(Δ)})
    Δs = [begin
      dim_xs = 1:ndims(xs)
      till_xs = ntuple((i -> i in dims ? (i in dim_xs ? size(xs,i) : 1) : 0), Val{ndims(Δ)})
      xs_in_Δ = ntuple(i -> till_xs[i] > 0 ? (start[i]+1:start[i]+till_xs[i]) : Colon(), Val{ndims(Δ)})
      d = reshape(Δ[xs_in_Δ...],size(xs))
      start = start .+ till_xs
      d
    end for xs in Xs]
    return (nothing, Δs...,)
  end
end

Base.reshape(xs::TrackedArray, dims::Union{Colon,Int64}...) = reshape(xs, dims)
Base.reshape(xs::TrackedArray, dims::Tuple{Vararg{Union{Int64,Colon}}}) = reshape(xs, Base._reshape_uncolon(xs, dims))
Base.reshape(xs::TrackedArray, dims::Tuple{Vararg{Int64}}) = track(reshape, xs, dims)

@grad reshape(xs, dims) = reshape(data(xs), dims), Δ -> (reshape(Δ, size(xs)),nothing)

Base.permutedims(xs::TrackedArray, dims) = track(permutedims, xs, dims)
@grad permutedims(xs, dims) = permutedims(data(xs), dims), Δ -> (permutedims(Δ, invperm(dims)),nothing)

function _kron(mat1::AbstractMatrix,mat2::AbstractMatrix)
    m1, n1 = size(mat1)
    mat1_rsh = reshape(mat1,(1,m1,1,n1))

    m2, n2 = size(mat2)
    mat2_rsh = reshape(mat2,(m2,1,n2,1))

    return reshape(mat1_rsh.*mat2_rsh, (m1*m2,n1*n2))
end

Base.kron(a::TrackedMatrix, b::TrackedMatrix)  = _kron(a, b)
Base.kron(a::TrackedMatrix, b::AbstractMatrix) = _kron(a, b)
Base.kron(a::AbstractMatrix, b::TrackedMatrix) = _kron(a, b)

# Reductions

Base.sum(xs::TrackedArray, dim) = track(sum, xs, dim)
Base.sum(xs::TrackedArray) = track(sum, xs)
Base.sum(f::Union{Function,Type},xs::TrackedArray) = sum(f.(xs))

@grad sum(xs, dim...) = sum(data(xs), dim...),
  Δ -> (zero(xs) .+ Δ, map(_->nothing,dim)...)

Base.prod(xs::TrackedArray, dim) = track(prod, xs, dim)
Base.prod(xs::TrackedArray) = track(prod, xs)
Base.prod(f::Union{Function, Type}, xs::TrackedArray) = prod(f.(xs))

@grad prod(xs) = prod(data(xs)), Δ -> (prod(xs) ./ xs .* Δ,)
@grad prod(xs, dim) = prod(data(xs), dim),
  Δ -> (nobacksies(:sum,
          reshape(.*(circshift.([reshape(data(xs), length(xs))], 1:length(xs)-1)...), size(xs)) .* Δ),
        nothing)

Base.findfirst(xs::TrackedArray, args...) = findfirst(xs.data, args...)

Base.mean(xs::TrackedArray) = track(mean, xs)
Base.mean(xs::TrackedArray, region) = track(mean, xs, region)

Base.maximum(xs::TrackedArray) = track(maximum, xs)
Base.maximum(xs::TrackedArray, region) = track(maximum, xs, region)
Base.minimum(xs::TrackedArray) = track(minimum, xs)
Base.minimum(xs::TrackedArray, region) = track(minimum, xs, region)

LinAlg.dot(xs::TrackedVector, ys::TrackedVector) = track(dot, xs, ys)
LinAlg.dot(xs::AbstractVector, ys::TrackedVector) = track(dot, xs, ys)
LinAlg.dot(xs::TrackedVector, ys::AbstractVector) = track(dot, xs, ys)

@grad dot(xs, ys) = dot(data(xs), data(ys)), Δ -> (Δ .* ys, Δ .* xs)

# Hacks to get std working
Base.std(x::TrackedArray; mean = Base.mean(x)) =
  sqrt.(sum((x .- mean).^2) ./ (length(x)-1))
Base.std(x::TrackedArray, dim; mean = Base.mean(x, dim)) =
  sqrt.(sum((x .- mean).^2, dim) ./ (size(x, dim)-1))

Base.vecnorm(x::TrackedArray, p::Real = 2) =
  sum(abs.(x).^p .+ eps(0f0))^(1/p) # avoid d(sqrt(x))/dx == Inf at 0

@grad mean(xs) = mean(data(xs)), Δ -> (Δ / length(xs),)
@grad mean(xs, region) = mean(data(xs), region), Δ -> (zero(xs) .+ Δ ./ prod(size(xs, region...)),nothing)

@grad function maximum(xs, r...)
  maximum(data(xs), r...), function (Δ)
    Δ′ = zero(xs)
    _, i = findmax(data(xs), r...)
    Δ′[i] = data(Δ)
    return (nobacksies(:maximum, Δ′),map(_->nothing,r)...)
  end
end
@grad function minimum(xs, r...)
  minimum(data(xs), r...), function (Δ)
    Δ′ = zero(xs)
    _, i = findmin(data(xs), r...)
    Δ′[i] = data(Δ)
    return (nobacksies(:minimum, Δ′),map(_->nothing,r)...)
  end
end

# BLAS

Base.diagm(x::TrackedVector) = track(diagm, x)
@grad diagm(x) = diagm(data(x)), Δ -> (diag(Δ),)

for f in :[*, Ac_mul_B, A_mul_Bc, A_mul_Bt, At_mul_B].args
  @eval begin
    import Base.$f
    $f(a::TrackedMatrix, b::TrackedMatrix)  = track($f, a, b)
    $f(a::TrackedMatrix, b::AbstractMatrix) = track($f, a, b)
    $f(a::AbstractMatrix, b::TrackedMatrix) = track($f, a, b)

    $f(a::TrackedMatrix, b::TrackedVector)  = track($f, a, b)
    $f(a::TrackedMatrix, b::AbstractVector) = track($f, a, b)
    $f(a::AbstractMatrix, b::TrackedVector) = track($f, a, b)

    $f(a::TrackedVector, b::TrackedVector)  = track($f, a, b)
    $f(a::TrackedVector, b::AbstractVector) = track($f, a, b)
    $f(a::AbstractVector, b::TrackedVector) = track($f, a, b)
  end
end

@grad a::AbstractMatrix * b::AbstractVecOrMat =
  data(a)*data(b), Δ -> (A_mul_Bt(Δ, b), At_mul_B(a, Δ))

@grad Ac_mul_B(a, b) = Ac_mul_B(data(a), data(b)), Δ -> (A_mul_Bt(Δ, b)', a*Δ)
@grad A_mul_Bc(a, b) = A_mul_Bc(data(a), data(b)), Δ -> (Δ * b, At_mul_B(a, Δ)')

@grad At_mul_B(a, b) = At_mul_B(data(a), data(b)), Δ -> (A_mul_Bt(Δ, b)', a*Δ)
@grad A_mul_Bt(a, b) = A_mul_Bt(data(a), data(b)), Δ -> (Δ * b, At_mul_B(a, Δ)')

# NNlib

using NNlib
import NNlib: softmax, ∇softmax, logsoftmax, ∇logsoftmax, conv, maxpool, meanpool

softmax(xs::TrackedArray) = track(softmax, xs)

@grad softmax(xs) = softmax(data(xs)), Δ -> (nobacksies(:softmax, ∇softmax(data(Δ), data(xs))),)

logsoftmax(xs::TrackedArray) = track(logsoftmax, xs)

@grad logsoftmax(xs) = logsoftmax(data(xs)), Δ -> (nobacksies(:logsoftmax, ∇logsoftmax(data(Δ), data(xs))),)

conv(x::TrackedArray,  w::TrackedArray;  kw...) = track_kw(conv, x, w; kw...)
conv(x::AbstractArray, w::TrackedArray;  kw...) = track_kw(conv, x, w; kw...)
conv(x::TrackedArray,  w::AbstractArray; kw...) = track_kw(conv, x, w; kw...)

@grad conv(x, w; kw...) =
  conv(data(x), data(w); kw...),
    Δ -> nobacksies(:conv,
      (NNlib.∇conv_data(data.((Δ, x, w))...; kw...),
       NNlib.∇conv_filter(data.((Δ, x, w))...; kw...)))

maxpool(x::TrackedArray, k; kw...) = track_kw(maxpool, x, k; kw...)

@grad function maxpool(x, k; kw...)
  y = maxpool(data(x), k; kw...)
  y, Δ -> (nobacksies(:maxpool, NNlib.∇maxpool(data.((Δ, y, x))..., k; kw...)), nothing)
end

meanpool(x::TrackedArray, k; kw...) = track_kw(meanpool, x, k; kw...)

@grad function meanpool(x, k; kw...)
  y = meanpool(data(x), k; kw...)
  y, Δ -> (nobacksies(:maxpool, NNlib.∇meanpool(data.((Δ, y, x))..., k; kw...)), nothing)
end

# Broadcasting

using ForwardDiff: Dual, partials, value

dualify(xs, n) = xs
dualify(xs::AbstractArray, ps) = map(x -> Dual(x, ps), xs)
dualify(xs::Real, ps) = Dual(xs, ps)

unbroadcast(x::Tuple, Δ) =
  x == size(Δ) ? Δ :
    reshape(sum(Δ, filter(n -> n > length(x) || x[n] == 1, 1:ndims(Δ))), x)

unbroadcast(x::Tuple{}, Δ) = sum(Δ)

function getpartial(Δ, x, i)
  @inbounds p = getindex(partials(x), i)
  return Δ * p
end

function ∇broadcast(f, args::Vararg{Any,N}) where N
  sizes = size.(args)
  dargs = map((x,i) -> dualify(data(x), ntuple(j -> i==j, Val{N})), args, ntuple(identity, Val{N}))
  out = broadcast(f, dargs...)
  eltype(out) <: Dual || return out
  y = value.(out)
  back = function (Δ_)
    Δ = data(Δ_)
    Δargs = ntuple(i -> getpartial.(Δ, out, i), Val{N})
    dxs = map((x, Δ) -> unbroadcast(x, Δ), sizes, Δargs)
    nobacksies(:broadcast, dxs)
  end
  # So we can return non-tracked arrays
  track(Call(back, tracker.(args)), y)
end

Base.Broadcast._containertype(::Type{<:TrackedReal}) = TrackedArray
Base.Broadcast._containertype(::Type{<:TrackedArray}) = TrackedArray
Base.Broadcast.promote_containertype(::Type{TrackedArray}, ::Type{TrackedArray}) = TrackedArray
Base.Broadcast.promote_containertype(::Type{Array}, ::Type{TrackedArray}) = TrackedArray
Base.Broadcast.promote_containertype(::Type{TrackedArray}, ::Type{Array}) = TrackedArray
Base.Broadcast.promote_containertype(::Type{TrackedArray}, ct) = TrackedArray
Base.Broadcast.promote_containertype(ct, ::Type{TrackedArray}) = TrackedArray
Base.Broadcast.broadcast_indices(::Type{TrackedArray}, A::Ref) = ()
Base.Broadcast.broadcast_indices(::Type{TrackedArray}, A) = indices(A)

Base.Broadcast.broadcast_c(f, ::Type{TrackedArray}, A, Bs...) = ∇broadcast(f, A, Bs...)


text:
using RDatasets, DataFrames, FixedEffectModels, Base.Test

x = dataset("plm", "Cigar")
x[:pState] = pool(x[:State])
x[:pYear] = pool(x[:Year])

result = reg(Sales ~ NDI, x)
show(result)
predict(result, x)
residuals(result, x)

reg(Sales ~ CPI + (Price = Pimin), x)

result = reg(Sales ~ CPI + (Price = Pimin), x)
predict(result, x)
residuals(result, x)
model_response(result, x)
@test  nobs(result) == 1380
@test_approx_eq  vcov(result)[1]  3.5384578251636785

show(reg(Sales ~ Price |> pState, x))
show(reg(Sales ~ CPI + (Price = Pimin) |> pState, x))

result = reg(Sales ~ Price, x)
@test maxabs(residuals(result, x)[1:10] .- [-39.2637, -37.48801, -34.38801, -36.09743, -36.97446, -43.15547, -41.22573, -40.83648, -34.52427, -28.91617]) <= 1e-4

result = reg(Sales ~ Price |> pState, x, save = true)
@test maxabs(result.augmentdf[:residuals][1:10] .- [-22.08499, -20.33318, -17.23318, -18.97645, -19.85547, -26.1161, -24.20627, -23.87674, -17.62624, -12.01018]) <= 1e-4

ext: jl


text:
# ---
# title: 636. Exclusive Time of Functions
# id: problem636
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Stack
# link: <https://leetcode.com/problems/exclusive-time-of-functions/description/>
# hidden: true
# ---
#
# On a **single-threaded** CPU, we execute a program containing `n` functions.
# Each function has a unique ID between `0` and `n-1`.
#
# Function calls are **stored in a[call
# stack](https://en.wikipedia.org/wiki/Call_stack)**: when a function call
# starts, its ID is pushed onto the stack, and when a function call ends, its ID
# is popped off the stack. The function whose ID is at the top of the stack is
# **the current function being executed**. Each time a function starts or ends,
# we write a log with the ID, whether it started or ended, and the timestamp.
#
# You are given a list `logs`, where `logs[i]` represents the `ith` log message
# formatted as a string `"{function_id}:{"start" | "end"}:{timestamp}"`. For
# example, `"0:start:3"` means a function call with function ID `0` **started at
# the beginning** of timestamp `3`, and `"1:end:2"` means a function call with
# function ID `1` **ended at the end** of timestamp `2`. Note that a function
# can be called **multiple times, possibly recursively**.
#
# A function's **exclusive time** is the sum of execution times for all function
# calls in the program. For example, if a function is called twice, one call
# executing for `2` time units and another call executing for `1` time unit, the
# **exclusive time** is `2 + 1 = 3`.
#
# Return _the **exclusive time** of each function in an array, where the value
# at the _`ith` _index represents the exclusive time for the function with
# ID_`i`.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2019/04/05/diag1b.png)
#
#
#
#     Input: n = 2, logs = ["0:start:0","1:start:2","1:end:5","0:end:6"]
#     Output: [3,4]
#     Explanation:
#     Function 0 starts at the beginning of time 0, then it executes 2 for units of time and reaches the end of time 1.
#     Function 1 starts at the beginning of time 2, executes for 4 units of time, and ends at the end of time 5.
#     Function 0 resumes execution at the beginning of time 6 and executes for 1 unit of time.
#     So function 0 spends 2 + 1 = 3 units of total time executing, and function 1 spends 4 units of total time executing.
#
#
# **Example 2:**
#
#
#
#     Input: n = 1, logs = ["0:start:0","0:start:2","0:end:5","0:start:6","0:end:6","0:end:7"]
#     Output: [8]
#     Explanation:
#     Function 0 starts at the beginning of time 0, executes for 2 units of time, and recursively calls itself.
#     Function 0 (recursive call) starts at the beginning of time 2 and executes for 4 units of time.
#     Function 0 (initial call) resumes execution then immediately calls itself again.
#     Function 0 (2nd recursive call) starts at the beginning of time 6 and executes for 1 unit of time.
#     Function 0 (initial call) resumes execution at the beginning of time 7 and executes for 1 unit of time.
#     So function 0 spends 2 + 4 + 1 + 1 = 8 units of total time executing.
#
#
# **Example 3:**
#
#
#
#     Input: n = 2, logs = ["0:start:0","0:start:2","0:end:5","1:start:6","1:end:6","0:end:7"]
#     Output: [7,1]
#     Explanation:
#     Function 0 starts at the beginning of time 0, executes for 2 units of time, and recursively calls itself.
#     Function 0 (recursive call) starts at the beginning of time 2 and executes for 4 units of time.
#     Function 0 (initial call) resumes execution then immediately calls function 1.
#     Function 1 starts at the beginning of time 6, executes 1 units of time, and ends at the end of time 6.
#     Function 0 resumes execution at the beginning of time 6 and executes for 2 units of time.
#     So function 0 spends 2 + 4 + 1 = 7 units of total time executing, and function 1 spends 1 unit of total time executing.
#
#
# **Example 4:**
#
#
#
#     Input: n = 2, logs = ["0:start:0","0:start:2","0:end:5","1:start:7","1:end:7","0:end:8"]
#     Output: [8,1]
#
#
# **Example 5:**
#
#
#
#     Input: n = 1, logs = ["0:start:0","0:end:0"]
#     Output: [1]
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 100`
#   * `1 <= logs.length <= 500`
#   * `0 <= function_id < n`
#   * `0 <= timestamp <= 109`
#   * No two start events will happen at the same timestamp.
#   * No two end events will happen at the same timestamp.
#   * Each function has an `"end"` log for each `"start"` log.
#
#
## @lc code=start
using LeetCode

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


text:
using DelimitedFiles;
using LinearAlgebra;

full = readdlm("data/y_o_full.txt");
sync = readdlm("data/y_o_sync.txt");
async = readdlm("data/y_o_async.txt");

ϵ = 1e-6;
ilist = 3:2:41;
tlist = full[ilist, 1];
is = 1;
ia = 1;
for i = ilist
    while abs(sync[is, 1] - tlist[i÷2]) > ϵ
        global is += 1
    end
    while abs(sync[ia, 1] - tlist[i÷2]) > ϵ
        global ia += 1
    end
    sync[is, i] = full[i, i];
    async[ia, i] = full[i, i];
end

ef = Vector{Bool}(undef, 41)
es = Vector{Bool}(undef, 41)
ea = Vector{Bool}(undef, 41)

for i = 1:41
    ef[i] = any(!isnan, full[:, i])
    es[i] = any(!isnan, sync[:, i])
    ea[i] = any(!isnan, async[:, i])
end

@show all(ef)
@show all(es)
@show all(ea)

open("data/y_o_sync_add.txt", "w") do ysa
    writedlm(ysa, sync)
end

open("data/y_o_async_add.txt", "w") do yaa
    writedlm(yaa, async)
end

max_stars_repo_name: Typas/Data-Assimilation-Project

max_issues_repo_name: Typas/Data-Assimilation-Project

max_forks_repo_name: Typas/Data-Assimilation-Project


text:
module Codim1
using LinearAlgebra
using Statistics: mean
include("solver.jl")
include("analysis.jl")
include("tools.jl")
include("resolve_point.jl")
include("display.jl")
include("diagnosis.jl")
end  # module


text:
#This opf model is borrowed from JuMP JuMPSupplement
#   reference: https://github.com/mlubin/JuMPSupplement/tree/master/acpower
#   modified to benchmark the Edge Pushing implementation

using JuMP
using Ipopt
using Compat

if parse(Int, ARGS[2]) == 1
    #use JuMP
elseif parse(Int, ARGS[2]) == 2
    include("interface_tape.jl")
    using TapeInterface
elseif parse(Int, ARGS[2]) == 3
    include("interface_check.jl")
    using TapeInterface
else
    @assert false
end

type Bus
    bustype::Int
    name::ASCIIString
    voltage0::Float64
    angle0::Float64
    p_gen::Float64
    q_gen::Float64
    q_min::Float64
    q_max::Float64
    p_load::Float64
    q_load::Float64
    g_shunt::Float64
    b_shunt0::Float64
    b_shunt_min::Float64
    b_shunt_max::Float64
    b_dispatch::Float64
    area::Float64
end

type Branch
    from::Int
    to::Int
    branchtype::Int
    r::Float64
    x::Float64
    c::Float64
    tap0::Float64
    tap_min0::Float64
    tap_max0::Float64
    def0::Float64
    def_min::Float64
    def_max::Float64
    g::Float64
    b::Float64
end

function opf(bus, branch,
    bus_voltage_min, bus_voltage_max,
    branch_tap_min, branch_tap_max,
    p_gen_upper, p_gen_lower)

    simplify = true
    @show simplify
    if parse(Int, ARGS[2]) == 1
        mod = Model(solver=IpoptSolver(max_iter=1), simplify_nonlinear_expressions=simplify)
    elseif parse(Int, ARGS[2]) == 2 || parse(Int, ARGS[2]) == 3
        mod = Model(solver=TapeSolver(IpoptSolver(max_iter=1)), simplify_nonlinear_expressions=simplify)
    else
        @assert false
    end

    nbus = length(bus)
    nbranch = length(branch)
    @variable(mod, bus_voltage_min[bus[i].bustype] <= bus_voltage[i=1:nbus] <= bus_voltage_max[bus[i].bustype], start=1)
    @variable(mod, bus[i].b_shunt_min <= bus_b_shunt[i=1:nbus] <= bus[i].b_shunt_max, start=bus[i].b_shunt0)
    @variable(mod, bus_angle[1:nbus], start=0)

    @variable(mod, branch_tap_min <= branch_tap[1:nbranch] <= branch_tap_max, start=1)
    @variable(mod, branch[i].def_min <= branch_def[i=1:nbranch] <= branch[i].def_max, start = branch[i].def0)

    in_lines = [Int[] for i in 1:nbus] # indices of incoming branches
    out_lines = [Int[] for i in 1:nbus] # indices of outgoing branches
    @show nbranch
    @show nbus
    for i in 1:nbranch
        b = branch[i]
        push!(out_lines[b.from],i)
        push!(in_lines[b.to],i)
        @assert 1 <= b.to <= nbus
    end

    @NLexpression(mod, Gself[k=1:nbus],
        bus[k].g_shunt +
        sum{branch[i].g*branch_tap[i]^2, i=out_lines[k]} +
        sum{ branch[i].g, i=in_lines[k]})

    @NLexpression(mod, Gout[i=1:nbranch],
        (0-branch[i].g*cos(branch_def[i])+branch[i].b*sin(branch_def[i]))*branch_tap[i])
    @NLexpression(mod, Gin[i=1:nbranch],
        (0-branch[i].g*cos(branch_def[i])-branch[i].b*sin(branch_def[i]))*branch_tap[i])

    @NLexpression(mod, Bself[k=1:nbus],
        bus_b_shunt[k] +
        sum{branch[i].b*branch_tap[i]^2 + branch[i].c/2, i=out_lines[k]} +
        sum{ branch[i].b + branch[i].c/2, i=in_lines[k]})

    @NLexpression(mod, Bin[i=1:nbranch],
        (branch[i].g*sin(branch_def[i])-branch[i].b*cos(branch_def[i]))*branch_tap[i])
    @NLexpression(mod, Bout[i=1:nbranch],
        (0-branch[i].g*sin(branch_def[i])-branch[i].b*cos(branch_def[i]))*branch_tap[i])

    # Minimize active power

    @NLobjective(mod, Min, sum{ (bus[k].p_load +
        sum{ bus_voltage[k] * bus_voltage[branch[i].from] *
            (Gin[i] * cos(bus_angle[k] - bus_angle[branch[i].from]) +
             Bin[i] * sin(bus_angle[k] - bus_angle[branch[i].from])), i=in_lines[k] } +
        sum{ bus_voltage[k] * bus_voltage[branch[i].to] *
            (Gout[i] * cos(bus_angle[k] - bus_angle[branch[i].to]) +
             Bout[i] * sin(bus_angle[k] - bus_angle[branch[i].to])), i=out_lines[k] }
        + bus_voltage[k]^2*Gself[k] )^2,
        k in 1:nbus; bus[k].bustype == 2 || bus[k].bustype == 3})
    for k in 1:nbus
        bus[k].bustype == 0 || continue
        # p_load
        @NLconstraint(mod, bus[k].p_gen - bus[k].p_load -
            sum{ bus_voltage[k] * bus_voltage[branch[i].from] *
                (Gin[i] * cos(bus_angle[k] - bus_angle[branch[i].from]) +
                 Bin[i] * sin(bus_angle[k] - bus_angle[branch[i].from])), i=in_lines[k] } -
            sum{ bus_voltage[k] * bus_voltage[branch[i].to] *
                (Gout[i] * cos(bus_angle[k] - bus_angle[branch[i].to]) +
                 Bout[i] * sin(bus_angle[k] - bus_angle[branch[i].to])), i=out_lines[k] }
            - bus_voltage[k]^2*Gself[k] == 0)
        # q_load
        @NLconstraint(mod, bus[k].q_gen - bus[k].q_load -
            sum{ bus_voltage[k] * bus_voltage[branch[i].from] *
                (Gin[i] * sin(bus_angle[k] - bus_angle[branch[i].from]) -
                 Bin[i] * cos(bus_angle[k] - bus_angle[branch[i].from])), i=in_lines[k] } -
            sum{ bus_voltage[k] * bus_voltage[branch[i].to] *
                (Gout[i] * sin(bus_angle[k] - bus_angle[branch[i].to]) -
                 Bout[i] * cos(bus_angle[k] - bus_angle[branch[i].to])), i=out_lines[k] }
             + bus_voltage[k]^2*Bself[k] == 0)
    end
    for k in 1:nbus
        (bus[k].bustype == 2 || bus[k].bustype == 3) || continue
        # q_inj
        @NLconstraint(mod, bus[k].q_min <= bus[k].q_load +
            sum{ bus_voltage[k] * bus_voltage[branch[i].from] *
                (Gin[i] * sin(bus_angle[k] - bus_angle[branch[i].from]) -
                 Bin[i] * cos(bus_angle[k] - bus_angle[branch[i].from])), i=in_lines[k] } +
            sum{ bus_voltage[k] * bus_voltage[branch[i].to] *
                (Gout[i] * sin(bus_angle[k] - bus_angle[branch[i].to]) -
                 Bout[i] * cos(bus_angle[k] - bus_angle[branch[i].to])), i=out_lines[k] } +
            - bus_voltage[k]^2*Bself[k] <= bus[k].q_max)
        # p_inj
        @NLconstraint(mod, 0 <= bus[k].p_load +
            sum{ bus_voltage[k] * bus_voltage[branch[i].from] *
                (Gin[i] * cos(bus_angle[k] - bus_angle[branch[i].from]) +
                 Bin[i] * sin(bus_angle[k] - bus_angle[branch[i].from])), i=in_lines[k] } +
            sum{ bus_voltage[k] * bus_voltage[branch[i].to] *
                (Gout[i] * cos(bus_angle[k] - bus_angle[branch[i].to]) +
                 Bout[i] * sin(bus_angle[k] - bus_angle[branch[i].to])), i=out_lines[k] } +
            + bus_voltage[k]^2*Gself[k] <= p_gen_upper*bus[k].p_gen)
    end

    for i in 1:nbus
        #---- FREEZE THE REFERENCE BUS ANGLE TO ZERO.
        if bus[i].bustype == 3
            setlowerbound(bus_angle[i], 0.0)
            setupperbound(bus_angle[i], 0.0)
        end
        #---- FREEZE ANY DISPATCHABLE SHUNTS.
        if bus[i].b_dispatch == 0
            setlowerbound(bus_b_shunt[i],bus[i].b_shunt0)
            setupperbound(bus_b_shunt[i],bus[i].b_shunt0)
        end
    end

    for i in 1:nbranch
        #---- FREEZE ANY BRANCH TAPS.
        if branch[i].branchtype == 0 || branch[i].branchtype == 3
            setlowerbound(branch_tap[i],1)
            setupperbound(branch_tap[i],1)
        end
        #---- FREEZE CERTAIN PHASE SHIFTERS.
        if branch[i].branchtype != 4
            setlowerbound(branch_def[i],getvalue(branch_def[i]))
            setupperbound(branch_def[i],getvalue(branch_def[i]))
        end
    end

    println("solving")
    solve(mod)

end

function prepdata(bus, branch)

    buses = Array(Bus,size(bus,1))
    refbus = 0
    shunt = 0
    busmap = Dict()
    for i in 1:size(bus,1)
        @assert !haskey(busmap, bus[i,1])
        busmap[bus[i,1]] = i
        buses[i] = Bus(bus[i,2:end]...)
        # rescale
        buses[i].p_gen /= 100
        buses[i].q_gen /= 100
        buses[i].q_min /= 100
        buses[i].q_max /= 100
        buses[i].p_load /= 100
        buses[i].q_load /= 100
        if buses[i].bustype == 3
            refbus += 1
        end
        if buses[i].b_dispatch == 0
            shunt += 1
        end
    end
    @show refbus
    @show shunt

    branches = Array(Branch,size(branch,1))
    branchtap = 0
    phaseshift = 0
    for i in 1:size(branch,1)
        branches[i] = Branch(branch[i,2:end]...,0,0)
        @assert branch[i,1] == i
        branches[i].to = busmap[branches[i].to]
        branches[i].from = busmap[branches[i].from]
        branches[i].g = branches[i].r /(branches[i].r^2 + branches[i].x^2)
        branches[i].b = -branches[i].x /(branches[i].r^2 + branches[i].x^2)
        # rescale
        branches[i].def_min *= 3.14159/180
        branches[i].def_max *= 3.14159/180
        branches[i].def0 *= -3.14159/180
        if branches[i].branchtype == 0 || branches[i].branchtype == 3
            branchtap += 1
        end
        if branches[i].branchtype != 4
            phaseshift += 1
        end
    end
    @show branchtap
    @show phaseshift
    return branches, buses
end

numbuses = ARGS[1]
branches, buses = prepdata(readdlm("IEEE$numbuses.bus"),readdlm("IEEE$numbuses.branch"))

bus_voltage_min = @compat Dict(0 => 0.85, 1 => 0.85, 2 => 0.92, 3 => 0.99)
bus_voltage_max = @compat Dict(0 => 1.15, 1 => 1.15, 2 => 1.08, 3 => 1.01)
branch_tap_min = 0.85
branch_tap_max = 1.15

p_gen_upper = 1.10
p_gen_lower = 0.90
opf(buses, branches,
    bus_voltage_min, bus_voltage_max,
    branch_tap_min, branch_tap_max,
    p_gen_upper, p_gen_lower)


text:
# Autogenerated wrapper script for p3a_jll for x86_64-w64-mingw32-cxx03
export libp3a

using Kokkos_jll
using MicrosoftMPI_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("p3a")
JLLWrappers.@declare_library_product(libp3a, "libp3a.dll")
function __init__()
    JLLWrappers.@generate_init_header(Kokkos_jll, MicrosoftMPI_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libp3a,
        "bin\\libp3a.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@testset "Turing independence" begin
    @model coinflip(y) = begin
        p ~ Beta(1, 1)
        N = length(y)
        for i = 1:N
            y[i] ~ Bernoulli(p)
        end
    end
    model = coinflip([1,1,0])
    model(SampleFromPrior(), LikelihoodContext())
end


text:
"""
    pred_mean(topredict::String, train::DataFrame, test::DataFrame)

Derives the RMSD between the test dataset and the mean of the train dataset.

# Arguments

- `topredict::String`: The name of the column of the dataframe of interest that is to be predicted.
- `train::DataFrame`: The dataset used to derive the mean.
- `test::DataFrame`: The dataset used to test the model.
"""

function pred_mean(topredict::String, train::DataFrame, test::DataFrame)
    predictions = Float64[]
    for i in 1:nrow(test)
        push!(predictions, mean(train[:, Symbol(topredict)]))
    end
    RMSD = rmsd(predictions, test[:, Symbol(topredict)])
    return RMSD
end

max_stars_repo_name: AlexLsn/CSO

max_issues_repo_name: AlexLsn/CSO

max_forks_repo_name: AlexLsn/CSO


text:
#MM dictionary

function MMwithRedirects(mm::Dict{Int64,Int64})
  mmr = Dict{Int64,Int64}()
  for m in mm
    finID = m[2]
    while true
      if haskey(mm, finID)
        finID = mm[finID]
      else
        break;
      end
    end
    mmr[m[1]] = deepcopy(finID)
  end
  return mmr
end

# point to the first instance in the map which may be the same feature
MM = Dict{Int64,Int64}()

MM[103] = 2
MM[273] = 2
MM[105] = 3
MM[274] = 3
MM[106] = 5
MM[270] = 5
MM[111] = 6
MM[275] = 6
MM[281] = 9
# early tests
MM[265] = 1
  MM[286] = 11 # driving data looks convincing
MM[279] = 18
      MM[284] = 20
MM[290] = 22
MM[291] = 23
MM[294] = 26 # repeated with different id as false association
MM[88] = 27
MM[285] = 27
MM[299] = 34
MM[317] = 34
MM[298] = 36
MM[299] = 43
MM[317] = 43

MM[289] = 48
      MM[310] = 54
MM[302] = 56
MM[308] = 59
MM[304] = 60
MM[301] = 61
  MM[287] = 65 # maybe 61, maybe more than one tree
      MM[334] = 68
MM[96] = 75
      MM[355] = 83
MM[107] = 86
MM[323] = 86
      MM[307] = 90
        MM[322] = 92
        MM[347] = 97
MM[261] = 116
      MM[267] = 121
MM[269] = 125
MM[249] = 126
MM[257] = 129
MM[196] = 132
  MM[251] = 132
MM[199] = 134
MM[252] = 134
MM[200] = 138
MM[258] = 138
  MM[182] = 139
  MM[253] = 139 #182
MM[208] = 139
  MM[247] = 144
MM[242] = 145
MM[240] = 149
MM[235] = 158
MM[238] = 164
MM[216] = 165
      MM[225] = 165
MM[227] = 167
  MM[209] = 197
    MM[201] = 193
MM[230] = 219
  MM[239] = 219
MM[226] = 220
#maybes
MM[244] = 135
MM[255] = 135
  MM[309] = 288
      MM[348] = 90#307
  MM[349] = 59# 337

# Beyond 140s
MM[408] = 6
MM[437] = 5
MM[438] = 3
MM[434] = 2
MM[409] = 22
MM[451] = 1
MM[411] = 55
MM[417] = 59
MM[416] = 56
  MM[412] = 20
  MM[409] = 22
  MM[418] = 44
  MM[415] = 48
    MM[378] = 23
      MM[454] = 125
      MM[459] = 121
      MM[457] = 116
      MM[466] = 126
        MM[330] = 91
        MM[472] = 135
        MM[474] = 139
        MM[478] = 144
          MM[605] = 387
          MM[620] = 395
          MM[601] = 381
          MM[606] = 377
            MM[495] = 20
            MM[575] = 22
            MM[553] = 45
            MM[571] = 23
            MM[570] = 26
            MM[568] = 43
              MM[481] = 254
              MM[480] = 167
              MM[231] = 167
              MM[419] = 27
              MM[447] = 9
              MM[421] = 18
              MM[373] = 59
              MM[550] = 90
              MM[560] = 90
              MM[535] = 83
              MM[512] = 99
  # laser to back axel lever arm not modeled in tracker and sighting constraints
  MM[748] = 629
  MM[755] = 631
  MM[741] = 633
  MM[739] = 654
  MM[744] = 654
    MM[761] = 622
    MM[765] = 368#626
    MM[770] = 647
    MM[769] = 653
    MM[771] = 664
# matches at 330s data
MM[828] = 48
MM[823] = 36
MM[817] = 56
MM[801] = 59
MM[835] = 22
MM[836] = 23
MM[831] = 6
MM[838] = 5
MM[846] = 3
MM[845] = 2
MM[832] = 18
MM[834] = 27
MM[826] = 44
MM[816] = 45
MM[818] = 51
MM[803] = 55
MM[847] = 11
MM[830] = 20

#Matches for 350
MM[855] = 125
MM[859] = 121
MM[861] = 116
MM[870] = 138
MM[865] = 126
  MM[867] = 132
MM[872] = 139
MM[883] = 189
MM[876] = 144
MM[882] = 193
MM[871] = 135
MM[879] = 162
MM[884] = 164
MM[877] = 145
MM[889] = 165
MM[881] = 149
  MM[880] = 197
    MM[852] = 1
    MM[863] = 464
      MM[445] = 11
    MM[840] = 26
      MM[504] = 6
      MM[497] = 5
      MM[868] = 467
      MM[486] = 467
    MM[869] = 129 #473
      MM[515] = 85
    MM[780] = 83
      MM[527] = 58
  MM[763] = 372
  MM[793] = 387
    MM[796] = 377
        MM[376] = 59
        MM[337] = 59
          MM[799] = 59
        MM[414] = 36
        MM[499] = 93
          MM[588] = 51
          MM[599] = 56
          MM[374] = 26
          MM[573] = 36
            MM[566] = 36#48
          MM[600] = 48
        MM[594] = 125
        MM[582] = 23
            MM[483] = 61
          # MM[556] = 45
      MM[875] = 134

  MM[848] = 597
  MM[588] = 51
      MM[502] = 463
      MM[562] = 75
        MM[513] = 27
        MM[490] = 27
          MM[470] = 134
          MM[468] = 132
      MM[469] = 138
      MM[473] = 129
    MM[887] = 167
    MM[890] = 167
    MM[886] = 158

      MM[404] = 26
      MM[406] = 26

# more risky associations
MM[297] = 44
MM[592] = 11
MM[589] = 6
MM[580] = 5
MM[627] = 397
MM[632] = 396
MM[617] = 391
MM[767] = 391
MM[613] = 372
  MM[862] = 85
    MM[558] = 54
  MM[496] = 1
  MM[487] = 1
    MM[494] = 48
    MM[574] = 48
    MM[567] = 48
    MM[525] = 48

    MM[488] = 2
    MM[503] = 2
      MM[538] = 43
      MM[533] = 44
        MM[569] = 44
      MM[525] = 48
        MM[556] = 56
      MM[484] = 91
      MM[551] = 82
        MM[562] = 75
        MM[563] = 59

  MM[630] = 371
  MM[626] = 368
    MM[577] = 20

# new for larger T=1400 map, submaps [60;1970;1490], spread=70
MM[3022] = 1489
MM[3021] = 1483
MM[3015] = 1484
MM[3024] = 1493
MM[3019] = 1497
MM[3013] = 1488
MM[3014] = 1494
MM[2995] = 1496
MM[3007] = 1509
MM[2987] = 6 # 1506
MM[2982] = 5 # 1502
MM[2979] = 3 # 1507
MM[2978] = 2 # 1508

MM[1905] = 6
MM[1913] = 5
MM[1911] = 3
MM[1912] = 2
MM[1956] = 2
#triplets
MM[1940] = 1438
MM[1487] = 1438
MM[1935] = 1434
MM[1935] = 1483
MM[1937] = 1436
MM[1489] = 1436
# MM[1929] = 143x
MM[1938] = 1435
MM[1493] = 1435
MM[1927] = 1428
MM[1488] = 1428
MM[1926] = 1427
MM[1503] = 1427
MM[1918] = 1509
MM[1921] = 1420
MM[1496] = 1420
MM[1920] = 1422
MM[1950] = 1494
MM[1946] = 1494
MM[1928] = 1494

MM[1506] = 6
MM[1502] = 5
MM[1507] = 3
MM[1508] = 2

## post 11AM
# from submaps [80;1490], spread=80
MM[1433] = 138
MM[1497] = 138
MM[1418] = 125
MM[1509] = 125
MM[1488] = 126
MM[1428] = 126
MM[1435] = 134
MM[1493] = 134
MM[1434] = 135
MM[1483] = 135
MM[1436] = 132
MM[1489] = 132
MM[1490] = 139
MM[1420] = 116
MM[1496] = 116
MM[1422] = 121

# from submaps [80;950], spread=80
MM[985] = 6
MM[972] = 5
MM[971] = 3
MM[970] = 2

MM[936] = 116
MM[995] = 125
MM[929] = 144
MM[924] = 139 # might be 923
MM[925] = 132
MM[932] = 138
MM[928] = 126
MM[927] = 135
MM[926] = 134
MM[931] = 145
MM[935] = 149

MM[1019] = 20
MM[1018] = 22

# submaps from [3080;2550], spread=100
# upper left corner loop closures
MM[3176] = 2645
MM[3164] = 2628
MM[3180] = 2644
MM[3169] = 2632
MM[3168] = 2638
MM[3176] = 2649 #?
MM[3167] = 2633
MM[3162] = 2625
MM[3170] = 2636
MM[2646] = 2636
MM[3179] = 2629
MM[2647] = 2629
MM[3158] = 2624
MM[3150] = 2613
MM[3141] = 2610
MM[3139] = 2608
MM[3153] = 2614
MM[3137] = 2606

## post 3:20PM, from submaps [100;2550;3080], spread=100
MM[2463] = 138
MM[2455] = 126
MM[2466] = 134
MM[2456] = 132
MM[2464] = 135
MM[2467] = 189#139
MM[2477] = 139# something wrong with (not found) 187
MM[2469] = 144
MM[2480] = 164
MM[2471] = 162
MM[2489] = 165
MM[2481] = 167

## from submaps [100;2900], spread=100
MM[2852] = 138
MM[2849] = 126
MM[2847] = 134
MM[2837] = 132
MM[2840] = 135
MM[2820] = 139
MM[2838] = 189
MM[2813] = 144
MM[2824] = 164
MM[2809] = 162
MM[2801] = 193
MM[2816] = 165
MM[2863] = 125
MM[2866] = 121
MM[2853] = 116
MM[2874] = 6
MM[2864] = 5
MM[2869] = 3
MM[2868] = 2

MM[2892] = 48
MM[2899] = 36
MM[2896] = 45
MM[2889] = 26

## post 5:30 solve, from submaps [700; 1200], spread=100
MM[1243] = 737
MM[1239] = 722
MM[1230] = 717
MM[1247] = 732
MM[1246] = 725
MM[1227] = 707
MM[1225] = 702
MM[1255] = 702
MM[1250] = 702

MM[1202] = 670
MM[1300] = 758
MM[1194] = 684
MM[1295] = 751
MM[1275] = 695

MM[1220] = 700
MM[1267] = 1220

MM[1211] = 691
MM[1260] = 1211
MM[1269] = 1208
MM[1273] = 1216

# from submaps [700;1800], spread=100
MM[1789] = 737
MM[1780] = 722
MM[1770] = 717
MM[1788] = 732
MM[1769] = 707
MM[1767] = 702
MM[1773] = 715
MM[1741] = 698
MM[1848] = 786
MM[1831] = 754
MM[1791] = 654
MM[1792] = 633
MM[1797] = 746
MM[1731] = 686
MM[1793] = 742
MM[1740] = 687
MM[1863] = 747
MM[1833] = 658

# from submaps [100;650;2300], spread=100
MM[2262] = 737
MM[2260] = 722
MM[2266] = 717
MM[2263] = 732
MM[2269] = 707

MM[2300] = 658
MM[747] = 658
MM[2303] = 629
MM[2306] = 631
MM[2297] = 633
MM[2277] = 654
MM[2295] = 746
MM[2278] = 742
MM[2296] = 743
MM[742] = 686
MM[743] = 676

# from submaps [100; 2950]
MM[2990] = 116
MM[2989] = 20
MM[2984] = 18
MM[2966] = 9
MM[2964] = 27
MM[2946] = 86
MM[2945] = 84

MM[2964] = 27
MM[2908] = 43
MM[2909] = 54

MM[2910] = 82
MM[2962] = 75

MM[2913] = 55

## after 8:50 PM solve
MM[3027] = 144
MM[3036] = 164
MM[3028] = 162

MM[2879] = 11
MM[2878] = 22
MM[2902] = 56
MM[2905] = 59
MM[2897] = 44
MM[2880] = 20
MM[2885] = 18
MM[2891] = 27
MM[2877] = 9

# more risky, submaps [650, 1120, 2300], spread=100
MM[2369] = 1147
MM[2310] = 1138
MM[2364] = 1137
MM[2326] = 1163
MM[2316] = 1141
MM[2319] = 1152

MM[1153] = 658
MM[1154] = 629
MM[1157] = 631

MM[1141] = 622
MM[1163] = 647

MM[2327] = 664
MM[2320] = 653

# Thu morn, from submaps [650;2300;1750], spread=100
MM[2308] = 652
MM[1734] = 683

MM[1727] = 676
MM[1795] = 676

MM[2314] = 1812
MM[1804] = 631
MM[1685] = 631
MM[1684] = 633
MM[1812] = 670
MM[2313] = 1732
MM[1805] = 1732
MM[1732] = 684
MM[1818] = 1734

MM[1733] = 681

MM[1803] = 695

#from submaps [400;600;800], spread=100
MM[800] = 381
MM[757] = 652
MM[788] = 652
MM[758] = 670

#from submaps [400;600;2100]
# from submaps [650;2150], spread=100
MM[2094] = 658
MM[2098] = 629
MM[2109] = 654
MM[2092] = 633

MM[2132] = 687
MM[2170] = 686
MM[2128] = 676
MM[2136] = 684
MM[2127] = 670
MM[2137] = 683

MM[2086] = 627

# Saturday, fromSubmaps=[1130,1350], spread=100
MM[1340] = 1134
MM[1321] = 1138
MM[1318] = 1137
MM[1341] = 1133
MM[1403] = 1120
MM[1388] = 1111
MM[1405] = 1113
MM[1384] = 1096
MM[1379] = 1088
MM[1382] = 1105
MM[1398] = 1098
MM[1392] = 1102

# Unsure associations
# MM[923] = 139 # ref 924
# MM[990] = 18
# MM[1298] = 654
# MM[1153] = 754
# MM[2934] = 68
# MM[2147] = 694

# purposeful bad associates, multimodal constraints
# MM[232] = 167
# MM[225] = 167
# MM[260] = 126
# MM[121] = 5
# MM[314] = 109
# MM[294] = 45 #  repeat of 26
# MM[344] = 44
# # slalums
# MM[284] = 105
# MM[556] = 45

#outrageous
# MM[210] = 100
# MM[2658] = 706
# MM[2616] = 693
# MM[2755] = 664
# MM[2720] = 773
# MM[2796] = 715
# MM[2607] = 633
# MM[2778] = 600

# randomly permute some fraction of initially correct loop closures
function MMRandErrs(MM::Dict{Int64,Int64};frac=0.1)
  MMr = Dict{Int64, Int64}()

  allkeys = collect(keys(MM))
  len = length(allkeys)
  @show clen = round(Int, frac*len)
  @show clen2 = round(Int, 0.5*clen)
  p1 = Dict{Int,Void}()
  p2 = Dict{Int,Void}()

  for m in MM    MMr[m[1]] = m[2]; end

  while length(p1) < clen2 && length(p2) < clen2
    prop1 = allkeys[rand(1:len)]
    if !haskey(p1, prop1) && !haskey(p2, prop1)
      for i in 1:100
        prop2 = allkeys[rand(1:len)]
        # want to try MMr[prop1] = MMr[prop2] & MMr[prop2] = MMr[prop1]
        mm1 = MMr[prop1]
        mm2 = MMr[prop2]
        # @show !haskey(p2, prop2), !haskey(p1, prop2), prop1 > mm2, prop2 > mm1
        if !haskey(p2, prop2) && !haskey(p1, prop2) && prop1 > mm2 && prop2 > mm1
          p1[prop1] = nothing
          p2[prop2] = nothing
          MMr[prop1] = mm2
          MMr[prop2] = mm1
          break
        end
      end
    end
  end

  # for i in 1:(clen2-1)
  #   temp = MMr[p[clen2+i]]
  #   MMr[p1[clen2+i]] = MMr[p1[i]]
  #   MMr[p1[i]] = temp
  # end
  return MMr
end

function showMMpermutes(MMr, MMe)
  for m in MMr  if MMe[m[1]] != m[2]   @show m[1], m[2], MMe[m[1]] end end
  nothing
end

function testMMRandErr(MM::Dict{Int64,Int64};frac::Float64=0.1)
  MMe = MMRandErrs(MM,frac=frac)
  fi = map(Int, [m[2] == MMe[m[1]] for m in MM])
  @show sum(fi.==0)/(sum(fi.==1)+0.0+sum(fi.==0))
end

#


text:
abstract type Asset end

abstract type Property <: Asset end
abstract type Investment <: Asset end
abstract type Cash <: Asset end

abstract type House <: Property end
abstract type Apartment <: Property end

abstract type FixedIncome <: Investment end
abstract type Equity <: Investment end

abstract type Art end

struct Painting <: Art
    artist::String
    title::String
end

"""
subtypetree(roottype, level=1, indent=4)
    Display the entire type hierarchy starting from the specified 'roottype'
"""
function subtypetree(roottype, level=1, indent=4)
    level == 1 && println(roottype)
    for s in subtypes(roottype)
        println(join(fill(" ", level * indent)) * string(s))
        subtypetree(s, level + 1, indent)
    end
end

# simple functions on abstract types
describe(a::Asset) = "Something valuable"
describe(e::Investment) = "Financial investment"
describe(e::Property) = "Physical property"

"""
    location(p::Property)
Returns the location of the property as a tuple of (latitude, longitude).
"""
location(p::Property) = error("Location is not defined int he concrete type")

function walking_distance(p1::Property, p2::Property)
    loc1 = location(p1)
    loc2 = location(p2)
    return abs(loc1.x - loc2.x) + abs(loc1.y - loc2.y)
end

mutable struct Stock <: Equity
    symbol::String
    name::String
end

const Thing = Union{Painting,Stock}

struct BasketOfThings
    stocks::Vector{Thing}
    reason::String
end

function describe(s::Stock)
    return s.symbol * "(" * s.name * ")"
end

struct StockHolding{T <: Real,P <: AbstractFloat}
    stock::Stock
    quantity::T
    price::P
    marketvale::P
end

abstract type Holding{P} end

struct StockHolding3{T,P} <: Holding{P}
    stock::Stock
    quantity::T
    price::P
    marketvalue::P
end

struct CashHolding{P} <: Holding{P}
    currency::String
    amount::P
    marketvalue::P
end


text:
using ModelingToolkit
using LinearAlgebra
using SparseArrays: sparse
using Test

@variables a,b,c,d,e,f,g,h,i

# test hashing
aa = a; # old a

@variables a

@test isequal(a, aa)
@test hash(a) == hash(aa)

@test isequal(get_variables(a+aa+1), [a])

@test hash(a+b ~ c+d) == hash(a+b ~ c+d)

# test some matrix operations don't throw errors
X = [0 b c; d e f; g h i]
@test iszero(simplify(det(X) - ((d * ((b * i) - (c * h))) + (g * ((b * f) - (c * e))))))
F = lu(X)
@test F.p == [2, 1, 3]
R = simplify.(F.L * F.U - X[F.p, :], polynorm=true)
@test iszero(R)
@test simplify.(F \ X) == I
@test ModelingToolkit._solve(X, X) == I
inv(X)
qr(X)

X2 = [0 b c; 0 0 0; 0 h 0]
@test_throws SingularException lu(X2)
F2 = lu(X2, check=false)
@test F2.info == 1

# test operations with sparse arrays and Operations
# note `isequal` instead of `==` because `==` would give another Operation

# test that we can create a sparse array of Operation
Oarray = zeros(Num, 2,2)
Oarray[2,2] = a
@test isequal(sparse(Oarray), sparse([2], [2], [a]))

# test Operation * sparse
@test isequal(a * sparse([2], [2], [1]), sparse([2], [2], [a * 1]))

# test sparse{Operation} + sparse{Operation}
A = sparse([2], [2], [a])
B = sparse([2], [2], [b])
@test isequal(A + B, sparse([2], [2], [a+b]))

# test sparse{Operation} * sparse{Operation}
C = sparse([1, 2], [2, 1], [c, c])
D = sparse([1, 2], [2, 1], [d, d])

@test isequal(C * D, sparse([1,2], [1,2], [c * d, c * d]))

@parameters t σ ρ β
@variables x(t) y(t) z(t)
D = Differential(t)
Dx = Differential(x)
Dy = Differential(y)
Dz = Differential(z)

eqs = [D(x) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]
J = Num[Dx(eqs[1].rhs) Dy(eqs[1].rhs) Dz(eqs[1].rhs)
 Dx(eqs[2].rhs) Dy(eqs[2].rhs) Dz(eqs[2].rhs)
 Dx(eqs[3].rhs) Dy(eqs[3].rhs) Dz(eqs[3].rhs)]

J = expand_derivatives.(J)
using LinearAlgebra
luJ = lu(J,Val(false))

using ModelingToolkit
@variables M[1:2,1:2]
inv(M)

@variables b[1:2]
M = [1 0; 0 2]
M \ b
M \ reshape(b,2,1)
M = [1 1; 0 2]
M \ reshape(b,2,1)

M = [1 a; 0 2]
M \ b
M \ [1, 2]

# test det
@variables X[1:4,1:4]
d1 = det(X, laplace=true)
d2 = det(X, laplace=false)
_det1 = eval(build_function(d1, X))
_det2 = eval(build_function(d2, X))
A = [1 1 1 1
     1 0 1 1
     1 1 0 1
     1 1 1 0]
@test _det1(A) == -1
@test _det2(A) == -1

@variables X[1:3,1:3]
d1 = det(X, laplace=true)
d2 = det(X, laplace=false)
_det1 = eval(build_function(d1, X))
_det2 = eval(build_function(d2, X))
A = [1 1 1
     1 0 1
     1 1 1]
@test _det1(A) == 0
@test _det2(A) == 0


text:
function make_model(in_path::AbstractString; out_path="./autogeneratedmodel",
                    host="bacteria", model="Kinetics", lang="julia")
  # default argument names
  InPath = in_path
  OutPath = out_path
  Host = host
  Model = model
  Lang = lang

  path_head = Base.@__DIR__
  # Error report: line number --> error message
  error_report = Dict{Int64, Array}()

  # include(joinpath(path_head, "preprocessor3.jl"))
  # include("preprocessor3.jl")
  println("\n---------loading sentences------------")
  # inputSentencesArray = Array of pairs of (line, line #)
  inputSentencesArray = getRidOfNewlineHashcommentBlankline(InPath)
  println("\n>>>>>tokenizing sentences" * ">"^10)
  # tokenizedInputSentencesArray = Array of pairs of (tokens, line #)
  tokenizedInputSentencesArray = sentenceTokenization2(inputSentencesArray)
  println("\n>>>>>normalizing and tagging" * ">"^10)
  reservedWordsPath = joinpath(path_head, "reservedWords.jl")
  # reservedWordsPath = "reservedWords.jl"
  # taggedSentencesArray = Array of pairs of (taggedTokens, line #)
  taggedSentencesArray = tokenClassification2(tokenizedInputSentencesArray, reservedWordsPath)
  # reshape for output observation
  printTagSen = [["$(y[1])/$(y[2])" for y in x] for (x, id) in taggedSentencesArray]
  # println(typeof(printTagSen))
  foreach(println, [join(ts, "  ") for ts in printTagSen])

  # include(joinpath(path_head, "biosymDecoder.jl"))
  # include("biosymDecoder.jl")
  println("\n------------information extraction------------")
  include(reservedWordsPath)
  # BioSymVerbInfo = Array of pairs of ((senType, line#), BioSym, paraSet)
  BioSymVerbInfo = extractBioSymVerbInformation(taggedSentencesArray,
                   reservedWords["SentenceType"], error_report)
  # println("print by foreach:")
  foreach(println, BioSymVerbInfo)
  println("\n>>>>>decoding bio symbols" * ">"^10)
  decodingBioSymGroups(BioSymVerbInfo, error_report)
  # BioSymVerbInfo = Array of pairs of ((senType, line#), BioSym, paraSet)
  println("\n>>>>>decoding bio symbols results" * ">"^10)
  printArrayOfTripleToken(BioSymVerbInfo, AbstractString)
  println("\n>>>>>type conversion dictionary" * ">"^10)
  typeConversionDict = setUpSymbolConvertionDict(BioSymVerbInfo, error_report)
  foreach(println, typeConversionDict)
  println("\n>>>>>replace each biosymbol string with generalBioSym
          composite" * ">"^10)
  # newVerbBiosymInfo = Array of ((senType, line#), BioSym, paraSet)
  newVerbBiosymInfo = replaceEachBioSymWithGeneralBioSym(BioSymVerbInfo,
                      typeConversionDict, error_report)
  printArrayOfTripleToken(newVerbBiosymInfo, generalBioSym)

  # include(joinpath(path_head, "semanticChecking.jl"))
  # include("semanticChecking.jl")
  println("\n------------right before IR generation-----------------")
  # preIR_DS = array of tuple ((verb, line #), reactant, product, catalyst)
  preIR_DS = semanticCheckingForEachTriplet(newVerbBiosymInfo,
             typeConversionDict, error_report)
  printArrayOfTripleToken(preIR_DS, generalBioSym)

  #######################
  # TO HERE, no more error report after this line
  #######################
  # print errors
  if length(error_report) != 0
    println("\n\n\n#################### Error report: Useful tips ########")
    for lineNum in sort(collect(keys(error_report)))
      println("in line $lineNum -->", )
      for err in error_report[lineNum]
        println("\t$(err)")
      end
    end
    exit()
  end

  println("\n------------model generation preparation------------------")
  # include(joinpath(path_head, "modelGeneration.jl"))
  # include("modelGeneration.jl")
  sys2userDict = Dict()
  for (key, val) in typeConversionDict
    sys2userDict[val] = key
  end
  # reform to match with NML_V1
  # preIR_DS = array of tuple ((verb, line #), [reactant, product, catalyst])
  (rnx, txtl, rnx_set, mRNA_set, m_protein_set) = preparing_rnxList_txtlDict(
    preIR_DS, sys2userDict, typeConversionDict)

  println(">>>>>rnx list" * ">"^10)
  foreach(x->println(x.rnxName), rnx)
  # Sorting, insert "BIOMASS" here
  (sorted_rnx_species_array, rnx_species2index_dict, rnx_index2species_dict,
   sorted_all_species_array, all_species2index_dict, all_index2species_dict, extra_species_num) =
    sorting_species_list(rnx_set, mRNA_set, m_protein_set, sys2userDict)
  println("------------")
  foreach(println, sorted_rnx_species_array)
  println("------------")
  foreach(println, sorted_all_species_array)

  # write program to disk
  output_file = Dict{String, String}()  # put all files in this dict
  # generate SBML
  XMLfile = generate_SBML_file(sys2userDict, rnx_set, rnx)
  output_file["network.xml"] = XMLfile
  # stoichiometric_matrix
  stoichiometric_matrix_buffer, stoichiometric_matrix = build_stoichiometric_matrix_buffer(rnx, rnx_species2index_dict)
  output_file["stoichiometry.dat"] = stoichiometric_matrix_buffer
  # generate MAT
  MATfile = generate_MAT_file(rnx, sorted_rnx_species_array, stoichiometric_matrix)

  println("\n------------$Model Modeling Framework-----")
  # mapping: Lang to (FBA/FVA distribution, Strategy, function prefix,
  #                   file_suffix)
  DISTRIBUTION = Dict(
      "julia"   => ("distribution/dis.jl", "JuliaStrategy.jl", "jl_", ".jl"),
      "python2" => ("distribution/dis.py2", "Python2Strategy.jl", "py2_", ".py"),
      "python"  => ("distribution/dis.py3", "Python3Strategy.jl", "py3_", ".py"),
      "python3" => ("distribution/dis.py3", "Python3Strategy.jl", "py3_", ".py"),
      "matlab"  => ("distribution/dis.m", "MATLABStrategy.jl", "m_", ".m")
  )
  target_lang = Lang
  model_frame = Model
  if target_lang in ["julia", "python3", "python2", "python", "matlab"]
    println("Generate model in $target_lang")
    # fn_prefix = DISTRIBUTION[target_lang][3]
    file_suffix = DISTRIBUTION[target_lang][4]
    dis_file_folder = DISTRIBUTION[target_lang][1]
  else
    println("Generate model in julia as default")
    # include(joinpath(path_head, DISTRIBUTION["julia"][2]))
    # fn_prefix = "jl_"
    target_lang = "julia"
    file_suffix = ".jl"
    dis_file_folder = DISTRIBUTION["julia"][1]
  end

  # mapping: Lang to functions

  if target_lang in ["python", "python3"]
    target_lang = "python"
  end
  fnmap = Dict(
      "julia"   => (jl_build_kinetics_buffer, jl_build_data_dictionary_buffer,
                    jl_build_simulation_buffer, jl_build_solveODEBalances_buffer,
                    jl_generate_FBA_data_dictionary, jl_build_FBA_solve_buffer,
                    jl_build_rRules_buffer, jl_build_rFBA_solve_buffer,
                    jl_build_rFBA_include_buffer,
                    jl_build_FBA_include_buffer, jl_build_kinetics_include_buffer),
      "python2" => (py2_build_kinetics_buffer, py2_build_data_dictionary_buffer,
                    py2_build_simulation_buffer, py2_build_solveODEBalances_buffer,
                    py2_generate_FBA_data_dictionary, py2_build_FBA_solve_buffer,
                    py2_build_rRules_buffer, py2_build_rFBA_solve_buffer),
      "python" => (py3_build_kinetics_buffer, py3_build_data_dictionary_buffer,
                    py3_build_simulation_buffer, py3_build_solveODEBalances_buffer,
                    py3_generate_FBA_data_dictionary, py3_build_FBA_solve_buffer,
                    py3_build_rRules_buffer, py3_build_rFBA_solve_buffer),
      "matlab"  => (m_build_kinetics_buffer, m_build_data_dictionary_buffer,
                    m_build_simulation_buffer, m_build_solveODEBalances_buffer,
                    m_generate_FBA_data_dictionary, m_build_FBA_solve_buffer,
                    m_build_rRules_buffer, m_build_rFBA_solve_buffer)
  )

  # Generate FBA model
  if model_frame in ["FBA", "FVA"]
    # generate DataDictionary.file_suffix
    if target_lang == "matlab"
      data_dictionary_buffer, maximize_product_buffer = fnmap[target_lang][5](
             rnx, sorted_rnx_species_array, extra_species_num)
      output_file["DataDictionary" * file_suffix] = data_dictionary_buffer
      output_file["maximizeProductDictionary" * file_suffix] = maximize_product_buffer
    else
      data_dictionary_buffer = fnmap[target_lang][5](rnx,
                                 sorted_rnx_species_array, extra_species_num)
      output_file["DataDictionary" * file_suffix] = data_dictionary_buffer
    end
    # generate Include.jl for Julia
    if target_lang == "julia"
      output_file["Include.jl"] = fnmap["julia"][length(fnmap["julia"])-1]()
    end
    # generate Solve.file_suffix
    output_file["Solve" * file_suffix] = fnmap[target_lang][6]()
    output_dir = OutPath
    # check for directory
    if !(isdir(output_dir))
      mkpath(output_dir)
    end
    # write to out_path
    for (file_name, file) in output_file
      write(joinpath(output_dir, file_name), file)
    end
    # copy to out_path
    for file in readdir(joinpath(path_head, dis_file_folder))
      if occursin(file_suffix, file)
        println(file)
        f = joinpath(path_head, dis_file_folder, file)
        cp(f, joinpath(output_dir, file); force=true)
      end
    end
    # copy model file (in SEML) to out_path
    cp(InPath, joinpath(output_dir, basename(InPath)); force=true)
    # mv MATfile to out_path
    mv(MATfile, joinpath(output_dir, MATfile); force=true)

  # Generate rFBA model
  elseif model_frame in ["rFBA"]
    # generate DataDictionary.file_suffix
    if target_lang == "matlab"
      data_dictionary_buffer, maximize_product_buffer = fnmap[target_lang][5](
             rnx, sorted_rnx_species_array, extra_species_num)
      output_file["DataDictionary" * file_suffix] = data_dictionary_buffer
      output_file["maximizeProductDictionary" * file_suffix] = maximize_product_buffer
    else
      data_dictionary_buffer = fnmap[target_lang][5](rnx,
                                 sorted_rnx_species_array, extra_species_num)
      output_file["DataDictionary" * file_suffix] = data_dictionary_buffer
    end
    # Include.jl for julia
    if target_lang == "julia"
      output_file["Include.jl"] = fnmap["julia"][length(fnmap["julia"])-2]()
    end
    # generate Solve.file_suffix
    output_file["rRules" * file_suffix] = fnmap[target_lang][7]()
    output_file["Solve" * file_suffix] = fnmap[target_lang][8]()
    output_dir = OutPath
    # check for directory
    if !(isdir(output_dir))
      mkpath(output_dir)
    end
    # write to out_path
    for (file_name, file) in output_file
      write(joinpath(output_dir, file_name), file)
    end
    # copy to out_path
    for file in readdir(joinpath(path_head, dis_file_folder))
      if occursin(file_suffix, file)
        println(file)
        f = joinpath(path_head, dis_file_folder, file)
        cp(f, joinpath(output_dir, file); force=true)
      end
    end
    # copy model file (in SEML) to out_path
    cp(InPath, joinpath(output_dir, basename(InPath)); force=true)
    # mv MATfile to out_path
    mv(MATfile, joinpath(output_dir, MATfile); force=true)

  # Generate Kinetic model
  else
    if model_frame != "Kinetics"
      println("Generate Kinetics model as default")
    end
    # generate Kinetics.suffix
    (kinetics_buffer, Monod_const, W_array, disassociation_const) =
      fnmap[target_lang][1](all_species2index_dict, rnx, txtl, sys2userDict)
    output_file["Kinetics" * file_suffix] = kinetics_buffer
    # generate DataDictionary.suffix
    host_type = Host
    data_dictionary_buffer = fnmap[target_lang][2](host_type, sorted_all_species_array,
      all_species2index_dict,
      sorted_rnx_species_array, rnx, txtl, Monod_const, W_array, disassociation_const,
      collect(mRNA_set), collect(m_protein_set))
    output_file["DataDictionary" * file_suffix] = data_dictionary_buffer
    # generate Balances.suffix
    ODE_simulation_buffer = fnmap[target_lang][3](extra_species_num)
    output_file["Balances" * file_suffix] = ODE_simulation_buffer
    # generate SolveBalances.suffix
    solveODEBalances_buffer = fnmap[target_lang][4](sorted_all_species_array,
      all_species2index_dict, collect(mRNA_set), collect(m_protein_set))
    output_file["SolveBalances" * file_suffix] = solveODEBalances_buffer
    # check for output folder
    output_dir = OutPath
    if !(isdir(output_dir))
      mkpath(output_dir)
    end
    if target_lang == "julia"
      # cp(joinpath(path_head, "Include.jl"), joinpath(output_dir, "Include.jl"); force=true)
      output_file["Include.jl"] = fnmap["julia"][length(fnmap["julia"])]()
    end
    # copy model (in SEML) to out_path
    cp(InPath, joinpath(output_dir, basename(InPath)); force=true)
    # mv MATfile to out_path
    mv(MATfile, joinpath(output_dir, MATfile); force=true)
    # write to output folder
    for (file_name, file) in output_file
      write(joinpath(output_dir, file_name), file)
    end
  end

end

# # test
# include("preprocessor3.jl")
# include("biosymDecoder.jl")
# include("semanticChecking.jl")
# include("modelGeneration.jl")
# include("JuliaStrategy.jl")
# include("Python2Strategy.jl")
# include("Python3Strategy.jl")
# include("MATLABStrategy.jl")
# include("SBMLgen.jl")
# include("MATgen.jl")
#
# make_model("../example/testcase/kinetcase.dat", model="Kinetics")

max_stars_repo_name: varnerlab/SEML

max_issues_repo_name: varnerlab/SEML

max_forks_repo_name: varnerlab/SEML


text:
struct GymEnv{T,Ta,To,P} <: AbstractEnv
    pyenv::P
    observation_space::To
    action_space::Ta
    state::P
end
export GymEnv

mutable struct AtariEnv{IsGrayScale,TerminalOnLifeLoss,N,S <: AbstractRNG} <: AbstractEnv
    ale::Ptr{Nothing}
    name::String
    screens::Tuple{Array{UInt8,N},Array{UInt8,N}}  # for max-pooling
    actions::Vector{Int}
    action_space::Base.OneTo{Int}
    observation_space::Space{Array{ClosedInterval{UInt8},N}}
    noopmax::Int
    frame_skip::Int
    reward::Float32
    lives::Int
    rng::S
end
export AtariEnv

mutable struct OpenSpielEnv{S,G} <: AbstractEnv
    state::S
    game::G
end
export OpenSpielEnv

# ??? can we safely ignore the `game` field here
Base.hash(e::OpenSpielEnv, h::UInt) = hash(e.state, h)
Base.:(==)(e::OpenSpielEnv, ee::OpenSpielEnv) = e.state == ee.state

mutable struct SnakeGameEnv{A,N,G} <: AbstractEnv
    game::G
    latest_snakes_length::Vector{Int}
    latest_actions::Vector{CartesianIndex{2}}
    is_terminated::Bool
end
export SnakeGameEnv

struct AcrobotEnvParams{T}
    link_length_a::T # [m]
    link_length_b::T # [m]
    link_mass_a::T # : [kg] mass of link 1
    link_mass_b::T # : [kg] mass of link 2
    # : [m] position of the center of mass of link 1
    link_com_pos_a::T
    # : [m] position of the center of mass of link 2
    link_com_pos_b::T
    # : Rotation related parameters
    link_moi::T
    max_torque_noise::T
    # : [m/s] maximum velocity of link 1
    max_vel_a::T
    # : [m/s] maximum velocity of link 2
    max_vel_b::T
    # : [m/s2] acceleration due to gravity
    g::T
    # : [s] timestep
    dt::T
    # : maximum steps in episode
    max_steps::Int
end

export AcrobotEnvParams

mutable struct AcrobotEnv{T,R <: AbstractRNG} <: AbstractEnv
    params::AcrobotEnvParams{T}
    state::Vector{T}
    action::Int
    done::Bool
    t::Int
    rng::R
    reward::T
    # difference in second link angular acceleration equation
    # as per python gym
    book_or_nips::String
    # array of available torques based on actions
    avail_torque::Vector{T}
end

export AcrobotEnv


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

"""
ACASX Study comparing performance of GP, MC, MCTS, GE Trees.
Single-threaded versions are used for more stable comparison.
combine_and_plot()
"""
module ACASX_Compare_Tree

export combine_and_plot, combine_mc_logs, combine_mcts_logs, combine_ge_logs, combine_gp_logs,
    master_log, master_plot, run_main

import Compat: ASCIIString, UTF8String
using GrammarExpts, GBDTs
using ExprSearch: GP, MC, MCTS, GE
using ACASX_GP_Tree, ACASX_GE_Tree, ACASX_MC_Tree, ACASX_MCTS_Tree
using LogJoiner
using RLESUtils, Loggers, MathUtils, LatexUtils, Sweeper
using DataFrames
using PGFPlots, TikzPictures

const CONFIG = "nvn_libcas098smallfilt_10K"
#const CONFIG = "nvn_dascfilt"
const STUDYNAME = "ACASX_Compare_Tree"
const MC_NAME = "ACASX_MC_Tree"
const MCTS_NAME = "ACASX_MCTS_Tree"
const GE_NAME = "ACASX_GE_Tree"
const GP_NAME = "ACASX_GP_Tree"

const RESULTDIR = joinpath(dirname(@__FILE__), "..", "..", "..", "results")

const MASTERLOG_FILE = joinpath(RESULTDIR, STUDYNAME, "masterlog.csv.gz")
const PLOTLOG_FILE =  joinpath(RESULTDIR, STUDYNAME, "plotlog.csv.gz")
const PLOTFILEROOT = joinpath(RESULTDIR, STUDYNAME, "plots")

resultpath(dir::ASCIIString="") = joinpath(RESULTDIR, dir)
studypath(dir::ASCIIString="") = joinpath(RESULTDIR, STUDYNAME, dir)

function combine_mc_logs(dir::AbstractString=studypath(MC_NAME))
    logjoin(dir, "acasx_mc_tree_log.txt", ["classifier_metrics",
        "interpretability_metrics", "computeinfo"], joinpath(dir, "subdirjoined"))
end
function combine_mcts_logs(dir::AbstractString=studypath(MCTS_NAME))
    logjoin(dir, "acasx_mcts_tree_log.txt", ["classifier_metrics",
        "interpretability_metrics", "computeinfo"], joinpath(dir, "subdirjoined"))
end
function combine_ge_logs(dir::AbstractString=studypath(GE_NAME))
    logjoin(dir, "acasx_ge_tree_log.txt", ["classifier_metrics",
        "interpretability_metrics", "computeinfo"], joinpath(dir, "subdirjoined"))
end
function combine_gp_logs(dir::AbstractString=studypath(GP_NAME))
    logjoin(dir, "acasx_gp_tree_log.txt", ["classifier_metrics",
        "interpretability_metrics", "computeinfo"], joinpath(dir, "subdirjoined"))
end

function master_log(namedirs::Vector{Pair{Symbol,ASCIIString}}=[:MC=>studypath(MC_NAME),
    :MCTS=>studypath(MCTS_NAME), :GE=>studypath(GE_NAME), :GP=>studypath(GP_NAME)])

    masterlog = nothing
    for (name, dir) in namedirs
        logs = load_log(TaggedDFLogger, joinpath(dir, "subdirjoined.txt"))

        #handle cpu time
        Dcpu = logs["computeinfo"]
        Dcpu = Dcpu[Dcpu[:parameter].=="cpu_time", [:value,:name]]
        rename!(Dcpu, :value, :cpu_time)
        Dcpu[:cpu_time] = map(x->parse(Float64,x), Dcpu[:cpu_time])

        D = join(logs["interpretability_metrics"], logs["classifier_metrics"]; on=[:name])
        D = join(D, Dcpu; on=[:name])
        D[:algorithm] = fill(name, nrow(D))
        if masterlog == nothing
            masterlog = D
        else
            append!(masterlog, D)
        end
    end
    writetable(MASTERLOG_FILE, masterlog)
    masterlog
end

master_plot(; kwargs...) = master_plot(readtable(MASTERLOG_FILE); kwargs...)

"""
Subsamples the collected data at 'subsample' rate to plot at a lower rate than collected
"""
function master_plot(masterlog::DataFrame; subsample::Int64=25000)
    D = masterlog

    delete!(D, :name)
    #aggregate over seed
    D = aggregate(D, [:algorithm], [mean, std, length, SEM])

    #workaround for naming in julia 0.4
    for s0 in names(D)
        str = string(s0)
        str = replace(str, "MathUtils.", "")
        s1 = Symbol(str)
        if s0 != s1
            rename!(D, s0, s1)
        end
    end

    writetable(PLOTLOG_FILE, D)
end

function combine_and_plot()
    combine_mc_logs()
    combine_mcts_logs()
    combine_ge_logs()
    combine_gp_logs()
    ml = master_log()
    master_plot(ml)
end

end #module


text:
import EnergyExpressions: QuantumOperator, NBodyTermFactor, NBodyTerm,
    OrbitalMatrixElement, NBodyMatrixElement
import AngularMomentumAlgebra: radial_integral

function ref_exchange_multipoles(a, b)
    # Table 1⁶ of Condon & Shortley. As noted by in ref [10] of
    #
    # - Cohl, H. S., Rau, A. R. P., Tohline, J. E., Browne, D. A., Cazes,
    #   J. E., & Barnes, E. I. (2001). Useful alternative to the multipole
    #   expansion of1/rpotentials. Physical Review A, 64(5),
    #   052509. http://dx.doi.org/10.1103/physreva.64.052509
    #
    # df, k = 5 is in error, and the numerator has to multiplied by 5
    # and the denominator replaced by 7623.
    reference = Dict((o"ks", o"kp") => ((1, 3, Dict((0,1) => 1,
                                                    (0,0) => 1)),),
                     (o"ks", o"kf") => ((3, 7, Dict((0,3) => 1,
                                                    (0,2) => 1,
                                                    (0,1) => 1,
                                                    (0,0) => 1)),),
                     (o"kp", o"kd") => ((1, 15, Dict((1,2) => 6,
                                                     (1,1) => 3,
                                                     (1,0) => 1,
                                                     (0,2) => 0,
                                                     (0,1) => 3,
                                                     (0,0) => 4,
                                                     (1,-2) => 0,
                                                     (1,-1) => 0)),
                                        (3, 245, Dict((1,2) => 3,
                                                      (1,1) => 9,
                                                      (1,0) => 18,
                                                      (0,2) => 15,
                                                      (0,1) => 24,
                                                      (0,0) => 27,
                                                      (1,-2) => 45,
                                                      (1,-1) => 30))),
                     (o"kd", o"kf") => ((1, 35, Dict((2,3) => 15,
                                                     (2,2) => 5,
                                                     (2,1) => 1,
                                                     (2,0) => 0,
                                                     (1,3) => 0,
                                                     (1,2) => 10,
                                                     (1,1) => 8,
                                                     (1,0) => 3,
                                                     (0,3) => 0,
                                                     (0,2) => 0,
                                                     (0,1) => 6,
                                                     (0,0) => 9,
                                                     (2,-3) => 0,
                                                     (2,-2) => 0,
                                                     (2,-1) => 0,
                                                     (1,-3) => 0,
                                                     (1,-2) => 0,
                                                     (1,-1) => 0)),
                                        (3, 315, Dict((2,3) => 10,
                                                      (2,2) => 20,
                                                      (2,1) => 24,
                                                      (2,0) => 20,
                                                      (1,3) => 25,
                                                      (1,2) => 15,
                                                      (1,1) => 2,
                                                      (1,0) => 2,
                                                      (0,3) => 25,
                                                      (0,2) => 0,
                                                      (0,1) => 9,
                                                      (0,0) => 16,
                                                      (2,-3) => 0,
                                                      (2,-2) => 0,
                                                      (2,-1) => 10,
                                                      (1,-3) => 0,
                                                      (1,-2) => 25,
                                                      (1,-1) => 15)),
                                        (5, 7623, Dict((2,3) => 5,
                                                       (2,2) => 25,
                                                       (2,1) => 75,
                                                       (2,0) => 175,
                                                       (1,3) => 35,
                                                       (1,2) => 120,
                                                       (1,1) => 250,
                                                       (1,0) => 400,
                                                       (0,3) => 140,
                                                       (0,2) => 315,
                                                       (0,1) => 450,
                                                       (0,0) => 500,
                                                       (2,-3) => 1050,
                                                       (2,-2) => 630,
                                                       (2,-1) => 350,
                                                       (1,-3) => 420,
                                                       (1,-2) => 560,
                                                       (1,-1) => 525))),
                     (o"ks", o"ks") => ((0, 1, Dict((0,0) => 1)),),
                     (o"ks", o"kd") => ((2, 5, Dict((0,2) => 1,
                                                    (0,1) => 1,
                                                    (0,0) => 1)),),
                     (o"kp", o"kp") => ((0, 1, Dict((1,1) => 1,
                                                    (1,0) => 0,
                                                    (0,0) => 1,
                                                    (1,-1) => 0)),
                                        (2, 25, Dict((1,1) => 1,
                                                     (1,0) => 3,
                                                     (0,0) => 4,
                                                     (1,-1) => 6))),
                     (o"kp", o"kf") => ((2, 175, Dict((1,3) => 45,
                                                      (1,2) => 30,
                                                      (1,1) => 18,
                                                      (1,0) => 9,
                                                      (0,3) => 0,
                                                      (0,2) => 15,
                                                      (0,1) => 24,
                                                      (0,0) => 27,
                                                      (1,-3) => 0,
                                                      (1,-2) => 0,
                                                      (1,-1) => 3)),
                                        (4, 189, Dict((1,3) => 1,
                                                      (1,2) => 3,
                                                      (1,1) => 6,
                                                      (1,0) => 10,
                                                      (0,3) => 7,
                                                      (0,2) => 12,
                                                      (0,1) => 15,
                                                      (0,0) => 16,
                                                      (1,-3) => 28,
                                                      (1,-2) => 21,
                                                      (1,-1) => 15))),
                     (o"kd", o"kd") => ((0, 1, Dict((2,2) => 1,
                                                    (2,1) => 0,
                                                    (2,0) => 0,
                                                    (1,1) => 1,
                                                    (1,0) => 0,
                                                    (0,0) => 1,
                                                    (2,-2) => 0,
                                                    (2,-1) => 0,
                                                    (1,-1) => 0)),
                                        (2, 49, Dict((2,2) => 4,
                                                     (2,1) => 6,
                                                     (2,0) => 4,
                                                     (1,1) => 1,
                                                     (1,0) => 1,
                                                     (0,0) => 4,
                                                     (2,-2) => 0,
                                                     (2,-1) => 0,
                                                     (1,-1) => 6)),
                                        (4, 441, Dict((2,2) => 1,
                                                      (2,1) => 5,
                                                      (2,0) => 15,
                                                      (1,1) => 16,
                                                      (1,0) => 30,
                                                      (0,0) => 36,
                                                      (2,-2) => 70,
                                                      (2,-1) => 35,
                                                      (1,-1) => 40))),
                     (o"kf", o"kf") => ((0, 1, Dict((3,3) => 1,
                                                    (3,2) => 0,
                                                    (3,1) => 0,
                                                    (3,0) => 0,
                                                    (2,2) => 1,
                                                    (2,1) => 0,
                                                    (2,0) => 0,
                                                    (1,1) => 1,
                                                    (1,0) => 0,
                                                    (0,0) => 1,
                                                    (3,-3) => 0,
                                                    (3,-2) => 0,
                                                    (3,-1) => 0,
                                                    (2,-2) => 0,
                                                    (2,-1) => 0,
                                                    (1,-1) => 0)),
                                        (2, 225, Dict((3,3) => 25,
                                                      (3,2) => 25,
                                                      (3,1) => 10,
                                                      (3,0) => 0,
                                                      (2,2) => 0,
                                                      (2,1) => 15,
                                                      (2,0) => 20,
                                                      (1,1) => 9,
                                                      (1,0) => 2,
                                                      (0,0) => 16,
                                                      (3,-3) => 0,
                                                      (3,-2) => 0,
                                                      (3,-1) => 0,
                                                      (2,-2) => 0,
                                                      (2,-1) => 0,
                                                      (1,-1) => 24)),
                                        (4, 1089, Dict((3,3) => 9,
                                                       (3,2) => 30,
                                                       (3,1) => 54,
                                                       (3,0) => 63,
                                                       (2,2) => 49,
                                                       (2,1) => 32,
                                                       (2,0) => 3,
                                                       (1,1) => 1,
                                                       (1,0) => 15,
                                                       (0,0) => 36,
                                                       (3,-3) => 0,
                                                       (3,-2) => 0,
                                                       (3,-1) => 42,
                                                       (2,-2) => 70,
                                                       (2,-1) => 14,
                                                       (1,-1) => 40)),
                                        (6, 7361*64, Dict((3,3) => 1,
                                                          (3,2) => 7,
                                                          (3,1) => 28,
                                                          (3,0) => 84,
                                                          (2,2) => 36,
                                                          (2,1) => 105,
                                                          (2,0) => 224,
                                                          (1,1) => 225,
                                                          (1,0) => 350,
                                                          (0,0) => 400,
                                                          (3,-3) => 924,
                                                          (3,-2) => 462,
                                                          (3,-1) => 210,
                                                          (2,-2) => 504,
                                                          (2,-1) => 378,
                                                          (1,-1) => 420))))

    g = CoulombInteraction()
    terms = NBodyTerm[]

    push_multipole!(k, coeff) =
        push!(terms, NBodyTerm([radial_integral([a,b], (k,g), [b,a])], coeff))

    ao,bo = Orbital(:k, a.orb.ℓ),Orbital(:k, b.orb.ℓ)
    i,j = ao ≤ bo ? (1,2) : (2,1)
    uv = (ao,bo)[[i,j]]

    m′,m = (a.m[1],b.m[1])[[i,j]]
    ao == bo && abs(m) > abs(m′) && ((m′,m) = (m,m′))
    # The table only stores the upper sign of the projection quantum
    # numbers, since the opposite signs have the same values.
    if iszero(m′)
        m = abs(m)
    else
        m *= sign(m′)
    end
    m′ = abs(m′)

    for (k,D,Ns) in reference[uv]
        N = Ns[(m′,m)]
        iszero(N) && continue
        push_multipole!(k, N/D)
    end

    NBodyMatrixElement(terms)
end

function ref_direct_multipoles(a, b)
    # Table 2⁶ of Condon & Shortley
    reference = Dict((o"kp", o"kp") => ((2, 25, Dict((1,1) => 1,
                                                     (1,0) => -2,
                                                     (0,0) => 4)),),
                     (o"kp", o"kd") => ((2, 35, Dict((1,2) => 2,
                                                     (1,1) => -1,
                                                     (1,0) => -2,
                                                     (0,2) => -4,
                                                     (0,1) => 2,
                                                     (0,0) => 4)),),
                     (o"kp", o"kf") => ((2, 75, Dict((1,3) => 5,
                                                     (1,2) => 0,
                                                     (1,1) => -3,
                                                     (1,0) => -4,
                                                     (0,3) => -10,
                                                     (0,2) => 0,
                                                     (0,1) => 6,
                                                     (0,0) => 8)),),
                     (o"kd", o"kd") => ((2, 49, Dict((2,2) => 4,
                                                     (2,1) => -2,
                                                     (2,0) => -4,
                                                     (1,1) => 1,
                                                     (1,0) => 2,
                                                     (0,0) => 4)),
                                        (4, 441, Dict((2,2) => 1,
                                                      (2,1) => -4,
                                                      (2,0) => 6,
                                                      (1,1) => 16,
                                                      (1,0) => -24,
                                                      (0,0) => 36))),
                     (o"kd", o"kf") => ((2, 105, Dict((2,3) => 10,
                                                      (2,2) => 0,
                                                      (2,1) => -6,
                                                      (2,0) => -8,
                                                      (1,3) => -5,
                                                      (1,2) => 0,
                                                      (1,1) => 3,
                                                      (1,0) => 4,
                                                      (0,3) => -10,
                                                      (0,2) => 0,
                                                      (0,1) => 6,
                                                      (0,0) => 8)),
                                        (4, 693, Dict((2,3) => 3,
                                                      (2,2) => -7,
                                                      (2,1) => 1,
                                                      (2,0) => 6,
                                                      (1,3) => -12,
                                                      (1,2) => 28,
                                                      (1,1) => -4,
                                                      (1,0) => -24,
                                                      (0,3) => 18,
                                                      (0,2) => -42,
                                                      (0,1) => 6,
                                                      (0,0) => 36))),
                     (o"kf", o"kf") => ((2, 225, Dict((3,3) => 25,
                                                      (3,2) => 0,
                                                      (3,1) => -15,
                                                      (3,0) => -20,
                                                      (2,2) => 0,
                                                      (2,1) => 0,
                                                      (2,0) => 0,
                                                      (1,1) => 9,
                                                      (1,0) => 12,
                                                      (0,0) => 16)),
                                        (4, 1089, Dict((3,3) => 9,
                                                       (3,2) => -21,
                                                       (3,1) => 3,
                                                       (3,0) => 18,
                                                       (2,2) => 49,
                                                       (2,1) => -7,
                                                       (2,0) => -42,
                                                       (1,1) => 1,
                                                       (1,0) => 6,
                                                       (0,0) => 36)),
                                        (6, 7361*64, Dict((3,3) => 1,
                                                          (3,2) => -6,
                                                          (3,1) => 15,
                                                          (3,0) => -20,
                                                          (2,2) => 36,
                                                          (2,1) => -90,
                                                          (2,0) => 120,
                                                          (1,1) => 225,
                                                          (1,0) => -300,
                                                          (0,0) => -400))))

    g = CoulombInteraction()
    terms = NBodyTerm[]

    push_multipole!(k, coeff) =
        push!(terms, NBodyTerm([radial_integral([a,b], (k,g), [a,b])], coeff))

    # The monopole is always present in the multipole expansion of the
    # direct interaction.
    push_multipole!(0, 1.0)

    ao,bo = Orbital(:k, a.orb.ℓ),Orbital(:k, b.orb.ℓ)
    i,j = ao ≤ bo ? (1,2) : (2,1)
    uv = (ao,bo)[[i,j]]
    if uv ∈ keys(reference)
        ref = reference[uv]

        m′m = (abs(a.m[1]), abs(b.m[1]))
        m′,m = uv[1] < uv[2] ? m′m[[i,j]] : reverse(minmax(m′m...))

        for (k,D,Ns) in ref
            N = Ns[(m′,m)]
            iszero(N) && continue
            push_multipole!(k, N/D)
        end
    end

    NBodyMatrixElement(terms)
end


text:
using SFML

const ipaddress = "127.0.0.1"
const port = 8080

socket = SocketTCP()

if connect(socket, ipaddress, port, 0) != SOCKET_DONE
    println("Connection to $ipaddress on port $port failed!")
    exit(1)
end

println("Connected")

packet = Packet()

while true
    packet = Packet()

    print("Message: ")
    msg = strip(readline(STDIN), '\n')

    write(packet, msg)
    send(socket, packet)

    in_packet = Packet()
    receive(socket, in_packet)
    println("Received: ", read_string(in_packet))
end


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

"""
    project_from_surface_to_plane(p, x0, n0)

Find the projection of vertex `p` to the auxiliary plane (x0, n0), where
`x0` is the origo of plane and `n0` is unit normal of plane.

# Example

```jldoctest
p = [1.0, 1.0, 3.0]
x0 = [0.0, 0.0, 0.0]
n0 = [0.0, 0.0, 1.0]
project_from_surface_to_plane(p, x0, n0)

# output
3-element Array{Float64,1}:
 1.0
 1.0
 0.0

```
"""
function project_from_surface_to_plane(p, x0, n0)
    return p - dot(p-x0, n0)*n0
end

function inv3(P::Matrix)
    n, m = size(P)
    @assert n == m == 3
    a, b, c, d, e_, f, g, h, i = P
    A = e_*i - f*h
    B = -d*i + f*g
    C = d*h - e_*g
    D = -b*i + c*h
    E = a*i - c*g
    F = -a*h + b*g
    G = b*f - c*e_
    H = -a*f + c*d
    I_ = a*e_ - b*d
    return 1/(a*A + b*B + c*C)*[A B C; D E F; G H I_]
end

"""
    project_from_plane_to_surface(Val{:Tri3}, p, n, x1, x2, x3)

Project vertex `p` to surface of Tri3 element.

# Example
```jldoctest
p = [1/3, 1/3, 0.0]
n = [0.0, 0.0, 1.0]
x1 = [0.0, 0.0, 1.0]
x2 = [1.0, 0.0, 1.0]
x3 = [0.0, 1.0, 1.0]
theta, d = project_from_plane_to_surface(Val{:Tri3}, p, n, x1, x2, x3)
round.(theta, 3), d

# output

([0.333, 0.333], 1.0)
```
"""
function project_from_plane_to_surface(::Type{Val{:Tri3}}, p, n, x1, x2, x3)
    A = [x1-x2 x1-x3 n]
    b = x1-p
    x = inv3(A)*b
    return x[1:2], x[3]
end


text:
function opstring(op::Integer)
    strings = [
        ">",
        "<",
        "+",
        "-",
        ".",
        ",",
        "[",
        "]",
    ]
    return strings[op]
end

function printbeforeop(
    op::Integer,
    memory::Vector{Uint8},
    memory_pointer::Integer,
)
    @printf(
        STDERR,
        "\n---\nBefore executing '%s': (%s, %d)\n",
        opstring(op),
        repr(memory),
        memory_pointer
    )
    return
end

function printafterop(
    op::Integer,
    memory::Vector{Uint8},
    memory_pointer::Integer,
)
    @printf(
        STDERR,
        "After executing '%s': (%s, %d)\n---\n\n",
        opstring(op),
        repr(memory),
        memory_pointer
    )
    return
end

function interpret(
    ops::Vector{Int};
    debug::Bool = false,
    io_in::IO = STDIN,
    io_out::IO = STDOUT
)
    memory_pointer = 1
    memory = zeros(Uint8, 1)
    opindex = 1

    n_ops = length(ops)

    while opindex <= n_ops
        op = ops[opindex]

        if debug
            printbeforeop(op, memory, memory_pointer)
        end

        # OP1: >
        #
        # Move the Memory Pointer to the next array cell.
        #
        if op == OP1
            memory_pointer += 1
            if memory_pointer > length(memory)
                push!(memory, 0x00)
            end
            opindex += 1

        # OP2: <
        #
        # Move the Memory Pointer to the previous array cell.
        #
        elseif op == OP2
            # TODO: Confirm that this is appropriate
            if memory_pointer > 1
                memory_pointer -= 1
            end
            opindex += 1

        # OP3: +
        #
        # Increment the array cell pointed at by the Memory Pointer.
        #
        elseif op == OP3
            memory[memory_pointer] += 1
            opindex += 1

        # OP4: -
        #
        # Decrement the array cell pointed at by the Memory Pointer.
        #
        elseif op == OP4
            memory[memory_pointer] -= 1
            opindex += 1

        # OP5: ,
        #
        # Read a character from STDIN and put its ASCII value into the cell
        # pointed at by the Memory Pointer.
        #
        elseif op == OP5
            memory[memory_pointer] = read(io_in, Uint8)
            opindex += 1

        # OP6: .
        #
        # Print the character with ASCII value equal to the value in the cell
        # pointed at by the Memory Pointer.
        #
        elseif op == OP6
            if debug
                @printf(STDERR, "Printing '%s'\n", char(memory[memory_pointer]))
            end
            print(io_out, char(memory[memory_pointer]))
            opindex += 1

        # OP7: [
        #
        # Move to the command following the matching Ook? Ook! if the value in
        # the cell pointed at by the Memory Pointer is zero. Note that
        # Ook! Ook? and Ook? Ook! commands nest like pairs of parentheses, and
        # matching pairs are defined in the same way as for parentheses.
        #
        elseif op == OP7
            if memory[memory_pointer] == 0x00
                level = 0
                opindex += 1
                complete = false
                while !complete && opindex <= n_ops
                    if ops[opindex] == OP7
                        level += 1
                        opindex += 1
                    elseif ops[opindex] == OP8
                        if level == 0
                            opindex += 1
                            complete = true
                        else
                            level -= 1
                            opindex += 1
                        end
                    else
                        opindex += 1
                    end
                end
            else
                opindex += 1
            end

        # OP8: ]
        #
        # Move to the command following the matching Ook! Ook? if the value in
        # the cell pointed at by the Memory Pointer is non-zero. This involves
        # traversing the ops array backwards.
        #
        elseif op == OP8
            if memory[memory_pointer] != 0x00
                level = 0
                opindex -= 1
                complete = false
                while !complete && opindex >= 1
                    if ops[opindex] == OP8
                        level += 1
                        opindex -= 1
                    elseif ops[opindex] == OP7
                        if level == 0
                            opindex += 1
                            complete = true
                        else
                            level -= 1
                            opindex -= 1
                        end
                    else
                        opindex -= 1
                    end
                end
            else
                opindex += 1
            end

        # Invalid operation
        else
            msg = @sprintf("Unknown op code: %d", ops[opindex])
            throw(ArgumentError(msg))
        end

        if debug
            printafterop(op, memory, memory_pointer)
        end
    end

    return
end

function brainfuck(
    s::String;
    debug::Bool = false,
    io_in::IO = STDIN,
    io_out::IO = STDOUT
)
    interpret(lexandparse(s), debug = debug, io_in = io_in, io_out = io_out)
end