text:
#######################################
###### Monomials and polynomials ######
#######################################

#  An SoSMonom represents a monomial: it maps each symbol to its exponent.
#  Exponents should always be strictly positive (never 0).
typealias SoSMonom Dict{Symbol,Int64}

#  An SoSPoly maps each monomial to its coefficient.
typealias SoSPoly Dict{SoSMonom,Float64}

#    there should be a faster option here that manages iterators, a la mergesort
macro dictplus(a, b)
    quote
        # A, including overlap with B
        map = Dict([ k => v + get(b,k,0) for (k,v) in a ])
        # B \ A
        for (k,v) in b
            if !haskey(a,k)
                map[k] = v
            end
        end
        map
    end
end

#  Multiplication of monomials.
*(a::SoSMonom,b::SoSMonom) = @dictplus(a,b)
+(a::SoSPoly,b::SoSPoly) = @dictplus(a,b)

#  Addition in-place of polynomials
function addpoly!(a::SoSPoly, b::SoSPoly)
    for (k,v) in b
        a[k] = get(a,k,0.0) + v
    end
end

#  Multiplication of polynomials -- just the naive method
function *(a::SoSPoly, b::SoSPoly)
    map = SoSPoly()
    for (k1,v1) in a
        for (k2,v2) in b
            kp = k1 * k2
            map[kp] = get(map,kp,0) + v1 * v2
        end
    end
    return map
end

#  Degree of monomial / polynomial
deg(m :: SoSMonom) = sum(values(m))
deg(p :: SoSPoly) = (length(p) == 0) ? 0 : maximum([deg(k) for (k,v) in p])

#  Enumerate all monomials of degree d.
#    Julia knows about combinations (without replacement), and we
#    want to choose variables with replacement. We simulate this
#    via (v+d-1 choose d).
monoms(varr :: Array{Symbol}, d :: Int64) =
    [ combToMonom(varr,co) for co in combinations(1:(length(varr)+d-1),d) ] :: Array{SoSMonom,1}
monoms(d :: Int64) = (d == 0) ? [SoSMonom()] : Array{SoSMonom,1}() # edge case

function combToMonom(varr :: Array{Symbol}, co :: Array{Int64})
    m = Dict{Symbol,Int64}() :: SoSMonom
    shift = 0
    for i in co
        k = varr[i+shift]
        m[k] = 1 + get(m,k,0);
        shift -= 1;
    end
    m
end

# Find all decompositions of mon into two monomials (m1,m2), each of degree
# at most degmax.
# XXX currently unused
function decomp(mon :: SoSMonom, degmax :: Int64)
    # maybe there's room for improvement here ...
    list = Tuple{SoSMonom,SoSMonom}[]
    varr = [k for (k,v) in mon]

    # for each candidate m1, see if its complement m2 is sensible
    mondeg = deg(mon)
    dmin = max(0,mondeg-degmax)
    dmax = min(mondeg,degmax)
    for deg in dmin:dmax
        for m1 in monoms(varr,deg)

            # construct the complement
            m2 = Dict([ k => (v-get(m1,k,0)) for (k,v) in mon ])

            # test for being sensible (ie non-negative exponents), and remove any zero entries
            if(length(m2) > 0 && minimum(values(m2)) < 0)
                continue
            end

            # remove any zero entries
            filter!((k,v)->(v!=0), m2)

            push!(list, (m1,m2))
        end
    end

    list
end

max_stars_repo_name: willperry/SumOfSquaresOptimization


text:
using Test
using SIMD
using ChipSort

@testset "Data buffers" begin

# types_to_test = [Int8, Int16, Int32, Int64, Float32, Float64]
types_to_test = [Int8, Float64]

function test_data_buffer(T, ::Val{N}) where N
    data = sort!(rand(T, N*N))
    db = DataBuffer(Val(N), data)
    @test all(all.([pop!(db) for k in 1:N] .== [Vec(tuple(data[k:k+N-1]...)) for k in 1:N:N*N]))
    @test pop!(db) == nothing
end

@testset for T in types_to_test
    for N in [4, 8]
        test_data_buffer(T, Val(N))
    end
end

function test_merge_node(T, ::Val{N}) where N
    data_a = sort!(rand(T, N*N))
    data_b = sort!(rand(T, N*N))
    alldata = sort!(vcat(data_a, data_b))

    left = DataBuffer(Val(N), data_a)
    right = DataBuffer(Val(N), data_b)
    mrg = MergeNode(left, right)

    @test all(all.([pop!(mrg) for k in 1:2*N] .== [Vec(tuple(alldata[k:k+N-1]...)) for k in 1:N:2*N*N]))
    @test pop!(mrg) == nothing
end

@testset for T in types_to_test
    for N in [4, 8]
        test_merge_node(T, Val(N))
    end
end

function test_multi_way_merge(T, ::Val{N}, ::Val{M}) where {N,M}
    data = [sort!(rand(T, N*N)) for _ in 1:M]
    alldata = sort!(vcat(data...))

    mrg = build_multi_merger(Val(N), data...)

    srt_hat = [pop!(mrg) for k in 1:M*N]
    srt_ref = [Vec(tuple(alldata[k:k+N-1]...)) for k in 1:N:M*N*N]

    @test all(all.(srt_hat .== srt_ref))
    @test pop!(mrg) == nothing
end

@testset for T in types_to_test
    for N in [4, 32]
        for M in [4, 32]
            test_multi_way_merge(T, Val(N), Val(M))
        end
    end
end

end

ext: jl


text:
using Test
using Mimi

@testset "Standard API" begin

# Test that the function does not error and returns a valid value
scc1 = MimiPAGE2009.compute_scc(year=2020)
@test scc1 isa Float64

# Test that a higher discount rate makes a lower scc value
scc2 = MimiPAGE2009.compute_scc(year=2020, eta=0., prtp=0.03)
@test scc2 < scc1

# Test with a modified model
m = MimiPAGE2009.get_model()
update_param!(m, :tcr_transientresponse, 3)
scc3 = MimiPAGE2009.compute_scc(m, year=2020)
@test scc3 > scc1

# Test get_marginal_model
mm = MimiPAGE2009.get_marginal_model(year = 2040)
mm[:ClimateTemperature, :rt_realizedtemperature]

# Test compute_scc_mm
result = MimiPAGE2009.compute_scc_mm(year=2050)
@test result.scc < scc1
@test result.mm isa Mimi.MarginalModel

end


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using SparseArrays, Base.Threads, LinearAlgebra, JuMP

"""
    _injection_shift_factors(; buses, lines)

Returns a (B-1)xL matrix M, where B is the number of buses and L is the number
of transmission lines. For a given bus b and transmission line l, the entry
M[l.offset, b.offset] indicates the amount of power (in MW) that flows through
transmission line l when 1 MW of power is injected at the slack bus (the bus
that has offset zero) and withdrawn from b.
"""
function _injection_shift_factors(;
    buses::Array{Bus},
    lines::Array{TransmissionLine},
)
    susceptance = _susceptance_matrix(lines)
    incidence = _reduced_incidence_matrix(lines = lines, buses = buses)
    laplacian = transpose(incidence) * susceptance * incidence
    isf = susceptance * incidence * inv(Array(laplacian))
    return isf
end

"""
    _reduced_incidence_matrix(; buses::Array{Bus}, lines::Array{TransmissionLine})

Returns the incidence matrix for the network, with the column corresponding to
the slack bus is removed. More precisely, returns a (B-1) x L matrix, where B
is the number of buses and L is the number of lines. For each row, there is a 1
element and a -1 element, indicating the source and target buses, respectively,
for that line.
"""
function _reduced_incidence_matrix(;
    buses::Array{Bus},
    lines::Array{TransmissionLine},
)
    matrix = spzeros(Float64, length(lines), length(buses) - 1)
    for line in lines
        if line.source.offset > 0
            matrix[line.offset, line.source.offset] = 1
        end
        if line.target.offset > 0
            matrix[line.offset, line.target.offset] = -1
        end
    end
    return matrix
end

"""
    _susceptance_matrix(lines::Array{TransmissionLine})

Returns a LxL diagonal matrix, where each diagonal entry is the susceptance of
the corresponding transmission line.
"""
function _susceptance_matrix(lines::Array{TransmissionLine})
    return Diagonal([l.susceptance for l in lines])
end

"""

    _line_outage_factors(; buses, lines, isf)

Returns a LxL matrix containing the Line Outage Distribution Factors (LODFs)
for the given network. This matrix how does the pre-contingency flow change
when each individual transmission line is removed.
"""
function _line_outage_factors(;
    buses::Array{Bus,1},
    lines::Array{TransmissionLine,1},
    isf::Array{Float64,2},
)::Array{Float64,2}
    incidence = Array(_reduced_incidence_matrix(lines = lines, buses = buses))
    lodf::Array{Float64,2} = isf * transpose(incidence)
    _, n = size(lodf)
    for i in 1:n
        lodf[:, i] *= 1.0 / (1.0 - lodf[i, i])
        lodf[i, i] = -1
    end
    return lodf
end


text:
function band(a,m,n,p)
r = (m+1)/2

for i=1:n
  for j=m+1:m+r-1
    a[i,j]=0.0
  end
end

for k=1:n
  max = 0.0
  i = k
  j = r
  while (i <= n) && (j >= 1)
    d = abs[i,j]
    if max < d
      max = d
      p[k] = i
    end
    i += 1
    j -= 1
  end
end

if max <= eps()
  m=0
  return
end

if p[k] != k
  i = r
  j = r+k-p[k]
  while (i <= m+r-1) && (i <= n-k+r)
    c=a[k,i]
    a[k,i]=a[p[k],j]
    a[p[k],j]=c
    i +=1
    j +=1
  end
end

a[k,r]=1.0/a[k,r]
h = r-1
i = k+1
while (h >= 1) && (i <= n)
end

return
end


text:
import BitIntegers: UInt256, @uint256_str
import SHA: sha256

Ast = Dict{String, Any}

abstract type AbstractTemplate end

Scope = Dict{String, AbstractTemplate}

struct Template <: AbstractTemplate
  scopes::Vector{Scope}
end

struct Selector
  id::Int64
  name::String
end

struct Signal
end

struct Constraint
  # = lc.QEQ|ERROR
end

struct Component
  signals::Vector{String}
  template::String
  params::Dict{String, UInt256}
end

struct Function
  # not used in exec
end

struct FunctionParams
  # copied directly from untyped ast
end

struct CtxErrorPos
  firstLine::Int64
  firstColumn::Int64
  lastLine::Int64
  lastColumn::Int64
end

struct CtxError
  pos::CtxErrorPos
  errStr::String
  errFile::String
  ast::Ast
  message::String
end

mutable struct Ctx
  scopes::Vector{Scope}
  signals::Dict{String, Signal}
  constraints::Vector{Constraint}
  components::Dict{String, Component}
  currentComponent::String
  functions::Dict{String, Function}
  functionParams::Dict{String, Vector{FunctionParams}}
  templates::Dict{String, Template}
  returnValue::UInt256
  fileName::String
  filePath::String
  includedFiles::Set{String}
  error:Union{Nothing, CtxError}
  Ctx() = new(Vector{Scope}(), Vector{Signal}(), Vector{Constraint}(), "", nothing)
end

## In executor
# ctx.constraints
# ctx.scopes
# ctx.signals
# ctx.error
# ctx.fileName
# ctx.filePath
# ctx.returnValue
# ctx.includedFiles
# ctx.components
# ctx.currentComponent
# ctx.functions
# ctx.functionParams

# ctx.templates

## In compiler/codegen/witness
# ctx.callFunction
# ctx.assert
# ctx.setPin
# ctx.setSignal
# ctx.setVar
# ctx.getPin
# ctx.getSignal
# ctx.getVar
# ctx.signalName2Idx
# ctx.signalNames
# ctx.totals

function exec(ctx::Ctx, ast::Ast)
  if (ast["type"] == "NUMBER") || (ast["type"] == "LINEARCOMBINATION") || (ast["type"] =="SIGNAL") || (ast["type"] == "QEQ")
      return ast
  elseif ast["type"] == "VARIABLE"
      return execVariable(ctx, ast)
  elseif ast["type"] == "PIN"
      return execPin(ctx, ast)
  elseif ast["type"] == "OP"
      if ast.op == "="
          return execVarAssignement(ctx, ast)
      elseif ast.op == "<--"
          return execSignalAssign(ctx, ast)
      elseif ast.op == "<=="
          return execSignalAssignConstrain(ctx, ast)
      elseif ast.op == "==="
          return execConstrain(ctx, ast)
      elseif ast.op == "+="
          return execVarAddAssignement(ctx, ast)
      elseif ast.op == "*="
          return execVarMulAssignement(ctx, ast)
      elseif ast.op == "+"
          return execAdd(ctx, ast)
      elseif ast.op == "-"
          return execSub(ctx, ast)
      elseif ast.op == "UMINUS"
          return execUMinus(ctx, ast)
      elseif ast.op == "*"
          return execMul(ctx, ast)
      elseif ast.op == "%"
          return execMod(ctx, ast)
      elseif ast.op == "PLUSPLUSRIGHT"
          return execPlusPlusRight(ctx, ast)
      elseif ast.op == "PLUSPLUSLEFT"
          return execPlusPlusLeft(ctx, ast)
      elseif ast.op == "MINUSMINUSRIGHT"
          return execMinusMinusRight(ctx, ast)
      elseif ast.op == "MINUSMINUSLEFT"
          return execMinusMinusLeft(ctx, ast)
      elseif ast.op == "/"
          return execDiv(ctx, ast)
      elseif ast.op == "\\"
          return execIDiv(ctx, ast)
      elseif ast.op == "**"
          return execExp(ctx, ast)
      elseif ast.op == "&"
          return execBAnd(ctx, ast)
      elseif ast.op == "&&"
          return execAnd(ctx, ast)
      elseif ast.op == "||"
          return execOr(ctx, ast)
      elseif ast.op == "<<"
          return execShl(ctx, ast)
      elseif ast.op == ">>"
          return execShr(ctx, ast)
      elseif ast.op == "<"
          return execLt(ctx, ast)
      elseif ast.op == ">"
          return execGt(ctx, ast)
      elseif ast.op == "<="
          return execLte(ctx, ast)
      elseif ast.op == ">="
          return execGte(ctx, ast)
      elseif ast.op == "=="
          return execEq(ctx, ast)
      elseif ast.op == "!="
          return execNeq(ctx, ast)
      elseif ast.op == "?"
          return execTerCon(ctx, ast)
      else
          error(ctx, ast, "Invalid operation: " + ast.op)
      end
  elseif ast["type"] == "DECLARE"
      if ast["declareType"] == "COMPONENT"
          return execDeclareComponent(ctx, ast)
      elseif (ast["declareType"] == "SIGNALIN")||(ast["declareType"] == "SIGNALOUT")||(ast["declareType"] == "SIGNAL")
          return execDeclareSignal(ctx, ast)
      elseif ast["declareType"] == "VARIABLE"
          return execDeclareVariable(ctx, ast)
      else
          error(ctx, ast, "Invalid declaration: " + ast["declareType"])
      end
  elseif ast["type"] == "FUNCTIONCALL"
      return execFunctionCall(ctx, ast)
  elseif ast["type"] == "BLOCK"
      return execBlock(ctx, ast)
  elseif ast["type"] == "FOR"
      return execFor(ctx, ast)
  elseif ast["type"] == "WHILE"
      return execWhile(ctx, ast)
  elseif ast["type"] == "IF"
      return execIf(ctx, ast)
  elseif ast["type"] == "RETURN"
      return execReturn(ctx, ast)
  elseif ast["type"] == "TEMPLATEDEF"
      return execTemplateDef(ctx, ast)
  elseif ast["type"] == "FUNCTIONDEF"
      return execFunctionDef(ctx, ast)
  elseif ast["type"] == "INCLUDE"
      return execInclude(ctx, ast)
  elseif ast["type"] == "ARRAY"
      return execArray(ctx, ast)
  else
      error(ctx, ast, "Invalid AST node type: " + ast["type"])
  end
end

function error(ctx::Ctx, ast::Ast, errStr::String)
  ctx.error = CtxError(
      CtxErrorPos(ast["first_line"], ast["first_column"], ast["last_line"], ast["last_column"]),
      errStr,
      ctx.fileName,
      ast: ast,
      errStr
  )
end

# function iterateSelectors(ctx, sizes, baseName, fn) {
#   if (sizes.length == 0) {
#       return fn(baseName)
#   }
#   const res = []
#   for (let i=0; i<sizes[0]; i++) {
#       res.push(iterateSelectors(ctx, sizes.slice(1), baseName+"["+i+"]", fn))
#       if (ctx.error) return null
#   }
#   return res
# }

# TODO - replace basename string to hashes
function iterateSelectors(ctx::Ctx, sizes::P, baseName::String, fn::T) where {T<: Function, P<: AbstractVector{Int64}}
  if length(sizes) == 0
      return fn(baseName)
  end
  res = []
  for i=1:sizes[1]
    push!(res, iterateSelectors(ctx, view(sizes, 2, length(sizes)), baseName+"["+i+"]", fn))
    if (ctx.error)
      return nothing
    end
  end
  return res
end

# function setScope(ctx, name, selectors, value) {
#   let l = getScopeLevel(ctx, name)
#   if (l==-1) l= ctx.scopes.length-1

#   if (selectors.length == 0) {
#       ctx.scopes[l][name] = value
#   } else {
#       setScopeArray(ctx.scopes[l][name], selectors)
#   }

#   function setScopeArray(a, sels) {
#       if (sels.length == 1) {
#           a[sels[0].value] = value
#       } else {
#           setScopeArray(a[sels[0]], sels.slice(1))
#       }
#   }
# }

# selectors struct changed
function setScope(ctx::Ctx, name::String, selectors::Vector{Selector}, value::Template)
  l = getScopeLevel(ctx, name)
  if l === nothing
    l = length(ctx.scopes)
  end

  if length(selectors) == 0
      ctx.scopes[l][name] = value
  else
    cursor = ctx.scopes[l][name]
    for s in view(selectors, 1, length(selectors) - 1)
      cursor = cursor.scopes[s.id][s.name]
    end
    s = selectors[end]
    cursor.scopes[s.id][s.name] = value
  end
end

# function getScope(ctx, name, selectors) {

#   const sels = []
#   if (selectors) {
#       for (let i=0; i< selectors.length; i++) {
#           const idx = exec(ctx, selectors[i])
#           if (ctx.error) return

#           if (idx.type != "NUMBER") return error(ctx, selectors[i], "expected a number")
#           sels.push( idx.value.toJSNumber() )
#       }
#   }

#   function select(v, s) {
#       s = s || []
#       if (s.length == 0)  return v
#       return select(v[s[0]], s.slice(1))
#   }

#   for (let i=ctx.scopes.length-1; i>=0; i--) {
#       if (ctx.scopes[i][name]) return select(ctx.scopes[i][name], sels)
#   }
#   return null
# }

# function getScope(ctx::Ctx, name::String, selectors::Vector{Selector})
#   const sels = []
#   if (selectors) {
#       for (let i=0; i< selectors.length; i++) {
#           const idx = exec(ctx, selectors[i])
#           if (ctx.error) return

#           if (idx.type != "NUMBER") return error(ctx, selectors[i], "expected a number")
#           sels.push( idx.value.toJSNumber() )
#       }
#   }

#   function select(v, s) {
#       s = s || []
#       if (s.length == 0)  return v
#       return select(v[s[0]], s.slice(1))
#   }

#   for (let i=ctx.scopes.length-1; i>=0; i--) {
#       if (ctx.scopes[i][name]) return select(ctx.scopes[i][name], sels)
#   }
#   return null
# }

function getScopeLevel(ctx::Ctx, name::String)
  for i=length(ctx.scopes):-1:1
    if haskey(ctx.scopes[i], name)
      return i
    end
  end
  return nothing
end

function execBlock(ctx::Ctx, ast::Ast)
    for i = 1:length(ast.statements)
        exec(ctx, ast.statements[i])
        if (ctx.returnValue || ctx.error)
            return
        end
    end
end

# function execTemplateDef(ctx::Ctx, ast::Ast)
#     const scope = ctx.scopes[0];  // Lets put templates always in top scope.
#     //    const scope = ctx.scopes[ctx.scopes.length-1]
#     if (getScope(ctx, ast.name)) {
#         return error(ctx, ast, "Name already exists: "+ast.name)
#     }
#     scope[ast.name] = {
#         type: "TEMPLATE",
#         params: ast.params,
#         block: ast.block,
#         fileName: ctx.fileName,
#         filePath: ctx.filePath,
#         scopes: copyScope(ctx.scopes)
#     }
# end
#
# function execFunctionDef(ctx::Ctx, ast::Ast)
#     const scope = ctx.scopes[0]; // Lets put functions always in top scope.
#     //    const scope = ctx.scopes[ctx.scopes.length-1]
#     if (getScope(ctx, ast.name)) {
#         return error(ctx, ast, "Name already exists: "+ast.name)
#     }
#     ctx.functionParams[ast.name] = ast.params
#     scope[ast.name] = {
#         type: "FUNCTION",
#         params: ast.params,
#         block: ast.block,
#         fileName: ctx.fileName,
#         filePath: ctx.filePath,
#         scopes: copyScope(ctx.scopes)
#     }
# end
#
# function execDeclareComponent(ctx::Ctx, ast::Ast)
#     const scope = ctx.scopes[ctx.scopes.length-1]
#
#     if (ast.name.type != "VARIABLE") return error(ctx, ast, "Invalid component name")
#     if (getScope(ctx, ast.name.name)) return error(ctx, ast, "Name already exists: "+ast.name.name)
#
#     const baseName = ctx.currentComponent ? ctx.currentComponent + "." + ast.name.name : ast.name.name
#
#     const sizes=[]
#     for (let i=0; i< ast.name.selectors.length; i++) {
#         const size = exec(ctx, ast.name.selectors[i])
#         if (ctx.error) return
#
#         if (size.type != "NUMBER") return error(ctx, ast.name.selectors[i], "expected a number")
#
#         sizes.push( size.value.toJSNumber() )
#     }
#
#
#     scope[ast.name.name] = iterateSelectors(ctx, sizes, baseName, function(fullName) {
#
#         ctx.components[fullName] = "UNINSTANTIATED"
#
#         return {
#             type: "COMPONENT",
#             fullName: fullName
#         }
#     })
#
#     return {
#         type: "VARIABLE",
#         name: ast.name.name,
#         selectors: []
#     }
# }
#
# function execInstantiateComponet(ctx, vr, fn) {
#
#     if (vr.type != "VARIABLE") return error(ctx, fn, "Left hand instatiate component must be a variable")
#     if (fn.type != "FUNCTIONCALL") return error(ctx, fn, "Right type of instantiate component must be a function call")
#
#     const componentName = vr.name
#     const templateName = fn.name
#
#     const scopeLevel = getScopeLevel(ctx, templateName)
#     if (scopeLevel == -1) return error(ctx,fn, "Invalid Template")
#     const template = getScope(ctx, templateName)
#
#     if (template.type != "TEMPLATE") return error(ctx, fn, "Invalid Template")
#
#
#     const paramValues = []
#     for (let i=0; i< fn.params.length; i++) {
#         const v = exec(ctx, fn.params[i])
#         if (ctx.error) return
#
#         paramValues.push(v)
#     }
#     if (template.params.length != paramValues.length) error(ctx, fn, "Invalid Number of parameters")
#
#     const vv = getScope(ctx, componentName, vr.selectors)
#
#     if (!vv) return error(ctx, vr, "Component not defined"+ componentName)
#
#     instantiateComponent(vv)
#
#     function instantiateComponent(varVal) {
#
#         function extractValue(v) {
#             if (Array.isArray(v)) {
#                 return v.map(extractValue)
#             } else {
#                 return v.value.toString()
#             }
#         }
#
#         if (Array.isArray(varVal)) {
#             for (let i =0; i<varVal.length; i++) {
#                 instantiateComponent(varVal[i])
#             }
#             return
#         }
#
#         if (ctx.components[varVal.fullName] != "UNINSTANTIATED") error(ctx, fn, "Component already instantiated")
#
#         const oldComponent = ctx.currentComponent
#         const oldFileName = ctx.fileName
#         const oldFilePath = ctx.filePath
#         ctx.currentComponent = varVal.fullName
#
#         ctx.components[ctx.currentComponent] = {
#             signals: [],
#             params: {}
#         }
#
#         const oldScopes = ctx.scopes
#
#         ctx.scopes = oldScopes.slice(0, scopeLevel+1)
#
#         if (template.params.length != paramValues.length) return error(ctx, fn, "Invalid number of parameters: " + templateName)
#
#         const scope = {}
#         for (let i=0; i< template.params.length; i++) {
#             scope[template.params[i]] = paramValues[i]
#             ctx.components[ctx.currentComponent].params[template.params[i]] = extractValue(paramValues[i])
#         }
#
#         ctx.components[ctx.currentComponent].template = templateName
#         ctx.fileName = template.fileName
#         ctx.filePath = template.filePath
#         ctx.scopes = copyScope( template.scopes )
#         ctx.scopes.push(scope)
#
#         execBlock(ctx, template.block)
#
#         ctx.fileName = oldFileName
#         ctx.filePath = oldFilePath
#         ctx.currentComponent = oldComponent
#         ctx.scopes = oldScopes
#     }
# end
#
# function execFunctionCall(ctx::Ctx, ast::Ast)
#
#     const scopeLevel = getScopeLevel(ctx, ast.name)
#     if (scopeLevel == -1) return error(ctx, ast, "Function not defined: " + ast.name)
#     const fnc = getScope(ctx, ast.name)
#
#     if (fnc.type != "FUNCTION") return error(ctx, ast, "Not a function: " + ast.name)
#
#     const paramValues = []
#     for (let i=0; i< ast.params.length; i++) {
#         const v = exec(ctx, ast.params[i])
#         if (ctx.error) return
#
#         paramValues.push(v)
#     }
#
#     if (ast.params.length != paramValues.length) error(ctx, ast, "Invalid Number of parameters")
#
#     const oldFileName = ctx.fileName
#     const oldFilePath = ctx.filePath
#
#     const oldScopes = ctx.scopes
#
#     ctx.scopes = oldScopes.slice(0, scopeLevel+1)
#
#     const scope = {}
#     for (let i=0; i< fnc.params.length; i++) {
#         scope[fnc.params[i]] = paramValues[i]
#     }
#
#     ctx.fileName = fnc.fileName
#     ctx.filePath = fnc.filePath
#     ctx.scopes = copyScope( fnc.scopes )
#     ctx.scopes.push(scope)
#
#     execBlock(ctx, fnc.block)
#
#     const res = ctx.returnValue
#     ctx.returnValue = null
#
#     ctx.fileName = oldFileName
#     ctx.filePath = oldFilePath
#     ctx.scopes = oldScopes
#
#     return res
# end

function execReturn(ctx::Ctx, ast::Ast)
    ctx.returnValue = exec(ctx, ast.value)
end

# function execDeclareSignal(ctx::Ctx, ast::Ast)
#     const scope = ctx.scopes[ctx.scopes.length-1]
#
#     if (ast.name.type != "VARIABLE") return error(ctx, ast, "Invalid component name")
#     if (getScope(ctx, ast.name.name)) return error(ctx, ast, "Name already exists: "+ast.name.name)
#
#     const baseName = ctx.currentComponent ? ctx.currentComponent + "." + ast.name.name : ast.name.name
#
#     const sizes=[]
#     for (let i=0; i< ast.name.selectors.length; i++) {
#         const size = exec(ctx, ast.name.selectors[i])
#         if (ctx.error) return
#
#         if (size.type != "NUMBER") return error(ctx, ast.name.selectors[i], "expected a number")
#         sizes.push( size.value.toJSNumber() )
#     }
#
#     scope[ast.name.name] = iterateSelectors(ctx, sizes, baseName, function(fullName) {
#         ctx.signals[fullName] = {
#             fullName: fullName,
#             direction: ast.declareType == "SIGNALIN" ? "IN" : (ast.declareType == "SIGNALOUT" ? "OUT" : ""),
#             private: ast.private,
#             component: ctx.currentComponent,
#             equivalence: "",
#             alias: [fullName]
#         }
#         ctx.components[ctx.currentComponent].signals.push(fullName)
#         return {
#             type: "SIGNAL",
#             fullName: fullName,
#         }
#     })
#     return {
#         type: "VARIABLE",
#         name: ast.name.name,
#         selectors: []
#     }
# end
#
# function execDeclareVariable(ctx::Ctx, ast::Ast)
#     const scope = ctx.scopes[ctx.scopes.length-1]
#
#     if (ast.name.type != "VARIABLE") return error(ctx, ast, "Invalid linear combination name")
#     if (getScope(ctx, ast.name.name)) return error(ctx, ast, "Name already exists: "+ast.name.name)
#
#     const sizes=[]
#     for (let i=0; i< ast.name.selectors.length; i++) {
#         const size = exec(ctx, ast.name.selectors[i])
#         if (ctx.error) return
#
#         if (size.type != "NUMBER") return error(ctx, ast.name.selectors[i], "expected a number")
#         sizes.push( size.value.toJSNumber() )
#     }
#
#     scope[ast.name.name] = iterateSelectors(ctx, sizes, "", function() {
#         return {
#             type: "NUMBER",
#             value: bigInt(0)
#         }
#     })
#
#     return {
#         type: "VARIABLE",
#         name: ast.name.name,
#         selectors: []
#     }
# end
#
# function execVariable(ctx::Ctx, ast::Ast)
#     let v
#     try {
#         v = getScope(ctx, ast.name, ast.selectors)
#     } catch(err) {
#         console.log(JSON.stringify(ast, null,1))
#     }
#     if (ctx.error) return
#
#     if (!v) return error(ctx, ast, "Variable not defined")
#
#     // If the signal has an assigned value (constant) just return the constant
#     if ((v.type == "SIGNAL") && (ctx.signals[v.fullName].value)) {
#         return {
#             type: "NUMBER",
#             value: ctx.signals[v.fullName].value
#         }
#     }
#     let res
#     res=v
#     return res
# end
#
# function execPin(ctx::Ctx, ast::Ast)
#     const component = getScope(ctx, ast.component.name, ast.component.selectors)
#     if (!component) return error(ctx, ast.component, "Component does not exists: "+ast.component.name)
#     if (ctx.error) return
#     let signalFullName = component.fullName + "." + ast.pin.name
#     for (let i=0; i< ast.pin.selectors.length; i++) {
#         const sel = exec(ctx, ast.pin.selectors[i])
#         if (ctx.error) return
#
#         if (sel.type != "NUMBER") return error(ctx, ast.pin.selectors[i], "expected a number")
#
#         signalFullName += "[" + sel.value.toJSNumber() + "]"
#     }
#     if (!ctx.signals[signalFullName]) error(ctx, ast, "Signal not defined:" + signalFullName)
#     return {
#         type: "SIGNAL",
#         fullName: signalFullName
#     }
# end
#
# function execFor(ctx::Ctx, ast::Ast)
#
#     ctx.scopes.push({})
#     exec(ctx, ast.init)
#     if (ctx.error) return
#
#     let v = exec(ctx, ast.condition)
#     if (ctx.error) return
#
#     if (typeof v.value != "undefined") {
#         while ((v.value.neq(0))&&(!ctx.returnValue)) {
#             exec(ctx, ast.body)
#             if (ctx.error) return
#
#             exec(ctx, ast.step)
#             if (ctx.error) return
#
#             v = exec(ctx, ast.condition)
#             if (ctx.error) return
#         }
#     }
#     ctx.scopes.pop()
# end
#
# function execWhile(ctx::Ctx, ast::Ast)
#     let v = exec(ctx, ast.condition)
#     if (ctx.error) return
#
#     if (typeof v.value != "undefined") {
#         while ((v.value.neq(0))&&(!ctx.returnValue)) {
#             exec(ctx, ast.body)
#             if (ctx.error) return
#
#             v = exec(ctx, ast.condition)
#             if (ctx.error) return
#         }
#     }
# end
#
# function execIf(ctx::Ctx, ast::Ast)
#     let v = exec(ctx, ast.condition)
#     if (ctx.error) return
#
#     if (typeof v.value != "undefined") {
#         if ((v.value.neq(0))&&(!ctx.returnValue)) {
#             exec(ctx, ast.then)
#             if (ctx.error) return
#         } else {
#             if (ast.else) {
#                 exec(ctx, ast.else)
#                 if (ctx.error) return
#             }
#         }
#     }
# end
#
# function execVarAssignement(ctx::Ctx, ast::Ast)
#     let v
#     if (ast.values[0].type == "DECLARE") {
#         v = exec(ctx, ast.values[0])
#         if (ctx.error) return
#     } else {
#         v = ast.values[0]
#     }
#     const num = getScope(ctx, v.name, v.selectors)
#     if (ctx.error) return
#
#     if ((typeof(num) != "object")||(num == null)) return  error(ctx, ast, "Variable not defined")
#
#     if (num.type == "COMPONENT") return execInstantiateComponet(ctx, v, ast.values[1])
#     if (ctx.error) return
# //    if (num.type == "SIGNAL") return error(ctx, ast, "Cannot assign to a signal with `=` use <-- or <== ops")
#
#     const res = exec(ctx, ast.values[1])
#     if (ctx.error) return
#
#     setScope(ctx, v.name, v.selectors, res)
#
#     return v
# end
#
# function execLt(ctx::Ctx, ast::Ast)
#     const a = exec(ctx, ast.values[0])
#     if (ctx.error) return
#     if (a.type != "NUMBER") return { type: "NUMBER" }
#     const b = exec(ctx, ast.values[1])
#     if (ctx.error) return
#     if (b.type != "NUMBER") return { type: "NUMBER" }
#     if (!a.value || !b.value) return { type: "NUMBER" }
#     return {
#         type: "NUMBER",
#         value: a.value.lt(b.value) ? bigInt(1) : bigInt(0)
#     }
# end
#
# function execGt(ctx::Ctx, ast::Ast)
#     const a = exec(ctx, ast.values[0])
#     if (ctx.error) return
#     if (a.type != "NUMBER") return { type: "NUMBER" }
#     const b = exec(ctx, ast.values[1])
#     if (ctx.error) return
#     if (b.type != "NUMBER") return { type: "NUMBER" }
#     if (!a.value || !b.value) return { type: "NUMBER" }
#     return {
#         type: "NUMBER",
#         value: a.value.gt(b.value) ? bigInt(1) : bigInt(0)
#     }
# end
#
# function execLte(ctx::Ctx, ast::Ast)
#     const a = exec(ctx, ast.values[0])
#     if (ctx.error) return
#     if (a.type != "NUMBER") return { type: "NUMBER" }
#     const b = exec(ctx, ast.values[1])
#     if (ctx.error) return
#     if (b.type != "NUMBER") return { type: "NUMBER" }
#     if (!a.value || !b.value) return { type: "NUMBER" }
#     return {
#         type: "NUMBER",
#         value: a.value.lesserOrEquals(b.value) ? bigInt(1) : bigInt(0)
#     }
# end
#
# function execGte(ctx::Ctx, ast::Ast)
#     const a = exec(ctx, ast.values[0])
#     if (ctx.error) return
#     if (a.type != "NUMBER") return { type: "NUMBER" }
#     const b = exec(ctx, ast.values[1])
#     if (ctx.error) return
#     if (b.type != "NUMBER") return { type: "NUMBER" }
#     if (!a.value || !b.value) return { type: "NUMBER" }
#     return {
#         type: "NUMBER",
#         value: a.value.greaterOrEquals(b.value) ? bigInt(1) : bigInt(0)
#     }
# end
#
# function execEq(ctx::Ctx, ast::Ast)
#     const a = exec(ctx, ast.values[0])
#     if (ctx.error) return
#     if (a.type != "NUMBER") return { type: "NUMBER" }
#     const b = exec(ctx, ast.values[1])
#     if (ctx.error) return
#     if (b.type != "NUMBER") return { type: "NUMBER" }
#     if (!a.value || !b.value) return { type: "NUMBER" }
#     return {
#         type: "NUMBER",
#         value: a.value.eq(b.value) ? bigInt(1) : bigInt(0)
#     }
# end
#
# function execNeq(ctx::Ctx, ast::Ast)
#     const a = exec(ctx, ast.values[0])
#     if (ctx.error) return
#     if (a.type != "NUMBER") return { type: "NUMBER" }
#     const b = exec(ctx, ast.values[1])
#     if (ctx.error) return
#     if (b.type != "NUMBER") return { type: "NUMBER" }
#     if (!a.value || !b.value) return { type: "NUMBER" }
#     return {
#         type: "NUMBER",
#         value: a.value.eq(b.value) ? bigInt(0) : bigInt(1)
#     }
# end
#
# function execBAnd(ctx::Ctx, ast::Ast)
#     const a = exec(ctx, ast.values[0])
#     if (ctx.error) return
#     if (a.type != "NUMBER") return { type: "NUMBER" }
#     const b = exec(ctx, ast.values[1])
#     if (ctx.error) return
#     if (b.type != "NUMBER") return { type: "NUMBER" }
#     if (!a.value || !b.value) return { type: "NUMBER" }
#     return {
#         type: "NUMBER",
#         value: a.value.and(b.value).and(__MASK__)
#     }
# end
#
# function execAnd(ctx::Ctx, ast::Ast)
#     const a = exec(ctx, ast.values[0])
#     if (ctx.error) return
#     if (a.type != "NUMBER") return { type: "NUMBER" }
#     const b = exec(ctx, ast.values[1])
#     if (ctx.error) return
#     if (b.type != "NUMBER") return { type: "NUMBER" }
#     if (!a.value || !b.value) return { type: "NUMBER" }
#     return {
#         type: "NUMBER",
#         value: (a.value.neq(0) && a.value.neq(0)) ? bigInt(1) : bigInt(0)
#     }
# end
#
# function execOr(ctx::Ctx, ast::Ast)
#     const a = exec(ctx, ast.values[0])
#     if (ctx.error) return
#     if (a.type != "NUMBER") return { type: "NUMBER" }
#     const b = exec(ctx, ast.values[1])
#     if (ctx.error) return
#     if (b.type != "NUMBER") return { type: "NUMBER" }
#     if (!a.value || !b.value) return { type: "NUMBER" }
#     return {
#         type: "NUMBER",
#         value: (a.value.neq(0) || a.value.neq(0)) ? bigInt(1) : bigInt(0)
#     }
# end
#
# function execShl(ctx::Ctx, ast::Ast)
#     const a = exec(ctx, ast.values[0])
#     if (ctx.error) return
#     if (a.type != "NUMBER") return { type: "NUMBER" }
#     const b = exec(ctx, ast.values[1])
#     if (ctx.error) return
#     if (b.type != "NUMBER") return { type: "NUMBER" }
#     if (!a.value || !b.value) return { type: "NUMBER" }
#     const v = b.value.greater(256) ? 256 : b.value.value
#     return {
#         type: "NUMBER",
#         value: a.value.shiftLeft(v).and(__MASK__)
#     }
# end
#
# function execShr(ctx::Ctx, ast::Ast)
#     const a = exec(ctx, ast.values[0])
#     if (ctx.error) return
#     if (a.type != "NUMBER") return { type: "NUMBER" }
#     const b = exec(ctx, ast.values[1])
#     if (ctx.error) return
#     if (b.type != "NUMBER") return { type: "NUMBER" }
#     if (!a.value || !b.value) return { type: "NUMBER" }
#     const v = b.value.greater(256) ? 256 : b.value.value
#     return {
#         type: "NUMBER",
#         value: a.value.shiftRight(v).and(__MASK__)
#     }
# end
#
# function execMod(ctx::Ctx, ast::Ast)
#     const a = exec(ctx, ast.values[0])
#     if (ctx.error) return
#     if (a.type != "NUMBER") return { type: "NUMBER" }
#     const b = exec(ctx, ast.values[1])
#     if (ctx.error) return
#     if (b.type != "NUMBER") return { type: "NUMBER" }
#     if (!a.value || !b.value) return { type: "NUMBER" }
#     return {
#         type: "NUMBER",
#         value: a.value.mod(b.value)
#     }
# end
#
# function execExp(ctx::Ctx, ast::Ast)
#     const a = exec(ctx, ast.values[0])
#     if (ctx.error) return
#     if (a.type != "NUMBER") return { type: "NUMBER" }
#     const b = exec(ctx, ast.values[1])
#     if (ctx.error) return
#     if (b.type != "NUMBER") return { type: "NUMBER" }
#     if (!a.value || !b.value) return { type: "NUMBER" }
#     return {
#         type: "NUMBER",
#         value: a.value.modPow(b.value, __P__)
#     }
# end
#
# function execDiv(ctx::Ctx, ast::Ast)
#     const a = exec(ctx, ast.values[0])
#     if (ctx.error) return
#     if (a.type != "NUMBER") return { type: "NUMBER" }
#     const b = exec(ctx, ast.values[1])
#     if (ctx.error) return
#     if (b.type != "NUMBER") return { type: "NUMBER" }
#     if (!a.value || !b.value) return { type: "NUMBER" }
#     if (b.value.isZero()) return error(ctx, ast, "Division by zero")
#     return {
#         type: "NUMBER",
#         value: a.value.times(b.value.modInv(__P__)).mod(__P__)
#     }
# end
#
# function execIDiv(ctx::Ctx, ast::Ast)
#     const a = exec(ctx, ast.values[0])
#     if (ctx.error) return
#     if (a.type != "NUMBER") return { type: "NUMBER" }
#     const b = exec(ctx, ast.values[1])
#     if (ctx.error) return
#     if (b.type != "NUMBER") return { type: "NUMBER" }
#     if (!a.value || !b.value) return { type: "NUMBER" }
#     if (b.value.isZero()) return error(ctx, ast, "Division by zero")
#     return {
#         type: "NUMBER",
#         value: a.value.divide(b.value)
#     }
# end
#
# function execAdd(ctx::Ctx, ast::Ast)
#     const a = exec(ctx, ast.values[0])
#     if (ctx.error) return
#     const b = exec(ctx, ast.values[1])
#     if (ctx.error) return
#
#     const res = lc.add(a,b)
#     if (res.type == "ERROR") return error(ctx, ast, res.errStr)
#
#     return res
# end
#
# function execSub(ctx::Ctx, ast::Ast)
#     const a = exec(ctx, ast.values[0])
#     if (ctx.error) return
#     const b = exec(ctx, ast.values[1])
#     if (ctx.error) return
#
#     const res = lc.sub(a,b)
#     if (res.type == "ERROR") return error(ctx, ast, res.errStr)
#
#     return res
# end
#
# function execUMinus(ctx::Ctx, ast::Ast)
#     const a = exec(ctx, ast.values[0])
#     if (ctx.error) return
#
#     const res = lc.negate(a)
#     if (res.type == "ERROR") return error(ctx, ast, res.errStr)
#
#     return res
# end
#
# function execMul(ctx::Ctx, ast::Ast)
#     const a = exec(ctx, ast.values[0])
#     if (ctx.error) return
#     const b = exec(ctx, ast.values[1])
#     if (ctx.error) return
#
#     const res = lc.mul(a,b)
#     if (res.type == "ERROR") return error(ctx, ast, res.errStr)
#
#     return res
# end
#
# function execVarAddAssignement(ctx::Ctx, ast::Ast)
#     const res = execAdd(ctx,{ values: [ast.values[0], ast.values[1]] } )
#     if (ctx.error) return
#     return execVarAssignement(ctx, { values: [ast.values[0], res] })
# end
#
# function execVarMulAssignement(ctx::Ctx, ast::Ast)
#     const res = execMul(ctx,{ values: [ast.values[0], ast.values[1]] } )
#     if (ctx.error) return
#     return execVarAssignement(ctx, { values: [ast.values[0], res] })
# end
#
# function execPlusPlusRight(ctx::Ctx, ast::Ast)
#     const resBefore = exec(ctx, ast.values[0])
#     if (ctx.error) return
#     const resAfter = execAdd(ctx,{ values: [ast.values[0], {type: "NUMBER", value: bigInt(1)}] } )
#     if (ctx.error) return
#     execVarAssignement(ctx, { values: [ast.values[0], resAfter] })
#     return resBefore
# end
#
# function execPlusPlusLeft(ctx::Ctx, ast::Ast)
#     if (ctx.error) return
#     const resAfter = execAdd(ctx,{ values: [ast.values[0], {type: "NUMBER", value: bigInt(1)}] } )
#     if (ctx.error) return
#     execVarAssignement(ctx, { values: [ast.values[0], resAfter] })
#     return resAfter
# end
#
# function execMinusMinusRight(ctx::Ctx, ast::Ast)
#     const resBefore = exec(ctx, ast.values[0])
#     if (ctx.error) return
#     const resAfter = execSub(ctx,{ values: [ast.values[0], {type: "NUMBER", value: bigInt(1)}] } )
#     if (ctx.error) return
#     execVarAssignement(ctx, { values: [ast.values[0], resAfter] })
#     return resBefore
# end
#
# function execMinusMinusLeft(ctx::Ctx, ast::Ast)
#     if (ctx.error) return
#     const resAfter = execSub(ctx,{ values: [ast.values[0], {type: "NUMBER", value: bigInt(1)}] } )
#     if (ctx.error) return
#     execVarAssignement(ctx, { values: [ast.values[0], resAfter] })
#     return resAfter
# end
#
# function execTerCon(ctx::Ctx, ast::Ast)
#     const cond = exec(ctx, ast.values[0])
#     if (ctx.error) return
#
#     if (!cond.value) return { type: "NUMBER" }
#
#     if (cond.value.neq(0)) {
#         return exec(ctx, ast.values[1])
#     } else {
#         return exec(ctx, ast.values[2])
#     }
# end
#
# function execSignalAssign(ctx::Ctx, ast::Ast)
#     let vDest
#     if (ast.values[0].type == "DECLARE") {
#         vDest = exec(ctx, ast.values[0])
#         if (ctx.error) return
#     } else {
#         vDest = ast.values[0]
#     }
#
#     let dst
#     if (vDest.type == "VARIABLE") {
#         dst = getScope(ctx, vDest.name, vDest.selectors)
#         if (ctx.error) return
#     } else if (vDest.type == "PIN") {
#         dst = execPin(ctx, vDest)
#         if (ctx.error) return
#     } else {
#         error(ctx, ast, "Bad assignement")
#     }
#
#     if (!dst) return  error(ctx, ast, "Signal not defined")
#     if (dst.type != "SIGNAL") return  error(ctx, ast, "Signal assigned to a non signal")
#
#     let sDest=ctx.signals[dst.fullName]
#     if (!sDest) return error(ctx, ast, "Invalid signal: "+dst.fullName)
#     while (sDest.equivalence) sDest=ctx.signals[sDest.equivalence]
#
#     if (sDest.value) return error(ctx, ast, "Signals cannot be assigned twice")
#
#     let src = exec(ctx, ast.values[1])
#     if (ctx.error) return
#
#
#     /*
#     let vSrc
#     if (ast.values[1].type == "DECLARE") {
#         vSrc = exec(ctx, ast.values[1])
#         if (ctx.error) return
#     } else {
#         vSrc = ast.values[1]
#     }
#
#     if (vSrc.type == "VARIABLE") {
#         src = getScope(ctx, vSrc.name, vSrc.selectors)
#         if (!src) error(ctx, ast, "Variable not defined: " + vSrc.name)
#         if (ctx.error) return
#     } else if (vSrc.type == "PIN") {
#         src = execPin(ctx, vSrc)
#     }
#     */
#
#     let assignValue = true
#     if (src.type == "SIGNAL") {
#         sDest.equivalence = src.fullName
#         sDest.alias = sDest.alias.concat(src.alias)
#         while (sDest.equivalence) sDest=ctx.signals[sDest.equivalence]
#         assignValue = false
#     }
#
#     if (assignValue) {
#         //        const resLC = exec(ctx, vSrc)
#         if (ctx.error) return
#
#         //        const v = lc.evaluate(ctx, resLC)
#         const v = lc.evaluate(ctx, src)
#
#         if (v.value) {
#             sDest.value = v.value
#         }
#     }
#
#     return vDest
# end
#
# function execConstrain(ctx::Ctx, ast::Ast)
#     const a = exec(ctx, ast.values[0])
#     if (ctx.error) return
#     const b = exec(ctx, ast.values[1])
#     if (ctx.error) return
#
#     const res = lc.sub(a,b)
#     if (res.type == "ERROR") return error(ctx, ast, res.errStr)
#
#     if (!lc.isZero(res)) {
#         ctx.constraints.push(lc.toQEQ(res))
#     }
#
#     return res
# end
#
# function execSignalAssignConstrain(ctx::Ctx, ast::Ast)
#     const v = execSignalAssign(ctx,endst)
#     if (ctx.error) return
#     execConstrain(ctx::Ctx, ast::Ast
#     if (ctx.error) return
#     return v
# end
#
# function execInclude(ctx::Ctx, ast::Ast)
#     const incFileName = path.resolve(ctx.filePath, ast.file)
#     const incFilePath = path.dirname(incFileName)
#
#     ctx.includedFiles = ctx.includedFiles || []
#     if (ctx.includedFiles[incFileName]) return
#
#     ctx.includedFiles[incFileName] = true
#
#     const src = fs.readFileSync(incFileName, "utf8")
#
#     if (!src) return error(ctx, ast, "Include file not found: "+incFileName)
#
#     const incAst = parser.parse(src)
#
#     const oldFilePath = ctx.filePath
#     const oldFileName = ctx.fileName
#     ctx.filePath = incFilePath
#     ctx.fileName = incFileName
#
#     exec(ctx, incAst)
#
#     ast.block = incAst
#
#     ctx.filePath = oldFilePath
#     ctx.fileName = oldFileName
# end
#
# function execArray(ctx::Ctx, ast::Ast)
#     const res = []
#
#     for (let i=0; i<ast.values.length; i++) {
#         res.push(exec(ctx, ast.values[i]))
#     }
#
#     return res
# end
#
# function copyScope(scope::Scope)
#     var scopesClone = []
#     for (let i=0; i<scope.length; i++) {
#         scopesClone.push(scope[i])
#     }
#     return scopesClone
# end
#


text:
# ~/.julia/config/startup.jl

let
    ESC = '\033'
    OSC = ']'    # Operating System Control code
    Ps  = 0      # Change Icon Name and Window Title to Pt
    BEL = '\007'
    f(Pt) = string(ESC, OSC, Ps, ';', Pt, BEL)
    commit_date = (first ∘ split)(Base.GIT_VERSION_INFO.date_string)
    (print ∘ f)("Julia $VERSION ($commit_date)")
end


text:
function _concatenated_seq_names(msas...; delim::String="_&_")
	concatenated_names = sequencenames(msas[1])
	for msa in msas[2:end]
		for (i, seq) in enumerate(sequencenames(msa))
			if seq == concatenated_names[i]
				continue
			end
			concatenated_names[i] *= delim * seq
		end
	end
	concatenated_names
end

function _concatenated_col_names(msas...)
	colnames = String[]
	msa_number = 0
	for msa in msas
		columns = columnnames(msa)
		first_col = first(columns)
		check_msa_change = '_' in first_col
		previous = ""
		msa_number += 1
		for col in columns
			if check_msa_change
				fields = split(col, '_')
				current = first(fields)
				if current != previous
					if previous != ""
						msa_number += 1
					end
					previous = string(current)
				end
				col = last(fields)
			end
			push!(colnames, "$(msa_number)_$col")
		end
	end
	colnames
end

_get_seq_lengths(msas...) = [ncolumns(msa) for msa in msas]

function _get_annot_types(fun, index, data::Annotations...)
	union(Set(k[index] for k in keys(fun(annot))) for annot in data)
end

function _concatenate_annotfile(data::Annotations...)
	annotfile = copy(getannotfile(data[1]))
	for ann in data[2:end]
		for (k, v) in getannotfile(ann)
			if haskey(annotfile, k)
				if k == "ColMap"
					annotfile[k] *= ',' * v
				else
					annotfile[k] *= "_&_" * v
				end
			else
				push!(annotfile, k => v)
			end
		end
	end
	annotfile
end

function _get_seqname_mapping(concatenated_seqnames, msas...)
	mapping = Dict{Tuple{Int, String}, String}()
	seq_names = hcat([sequencenames(msa) for msa in msas]...)
	nseq, nmsa = size(seq_names)
	@assert nseq == length(concatenated_seqnames)
	for i in 1:nseq
		for j in 1:nmsa
			mapping[(j, seq_names[i, j])] = concatenated_seqnames[i]
		end
	end
	mapping
end

function _concatenate_annotsequence(seqname_mapping, data::Annotations...)
	annotsequence = Dict{Tuple{String,String},String}()
	for (i, ann) in enumerate(data)
		for ((seqname, annot_name), value) in getannotsequence(ann)
			concatenated_seqname = get(seqname_mapping, (i, seqname), seqname)
			new_key = (concatenated_seqname, annot_name)
			if haskey(annotsequence, new_key)
				if annot_name == "SeqMap"
					annotsequence[new_key] *= ',' * value
				else
					annotsequence[new_key] *= "_&_" * value
				end
			else
				push!(annotsequence, new_key => value)
			end
		end
	end
	annotsequence
end

function _fill_and_update!(dict, last, key, i, value, seq_lengths)
	if haskey(dict, key)
		if last[key] == i - 1
			dict[key] *= value
		else
			previous = sum(seq_lengths[(last[key]+1):(i-1)])
			dict[key] *= " "^previous * value
		end
		last[key] = i
	else
		if i == 1
			push!(dict, key => value)
		else
			previous = sum(seq_lengths[1:(i-1)])
			push!(dict, key => " "^previous * value)
		end
		push!(last, key => i)
	end
end

function _fill_end!(dict, seq_lengths, entity)
	total_length = sum(seq_lengths)
	for (key, value) in dict
		current_length = length(value)
		if current_length < total_length
			dict[key] *= " "^(total_length - current_length)
		elseif current_length > total_length
			throw(ErrorException(
				"There are $current_length $entity annotated instead of $total_length."))
		end
	end
	dict
end

function _concatenate_annotcolumn(seq_lengths, data::Annotations...)
	annotcolumn = Dict{String,String}()
	last = Dict{String,Int}()
	for (i, ann) in enumerate(data)
		for (annot_name, value) in getannotcolumn(ann)
			_fill_and_update!(annotcolumn, last, annot_name, i, value, seq_lengths)
		end
	end
	_fill_end!(annotcolumn, seq_lengths, "columns")
end

function _concatenate_annotresidue(seq_lengths, seqname_mapping, data::Annotations...)
	annotresidue = Dict{Tuple{String,String},String}()
	last = Dict{Tuple{String,String},Int}()
	for (i, ann) in enumerate(data)
		for ((seqname, annot_name), value) in getannotresidue(ann)
			concatenated_seqname = get(seqname_mapping, (i, seqname), seqname)
			new_key = (concatenated_seqname, annot_name)
			_fill_and_update!(annotresidue, last, new_key, i, value, seq_lengths)
		end
	end
	_fill_end!(annotresidue, seq_lengths, "residues")
end

function Base.hcat(msa::T...) where T <: AnnotatedAlignedObject
	concatenated_msa = NamedArray(hcat(getresidues.(msa)...))
	seqnames = _concatenated_seq_names(msa...)
	colnames = _concatenated_col_names(msa...)
	setnames!(concatenated_msa, seqnames, 1)
	setnames!(concatenated_msa, colnames, 2)
	seqname_mapping = _get_seqname_mapping(seqnames, msa...)
	seq_lengths = _get_seq_lengths(msa...)
	old_annot = annotations.([msa...])
	new_annot = Annotations(
		_concatenate_annotfile(old_annot...),
		_concatenate_annotsequence(seqname_mapping, old_annot...),
		_concatenate_annotcolumn(seq_lengths, old_annot...),
		_concatenate_annotresidue(seq_lengths, seqname_mapping, old_annot...)
	)
	if haskey(new_annot.file, "HCat")
		delete!(new_annot.file, "HCat")
	end
	setannotfile!(
		new_annot,
		"HCat",
		join((replace(col, r"_[0-9]+$" => "") for col in colnames), ',')
	)
	T(concatenated_msa, new_annot)
end

function Base.hcat(msa::T...) where T <: UnannotatedAlignedObject
	concatenated_msa = NamedArray(hcat(getresidues.(msa)...))
	seqnames = _concatenated_seq_names(msa...)
	colnames = _concatenated_col_names(msa...)
	setnames!(concatenated_msa, seqnames, 1)
	setnames!(concatenated_msa, colnames, 2)
	T(concatenated_msa)
end

"""
It returns a vector of numbers from `1` to N for each column that indicates the source MSA.
The mapping is annotated in the `"HCat"` file annotation of an
`AnnotatedMultipleSequenceAlignment` or in the column names of an `NamedArray` or
`MultipleSequenceAlignment`.
"""
function gethcatmapping(msa::AnnotatedMultipleSequenceAlignment)
    annot = getannotfile(msa)
    if haskey(annot, "HCat")
        return _str2int_mapping(annot["HCat"])
    else
        return gethcatmapping(namedmatrix(msa))
    end
end

function gethcatmapping(msa::NamedResidueMatrix{AT}) where AT <: AbstractMatrix
	colnames = names(msa, 2)
	if !isempty(colnames)
		if !occursin('_', colnames[1])
			throw(ErrorException(
				"The column names have not generated by `hcat` on an `AnnotatedMultipleSequenceAlignment` or `MultipleSequenceAlignment`."
				))
		end
    	Int[ parse(Int, replace(col, r"_[0-9]+$" => "")) for col in colnames]
	else
		throw(ErrorException("There are not column names!"))
	end
end

gethcatmapping(msa::MultipleSequenceAlignment) = gethcatmapping(namedmatrix(msa))


text:
# Reset the search direction if it becomes corrupted
# return true if the direction was changed
reset_search_direction!(state, d, method) = false # no-op

_alphaguess(a) = a
_alphaguess(a::Number) = LineSearches.InitialStatic(alpha=a)

# Note that for these resets we're using `gradient(d)` but we don't need to use
# project_tangent! here, because we already did that inplace on gradient(d) after
# the last evaluation (we basically just always do it)
function reset_search_direction!(state, d, method::BFGS)
    n = length(state.x)
    T = eltype(state.x)

    if method.initial_invH == nothing
        if method.initial_stepnorm == nothing
            state.invH .= Matrix{T}(I, n, n)
        else
            initial_scale = method.initial_stepnorm * inv(norm(gradient(d), Inf))
            state.invH.= Matrix{T}(initial_scale*I, n, n)
        end
    else
        state.invH .= method.initial_invH(state.x)
    end
#    copyto!(state.invH, method.initial_invH(state.x))
    state.s .= .-gradient(d)
    return true
end

function reset_search_direction!(state, d, method::LBFGS)
    state.pseudo_iteration = 1
    state.s .= .-gradient(d)
    return true
end

function reset_search_direction!(state, d, method::ConjugateGradient)
    state.s .= .-state.pg
    return true
end

function perform_linesearch!(state, method, d)
    # Calculate search direction dphi0
    dphi_0 = real(dot(gradient(d), state.s))
    # reset the direction if it becomes corrupted
    if dphi_0 >= zero(dphi_0) && reset_search_direction!(state, d, method)
        dphi_0 = real(dot(gradient(d), state.s)) # update after direction reset
    end
    phi_0  = value(d)

    # Guess an alpha
    method.alphaguess!(method.linesearch!, state, phi_0, dphi_0, d)

    # Store current x and f(x) for next iteration
    state.f_x_previous = phi_0
    copyto!(state.x_previous, state.x)

    # Perform line search; catch LineSearchException to allow graceful exit
    try
        state.alpha, ϕalpha =
            method.linesearch!(d, state.x, state.s, state.alpha,
                               state.x_ls, phi_0, dphi_0)
        return true # lssuccess = true
    catch ex
        if isa(ex, LineSearches.LineSearchException)
            state.alpha = ex.alpha
            # We shouldn't warn here, we should just carry it to the output
            # @warn("Linesearch failed, using alpha = $(state.alpha) and
            # exiting optimization.\nThe linesearch exited with message:\n$(ex.message)")
            return false # lssuccess = false
        else
            rethrow(ex)
        end
    end
end


text:
baremodule ThreadedScans

function scan! end

function simple! end
function dac! end
function linear! end
function partitioned_hillis_steele! end

module Internal

using ArgCheck: ArgCheck, @argcheck, @check

using ..ThreadedScans: ThreadedScans

include("utils.jl")
include("promises.jl")
include("basecases.jl")
include("simple.jl")
include("dac.jl")
include("linear.jl")
include("hillis_steele.jl")
include("api.jl")

# Use README as the docstring of the module:
@doc let path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    replace(read(path, String), r"^```julia"m => "```jldoctest README")
end ThreadedScans

end  # module Internal

end  # baremodule ThreadedScans


text:
using RoughSurface
using Documenter
DocMeta.setdocmeta!(RoughSurface, :DocTestSetup, :(using RoughSurface); recursive=true)

makedocs(;
    modules=[RoughSurface],
    authors="nil-is-lin",
    repo="https://github.com/nil-is-lin/RoughSurface.jl/blob/{commit}{path}#L{line}",
    sitename="RoughSurface.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://nil-is-lin.github.io/RoughSurface.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Test" => "test.md"
    ],
)

deploydocs(;
    repo="github.com/nil-is-lin/RoughSurface.jl",
)


text:
# # [Open Circular Port Graph Examples](@id cpg_examples)
#
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/examples/CPG_examples.ipynb)

using AlgebraicDynamics.DWDDynam
using AlgebraicDynamics.CPortGraphDynam
using AlgebraicDynamics.CPortGraphDynam: draw, barbell, gridpath, grid, meshpath

using Catlab
using Catlab.WiringDiagrams
using Catlab.WiringDiagrams.CPortGraphs
using Catlab.Theories
using Catlab.CategoricalAlgebra

using OrdinaryDiffEq
using DynamicalSystems
using Plots, Plots.PlotMeasures

using PrettyTables

# ## SIR Epidemiology Model

# An SIR epidemiology model has three types of people: susceptible, infected, and recovered. When a susceptible person interacts with infected person, the susceptible person also becomes infected. Over time infected people  recover. Transition rates determine the how frequently susceptible people come into contact with infected people and how fast infected people recover. The system evolves according to the law of mass action.

# In a multi-city SIR model, each city has susceptible, infected, and recovered populations. To see the spread of the disease we will consider both susceptible and infected people moving between cities. To define a multi-city SIR model, we can compose multiple single-city SIR models using the composition syntax of open CPGs. The composition pattern will consist of three boxes each of which will be filled by a single-city SIR model. Ports expose the susceptible and infected populations of each city. One set of wires connect the susceptible and infected populations of cities 1 and 2. A second set of wires connect the susceptible and infected popuation of cities 2 and 3.

## Define the composition pattern
d₀ = OpenCPortGraph()
add_parts!(d₀, :Box, 1)
d₁ = barbell(2)
F = ACSetTransformation((Box=[2],), d₀, d₁)
G = ACSetTransformation((Box=[1],), d₀, d₁)
d₂ = apex(pushout(F,G))

## Define the primitive systems
β, μ, α₁, α₂ = 0.4, 0.4, 0.01, 0.01

sirfuncb = (u,x,p,t)->[-β*u[1]*u[2] - α₁*(u[1]-x[1]), # Ṡ
                        β*u[1]*u[2] - μ*u[2] - α₂*(u[2]-x[2]), #İ
                        μ*u[2] # Ṙ
                        ]
sirfuncm = (u,x,p,t)->[-β*u[1]*u[2] - α₁*(u[1]-(x[1]+x[3])/2),
                        β*u[1]*u[2] - μ*u[2] - α₂*(u[2]-(x[2]+x[4])/2),
                        μ*u[2]
                       ]

boundary  = ContinuousMachine{Float64}(2,3,sirfuncb, u->u[1:2])
middle    = ContinuousMachine{Float64}(4,3, sirfuncm, u->u[[1,2,1,2]])

## Compose
threecity = oapply(d₂, [boundary,middle,boundary])

# First, we will approximate the solution to the three city SIR model using Euler's method. The initial condition has 100 susceptible people in each city a single infected person in the first city. We show the infected populations in each city over time.

u0 = [100,1,0,100,0,0,100,0,0.0]

h = 0.01
nsteps = 100

threecity_approx = euler_approx(threecity, h)
traj = trajectory(threecity_approx, u0, [], nothing, nsteps)

map(traj) do u
    return (i1=u[2], i2=u[5], i3=u[8])
end |> pretty_table

# Next, we will solve the continuous system and plot the results. Over time the infected populations increase and the susceptible populations decrease. The delays in the plots illustrate how the disease spreads from city 1 to city 2 and then from city 2 to city 3.

## Solve and plot
tspan = (0.0, 1.0)

prob = ODEProblem(threecity, u0, tspan)
sol = solve(prob, Tsit5(); dtmax = 0.01)

plot(sol, lw=2, title = "SIR Epidemiology Model", bottom_margin=10mm, left_margin=10mm,
    label=["S" "I" "R"])

# ## Cellular automata

#  Circular port graphs are particularly useful for modeling systems where the composition pattern is given by a grid and where the dynamics are repetative. In the case of cellular automata the composition pattern is a row of $n$ cells each of which is connected to its two neighbors. The primitive systems are identical machines whose discrete dynamics are a specified rule. See [here](https://mathworld.wolfram.com/ElementaryCellularAutomaton.html) for a complete set of rules and the patterns they generate.

function Rule(k::Int)
    (left_neighbor, x, right_neighbor) ->
    Bool(digits(k, base=2, pad=8)[1 + right_neighbor + 2*x + 4*left_neighbor])
end

## Define the composition pattern
n = 100
row = apex(gridpath(n, 1))

## Define the primitive system which will be repeated for each cell
rule = DiscreteMachine{Bool}(2, 1, 2, (u, x, p, t)->Rule(p)(x[2], u[1], x[1]),
            u->[u[1], u[1]])

## Compose
automaton = oapply(row, rule)

## Solve and plot
u0 = zeros(Int, n); u0[Int(n/2)] = 1

rule_number = 126
traj = trajectory(automaton, u0, [0,0], rule_number, 100)
spy(Matrix(traj))


text:
loss1 = RobustModels.L2Loss()
loss2 = RobustModels.TukeyLoss()

@testset "linear: Expectile estimators" begin
    m  = fit(RobustLinearModel, form, data, MEstimator(loss1); method=:cg, maxiter=50, verbose=false, initial_scale=:mad)
    m1 = fit(RobustLinearModel, form, data, GeneralizedQuantileEstimator(loss1, 0.5); method=:cg, maxiter=50, verbose=false, initial_scale=:mad)

    # refit
    # cannot change type to GeneralQuantileEstimator
    @test_throws TypeError refit!(m; quantile=0.5)

    @testset "$(τ) $(name)" for τ in range(0.1, 0.9, step=0.1), name in ("Expectile", "Quantile")
        println("\n\t\u25CF Estimator: $name($τ)")
        est = getproperty(RobustModels, Symbol(name * "Estimator"))(τ)
        m2 = fit(RobustLinearModel, form, data, est; method=:cg, maxiter=50, rtol=1e-4, initial_scale=:mad)
        m3 = fit(RobustLinearModel, form, data, est; method=:chol, maxiter=50, rtol=1e-4, initial_scale=:mad)
        println("rlm(cg)  : ", coef(m2))
        println("rlm(chol): ", coef(m3))
        println(m2)
        if name == "Expectile"
            @test isapprox(coef(m2), coef(m3); rtol=1e-4)
            if τ==0.5
                @test isapprox(coef(m), coef(m2); rtol=1e-4)
            end

            # refit with new quantile
            refit!(m1; quantile=τ)
            @test isapprox(coef(m1), coef(m2); rtol=1e-4)
        else
            @test isapprox(coef(m2), coef(m3); rtol=1e-2)
        end
    end
end

@testset "linear: M-Quantile estimators" begin
    @testset "$(name) estimator" for name in ("Huber", "L1L2", "Fair", "Logcosh", "Arctan") #, "Cauchy", "Geman", "Welsch", "Tukey", "YohaiZamar")
        println("\n\t\u25CF M-Quantile Estimator: $name")
        loss = getproperty(RobustModels, Symbol(name * "Loss"))()
        est = MEstimator(loss)
        m1 = fit(RobustLinearModel, form, data, est; method=:cg, maxiter=50, verbose=false, initial_scale=:mad)
        println("no quant, rlm(cg)  : ", coef(m1))
        est1 = GeneralizedQuantileEstimator(loss)
        @testset "$(τ) quantile" for τ in range(0.1, 0.9, step=0.1)
            est = GeneralizedQuantileEstimator(loss, τ)
            m2 = fit(RobustLinearModel, form, data, est; method=:cg, maxiter=50, verbose=false, initial_scale=:mad)
            m3 = fit(RobustLinearModel, form, data, est; method=:chol, maxiter=50, initial_scale=:mad)
            m4 = fit(RobustLinearModel, form, data, est1; quantile=τ, method=:chol, maxiter=50, initial_scale=:mad)
            println("quant $(τ) rlm(cg)   : ", coef(m2))
            println("quant $(τ) rlm(chol) : ", coef(m3))
            println("quant $(τ) rlm2(chol): ", coef(m4))
            if τ==0.5
                @test isapprox(coef(m1), coef(m2); rtol=1e-5)
            end
            if name != "L1"
                @test isapprox(coef(m2), coef(m3); rtol=1e-2)
            end
            @test isapprox(coef(m3), coef(m4); rtol=1e-5)
        end
    end
end


text:
#=
Square-function with gradient
=#
function f_square(x)
    return x^2
end

function g_square(x)
    return 2*x
end

#=
Logistic-regression loss-function, gradient and sub-gradient
=#

function f_logreg()
    return function(x_i,y_i,θ,λ)
        return log(1+exp(-y_i*dot(θ,x_i)))
    end
end

function f_logreg(X,y,θ,λ)
    return function(θ)
        data = X
        labels = y
        λ = λ
        funcVal = 0
        # size(data,1) = 1000 datapoints
        for i in 1:size(data,1)
            # x = 784*1 column vector
            x_i = data[i:i,:]'
            funcVal += log( 1+ exp(-labels[i] * dot(θ,x_i) ) )
        end
        funcVal += λ/2*(norm(θ,2)^2)
        return funcVal
    end
end

function g_logreg(X,y,θ,λ)
    return function(θ)
        data = X
        labels = y
        λ = λ
        # n-dim-column-vector
        array = zeros(size(data,2))
        # add l2 norm
        array += λ * θ
        # iterate datapoints
        for i in 1:size(data,1)
            # pick i-th datapoint
        	x_i = data[i:i,:]'
            exponential = exp(-labels[i]*dot(θ,x_i))
            array += (-labels[i]*x_i * exponential)/(1+exponential)
        end
        return array
    end
end

function sub_g_logreg()
        return function(x_i,y_i,θ,λ)
              exponential = exp(-y_i*dot(θ,x_i))
              return (-y_i * x_i * exponential) / ( 1 + exponential ) + λ * θ
        end
end

#=
hinge-loss - Loss-function and gradient
=#
function f_svm(X,y,θ,λ)
    return function(θ)
        X = X
        y = y
        λ = λ
        fval = 0
        # datapoints
        for i in 1:size(X,1)
            x_i = X[i:i,:]'
            fval += max(0, (1 - y[i]*dot(θ,x_i)))
        end
        fval = fval / size(X,1)
        fval += λ/2 * norm(θ,2)
        return fval
    end
end

#=
Compute the gradient for one datapoint for hinge-loss.
=#
function sub_g_svm()
    return function(x_i,y_i,θ,λ)
      val = 1 - (y_i * dot(θ,x_i))
      if val < 0
          return λ * θ;
      else
          return -y_i * x_i + λ * θ;
      end
  end
end

#=
Rosenbrock-function and its gradient.
=#
function f_rosenbrock(x_k)
    return 100*(x_k[2]-x_k[1]^2)^2+(1-x_k[1])^2
end

function g_rosenbrock(x_k)
    return [ -400*x_k[1]*(x_k[2]-x_k[1]^2)-2(1-x_k[1]) ; 200*x_k[1]*(x_k[2]-x_k[1]^2) ]
end

max_stars_repo_name: zagazao/num-opt

max_issues_repo_name: zagazao/num-opt

max_forks_repo_name: zagazao/num-opt


text:
__precompile__()

#===============================================================================
 This package contains the metaheuristic the feasibility pump using CPLEX
===============================================================================#

module FeasibilityPump

using CPLEX
using LinearAlgebra
using Random
using SparseArrays

include("solver.jl")
include("rounding.jl")
include("projection.jl")
include("perturbation.jl")
include("constants.jl")

#export ParametersFP
export feasibilitypump
export setparamfp!
export getparamnamefp
export getparamfp
export initializeparametersfp
export __compilefp__

end # module

max_stars_repo_name: remi-garcia/FeasibilityPump

max_issues_repo_name: remi-garcia/FeasibilityPump

max_forks_repo_name: remi-garcia/FeasibilityPump


text:
module Day10
export solvePart1, solvePart2

struct Point
    position::Vector{Int64}
    velocity::Vector{Int64}
end

function getInput()
    open("day10input.txt") do f
        posvel = []
        for line in eachline(f)
            xLoc = findnext(r"[-0-9]+", line, 1)
            x = parse(Int64, line[xLoc])
            yLoc = findnext(r"[-0-9]+", line, xLoc[end]+1)
            y = parse(Int64, line[yLoc])
            xVelLoc = findnext(r"[-0-9]+", line, yLoc[end]+1)
            xVel = parse(Int64, line[xVelLoc])
            yVelLoc = findnext(r"[-0-9]+", line, xVelLoc[end]+1)
            yVel = parse(Int64, line[yVelLoc])
            push!(posvel, Point([x, y], [xVel, yVel]))
        end
        return posvel
    end
end

function getBoundingBox(points)
    xs = map(point -> point.position[1], points)
    ys = map(point -> point.position[2], points)
    return min(xs...), max(xs...), min(ys...), max(ys...)
end

function getBoundingBoxMagnitude(points)
    xMin, xMax, yMin, yMax = getBoundingBox(points)
    (xMax - xMin) + (yMax - yMin)
end

function getSuccessorPoints(points)
    next = []
    for point in points
        newPosition = point.position + point.velocity
        push!(next, Point(newPosition, point.velocity))
    end
    return next
end

function sPrintPoints(points)
    xMin, xMax, yMin, yMax = getBoundingBox(points)

    locContainsPoint = Dict()
    for point in points
        locContainsPoint[point.position] = true
    end

    s = []

    for y in yMin:yMax
        for x in xMin:xMax
            if get(locContainsPoint, [x, y], false)
                push!(s, "#")
            else
                push!(s, " ")
            end
        end
        push!(s, "\n")
    end
    return join(s, "")
end

function solve()
    points = getInput()
    size = getBoundingBoxMagnitude(points)
    next = getSuccessorPoints(points)
    nextSize = getBoundingBoxMagnitude(next)
    seconds = 0
    while nextSize < size
        points = next
        size = nextSize
        next = getSuccessorPoints(points)
        nextSize = getBoundingBoxMagnitude(next)
        seconds += 1
    end
    return (sPrintPoints(points), seconds)
end

function solvePart1()
    return solve()[1]
end

function solvePart2()
    return solve()[2]
end
end


text:
# candlesticks

# Makie
# Pkg.add("MarketData")
using CairoMakie, Random, TimeSeries, MarketData, Dates
let
    tempo = string.(timestamp(MarketData.ohlc))
    lentime = length(tempo)
    slice_dates = range(1,lentime, step = lentime÷20)

    colors = values(MarketData.ohlc.Close) .> values(MarketData.ohlc.Open)
    lowV = values(MarketData.ohlc.Low)
    highV = values(MarketData.ohlc.High)
    linesegs = []
    for i in 1:lentime
        push!(linesegs, CairoMakie.Point2f0(i, lowV[i]))
        push!(linesegs, CairoMakie.Point2f0(i, highV[i]))
    end
    linesegs = CairoMakie.Point2f0.(linesegs)
    cmap = ("#f64325", "#78f518")

    fig = Figure(resolution = (1200, 700), font = "sans", fontsize = 20)
    ax = Axis(fig, backgroundcolor = :black, ygridcolor = "#65866b",
        xgridcolor = "#65866b", xgridstyle=:dash, ygridstyle=:dash)
    barplot!(ax, 1:lentime, values(MarketData.ohlc.Open), fillto = values(MarketData.ohlc.Close),
        color = colors,strokewidth = 0.5, strokecolor = colors, colormap = cmap)
    linesegments!(ax, linesegs, color = colors, colormap = cmap)
    CairoMakie.xlims!(ax, 200, 300)
    CairoMakie.ylims!(ax, 13,24)
    ax.xticks = (slice_dates, tempo[slice_dates])
    ax.xticklabelrotation = π/4
    fig[1,1] = ax
    fig
end

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
using DrWatson
@quickactivate "GCN_HM_GRN-Integration"

exp = "inner_batch"

using DataFrames, Latexify
df = collect_results(datadir(exp));

using Plots
theme(:juno)
gr()

plts = []
for i ∈ 1:size(df, 1)
	plt = plot([df.auc_list[i], df.losses_train[i], df.losses_test[i]],
	        labels = ["auc" "train" "test"], legend = :right,
	        xlabel = "Epochs", title = "Training for Cell Line $(df.cl[i]) with batch = $(df.batch[i]) and $(df.layer[i]) layers (Best AUC: $(DrWatson.roundval(df.auc[i]; digits = 4, scientific = 4)))", size = (800, 500), dpi = 100)
	push!(plts, plt)
end
p = plot(plts..., layout=(4, 3), size = (3200, 2000), dpi = 200)
savefig(p, plotsdir(exp, "summary.png"))


text:
"""
    rr_vle_vapor_fraction(K,z)

Given a vector of K values and a vector of compositions, calculates the vapor fraction `β`.
the algorithm is a modification of _(1)_ , with safeguards for extreme cases.

If the algorithm fails to converge, returns `NaN`. if it converges to a value `β ∉ [0,1]`, returns `-Inf` or `Inf`, depending on the case.

1. Vassilis Gaganis, "Solution of the Rachford Rice equation using perturbation analysis", Fluid Phase Equilibria, Volume 536,2021,112981
"""
function rr_vle_vapor_fraction(K,z)
    #Solution of the Rachford Rice equation using perturbation analysis
    #https://doi.org/10.1016/j.fluid.2021.112981
    #strategy: find 4 strongest compositions
    #calculate exact β for those 4
    #real_β = β_ℵ1 + ε*β_ℵ2(ε), where ℵ1 is the strong group of component and
    #ℵ2 is the weak group of components.
    #on ε = 0, we have the aprox solution
    #on ε = 1 we have the real solution
    #ε*β_ℵ2(ε) ≈ ∑ci*ε^i
    #we obtain the coefficients and evaluate ε = 1.
    #for n = length(z) <= 4, we return exact solutions.
    n_z = length(z)
    _0 = zero(first(z)+first(K))
    _1 = one(_0)
    kmin = minimum(K)
    kmax = maximum(K)
    βmax = 1/(1-kmin)
    βmin = 1/(1-kmax)

    kmax < 1 && return (-_1/_0)
    kmin > 1 && return (_1/_0)
    if n_z <= 4
        return rr_vle_vapor_fraction_exact(K,z)
    end

    #z_ℵ1 is normalized
    ℵ1,k_ℵ1,z_ℵ1 = Clapeyron.rr_find_strongest(K,z)
    β_near0,β_near1 = extrema(((_1/(_1 -k_ℵ1[1])),(_1/(_1 -k_ℵ1[2]))))
    near_mean = sqrt(abs(β_near0)*abs(_1-β_near1))

    β0 = rr_vle_vapor_fraction_exact(k_ℵ1,z_ℵ1)
    b0 = β0
    M1,M2 = _0,_0
    N1,N2 = _0,_0
    Ξ1,Ξ2 = _0,_0
    P1 = _0
    Λ2 = _0
    for i in 1:length(z)
        βi = _1/(_1-K[i])
        γi = βi - b0
        γiinv = _1/γi
        λi = z[i]*γiinv
        μi = λi*γiinv
        νi = μi*γiinv
        ξi = νi*γiinv
        if i in ℵ1
            M1 += μi
            N1 += νi
            Ξ1 += ξi
            P1 += ξi*γiinv
        else
            M2 += μi
            N2 += νi
            Ξ2 += ξi
            Λ2 += λi
        end
    end
    b1 = -Λ2/M1
    b2 = -(b1*b1*N1 + b1*M2)/M1
    b11 = b1*b1
    b111 = b11*b1
    b12 = b1*b2
    b22 = b2*b2
    b3 = -((b111*Ξ1+2*b12*N1) + (b11*N2+b2*M2))/M1
    t1 = b11*b11*P1 + 3*b11*b2*Ξ1 + (b22+2*b1*b3)*N1
    t2 = b111*Ξ2 + 2*b12*N2 + b3*M2
    b4 = -(t1+t2)/M1
    βaprox = (b0+b1+b2+b3+b4) #β(ε=1)
    βsol1 =  Solvers.ad_newton(x->rr_flash_eval(K,z,x),βaprox)

    #converged to -2.2e-16 or lower, practically 0
    if abs(βsol1) < eps(_1)
        βsol1 =_0
    end
    if 0 <= βsol1 <= 1 #converged as expected inside range [0,1]
        βsol = βsol1
    elseif iszero(near_mean) #one βi is exacly zero or exactly one
        (rr_flash_eval(K,z,_1) < eps(_1)) && (βsol = _1)
        (rr_flash_eval(K,z,_0) < eps(_1)) && (βsol = _0)
    elseif (βsol1 < β_near0) | (βsol1 > β_near1)     #in this case, we have two asymptotes really, really near one and zero
        βaprox_near0 = near_mean
        βaprox_near1 = _1 - near_mean
        βsol_near0 =  Solvers.ad_newton(x->rr_flash_eval(K,z,x),βaprox_near0)
        βsol_near1 =  Solvers.ad_newton(x->rr_flash_eval(K,z,x),βaprox_near1)
        βsol_near0,βsol_near1
        if 0 <= βsol_near0 <= 1
            βsol = βsol_near0
        elseif 0 <= βsol_near1 <= 1
            βsol = βsol_near1
        end
    else
        βsol = _0/_0
    end
    return βsol

end
function rr_find_strongest(K,z)
    _0 = zero(first(z)+first(K))
    _1 = one(_0)
    (kmin,idmin) =  findmin(K)
    (kmax,idmax) =  findmax(K)
    βmin = _1/(_1-kmax)
    βmax = _1/(_1-kmin)
    idx = 0
    for i in 1:length(z)
        idx +=1
        Ki = K[i]
        idx in (idmin,idmax) && continue
        βi = _1/(_1-Ki)
        if βmin <= βi <= 0
            idmin = idx
            βmin = βi
        elseif 1 <= βi <=βmax
            idmax = idx
            βmax = βi
        end
    end

    #(4, 2, 1, 5)
    id1 = 0
    id2 = 0

    idx = 0
    strong_zβ = (-_1,-_1)
    #find the strongest components, by zi/|βi|, ignoring βmin and βmax
    for i in 1:length(z)
        Ki = K[i]
        zi = z[i]
        idx+=1
        idx in (idmin,idmax) && continue
        zβi =abs(zi*(_1-Ki))
        zβ1,zβ2 = strong_zβ
        if zβ2 <= zβi
            strong_zβ = (zβ2,zβi)
            id1 = id2
            id2 = idx
        elseif zβ1 <= zβi
            strong_zβ = (zβi,zβ2)
            id1 = idx
        end
    end

    #indices of the strongest components
    idxs = (idmin,idmax,id1,id2)
    sumz = z[idmin] + z[idmax] +z[id1] +z[id2]
    invsumz = _1/sumz
    #normalize z
    zs = (invsumz*z[idmin],invsumz*z[idmax],invsumz*z[id1],invsumz*z[id2])
    ks = (K[idmin],K[idmax],K[id1],K[id2])
    return idxs,ks,zs
    end

    function rr_vle_vapor_fraction_exact(K,z)
    #if this function is called, then we are sure that there is a solution in the interval [0,1]
    _0 = zero(first(K)+first(z))
    n = length(z)
    _1 = one(_0)
    if  n == 2
        z1,z2 = z
        k1,k2 = K
        b1,b2 = 1/(1-k1),1/(1-k2)
        a1 = (z1*b2 + z2*b1)
        a0 = z1+z2
        return a1/a0
    elseif n == 3
        #0 = a0 + a1β + a2β^2
        z1,z2,z3 = z
        k1,k2,k3 = K
        b1,b2,b3 = 1/(1-k1),1/(1-k2),1/(1-k3)
        a2 =(z1 + z2 + z3)
        a1 = -b1*(z2 + z3) - b2*(z1 + z3) - b3*(z1 + z2)
        a0 = b1*b2*z3 + b1*b3*z2 + b2*b3*z1
        Δ = a1*a1 - 4*a0*a2
        inva2 = 1/(2*a2)
        Δ2 = sqrt(Δ)*inva2
        x = -a1*inva2
        β1 = x-Δ2
        β2 = x+Δ2
        βmax = max(β1,β2)
        βmin = min(β1,β2)
        if 0 < β1 < 1
            return β1
        elseif 0 < β2 < 1
            return β2
        elseif isfinite(β1+β2)
            kmin = min(k1,k2,k3)
            kmax = max(k1,k2,k3)
            kmin > 1 && return βmax
            kmax < 1 && return βmin
        else
            return _0/_0
        end
    elseif n == 4
        #0 = a0 + a1β + a2β^2 + a3β^3
        z1,z2,z3,z4 = z
        k1,k2,k3,k4 = K
        b1,b2,b3,b4 = 1/(1-k1),1/(1-k2),1/(1-k3),1/(1-k4)
        a3 =(z1 + z2 + z3 + z4)
        a2 = -b1*(z2 + z3 + z4) - b2*(z1 + z3 + z4) - b3*(z1 + z2 + z4)
        a1 = b1*b2*(z3+z4) + b1*b3*(z2+z4) + b2*b3*(z1+z4) + b1*b4*(z2+z3) + b2*b4*(z1+z3) + b3*b4*(z1+z2)
        a0 = -b1*b2*b3*z4 - b1*b2*b4*z3 - b1*b3*b4*z2 - b2*b3*b4*z1
        res1 =  Solvers.roots3(a0,a1,a2,a3)
        βmax = max(b1,b2,b3,b4)
        βmin = min(b1,b2,b3,b4)
        clx1,clx2,clx3 = res1
        r1,r2,r3 = real(clx1),real(clx2),real(clx3)
        rsum = r1+r2+r3
        rmax,rmin = extrema((r1,r2,r3))
        rmid = rsum - rmax - rmin
        return rmid
    else
        return _0/_0
    end
    end

function rr_flash_eval(K,z,β,normalize=true)
        _0 = zero(first(z)+first(K)+first(β))
        _1 = one(_0)
        if normalize
            sumz = sum(z)
            invsumz = _1/sumz
        else
            invsumz = _1
        end
        res = _0
        for i in 1:length(z)
            Kim1 = K[i] - _1
            res += invsumz*z[i]*Kim1/(1+β*Kim1)
        end
    return res
end

"""
    rr_flash_vapor(k,z,β)
Returns the gas phase composition, given k-values `k`, the initial molar composition `z` and the molar vapor fraction ´β´.

Each gas phase composition is calculated acording to:

    xvᵢ = kᵢzᵢ/(1+ β(kᵢ-1))
"""
function rr_flash_vapor(k,z,β)
    function f(ki,zi)
    _1 = one(ki)
        return ki*zi/(_1+β*(ki-_1))
    end

    return map(f,k,z)
end

"""
    rr_flash_liquid(k,z,β)
Returns the liquid phase composition, given k-values `k`, the initial molar composition `z` and the molar vapor fraction ´β´.

Each gas phase composition is calculated acording to:

    xlᵢ = zᵢ/(1+ β(kᵢ-1))
"""
function rr_flash_liquid(k,z,β)
    function f(ki,zi)
        _1 = one(ki)
        return zi/(_1+β*(ki-_1))
    end
    return map(f,k,z)
end


text:
module IonCLI

using Comonicon
using IonBase

const CASTED_COMMANDS = IonBase.CASTED_COMMANDS

"The Ion manager."
@main

end


text:
using Test
using BetweenFlags

code_dir = joinpath(dat_dir, "julia", "simple_func")

"""
    post_process

Remove end of line ("\r\n") (Windows)
Remove carriage return ("\r") Mac OS (before X)
"""
function post_process(s)
  s = replace(s, "\r\n" => "\n")
  s = replace(s, "\r" => "\n")
  return s
end

@testset "julia/simple_func" begin

  filename = joinpath(code_dir, "input.jl")
  code = open(f->read(f, String), filename)

  flag_set = FlagSet([
  FlagPair{ScopeType}(
    StartFlag("function", ["\n",""], [" "]),
    StopFlag( "end",      ["\n","\r"], ["\n","\r"];)
  ),
  FlagPair{ScopeType}(
    StartFlag("if",       ["\n"], [" "]),
    StopFlag( "end",      ["\n","\r"], ["\n","\r"])
  ),
  FlagPair{ScopeType}(
    StartFlag("for",      ["\n"], [" "]),
    StopFlag( "end",      ["\n","\r"], ["\n","\r"])
  )])

  token_stream = TokenStream(code, flag_set)
  export_results && export_plot(token_stream,code; path=code_dir)

  ####
  #### for-end
  ####
  section = token_stream("for-end")
  section_expected = open(f->read(f, String), joinpath(code_dir, "expected_for.jl"))

  # TODO: account for return carriage
  @test post_process(section)==post_process(section_expected)

  # To remove:
  # filename = joinpath(code_dir, "output_for.jl")
  # open(filename,"w") do io
  #   print(io, section)
  # end

  ####
  #### if-end
  ####
  section = token_stream("if-end")
  section_expected = open(f->read(f, String), joinpath(code_dir, "expected_if.jl"))
  # # TODO: account for return carriage
  @test post_process(section)==post_process(section_expected)

  # To remove:
  # filename = joinpath(code_dir, "output_if.jl")
  # open(filename,"w") do io
  #   print(io, section)
  # end

end

@testset "Julia: if/function code block" begin

  text = "Foo
if cond
    function myfunc()
        more stuff
        if cond
            print('foobar')
        else
            print('foobar')
        end
        for cond2
            print('foobar')
        else
            print('foobar')
        end
        foobaz
    end
end
baz"

  text_expected = "
if cond
    function myfunc()
        more stuff
        if cond
            print('foobar')
        else
            print('foobar')
        end
        for cond2
            print('foobar')
        else
            print('foobar')
        end
        foobaz
    end
end
"
  flag_set = FlagSet([
  FlagPair{ScopeType}(
    StartFlag("function", ["\n"," "], [" "]),
    StopFlag( "end",      ["\n","\r"], ["\n","\r"])
  ),
  FlagPair{ScopeType}(
    StartFlag("if",       ["\n", " "], [" "]),
    StopFlag( "end",      ["\n","\r", " "], ["\n","\r"])
  ),
  FlagPair{ScopeType}(
    StartFlag("for",      ["\n", " "], [" "]),
    StopFlag( "end",      ["\n","\r", " "], ["\n","\r"])
  )])

  token_stream = TokenStream(text, flag_set)

  @test token_stream("if-end") == text_expected
  export_results && export_plot(token_stream, text; path=".")

end


text:
using GeneralizedGenerated
using TupleVectors: chainvec
import MeasureTheory: testvalue

export testvalue
EmptyNTtype = NamedTuple{(),Tuple{}} where T<:Tuple

# function testvalue(d::ConditionalModel, N::Int)
#     r = chainvec(testvalue(d), N)
#     for j in 2:N
#         @inbounds r[j] = testvalue(d)
#     end
#     return r
# end

# testvalue(d::ConditionalModel, N::Int) = testvalue(d, N)

@inline function testvalue(c::ConditionalModel)
    m = Model(c)
    return _testvalue(getmoduletypencoding(m), m, argvals(c))
end

@inline function testvalue(m::Model)
    return _testvalue(getmoduletypencoding(m), m, NamedTuple())
end

sourceTestvalue(m::Model) = sourceTestvalue()(m)
sourceTestvalue(jd::ConditionalModel) = sourceTestvalue(jd.model)

export sourceTestvalue
function sourceTestvalue()
    function(_m::Model)
        proc(_m, st::Assign)  = :($(st.x) = $(st.rhs))
        proc(_m, st::Sample)  = :($(st.x) = testvalue($(st.rhs)))
        proc(_m, st::Return)  = :(return $(st.rhs))
        proc(_m, st::LineNumber) = nothing

        vals = map(x -> Expr(:(=), x,x),parameters(_m))

        wrap(kernel) = @q begin
            $kernel
            $(Expr(:tuple, vals...))
        end

        buildSource(_m, proc, wrap) |> MacroTools.flatten
    end
end

@gg function _testvalue(M::Type{<:TypeLevel}, _m::Model, _args)
    body = type2model(_m) |> sourceTestvalue() |> loadvals(_args, NamedTuple())
    @under_global from_type(_unwrap_type(M)) @q let M
        $body
    end
end

@gg function _testvalue(M::Type{<:TypeLevel}, _m::Model, _args::NamedTuple{()})
    body = type2model(_m) |> sourceTestvalue()
    @under_global from_type(_unwrap_type(M)) @q let M
        $body
    end
end


text:
using LinearAlgebra
using GenericSchur

include("interface.jl")
# abstract type AbstractDiffEqOperator{T} end
# abstract type AbstractDiffEqLinearOperator{T} <: AbstractDiffEqOperator{T} end

struct LinearHOODEOperator{T} <: DiffEqBase.AbstractDiffEqLinearOperator{T}
    epsilon::T
    A::AbstractMatrix
end
DiffEqBase.isinplace(linop::LinearHOODEOperator, n::Int) = false

function LinearHOODEOperator(mat::AbstractMatrix{T}) where {T}
    # tab = zeros(Int,size(mat,1))
    # sz = 0
    # testlgn(x)=count(y->y!=0,x)!=0
    # for i=1:size(mat,1)
    #     if testlgn(mat[i,:]) && testlgn(mat[:,i])
    #         sz += 1
    #         tab[sz] = i
    #     end
    # end
    # newmat = zeros(T,sz,sz)
    # for i=1:sz, j=1:sz
    #     newmat[i,j] = mat[tab[i],tab[j]]
    # end
    P = eigvecs(mat .+ 0im) # ca bug avec les reels
    Ad = inv(P) * mat * P

    #    epsilon = 1/norm(Ad, Inf) ca bug avec Double64
    epsilon = 1 / maximum(abs.(Ad))
    A = epsilon * mat
    Aint = round.(A)
    @assert norm(A - Aint, Inf) <= eps(T) "The A matrix must be with integers A=$A"
    epsilon = norm(Aint, Inf) / norm(mat, Inf)
    LinearHOODEOperator{T}(epsilon, Aint)
end

*(v::T, L::LinearHOODEOperator{T}) where {T<:Number} =
    LinearHOODEOperator(L.epsilon / v, L.A)
Base.@propagate_inbounds Base.convert(::Type{AbstractMatrix}, L::LinearHOODEOperator) =
    (1 / L.epsilon) * L.A
isconstant(_::LinearHOODEOperator) = true

function LinearHOODEOperator(linop::DiffEqArrayOperator)
    linop.update_func != DiffEqBase.DEFAULT_UPDATE_FUNC &&
        error("no update operator function for HOODEAB Alg")
    LinearHOODEOperator(linop.A)
end

LinearHOODEOperator(linop::LinearHOODEOperator) = linop

function LinearHOODEOperator(
    odefct::ODEFunction{iip,LinOp},
) where {iip,LinOp<:DiffEqBase.AbstractDiffEqLinearOperator}
    return LinearHOODEOperator(odefct.f)
end

isSplitODEProblem(probtype::Any) = false
isSplitODEProblem(probtype::SplitODEProblem) = true

"""
    HOODEAB( order=4, ntau=32)

Algorithm for High-Oscillatory equation
"""
struct HOODEAB{order,ntau} <: DiffEqBase.AbstractODEAlgorithm
    HOODEAB(order::Int = 4; ntau = 32) = new{order,ntau}()
end
export HOODEAB

# OtherHOODE : structure used to store HOODESolution data in the dstats fields
struct OtherHOODE
    par_u0::PrepareU0
    p_coef::CoefExpAB
    absprec::Union{Nothing,Float64}
    relprec::Union{Nothing,Float64}
end

"""
    solve(prob::ODEProblem, alg::HOODEAB{order, ntau};
    dt=nothing,
    kwargs...
    ) where {order,ntau}

common interface solver for Highly oscillatory problems, that an ODE of this form

```math
\\frac{\\delta u(t)}{\\delta t} = \\frac{1}{\\varepsilon} A + F(u(t), t)
```

where ``u \\in \\R^n`` and  ``0 < \\varepsilon < 1``
``A`` must be a **periodic matrix** i.e. ``e^{t A} = e^{(t+\\pi) A}`` for any ``t \\in \\R``

## Argument :
- `prob::ODEProblem` : The problem to solve
- `alg::HOODEAB{order, ntau}` : the Adams-Bashforth HOODE algorithm

## Keywords :
- `dt` : duration of a time interval
- `kwargs...` : other keywords

"""
function DiffEqBase.solve(
    prob::ODEProblem,
    alg::HOODEAB{order,ntau};
    dt = nothing,
    kwargs...,
) where {order,ntau}
    isSplitODEProblem(prob.problem_type) || error("HOODEAB alg need SplitODEProblem type")
    (!isnothing(dt) && haskey(kwargs, :nb_t)) &&
        error("Only one of dt and nb_t must be defined")
    haskey(kwargs, :order) &&
        error("order must defined as parameter of algorithm : HOODEAB(order)")
    linop = LinearHOODEOperator(prob.f.f1)
    fct = prob.f.f2.f
    p = typeof(prob.p) == DiffEqBase.NullParameters ? missing : prob.p
    ho_prob = if haskey(prob.kwargs, :B)
        HOODEProblem(fct, prob.u0, prob.tspan, p, linop.A, linop.epsilon, prob.kwargs[:B])
    else
        HOODEProblem(fct, prob.u0, prob.tspan, p, linop.A, linop.epsilon)
    end
    if !haskey(kwargs, :nb_t) && !isnothing(dt)
        kwargs = (kwargs..., nb_t = Int64(round((prob.tspan[2] - prob.tspan[1]) / dt)))
    end

    sol = DiffEqBase.solve(
        ho_prob,
        HOODETwoScalesAB();
        nb_tau = ntau,
        order = order,
        kwargs...,
    )
    other = OtherHOODE(sol.par_u0, sol.p_coef, sol.absprec, sol.relprec)
    return DiffEqBase.build_solution(
        prob, # prob::Union{AbstractODEProblem,AbstractDDEProblem},
        HOODETwoScalesAB(), # alg,
        sol.t, # t,
        sol.u, # u,
        dense = sol.dense,
        interp = sol.interp,
        retcode = sol.retcode,
        destats = other,
    )
end


text:
using LinearAlgebra
using Plots
using IJulia

#Where to put these files?
include("cheb.jl")
include("parameters.jl")
include("geometry.jl")
include("profile.jl")
include("build_matrices.jl")
include("compute_spectrum.jl")
include("mesh.jl")
include("plot_fields.jl")

params=parameters(
     H = 1,
    Ly = 20,
    f₀ = 1,
     g = 10,
    Lj = 1,
    Uj = 0.1,
    Ny = 256,
    dk = 5e-2,
    kₘ = 2,
    );

params.Fr = (params.f₀ * params.Lj)^2/(params.g * params.H);

# Wavenumbers
ks = collect(params.dk:params.dk:params.kₘ) / params.Lj;
Nk = length(ks);

# Add FD option
# Grid and basic state
         Ny = params.Ny;
Dy,  Dy2, y = geometry(Ny, params.Ly, "cheb");
dQdy,  U, Ψ = profile(  y, params);

plot_basic_state(y, dQdy, U, Ψ, params.Ly, "basic_state.png")

# initialize fields to store
Nmodes = 2;
     σ = zeros(Nmodes, Nk);
     ω = zeros(Nmodes, Nk);
σmodes = zeros(ComplexF64, Ny+1, Nmodes, Nk);

# Do we want a structure with soln.U etc?
# Compute growth rates
for cnt in 1:length(ks)
    local k = ks[cnt]
    local A, B = build_AB(k, dQdy, U, Ψ, Dy, Dy2, y, params);
    local σ[:,cnt], ω[:,cnt], σmodes[2:end-1,:,cnt] = compute_spectrum(A, B, k, params, Nmodes);
end

plot_growth_rates(ks, σ, Nmodes, "growth_rates.png")

mode_number = 1;
      σ_max = maximum(σ[mode_number,:]);
    k_index = sortperm(σ[mode_number,:],rev=true)[1];
          k = ks[k_index]; # pick wavenumber

plot_1D_streamfunction(k_index, k, Ny, y, σmodes, mode_number, "modes_1D_streamfunction.png")
plot_2D_streamfunction(k_index, k, Ny, y, σmodes, mode_number, "modes_2D_streamfunction.png")
plot_2D_vorticity(k_index, k, Ny, y, Dy2, σmodes, mode_number, "modes_2D_vorticity.png")


text:
using BenchmarkTools
SUITE = BenchmarkGroup()
SUITE["filter_map_reduce"] = include("bench_filter_map_reduce.jl")
SUITE["filter_map_map!"] = include("bench_filter_map_map!.jl")
SUITE["partition_by"] = include("bench_partition_by.jl")


text:
using Test
using JotTest2

input = Dict("double" => 0.1)
input_doubled = JotTest2.response_func(input)
@test isa(input_doubled, Float64)
@test input_doubled == 0.2

input = Dict("add suffix" => "suffixme_")
input_suffixed = JotTest2.response_func(input)
@test isa(input_suffixed, String)


text:
using Test
using MIPVerify
using MIPVerify: LInfNormBoundedPerturbationFamily
@isdefined(TestHelpers) || include("../TestHelpers.jl")

@timed_testset "integration.jl" begin
    mnist = read_datasets("MNIST")
    nn_wk17a = get_example_network_params("MNIST.WK17a_linf0.1_authors")
    @testset "batch_find_untargeted_attack" begin
        mktempdir() do dir
            MIPVerify.batch_find_untargeted_attack(
                nn_wk17a,
                mnist.test,
                [1], # robust sample
                TestHelpers.get_main_solver(),
                solve_rerun_option = MIPVerify.never,
                pp = MIPVerify.LInfNormBoundedPerturbationFamily(0.1),
                norm_order = Inf,
                rebuild = true,
                tightening_algorithm = lp,
                tightening_solver = TestHelpers.get_tightening_solver(),
                cache_model = false,
                solve_if_predicted_in_targeted = false,
                save_path = dir,
            )
        end

        mktempdir() do dir
            MIPVerify.batch_find_untargeted_attack(
                nn_wk17a,
                mnist.test,
                [9, 248], # non-robust and misclassified sample
                TestHelpers.get_main_solver(),
                solve_rerun_option = MIPVerify.never,
                pp = MIPVerify.LInfNormBoundedPerturbationFamily(0.1),
                norm_order = Inf,
                rebuild = true,
                tightening_algorithm = interval_arithmetic,
                cache_model = false,
                solve_if_predicted_in_targeted = false,
                save_path = dir,
            )
        end
    end

    @testset "batch_find_targeted_attack" begin
        mktempdir() do dir
            MIPVerify.batch_find_targeted_attack(
                nn_wk17a,
                mnist.test,
                [1],
                TestHelpers.get_main_solver(),
                solve_rerun_option = MIPVerify.never,
                pp = MIPVerify.LInfNormBoundedPerturbationFamily(0.1),
                norm_order = Inf,
                tightening_algorithm = interval_arithmetic,
                tightening_solver = TestHelpers.get_tightening_solver(),
                cache_model = false,
                solve_if_predicted_in_targeted = false,
                target_labels = [1, 8],
                save_path = dir,
            )
        end
    end

    @testset "batch_build_model" begin
        mktempdir() do dir
            MIPVerify.batch_build_model(
                nn_wk17a,
                mnist.test,
                [1],
                TestHelpers.get_tightening_solver(),
                pp = MIPVerify.LInfNormBoundedPerturbationFamily(0.1),
                tightening_algorithm = interval_arithmetic,
            )
        end
    end
end


text:
using MLDatasets
using ImageView
using CuArrays
# using VariationalAE

function show_pic(x)
    reshape(x, 28, 28)' |> imshow
end

function show_recon(m, x)
    x_ = reconstruct(m, x)
    show_pic(x)
    show_pic(x_)
end

function run()
    m = VAE{Float32}(784, 500, 500, 20, 500, 500, 784) |> to_cuda

    X, _ = MNIST.traindata()
    X = convert(Matrix{Float32}, reshape(X, 784, 60000))
    @time m = fit!(m, X; n_epochs=2, cuda=true)

    # check reconstructed image
    for i=1:2:10
        show_recon(m, X[:, i])
    end

    # check learned features
    for i=1:2:10
        generate(m, i) |> show_pic
    end
end

function main()
    ex = quote
        (mu, log_sigma2) = encode(m, x)
        z = mu .+ sqrt.(exp.(log_sigma2)) .* eps
        x_rec = decode(m, z)
        rec_loss = -(sum(x .* log.(1.0e-10 + x_rec) + (1 - x) .* log.((1.0e-10 + 1) - x_rec), 1))
        KLD = -0.5 * sum(((1 + log_sigma2) .- mu .^ 2) - exp.(log_sigma2), 1)
        cost = mean(rec_loss .+ KLD)
    end
    m = VAE{Float32}(784, 500, 500, 20, 500, 500, 784) |> to_cuda
    x = cu(rand(Float32, 784, 100))
    inputs = [:m => m, :eps => cu(rand(Float32, 20, 100)), :x => x]
    g = ExGraph(ex; inputs...)

    encode(m, x)
end


text:
using
    Test,
    Distributions,
    ImageQualityIndexes,
    Random,
    Plots,
    Statistics,
    Wavelets,
    WaveletsExt

@testset "Utils" begin include("utils.jl") end
@testset "Visualizations" begin include("visualizations.jl") end
@testset "Transforms" begin include("transforms.jl") end
@testset "Denoising" begin include("denoising.jl") end
@testset "Best Basis" begin include("bestbasis.jl") end
@testset "LDB" begin include("ldb.jl") end