text:
abstract AbstractCTree
typealias IntArrays Array{Array{Int, 1}, 1}

################
## CTreePaths ##
################

@doc doc"""
CTree representation as array of individual paths to leaf nodes. Each
possible path to a final node is stored as an individual `Array{Int,
1}`, while the root node is stored separately in `root`. Paths
generally have different lengths, so that they are not stored as
`Array{Int, 2}`.

- `root`::Int

- `paths`::Array{Array{Int, 1}, 1}

Paths will be stored in lexicographical order.

## Constructor

````julia
CTreePaths(3, [1, 2], [1, 4], [5], [6, 7])
````
"""->
type CTreePaths <: AbstractCTree
    root::Int
    paths::IntArrays

    function CTreePaths(rt::Int, paths::IntArrays)
        ## constraints: paths need to be sorted lexicographically
        pathsSorted = sortPaths(paths)
        return new(rt, pathsSorted)
    end
end

##################
## constructors ##
##################

## construction with individual paths as variable arguments
function CTreePaths(rt::Int, paths...)
    ## paths need to be Array{Int, 1} or Array{Int, 2}
    nPaths = length(paths)
    ps = Array(Array{Int, 1}, nPaths)
    for ii=1:nPaths
        ps[ii] = paths[ii][:]
    end
    return CTreePaths(rt, ps)
end

## dg function
##------------

@doc doc"""
Delegate functions to single field of composite type.
"""->
function dg(rt::CTreePaths)
    return rt.paths
end

####################
## CondTreeParRef ##
####################
@doc doc"""
CTree representation in parent reference notation. `tree` is an
`Array{Int, 1}`, where entry ii denotes the single parent of node ii.
The root node is indicated with parent 0, while not yet connected
nodes have parent -1.
"""->
type CTreeParRef <: AbstractCTree
    tree::Array{Int, 1}
end

## dg function
##------------

@doc doc"""
Delegate functions to single field of composite type.
"""->
function dg(rt::CTreeParRef)
    return rt.tree
end

#################
## conversions ##
#################

import Base.convert

## convert to CTreePaths
##----------------------

@doc doc"""
Convert parent referencing notation to tree paths notation.

The algorithm starts by finding all final nodes as those that are not
listed as parent node for any other variable. For the case of
incomplete trees, nodes not yet connected are encoded with a value of
-1.
"""->
function convert(::Type{CTreePaths}, tr::CTreeParRef)
    parNot = tr.tree

    nVars = length(parNot)

    ## find notes without children
    childlessNodes = setdiff([1:nVars], parNot)

    ## exclude nodes that do not appear at all
    notAppearingNodes = find(parNot .== -1)

    endNodes = setdiff(childlessNodes, notAppearingNodes)

    nPaths = length(endNodes)

    ## preallocate array of arrays
    treePaths = Array(Array{Int, 1}, nPaths)

    rootNode = -1
    for ii=1:nPaths
        currPath = Array(Int, 1)
        currPath[1] = endNodes[ii]
        belowRoot = true

        ## walk along path
        while belowRoot
            parentNode = parNot[currPath[end]]
            if parentNode == 0
                belowRoot = false
            else
                append!(currPath, [parentNode])
            end
        end
        rootNode = currPath[end]
        treePaths[ii] = flipud(currPath)[2:end]
    end

    return CTreePaths(rootNode, treePaths)
end

function convert(::Type{CTreePaths}, tr::CTreePaths)
    return tr
end

function convert(::Type{CTreePaths}, tr::Array{Int, 1})
    trPar = CTreeParRef(tr)
    return convert(CTreePaths, trPar)
end

function convert(::Type{CTreeParRef}, tr::Array{Int, 1})
    return CTreeParRef(tr)
end

## convert to CTreeParRef
##-----------------------

@doc doc"""
Transform CTree given as list of paths to final nodes to parent
referencing notation. For parent referencing, generally all variables
from 1 to N should occur in the tree. However, even unfinished trees
are supported, and missing variables are denoted with a parent of -1.
To determine the last variable `convert` simply takes the maximum
value of occurring nodes in all paths.
"""->
function convert(::Type{CTreeParRef}, tr::CTreePaths, nVars::Int)
    ## transform single tree to parent notation vector
    parNot = -1*ones(Int, nVars)
    parNot[tr.root] = 0

    nPaths = length(tr.paths)

    for ii=1:nPaths
        currPath = tr.paths[ii]
        nNodes = length(currPath)
        parNot[currPath[1]] = tr.root

        for jj=1:(nNodes-1)
            ## get parent for each node
            parNot[currPath[end+1-jj]] = currPath[end-jj]
        end
    end
    return CTreeParRef(parNot)
end

function convert(::Type{CTreeParRef}, tr::CTreePaths)
    nPaths = width(tr)
    nVars = maximum(Int[maximum(tr.paths[ii]) for ii=1:nPaths])
    nVars = maximum([tr.root, nVars])
    return convert(CTreeParRef, tr, nVars)
end

function convert(::Type{CTreeParRef}, tr::CTreeParRef)
    return tr
end

#######################
## display functions ##
#######################

import Base.Multimedia.display
@doc doc"""
All AbstractCTree subtypes will basically display identically, showing
all paths of the tree to the leaf nodes. However, they additionally
show their type first, so that different types still can be
distinguised.
"""->
function displayInner(trPaths::CTreePaths)
    for ii=1:width(trPaths)
        println(join(["   ", string(trPaths.paths[ii])]))
    end
    return Nothing
end

## display AbstractCTree
function display(tr::AbstractCTree)
    println("type: $(typeof(tr))")
    trPaths = convert(CTreePaths, tr)
    displayInner(trPaths)
end

## ## CTree arrays
## ##-------------
## function display(trArr::Array{AbstractCTree, 1})
##     ## convert to array of CTreePaths and display that
##     nTrees = length(trArr)
##     for ii=1:nTrees
##         println("type: $(typeof(trArr[ii]))")
##         currTree = convert(CTreePaths, trArr[ii])
##         println("$(currTree.root):")
##         displayInner(currTree)
##     end
## end

function display(trArr::Array{CTreePaths, 1})
    ## convert to array of CTreePaths and display that
    nTrees = length(trArr)
    for ii=1:nTrees
        currTree = trArr[ii]
        println("type: $(typeof(currTree))")
        println("$(currTree.root):")
        displayInner(currTree)
    end
end

function display(trArr::Array{CTreeParRef, 1})
    ## convert to array of CTreePaths and display that
    nTrees = length(trArr)
    for ii=1:nTrees
        println("type: $(typeof(trArr[ii]))")
        currTree = convert(CTreePaths, trArr[ii])
        println("$(currTree.root):")
        displayInner(currTree)
    end
end

#########################
## dimension functions ##
#########################

@doc doc"""
`Int` length of the longest path of an `Array` of paths.
"""->
function maxDepth(paths::IntArrays)
    maxD = 0
    for ii=1:length(paths)
        currLength = length(paths[ii])
        if currLength > maxD
            maxD = currLength
        end
    end
    return maxD
end

@doc doc"""
`Int` length of the longest path of a `CTree` object.
"""->
function maxDepth(tr::CTreePaths)
    return maxDepth(tr.paths)
end

function maxDepth(tr::AbstractCTree)
    trPaths = convert(CTreePaths, tr)
    return maxDepth(trPaths)
end

@doc doc"""
`Array{Int, 1}` with depths / lengths to each final node in a tree.
"""->
function getDepths(tr::CTreePaths)
    nPaths = width(tr)
    depths = Array(Int, nPaths)
    for ii=1:nPaths
        depths[ii] = length(tr.paths[ii])
    end
    return depths
end

function getDepths(tr::AbstractCTree)
    trPaths = convert(CTreePaths, tr)
    return getDepths(trPaths)
end

@doc doc"""
`Int` number of different paths / leaf nodes of a given CTree.
"""->
function width(tr::CTreePaths)
    return length(tr.paths)
end

function width(tr::AbstractCTree)
    trPaths = convert(CTreePaths, tr)
    return width(trPaths)
end

########################
## check for equality ##
########################
@doc doc"""
For equality to hold both CTrees must have equal representation
(type). Besides, root node and paths must be identical.
"""->
function ==(tr1::CTreePaths, tr2::CTreePaths)
    return (tr1.root == tr2.root) & (tr1.paths == tr2.paths)
end

function ==(tr1::AbstractCTree, tr2::AbstractCTree)
    if typeof(tr1) != typeof(tr2)
        return false
    else
        trPaths1 = convert(CTreePaths, tr1)
        trPaths2 = convert(CTreePaths, tr2)
        return ==(trPaths1, trPaths2)
    end
end

##############
## getindex ##
##############

import Base.getindex

function getindex(tr::CTreePaths, pathInd::Int, pathNode::Int)
    ## single entry
    return tr.paths[pathInd][pathNode]
end

function getindex(tr::CTreePaths, pathInd::Int)
    ## single path
    return tr.paths[pathInd]
end

function getindex(tr::CTreePaths, pathInds::Array{Int, 1},
                  pathNodes::Array{Int, 1})
    ## multiple values, sorted, only unique entries
    nVals = length(pathInds)
    if nVals != length(pathNodes)
        error("index arrays must be of equal length")
    end
    vals = Array(Int, nVals)
    for ii=1:nVals
        vals[ii] = tr[pathInds[ii], pathNodes[ii]]
    end
    return sort(unique(vals))
end

#############
## allVals ##
#############

@doc doc"""
`Array{Int, 1}` of all occurring values in a tree, including the root
node.
"""->
function allNodes(tr::CTreePaths)
    vals = [tr.root, [tr.paths...]]
    return sort(unique(vals))
end

function allNodes(tr::CTreeParRef)
    return find(tr.tree .>= 0)
end

## function allNodes(tr::AbstractCTree)
##     trPar = convert(CTreeParRef, tr)
##     return allNodes(trPar)
## end

@doc doc"""
`Array{Int, 1}` of all occurring values in a tree, excluding the root
node.
"""->
function allPathNodes(tr::CTreePaths)
    vals = [[tr.paths...]]
    return sort(unique(vals))
end

function allPathNodes(tr::CTreeParRef)
    return find(tr.tree .> 0)
end

## function allPathNodes(tr::AbstractCTree)
##     trPar = convert(CTreeParRef, tr)
##     return allPathNodes(trPar)
## end

## getRoot
##--------

function getRoot(tr::CTreeParRef)
    return find(dg(tr) .== 0)[1]
end

function getRoot(tr::CTreePaths)
    return tr.root
end

#############
## condSet ##
#############

@doc doc"""
`Bool` indicating whether a given conditional set `Array{Int, 1}`
occurs at the beginning of any path of a `CTreePaths` instance. Note that
the exact order of numbers within the conditional set does not matter.

## Example

tr = Copulas.CTreePaths(6, [1 2 3], [1 2 4], [5 3])
@test Copulas.condSetChk([2, 1], tr)
@test !(Copulas.condSetChk([2, 4], tr))

"""->
function condSetChk(arr::Array{Int, 1}, tr::CTreePaths)
    ## does conditional set occur at beginning of some path?
    nNodes = length(arr)

    nPaths = width(tr)
    for ii=1:nPaths
        if length(tr[ii]) >= nNodes
            if issubset(arr, tr[ii][1:nNodes])
                return true
            end
        end
    end
    return false
end

function condSetChk(arr::Array{Int, 1}, tr::AbstractCTree)
    trPaths = convert(CTreePaths, tr)
    return condSetChk(arr, trPaths)
end

## findAndSortCondSet
##-------------------

@doc doc"""
Find a given conditioning set in the tree of a variable given in
parent notation. If the conditioning set occurs, it will be returned
with the same sorting as it appears in the tree: the first entry is
directly connected to the root node. If the conditioning set is not
found, the function throws an error.

The algorithm first translates the tree in parent notation into a
`Tree`, which might be unnecessary and time consuming.

And alternative would be to search for the conditioning set directly
in parent notation. A conditioning set of length 3 could be found as a
sequence that arrives at the root node in the 4th step, with all steps
including nodes that are part of the conditioning set.
"""->
function findAndSortCondSet(tr::CTreeParRef,
                     condSet::Array{Int, 1})
    trPar = convert(CTreePaths, tr)

    nNodes = length(condSet)

    nPaths = width(trPar)
    for ii=1:nPaths
        if length(trPar[ii]) >= nNodes
            if issubset(condSet, trPar[ii][1:nNodes])
                return trPar[ii][1:nNodes]
            end
        end
    end
    error("Conditional set does not occur.")
end

## getPathToRoot
##--------------

@doc doc"""
For tree given in parent notation, get full path from variable `var`
to root node. Path does neither include `var` nor `root`. The first
element is connected to `root`, the last one to `var`.
"""->
function getPathToRoot(tr::CTreeParRef, var::Int)
    parNot = tr.tree
    rootNode = find(parNot .== 0)[1]
    notAtRoot = true
    path = Int[]
    currVar = var
    while notAtRoot
        currVar = parNot[currVar]
        if currVar == rootNode
            notAtRoot = false
        else
            push!(path, currVar)
        end
    end
    return flipud(path)
end

function getPathToRoot(tr::AbstractCTree, var::Int)
    trPar = convert(CTreeParRef, tr)
    return getPathToRoot(trPar, var)
end

#############
## attach! ##
#############

@doc doc"""
Attach new node to unfinished conditioning tree of type CTreePaths
with given conditioning set.
"""->
function attach!(tr::CTreePaths, newNode::Int, condSet::Array{Int, 1})
    nPaths = width(tr)
    nNodes = length(condSet)
    for ii=1:nPaths
        path = tr.paths[ii]
        if length(path) >= nNodes
            if issubset(condSet, path[1:nNodes])
                if length(path) == nNodes
                    ## attach to path
                    push!(path, newNode)
                    break
                else
                    ## new path
                    condSetSorted = path[1:nNodes]
                    newPath = [condSetSorted, newNode]
                    push!(tr.paths, newPath)
                    break
                end
            end
        end
    end
    ## resort arrays
    nPaths = length(tr.paths)
    pathsSorted = sortPaths(tr.paths)
    for ii=1:nPaths
        tr.paths[ii] = pathsSorted[ii]
    end
end

function attach!(tr::CTreePaths, newNode::Int, xx::Array{None,1})
    attach!(tr, newNode, Array(Int, 0))
end


text:
#######################################################################################
#                                   Load Packages
#######################################################################################
cd(@__DIR__)
using Pkg
Pkg.activate("../..")
using StatsPlots, Revise, ACTRModels, Distributions, Turing, DataFrames
include("Siegler_Model.jl")
include("../../utilities/plotting.jl")
Random.seed!(82210)
#######################################################################################
#                                   Generate Data
#######################################################################################
# mismatch penalty parameter
δ = 16.0
# retrieval threshold parameter
τ = -0.45
# logistic scalar for activation noise
s = 0.5
# fixed parameters
parms = (mmp = true,noise = true,mmp_fun = sim_fun,ter = 2.05)
# number of blocks
n_blocks = 5
# stimuli for simulation
stimuli = [(num1 = 1,num2 = 1), (num1 = 1,num2 = 2), (num1 = 1,num2 = 3), (num1 = 2,num2 = 2),
    (num1 = 2,num2 = 3), (num1 = 3,num2 = 3)]
# generate data for each block
temp = map(x -> simulate(stimuli, parms; δ, τ, s), 1:n_blocks)
# flatten data
data = vcat(vcat(temp...)...)
data[1]
#######################################################################################
#                                    Define Model
#######################################################################################
@model model(data, parms) = begin
    δ ~ Normal(16, 8)
    τ ~ Normal(-0.45, 1)
    s ~ truncated(Normal(0.5, 0.5), 0.0, Inf)
    data ~ Siegler(δ, τ, s, parms)
end
#######################################################################################
#                                 Estimate Parameters
#######################################################################################
# Settings of the NUTS sampler.
n_samples = 1000
delta = 0.85
n_adapt = 1000
n_chains = 4
specs = NUTS(n_adapt, delta)
# Start sampling.
chain = sample(model(data, parms), specs, MCMCThreads(), n_samples, n_chains, progress=true)
#######################################################################################
#                                  Posterior Correlations
#######################################################################################
cors = cor(chain)
#######################################################################################
#                                      Plot Chains
#######################################################################################
pyplot()
posteriors = plot(chain, seriestype=:density, grid=false, titlefont=font(10), layout=(3,1),
    size=(250,350), xaxis=font(8), yaxis=font(8))
plot!(title="", xlabel="δ", subplot=1)
plot!(title="", xlabel="τ", subplot=2)
plot!(title="", xlabel="s", subplot=3)
#######################################################################################
#                                  Posterior Predictive
#######################################################################################
preds = posterior_predictive(x -> simulate(stimuli, parms; x...), chain, 1000)
preds = vcat(vcat(preds...)...)
df = DataFrame(preds)
sort!(stimuli)
p4 = rt_histogram(df, stimuli; size = (300,300))
p5 = response_histogram(df, stimuli; size = (300,300))
savefig(p4, "Siegler_RT_Predictions.eps")
savefig(p5, "Siegler_Response_Predictions.eps")
savefig(posteriors, "Siegler_Posteriors.eps")

ext: jl


text:
# import Iterators.flatten
function textDocument_semanticTokens_request(params::SemanticTokensParams, server::LanguageServerInstance, conn)
    doc = getdocument(server, URI2(params.textDocument.uri))
    offset = get_offset(doc, params.position)
    identifier = get_identifier(getcst(doc), offset)
    identifier !== nothing || return nothing
    highlights = DocumentHighlight[]
    for_each_ref(identifier) do ref, doc1, o
        if doc1._uri == doc._uri
            kind = StaticLint.hasbinding(ref) ? DocumentHighlightKinds.Write : DocumentHighlightKinds.Read
            push!(highlights, DocumentHighlight(Range(doc, o .+ (0:ref.span)), kind))
        end
    end
    return isempty(highlights) ? nothing : highlights
end

struct SemanticToken
    # token line number, relative to the previous token
    deltaLine::UInt32
    # token start character, relative to the previous token
    # (relative to 0 or the previous token’s start if they are on the same line)
    deltaStart::UInt32
    # the length of the token.
    length::UInt32
    # will be looked up in SemanticTokensLegend.tokenTypes
    tokenType::UInt32
    # each set bit will be looked up in SemanticTokensLegend.tokenModifiers
    tokenModifiers::UInt32
end

function SemanticToken(deltaLine::UInt32,
                       deltaStart::UInt32,
                       length::UInt32,
                       tokenType::String,
                       tokenModifiers::String)
    # TODO look up int encodings for tokenType and tokenModifiers
    SemanticToken(
        deltaLine,
        deltaStart,
        length,
        semantic_token_encoding(tokenType),
        0 # TODO
    )
end

# function SemanticToken(ex::EXPR)

# end

function SemanticTokens(tokens::Vector{SemanticToken}) :: SemanticTokens
    token_vectors = map(tokens) do token
        # token_index = i - 1
        [
            token.deltaLine,
            token.deltaStart,
            token.length,
            token.tokenType,
            token.tokenModifiers
        ]
    end
    SemanticTokens(Iterators.flatten(token_vectors) |> collect)
end

function textDocument_semanticTokens_full_request(params::SemanticTokensParams,
                                                  server::LanguageServerInstance, conn) :: Union{SemanticTokens,Nothing}
    uri = params.textDocument.uri
    doc = getdocument(server, URI2(uri))
    # return nothing
    ts = collect_semantic_tokens(getcst(doc), doc)
    return SemanticTokens(ts)
    # return collect_document_symbols(getcst(doc), server, doc)

    # doc = getdocument(server, URI2(params.textDocument.uri))
    # offset = get_offset(doc, params.position)
    # identifier = get_identifier(getcst(doc), offset)
    # identifier !== nothing || return nothing
    # highlights = DocumentHighlight[]
    # for_each_ref(identifier) do ref, doc1, o
    #     if doc1._uri == doc._uri
    #         kind = StaticLint.hasbinding(ref) ? DocumentHighlightKinds.Write : DocumentHighlightKinds.Read
    #         push!(highlights, DocumentHighlight(Range(doc, o .+ (0:ref.span)), kind))
    #     end
    # end
    # return isempty(highlights) ? nothing : highlights
end
# import CSTParser:
"""

parse applies these types
 Parser.SyntaxNode
  Parser.EXPR
  Parser.INSTANCE
    Parser.HEAD{K}
    Parser.IDENTIFIER
    Parser.KEYWORD{K}
    Parser.LITERAL{K}
    Parser.OPERATOR{P,K,dot}
    Parser.PUNCTUATION
  Parser.QUOTENODE

┌ Info:   1:60  file(  new scope lint )
│   1:60   function(  Binding(main:: (1 refs)) new scope)
│   1:8     call( )
│   1:4      main *
│   9:48    block( )
│   9:28       1:2   OP: =( )
│   9:10      s Binding(s:: (3 refs)) *
│  11:26      STRING: hello world!
│  29:40     call( )
│  29:35      println *
│  36:36      s *
│  41:48     macrocall( )
│  41:46      @show *
│  47:46      NOTHING: nothing
└  47:48      s *
"""
function collect_semantic_tokens(expr::EXPR, doc :: Document, pos=0):: Vector{SemanticToken}
    tokens = SemanticToken[]

    offset = pos
    for ex in expr
        # if isempty(ex)
            # leaf of parse tree

		# C.isidentifier(ex)
		kind = semantic_token_kind(ex)
		if kind !== nothing
            @info "Expression" ex kind
            # identifier, identifier_pos = get_identifier_pos(ex, 0, 0)
            name = C.get_name(ex)
            name_offset = 0
            # get the offset of the name expr
            if name !== nothing
                found = false
                for x in ex
                    if x == name
                        found = true
                        break
                    end
                    name_offset += x.fullspan
                end
                if !found
                    name_offset = -1
                end
            end
            @info "name" name name_offset
			line, char = get_position_at(doc, offset)
			# add this token
			token = SemanticToken(
				line,
				char,
				ex.span,
				semantic_token_encoding(kind),
				0
			)
			push!(tokens, token)
		end
        if !isempty(ex)
            # there are more nodes on this tree
            sub_tokens = collect_semantic_tokens(ex, doc, offset)
            push!(tokens, sub_tokens...)
        end
        offset += ex.fullspan
    end
    return tokens
    # return collect_semantic_tokens(ex)
    # @info get_toks(doc, 0)
end

"""
Get the semantic token kind for `expr`, which is assumed to be an identifier
"""
function semantic_token_kind(expr::EXPR) :: Union{String, Nothing}
    # C.isidentifier(expr) || return nothing
    return if C.defines_function(expr)# || C.is_func_call(expr)
        SemanticTokenKinds.Function
    elseif C.defines_struct(expr)
        SemanticTokenKinds.Struct
    elseif C.defines_macro(expr)
        SemanticTokenKinds.Macro
    # elseif C.isoperator(expr)
        # SemanticTokenKinds.Operator
    end
end
const C = CSTParser
# function collect_semantic_tokens(ex::EXPR, pos=0) :: Vector{SemanticToken}
#     tokens = SemanticToken[]

#     @info "expr.head = $(ex.head)"
#     token_kind = if C.defines_function(ex) || C.is_func_call(ex)
#         SemanticTokenKinds.Function
#     elseif C.defines_struct(ex)
#         SemanticTokenKinds.Struct
#     elseif C.defines_macro(ex)
#         SemanticTokenKinds.Macro
#     elseif C.isoperator(ex)
#         SemanticTokenKinds.Operator
#     end
#     @info token_kind

#     if token_kind !== nothing
#         token = SemanticToken(
#             0,
#             pos,
#             ex.span,
#             semantic_token_encoding(token_kind),
#             0
#         )
#         push!(tokens, token)
#     end
#     pos = pos + span(tokens)

#     if ex.args !== nothing
         # sub_tokens = map(ex.args) do arg
         #     collect_semantic_tokens(arg, pos)
         # end |> Iterators.flatten |> collect
         # pos = pos + span(sub_tokens)
         # push!(tokens, sub_tokens...)
#     end

#     return tokens
# end

span(token :: SemanticToken) = token.length
function span(tokens :: Vector{SemanticToken})
    lengths = map(tokens) do token
        token.length
    end

    reduce(+, lengths)
end


text:
module Util

using Jld2Mat.Params

export printlog

const printlog = verbose ? println : (args...) -> nothing

end # module


text:
######################################################################
# Event emitter & listener system leveraging dynamic multiple dispatch.
# -----
# Licensed under MIT License
module Events
import ..ExtraFun
import ExtraFun: Ident, indexof

export @event
macro event(fn::Expr)
    @assert fn.head === :function || fn.head === :(=) "expected function declaration"
    @assert length(fn.args) == 2 "unexpected function data"

    fn_head, fn_body = fn.args

    @assert fn_head.head === :call "expected function declaration"

    curly  = fn_head.args[1]
    params = nothing
    args   = fn_head.args[2:length(fn_head.args)]
    if args[1].head === :parameters
        params = args[1]
        args = args[2:length(args)]
    end

    @assert curly.head === :curly "expected curly function with event name symbol"

    listener_name = curly.args[1]
    event_name    = curly.args[2].value::Symbol

    if !haskey(_static_events, event_name)
        _static_events[event_name] = Symbol[]
    end
    push!(_static_events[event_name], listener_name)

    result = :(Events.listen(::ExtraFun.Ident{$(QuoteNode(event_name))}, ::ExtraFun.Ident{$(QuoteNode(listener_name))}) = $fn_body)
    if params !== nothing
        insert!(result.args[1].args, 2, params)
    end
    append!(result.args[1].args, args)
    esc(result)
end
isdotsequence(::Symbol) = true
isdotsequence(expr::Expr) = expr.head === :. && isdotsequence(args[2])

export emit
function emit(evt_name::Symbol, args...; kwargs...)
    errors = Exception[]

    if haskey(_static_events, evt_name)
        for fn_name ∈ _static_events[evt_name]
            try
                listen(Ident{evt_name}(), Ident{fn_name}(), args...; kwargs...)
            catch ex
                push!(errors, ex)
            end
        end
    end

    if haskey(_callback_events, evt_name)
        for cb ∈ _callback_events[evt_name]
            try
                cb(args...; kwargs...)
            catch ex
                push!(errors, ex)
            end
        end
    end

    if !isempty(errors)
        throw(CompositeException(errors))
    end
    nothing
end

export listen
function listen() end
function listen(callback, evt_name::Symbol)
    if !haskey(_callback_events, evt_name)
        _callback_events[evt_name] = []
    end
    push!(_callback_events[evt_name], callback)
    nothing
end

export listen_once
function listen_once(callback, evt_name::Symbol)
    function wrapper(args...; kwargs...)
        callback(args...; kwargs...)
        deafen(evt_name, wrapper)
    end
    listen(wrapper, evt_name)
end

export deafen
function deafen(evt_name::Symbol, callback)
    if haskey(_callback_events, evt_name)
        idx = indexof(_callback_events[evt_name], callback; strict=true)
        if idx !== nothing
            deleteat!(_callback_events[evt_name], idx)
        end
    end
    nothing
end

const _static_events = Dict{Symbol, Vector{Symbol}}()
const _callback_events = Dict{Symbol, Vector}()
end # module Events

max_stars_repo_name: Kirusifix/ExtraFun

max_issues_repo_name: Kirusifix/ExtraFun

max_forks_repo_name: Kirusifix/ExtraFun


text:
cd(@__DIR__)
using Pkg
pkg"activate ."

push!(LOAD_PATH, pwd(), "src")

function main()
  include(joinpath("src", "Static.jl"))
end; main()


text:
"""
Writes out the alpha vectors in the `.alpha` file format
"""
function Base.write(io::IO, alphas::POMDPAlphas)

	vector_length = size(alphas.alpha_vectors, 1)

	for i in 1 : length(alphas.alpha_actions)
		println(io, alphas.alpha_actions[i])

		for j in 1 : vector_length
			@printf(io, "%.25f", alphas.alpha_vectors[j,i])
			if j != vector_length
				print(io, " ")
			else
				print(io, "\n")
			end
		end

		println(io, "")
	end
end

"""
Write out a `.pomdp` file using the POMDPs.jl interface
Specification: http://cs.brown.edu/research/ai/pomdp/examples/pomdp-file-spec.html
"""
function Base.write(io::IO, pomdp::POMDP)

    ss = states(pomdp)
    as = actions(pomdp)
    os = observations(pomdp)
	println(io, "discount: ", discount(pomdp))
	println(io, "values: reward") # NOTE(tim): POMDPs.jl assumes rewards rather than costs
    println(io, "states: ", length(ss))
    println(io, "actions: ", length(as))
    println(io, "observations: ", length(os), "\n")

	# --------------------------------------------------------------------------
	# TRANSITION
	#
	# T: <action>
	# %f %f ... %f
	# %f %f ... %f
	# ...
	# %f %f ... %f
	#
	# each row corresponds to one of the start states and
	# each column specifies one of the ending states

    pomdp_states = ordered_states(pomdp)
    pomdp_actions = ordered_actions(pomdp)

	for (action_index, a) in enumerate(pomdp_actions)
		println(io, "T:", action_index-1)
		for s in pomdp_states
			T = transition(pomdp, s, a)
			for (end_state_index, s2) in enumerate(pomdp_states)
				print(io, pdf(T, s2))
				if end_state_index != length(pomdp_states)
					print(io, " ")
				else
					print(io, "\n")
				end
			end
		end
		print(io, "\n")
	end

	# --------------------------------------------------------------------------
	# OBSERVATION
	#
	# O: <action>
	# %f %f ... %f
	# %f %f ... %f
	# ...
	# %f %f ... %f
	#
	# each row corresponds to one of the states and
	# each column specifies one of the observations

    pomdp_observations = ordered_observations(pomdp)

	for (action_index, a) in enumerate(pomdp_actions)
		println(io, "O:", action_index-1)
		for (state_index, s) in enumerate(pomdp_states)
			O = observation(pomdp, a, s)
			for (obs_index, o) in enumerate(pomdp_observations)
				print(io, pdf(O, o))
				if obs_index != length(pomdp_observations)
					print(io, " ")
				else
					print(io, "\n")
				end
			end
		end
		print(io, "\n")
	end

	# --------------------------------------------------------------------------
	# REWARD
	#
	# R: <action> : <start-state> : <end-state> : <observation> %f
	# where * is used for <end-state> and <observation> to indicate a wildcard
	# that is expanded to all existing entities

	for (action_index, a) in enumerate(pomdp_actions)
		for (start_state_index, s) in enumerate(pomdp_states)
            for (end_state_index, sp) in enumerate(pomdp_states)
                for (obs_index, o) in enumerate(pomdp_observations)
                    r = reward(pomdp, s, a, sp, o)
                    println(io, "R:", action_index-1, ":", start_state_index-1,
                            ":", end_state_index-1, ":", obs_index-1, " ", r)
                end
            end
		end
	end
end


text:
@testset "792.number-of-matching-subsequences.jl" begin
    @test num_matching_subseq("dcaog", ["cat", "cop", "dog", "cao"]) == 2
    @test num_matching_subseq("abcde", ["a", "bb", "acd", "ace"]) == 3
end


text:
import LinearAlgebra: norm_sqr

@inline lj_potential_uij(
    r::Float64; ϵ = 1.0, σ = 1.0
) = let r² = r ^ 2, σ² = σ ^ 2
    @fastmath 4 * ϵ * ((σ² ^ 6) / (r² ^ 6) - (σ² ^ 3) / (r² ^ 3))
end

@inline lj_potential_uij(
    r; ϵ = 1.0, σ = 1.0
) = let r² = norm_sqr(r), σ² = σ ^ 2
    @fastmath 4 * ϵ * ((σ² ^ 6) / (r² ^ 6) - (σ² ^ 3) / (r² ^ 3))
end

@inline lj_potential_uij(
    r₁, r₂;
    ϵ = 1.0, σ = 1.0,
    dist = default_distance
) = lj_potential_uij(dist(r₁, r₂), ϵ = ϵ, σ = σ)

_lj_potential_wij(
    r², ϵ = 1.0, σ² = 1.0
) = @fastmath 24 * ϵ * (2 * (σ² ^ 6) / (r² ^ 6) - (σ² ^ 3) / (r² ^ 3))

@inline function lj_potential_fij(
    r; ϵ = 1.0, σ = 1.0
)
    r² = norm_sqr(r)
    wij = _lj_potential_wij(r², ϵ, σ ^ 2)
    fij = r * wij / r²
end

@inline lj_potential_fij(
    r₁, r₂;
    ϵ = 1.0, σ = 1.0,
    dist = default_distance
) = lj_potential_fij(dist(r₁, r₂), ϵ = ϵ, σ = σ)

@inline lj_potential_uij_cutoff(
    r;
    ϵ = 1.0, σ = 1.0,
    r_cutoff = 2.5σ
) = norm_sqr(r) ≥ r_cutoff ^ 2 ? 0.0 :
    lj_potential_uij(r, ϵ = ϵ, σ = σ) - lj_potential_uij(r_cutoff, ϵ = ϵ, σ = σ)

@inline lj_potential_uij_cutoff(
    r₁, r₂;
    ϵ = 1.0, σ = 1.0,
    dist = default_distance,
    r_cutoff = 2.5σ
) = lj_potential_uij_cutoff(dist(r₁, r₂), ϵ = ϵ, σ = σ, r_cutoff = r_cutoff)

@inline lj_potential_fij_cutoff(
    r;
    ϵ = 1.0, σ = 1.0,
    r_cutoff = 2.5σ
) = norm_sqr(r) ≥ r_cutoff ^ 2 ? zero(Vector3) :
    lj_potential_fij(r, ϵ = ϵ, σ = σ)

@inline lj_potential_fij_cutoff(
    r₁, r₂;
    ϵ = 1.0, σ = 1.0,
    dist = default_distance,
    r_cutoff = 2.5σ
) = lj_potential_fij_cutoff(dist(r₁, r₂), ϵ = ϵ, σ = σ, r_cutoff = r_cutoff)


text:
#
# Checks if the data is ok
#

function check_size(X,c)
  if (length(size(X)) > 2) || (length(size(X)) == 2 && size(X,2) != 1)
    error(" Only 1D arrays are accepted, and got $c with dimensions = ",size(X))
  end
  return vec(copy(X))
end

function checkdata(X :: AbstractArray{<:Real}, Y :: AbstractArray{<:Real}, options :: Options)
  if length(X) != length(Y)
    error(" Input x and y vectors must have the same length.")
  end
  X = check_size(X,"X")
  Y = check_size(Y,"Y")
  # Set some reasonable ranges for the initial guesses
  if options.p0_range[1] ≈ 0. && options.p0_range[2] ≈ 0.
    range = maximum(Y) - minimum(Y)
    options.p0_range[1] = -100*range
    options.p0_range[2] = +100*range
  end
  return X, Y
end


text:
using GeometryBasics, Makie
using FileIO, Colors
using RPRMakie
RPRMakie.activate!(iterations=200)
earth = load(Makie.assetpath("earth.png"))
m = uv_mesh(Tesselation(Sphere(Point3f(0), 1.0f0), 60))
f, ax, mplot = Makie.mesh(m; color=earth)
Makie.mesh!(ax, Sphere(Point3f(2, 0, 0), 0.1f0); color=:red)
x, y = collect(-8:0.5:8), collect(-8:0.5:8)
z = [sinc(√(X^2 + Y^2) / π) for X in x, Y in y]
wireframe!(ax, -2 .. 2, -2 .. 2, z)
display(ax.scene)


text:
using LinearAlgebra

using ElectromagneticFields: crossproduct

export InitialConditions, InitialConditionsGC
export charged_particle, guiding_center, pauli_particle
export e, me, mp, md, mα

const e  = 1.602176634E-19    # electron charge
const me = 9.1093837015E-31   # electron mass
const mp = 1.6726219237E-27   # proton mass
const md = 3.3435837724E-27   # deuteron mass
const mα = 6.6446573357E-27   # α-particle mass

"""
Store initial conditions for charged particle, Pauli particle and guiding center models.

Fields:
* `x`: particle position
* `X`: gyro center position
* `ρ`: gyro radius vector
* `vvec`: velocity vector
* `vpar`: parallel velocity vector
* `vper`: perpendicular velocity vector
* `v`: absolute value of velocity
* `u`: absolute value of parallel velocity
* `μ`: magnetic moment
* `θ`: gyro angle ∈ [0,2π]
* `α`: pitch angle ∈ [0,π/2]
* `ω`: gyro frequency
* `mass`: mass in kg
* `energy`: energy in eV
* `charge`: charge number
"""
struct InitialConditions{T <: Real}
    x::Vector{T}
    X::Vector{T}
    ρ::Vector{T}

    vvec::Vector{T}
    vpar::Vector{T}
    vper::Vector{T}

    v::T
    u::T
    μ::T

    θ::T
    α::T
    ω::T

    mass::T
    energy::T
    charge::Int
end

"""
Compute initial conditions from the following arguments:
* `X`: gyro center position
* `θ`: gyro angle
* `α`: pitch angle
* `E`: energy
* `M`: mass
* `C`: charge number
* `a`, b, c: magnetic field unit vectors in physical coordinates
* `b`: magnetic field unit vector in contravariant coordinates
* `B`: amplitude of magnetic field
* `g̅`: inverse metric coefficients
* `DF̄`: inverse Jacobian matrix
* `J`: Jacobian determinant
* `l=1`: length normalization
"""
function InitialConditions(X::AbstractVector{T}, θ::T, α::T, Etot::T, M::T, C, â::Function, b̂::Function, ĉ::Function, b::Function, B::Function, g̅::Function, DF̄::Function, J::Function; l₀=1) where {T}
    # computet parallel and perpedicular energy according to pitch angle
    Epar = Etot * (1 - sin(α))
    Eper = Etot * sin(α)

    # compute absolute value of velocity v from energy Etot,
    # parallel velocity u, perpendicular velocity w, and
    # magnetic moment μ
    v = sqrt(2 * Etot * M / e) / l₀
    u = sqrt(2 * Epar * M / e) / l₀
    w = sqrt(2 * Eper * M / e) / l₀
    μ = w^2 / 2B(0,X)

    # compute v̂par and v̂per in physical coordinates
    v̂par = u * b̂(0,X)
    v̂per = w * ( - â(0,X) * sin(θ) - ĉ(0,X) * cos(θ) )

    # transform vpar and vper to contravariant coordinates
    vpar = DF̄(0,X) * v̂par
    vper = DF̄(0,X) * v̂per

    # compute total velocity vector field
    vvec = vpar .+ vper

    # compute gyro radius vector ρ and particle position x
    ω = e * B(0,X) / M
    ρ = [crossproduct(b(0,X), vvec, g̅(0,X), J(0,X), l) for l in 1:3] / B(0,X)
    x = X .+ ρ

    InitialConditions{T}(x, X, ρ, vvec, vpar, vper, v, u, μ, θ, α, ω, M, Etot, C)
end

"""
Compute initial conditions from the following arguments:
* `X`: gyro center position
* `θ`: gyro angle
* `u`: absolute value of parallel velocity
* `μ`: magnetic moment
* `M`: mass
* `C`: charge number
* `a`, b, c: magnetic field unit vectors in physical coordinates
* `b`: magnetic field unit vector in contravariant coordinates
* `B`: amplitude of magnetic field
* `g̅`: inverse metric coefficients
* `DF̄`: inverse Jacobian matrix
* `J`: Jacobian determinant
* `l=1`: length normalization
"""
function InitialConditionsGC(X::AbstractVector{T}, θ::T, u::T, μ::T, M::T, C, â::Function, b̂::Function, ĉ::Function, b::Function, B::Function, g̅::Function, DF̄::Function, J::Function; l₀=1) where {T}
    # compute absolute value of perpendicular and total velocity
    w = sqrt(2B(0,X) * μ)
    v = sqrt(u^2 + w^2)

    # compute perpendicular and total energy
    Eper = e * l₀^2 * w^2 / 2 / M
    Etot = e * l₀^2 * v^2 / 2 / M

    # compute pitch angle
    α = asin(Eper / Etot)

    # compute v̂par and v̂per in physical coordinates
    v̂par = u * b̂(0,X)
    v̂per = w * ( - â(0,X) * sin(θ) - ĉ(0,X) * cos(θ) )

    # transform vpar and vper to contravariant coordinates
    vpar = DF̄(0,X) * v̂par
    vper = DF̄(0,X) * v̂per

    # compute total velocity vector field
    vvec = vpar .+ vper

    # compute gyro radius vector ρ and particle position x
    ω = e * B(0,X) / M
    ρ = [crossproduct(b(0,X), vvec, g̅(0,X), J(0,X), l) for l in 1:3] / B(0,X)
    x = X .+ ρ

    InitialConditions{T}(x, X, ρ, vvec, vpar, vper, v, u, μ, θ, α, ω, M, Etot, C)
end

function Base.show(io::IO, ics::InitialConditions)
    print(io, "Charged Particle Initial Conditions with\n")
    print(io, "  x   = ", ics.x, "\n")
    print(io, "  X   = ", ics.X, "\n")
    print(io, "  ρ   = ", ics.ρ, "\n")
    print(io, "  v   = ", ics.vvec, "\n")
    print(io, "  v∥  = ", ics.vpar, "\n")
    print(io, "  v⟂  = ", ics.vper, "\n")
    print(io, "  |v| = ", ics.v, "\n")
    print(io, "  u   = ", ics.u, "\n")
    print(io, "  μ   = ", ics.μ, "\n")
    print(io, "  θ   = ", ics.θ, "\n")
    print(io, "  α   = ", ics.α, "\n")
    print(io, "  ω   = ", ics.ω, "\n")
    print(io, "  M   = ", ics.mass, "\n")
    print(io, "  E   = ", ics.energy, "\n")
    print(io, "  C   = ", ics.charge)
end

"""
Extracts the charged particle initial conditions and returns the tuple `(x,v)`.
If the keyword argument `noncanonical` is set to `true`, the functions returns
the vector `vcat(x,v)`.
"""
function charged_particle(ics::InitialConditions{T}; noncanonical=false) where {T}
    if noncanonical
        return vcat(ics.x, ics.v)
    else
        return (ics.x, ics.vvec)
    end
end

"""
Extracts the guiding center initial conditions and returns the tuple `(vcat(X,u),μ)`.
"""
function guiding_center(ics::InitialConditions{T}) where {T}
    (vcat(ics.X, ics.u), ics.μ)
end

"""
Extracts the Pauli particle initial conditions and returns the tuple `(X,vpar,μ)`.
If the keyword argument `noncanonical` is set to `true`, the functions returns
the tuple `(vcat(X,vpar),μ)`.
"""
function pauli_particle(ics::InitialConditions{T}; noncanonical=false) where {T}
    if noncanonical
        return (vcat(ics.X, ics.vpar), ics.μ)
    else
        return (ics.X, ics.vpar, ics.μ)
    end
end


text:
module EscherTodo

using FunctionalCollections
using Reactive
using Escher

# This is an example app to demonstrate Components

include("item.jl")
include("list.jl")

function load_list()
    TodoList("Todo",
        pvec([TodoItem(false, "Make presentation"),
              TodoItem(false, "Don't lose it")]))
end

function app(window)
    input = Signal(Any, nothing)

    list = foldp(load_list(), input) do lst, action
        @show update(lst, action)
    end

    map(list) do l
        view(l) >>> input
    end
end

function debug_app(window)
    input = Signal(Any, nothing)
    idx_signal = Signal(Int, 1)

    list = foldp(load_list(), input) do lst, action
        @show update(lst, action)
    end

    history = foldp(push, pvec(Any[value(list)]), list)
    current_index = merge(map(length, history), idx_signal)

    map(current_index, history) do i, hist
        vbox(
            hbox("Recorded states: ",
                 slider(1:length(hist),
                        value=length(hist)) >>> idx_signal
            ) |> packacross(center),
            view(hist[i]) >>> input
        )
    end
end

end # module


text:
#################################################################################################
# AIM 4: group comparison performance
#################################################################################################
## Deps
using StatsPlots
using Distributed
using Distributions
using DelimitedFiles
@everywhere using Pkg
@everywhere Pkg.activate("./") # <=
@everywhere using CpelNano

## Source code
root_dir = dirname(dirname(pathof(CpelNano)))
include("$(root_dir)/calls/PaperSimulations/Aim-4/Aim-4-Source.jl")

## Constants

# IO
const aim_dir = "/Users/jordiabante/OneDrive - Johns Hopkins/CpelNano/Data/Simulations/Aim-4/" # <=
const blind_friend_col = ["#999999","#E69F00","#56B4E9","#009E73","#F0E442","#0072B2","#D55E00","#CC79A7"]

# Set up
@everywhere const M = 5
@everywhere const POBS = 1.0
@everywhere const Δ_LEVEL_MEAN = 0.13
@everywhere const LEVEL_STD_M = 3.0 / sqrt(4000 / 450)
@everywhere const LEVEL_STD_U = 3.0 / sqrt(4000 / 450)
@everywhere const LEVEL_MEAN_M = Δ_LEVEL_MEAN / 2.0
@everywhere const LEVEL_MEAN_U = -Δ_LEVEL_MEAN / 2.0
@everywhere const EF_PERC_RNG = 0.05:0.05:0.25

## Default attributes
default(titlefont=(14, "arial"),guidefont=(16, "arial"),tickfont=(12, "arial"))

#################################################################################################
# Calls
#################################################################################################

# Arguments
eff_type = ARGS[1]
matched = parse(Bool, ARGS[2])
n_sim_reps = parse(Int64, ARGS[3])
CpelNano.print_log("Running ROC simulations for $(matched) $(eff_type) w/ $(n_sim_reps) reps")

# Run
if matched

    ### Matched

    # Vars
    n_grp_reps = 6

    # Simulations
    eff_size_arr, tpr, fpr = run_roc_sim(n_sim_reps, n_grp_reps, eff_type, matched)

    # Plot
    plt = plt_roc(eff_size_arr, tpr, fpr, "ROC - $(eff_type) - matched ")
    savefig(plt, "$(aim_dir)/ROC/ROC-Matched-$(eff_type).pdf")

else

    ### Unmatched

    # Vars
    n_grp_reps = 5

    # Simulations
    eff_size_arr, tpr, fpr = run_roc_sim(n_sim_reps, n_grp_reps, eff_type, matched)

    # Plot
    plt = plt_roc(eff_size_arr, tpr, fpr, "ROC - $(eff_type) - unmatched ")
    savefig(plt, "$(aim_dir)/ROC/ROC-Unmatched-$(eff_type).pdf")

end


text:
ts = ["the three truths", "ababababab"]
tsub = ["th", "abab"]

println("Test of non-overlapping substring counts.")
for i in 1:length(ts)
    print(ts[i], " (", tsub[i], ") => ")
    println(length(matchall(Regex(tsub[i]), ts[i])))
end
println()
println("Test of overlapping substring counts.")
for i in 1:length(ts)
    print(ts[i], " (", tsub[i], ") => ")
    println(length(matchall(Regex(tsub[i]), ts[i], true)))
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


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

using Markdown
using InteractiveUtils

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

# ╔═╡ 6df3beed-24a7-4b26-a315-0520f4863190
develop=true

# ╔═╡ 9701cbe0-d048-11eb-151b-67dda7b72b71
begin
   using Pkg
   Pkg.activate(mktempdir())
   Pkg.add("Revise")
   using Revise
   Pkg.add(["PyPlot","PlutoUI","Plots"])
	if develop
		Pkg.add("GLMakie")
	   Pkg.develop(["PlutoVista","GridVisualize", "ExtendableGrids"])
	else
	   Pkg.add(["PlutoVista","GridVisualize","ExtendableGrids"])
	end
   using PlutoVista,GridVisualize,ExtendableGrids,PlutoUI
   import Plots,PyPlot
	if develop
		import GLMakie
	else
		GLMakie=nothing
	end
end

# ╔═╡ 62c6c85c-e374-4315-b737-b7a0005d2985
Plots.gr()

# ╔═╡ 68e2c958-b417-4ba1-9577-697304fe140a
TableOfContents()

# ╔═╡ b35d982d-1fa9-413d-b008-892b4f241097
md"""
# Test notebook for PlutoVista backend of GridVisualize
"""

# ╔═╡ 00b04f6b-34a6-4e30-9864-d273305281d4
md"""
## 1D scalar plot
"""

# ╔═╡ a20d74c9-16da-408a-b247-0c17321888f9
function testplot1(Plotter)
	grid=simplexgrid(0:0.01:10)
	scalarplot(grid,map(sin,grid),Plotter=Plotter,resolution=(600,200),markershape=:star5,markevery=20, xlabel="x",ylabel="z",legend=:rt,label="sin")
end

# ╔═╡ cc17187f-404c-4c31-8625-fa067eea7273
testplot1(PyPlot)

# ╔═╡ ad1c4dd6-7c6b-4433-a1ac-b2f817ba5d81
testplot1(Plots)

# ╔═╡ 33482af8-3542-4723-ae43-770a789b69b3
testplot1(PlutoVista)

# ╔═╡ d21982eb-aa7e-424c-932e-f0b1460f19b9
testplot1(GLMakie)

# ╔═╡ c4eeb06f-932e-4acc-8e5b-f2a7f9242a42
function testplot2(Plotter;t=0)
	p=GridVisualizer(Plotter=Plotter, resolution=(500,200),legend=:rt,xlabel="x")
	grid=simplexgrid(0:0.01:10)
	scalarplot!(p,grid,map(x->sin(x-t),grid),Plotter=Plotter,color=:red,label="sin(x-$(t))",linestyle=:dash)
	scalarplot!(p,grid,map(cos,grid),Plotter=Plotter,color=:green,clear=false,label="cos",linestyle=:dashdot,linewidth=3)
	reveal(p)
end

# ╔═╡ 29ca4775-6ba5-474c-bd2c-8f770b09addd
testplot2(PyPlot)

# ╔═╡ c0ab77e8-01ea-436d-85f2-34e253944f11
testplot2(Plots)

# ╔═╡ 4f73fb99-fab9-4ca4-addb-7d12d16ceba0
testplot2(GLMakie)

# ╔═╡ 84192945-d4b6-4949-8f06-d94e04a7a56d
testplot2(PlutoVista)

# ╔═╡ 7fbaf93f-3cfb-47d0-8252-487e60ba3e54
md"""
### Changing data by re-creating a plot
"""

# ╔═╡ 63fe3259-7d79-40ec-98be-e0592e40ee6b
@bind t2 PlutoUI.Slider(0:0.1:5,show_value=true)

# ╔═╡ 4de6b5c9-4d2d-4bcb-bc88-c6f50a23f9b6
testplot2(PlutoVista,t=t2)

# ╔═╡ 2061e7fd-c740-4d4b-af5b-7a3a9444aafd
md"""
### Changing data by updating the plot

For this pattern, we observe a notable difference to other backends:  the plot with PlutoVista appears above the cell where the GridVisualizer is created instead of the cell where data are plotted. The reason is that for updating data, we need to have a visualization context which stays the same.
"""

# ╔═╡ f84beb4f-4136-4e5a-ba43-279b703fc75f
begin
	X2=0:0.001:10
	grid2=simplexgrid(collect(X2))
	f2(t)=map( (x)->sin(x^2-t),grid2)
end

# ╔═╡ 29fa4467-65ee-4dad-a660-5197864ddbdc
md"""
t4: $(@bind t4 PlutoUI.Slider(-10:0.1:10, default=0, show_value=true))
"""

# ╔═╡ c1278fb2-3e75-445f-893a-b8b8a7e931d3
p=GridVisualizer(Plotter=PlutoVista,resolution=(600,200),dim=1,legend=:lt);p

# ╔═╡ 661531f7-f740-4dd4-9a59-89ddff06ba5c
scalarplot!(p,X2,f2(t4),show=true,clear=true,color=:red,label="t=$(t4)")

# ╔═╡ ed9b80e5-9678-4ba6-bb36-c2e0674ed9ba
md"""
## 1D grid plot
"""

# ╔═╡ 9ce4f63d-cd96-48d7-a637-07cb84fa88ab
function testgridplot(Plotter)
	grid=simplexgrid(0:1:10)
	cellmask!(grid,[0.0],[5],2)
	bfacemask!(grid,[5],[5],3)
	gridplot(grid, Plotter=Plotter,resolution=(600,200),legend=:rt)
end

# ╔═╡ 77eeefc7-e416-426b-8f87-1bc8439dae6d
testgridplot(PyPlot)

# ╔═╡ d503ee1e-1e1f-4235-b286-dc3137a2c96a
testgridplot(PlutoVista)

# ╔═╡ 395d05f4-46a9-4e85-804c-35c79acc96e3
testgridplot(GLMakie)

# ╔═╡ ae1fe1ab-4a0e-4c80-bd6f-912201fb4bb4
md"""
## 2D Scalar plot
"""

# ╔═╡ d5258595-60e4-406f-a71e-69111cdad8b9
function testplot3(Plotter)
	X=0:0.1:10
	grid=simplexgrid(X,X)
	f=map( (x,y)->sin(x)*atan(y),grid)
	scalarplot(grid,f,Plotter=Plotter,
		resolution=(300,300),limits=(-π/2,π/2))
end

# ╔═╡ c98a90bf-1a3e-4681-a3b0-663c6844df6b
testplot3(PyPlot)

# ╔═╡ a0c3067b-3aa5-493e-b132-89746483b5ce
testplot3(Plots)

# ╔═╡ 0998a9a7-d57a-476e-aacd-bee9396e9b8f
testplot3(PlutoVista)

# ╔═╡ cefb38c1-159e-42db-8088-294573fcece2
md"""
### Changing data

Generally, as above, with Plutovista we need two cells - one with the graph shown, and a second one which triggers the modification.
"""

# ╔═╡ a9f4f98f-ec2f-42d6-88da-4a8a6f727e93
begin
	X=0:0.05:10
	grid=simplexgrid(X,X)
	f(t)=map( (x,y)->sin(x-t)*atan(y)*cos((y-t)),grid)
end

# ╔═╡ faa59bbd-df1f-4c62-9a77-c4752c6a6df4
vis=GridVisualizer(resolution=(300,300),Plotter=PlutoVista,dim=2);vis

# ╔═╡ 6f1707ed-79ab-42dc-8ad8-d66a9e1a65b3
md"""
t= $(@bind t PlutoUI.Slider(-10:0.1:10, default=0, show_value=true))
"""

# ╔═╡ 412c905f-050c-4b78-a66f-0d03978e7edf
scalarplot!(vis,grid,f(t),limits=(-π/2,π/2),show=true,levels=4)

# ╔═╡ e9bc2dae-c303-4063-9ea9-36f95f93371c
md"""
## 2D Grid plot
"""

# ╔═╡ 2b3cb0f4-0656-4981-bec6-48785caf2994
function testgridplot2d(Plotter)
	X=-1:0.2:1
	grid=simplexgrid(X,X)
	gridplot(grid,Plotter=Plotter,resolution=(300,300))
end

# ╔═╡ db488643-2e6b-40d6-ba81-126c752953c5
testgridplot2d(PyPlot)

# ╔═╡ 1388c246-be49-4757-a2cc-a685642b6b37
testgridplot2d(PlutoVista)

# ╔═╡ 3f807dfe-7ad3-4740-9e4f-23a39c15603c
testgridplot2d(GLMakie)

# ╔═╡ 5eee8f1d-49ca-4e95-bd14-fe415b0c15e5
md"""
## 3D Scalar plot

Here we use the possibility to update plots to allow moving isosurfaces and plane cuts.
"""

# ╔═╡ 0c99daca-f9a8-4116-867b-e13461c3e754
function grid3d(;n=15)
    X=collect(0:1/n:1)
    g=simplexgrid(X,X,X)
end

# ╔═╡ 82ccfd24-0053-4399-9bc8-b2e4010bbc92
function func3d(;n=15)
    g=grid3d(n=n)
    g, map((x,y,z)->sinpi(2*x)*sinpi(3.5*y)*sinpi(1.5*z),g)
end

# ╔═╡ 8b20f720-5470-4da7-bbb6-b746e887046e
g3,f3=func3d(n=19)

# ╔═╡ c0a0ea34-6fc3-4409-934e-086a1a36f94e
p3d=GridVisualizer(Plotter=PlutoVista,resolution=(300,300),dim=3);p3d

# ╔═╡ 35be5ef4-0664-4196-8f10-cf71ec7cb371
md"""
f: $(@bind flevel Slider(0:0.01:1,show_value=true,default=0.45))

x: $(@bind xplane Slider(0:0.01:1,show_value=true,default=0.45))
y: $(@bind yplane Slider(0:0.01:1,show_value=true,default=0.45))
z: $(@bind zplane Slider(0:0.01:1,show_value=true,default=0.45))
"""

# ╔═╡ ecd941a0-85b7-4bb7-a903-b19a500198e1
scalarplot!(p3d,g3,f3;levels=[flevel],xplanes=[xplane],yplanes=[yplane],zplanes=[zplane],colormap=:hot,outlinealpha=0.05,show=true,levelalpha=0.5)

# ╔═╡ d924d90d-4102-4ae8-b8de-254a17a5d4df
X4=-1:0.1:1; g4=simplexgrid(X4,X4,X4)

# ╔═╡ 57ed5eea-bc1c-45eb-b4d3-dc63088db21a
scalarplot(g4,map( (x,y,z)-> 0.01*exp(-10*(x^2+y^2+z^2)),g4),levels=3,Plotter=PlutoVista)

# ╔═╡ 597849e9-b9a7-4728-a278-7571d7c1a625
scalarplot(Plotter=PyPlot,g3,f3;resolution=(300,300),levels=[0.5],
	xplanes=[0.4],yplanes=[0.4],zplanes=[0.4],show=true,colormap=:viridis)

# ╔═╡ 4b9113d2-10bd-4f7a-a2b8-22092656c6b3
md"""
## 3D grid plot
"""

# ╔═╡ 81f0a07d-3d0c-4e7a-9684-1ca4d584b210
gridplot(Plotter=PyPlot,g3; resolution=(300,300),
	xplanes=[1.0],yplanes=[1.0],zplanes=[0.4],show=true)

# ╔═╡ f78196ca-d972-4fa6-bdc2-e76eba7ca5a1
p3dgrid=GridVisualizer(Plotter=PlutoVista,resolution=(300,300),dim=3)

# ╔═╡ 7dd92757-c100-4158-baa8-1d9218c39aa7
md"""
x: $(@bind gxplane Slider(0:0.01:1,show_value=true,default=1.0))
y: $(@bind gyplane Slider(0:0.01:1,show_value=true,default=1.0))
z: $(@bind gzplane Slider(0:0.01:1,show_value=true,default=0.45))
"""

# ╔═╡ 840c80b7-5186-45a5-a8df-ec4fb50a5dbb
gridplot!(p3dgrid,g3; xplane=gxplane,yplane=gyplane,zplane=gzplane,show=true)

# ╔═╡ e656dac5-466e-4c07-acfa-0478ad000cb2
function qv2d(;Plotter=nothing,n=20,stream=false,kwargs...)

    X=0.0:10/(n-1):10

    f(x,y)=sin(x)*cos(y)+0.05*x*y
    fx(x,y)=-cos(x)*cos(y)-0.05*y
    fy(x,y)=sin(x)*sin(y)-0.05*x

    grid=simplexgrid(X,X)

    v=map(f,grid)
    ∇v=vcat(map(fx,grid)',map(fy,grid)')

    gvis=GridVisualizer(Plotter=Plotter,resolution=(600,600))
    scalarplot!(gvis,grid,v,colormap=:summer,levels=7)
	if stream
    	streamplot!(gvis,grid,∇v;clear=false, show=true,kwargs...)
	else
  	 vectorplot!(gvis,grid,∇v;clear=false, show=true,kwargs...)
	end
	reveal(gvis)
end

# ╔═╡ 812af347-7606-4c54-b155-88322d20d921
qv2d(Plotter=PlutoVista,n=50,spacing=0.35)

# ╔═╡ 9462f2a4-e7ef-469f-85d4-12f3fa411091
qv2d(Plotter=PyPlot,n=100,spacing=0.35)

# ╔═╡ 8293e0a4-ae15-4929-b20d-95369788f3c9
qv2d(Plotter=PyPlot,n=100,spacing=0.06,stream=true,color=:darkred)

# ╔═╡ 7c66ca9f-7365-489e-ae50-f8c3f9eb145b
qv2d(Plotter=GLMakie)

# ╔═╡ f03cb172-82c7-4187-be43-01d307977713
qv2d(Plotter=Plots);

# ╔═╡ Cell order:
# ╠═6df3beed-24a7-4b26-a315-0520f4863190
# ╠═9701cbe0-d048-11eb-151b-67dda7b72b71
# ╠═62c6c85c-e374-4315-b737-b7a0005d2985
# ╠═68e2c958-b417-4ba1-9577-697304fe140a
# ╟─b35d982d-1fa9-413d-b008-892b4f241097
# ╟─00b04f6b-34a6-4e30-9864-d273305281d4
# ╠═a20d74c9-16da-408a-b247-0c17321888f9
# ╠═cc17187f-404c-4c31-8625-fa067eea7273
# ╠═ad1c4dd6-7c6b-4433-a1ac-b2f817ba5d81
# ╠═33482af8-3542-4723-ae43-770a789b69b3
# ╠═d21982eb-aa7e-424c-932e-f0b1460f19b9
# ╠═c4eeb06f-932e-4acc-8e5b-f2a7f9242a42
# ╠═29ca4775-6ba5-474c-bd2c-8f770b09addd
# ╠═c0ab77e8-01ea-436d-85f2-34e253944f11
# ╠═4f73fb99-fab9-4ca4-addb-7d12d16ceba0
# ╠═84192945-d4b6-4949-8f06-d94e04a7a56d
# ╟─7fbaf93f-3cfb-47d0-8252-487e60ba3e54
# ╟─63fe3259-7d79-40ec-98be-e0592e40ee6b
# ╠═4de6b5c9-4d2d-4bcb-bc88-c6f50a23f9b6
# ╟─2061e7fd-c740-4d4b-af5b-7a3a9444aafd
# ╠═f84beb4f-4136-4e5a-ba43-279b703fc75f
# ╟─29fa4467-65ee-4dad-a660-5197864ddbdc
# ╠═c1278fb2-3e75-445f-893a-b8b8a7e931d3
# ╠═661531f7-f740-4dd4-9a59-89ddff06ba5c
# ╠═ed9b80e5-9678-4ba6-bb36-c2e0674ed9ba
# ╠═9ce4f63d-cd96-48d7-a637-07cb84fa88ab
# ╠═77eeefc7-e416-426b-8f87-1bc8439dae6d
# ╠═d503ee1e-1e1f-4235-b286-dc3137a2c96a
# ╠═395d05f4-46a9-4e85-804c-35c79acc96e3
# ╟─ae1fe1ab-4a0e-4c80-bd6f-912201fb4bb4
# ╠═d5258595-60e4-406f-a71e-69111cdad8b9
# ╠═c98a90bf-1a3e-4681-a3b0-663c6844df6b
# ╠═a0c3067b-3aa5-493e-b132-89746483b5ce
# ╠═0998a9a7-d57a-476e-aacd-bee9396e9b8f
# ╟─cefb38c1-159e-42db-8088-294573fcece2
# ╠═a9f4f98f-ec2f-42d6-88da-4a8a6f727e93
# ╠═faa59bbd-df1f-4c62-9a77-c4752c6a6df4
# ╠═412c905f-050c-4b78-a66f-0d03978e7edf
# ╟─6f1707ed-79ab-42dc-8ad8-d66a9e1a65b3
# ╟─e9bc2dae-c303-4063-9ea9-36f95f93371c
# ╠═2b3cb0f4-0656-4981-bec6-48785caf2994
# ╠═db488643-2e6b-40d6-ba81-126c752953c5
# ╠═1388c246-be49-4757-a2cc-a685642b6b37
# ╠═3f807dfe-7ad3-4740-9e4f-23a39c15603c
# ╟─5eee8f1d-49ca-4e95-bd14-fe415b0c15e5
# ╠═0c99daca-f9a8-4116-867b-e13461c3e754
# ╠═82ccfd24-0053-4399-9bc8-b2e4010bbc92
# ╠═8b20f720-5470-4da7-bbb6-b746e887046e
# ╠═c0a0ea34-6fc3-4409-934e-086a1a36f94e
# ╟─35be5ef4-0664-4196-8f10-cf71ec7cb371
# ╠═ecd941a0-85b7-4bb7-a903-b19a500198e1
# ╠═d924d90d-4102-4ae8-b8de-254a17a5d4df
# ╠═57ed5eea-bc1c-45eb-b4d3-dc63088db21a
# ╠═597849e9-b9a7-4728-a278-7571d7c1a625
# ╟─4b9113d2-10bd-4f7a-a2b8-22092656c6b3
# ╠═81f0a07d-3d0c-4e7a-9684-1ca4d584b210
# ╠═f78196ca-d972-4fa6-bdc2-e76eba7ca5a1
# ╟─7dd92757-c100-4158-baa8-1d9218c39aa7
# ╠═840c80b7-5186-45a5-a8df-ec4fb50a5dbb
# ╠═e656dac5-466e-4c07-acfa-0478ad000cb2
# ╠═812af347-7606-4c54-b155-88322d20d921
# ╠═9462f2a4-e7ef-469f-85d4-12f3fa411091
# ╠═8293e0a4-ae15-4929-b20d-95369788f3c9
# ╠═7c66ca9f-7365-489e-ae50-f8c3f9eb145b
# ╠═f03cb172-82c7-4187-be43-01d307977713


text:
import ...BalanceLaws:
    vars_state,
    init_state_prognostic!,
    init_state_auxiliary!,
    update_auxiliary_state!,
    integral_load_auxiliary_state!,
    integral_set_auxiliary_state!

struct VerticalIntegralModel{M} <: BalanceLaw
    ocean::M
    function VerticalIntegralModel(ocean::M) where {M}
        return new{M}(ocean)
    end
end

vars_state(tm::VerticalIntegralModel, st::Prognostic, FT) =
    vars_state(tm.ocean, st, FT)

function vars_state(m::VerticalIntegralModel, ::Auxiliary, T)
    @vars begin
        ∫x::SVector{2, T}
    end
end

init_state_auxiliary!(
    tm::VerticalIntegralModel,
    A::MPIStateArray,
    grid,
    direction,
) = nothing

function vars_state(m::VerticalIntegralModel, ::UpwardIntegrals, T)
    @vars begin
        ∫x::SVector{2, T}
    end
end

@inline function integral_load_auxiliary_state!(
    m::VerticalIntegralModel,
    I::Vars,
    Q::Vars,
    A::Vars,
)
    I.∫x = A.∫x

    return nothing
end

@inline function integral_set_auxiliary_state!(
    m::VerticalIntegralModel,
    A::Vars,
    I::Vars,
)
    A.∫x = I.∫x

    return nothing
end

function update_auxiliary_state!(
    dg::DGModel,
    tm::VerticalIntegralModel,
    x::MPIStateArray,
    t::Real,
)
    A = dg.state_auxiliary

    # copy tendency vector to aux state for integration
    function f!(::VerticalIntegralModel, x, A, t)
        @inbounds begin
            A.∫x = @SVector [x.u[1], x.u[2]]
        end

        return nothing
    end
    update_auxiliary_state!(f!, dg, tm, x, t)

    # compute integral for Gᵁ
    indefinite_stack_integral!(dg, tm, x, A, t) # bottom -> top

    return true
end


text:
function lab_18(r::Robot)

    putmarker!(r)

    path = go_to_west_south_corner_and_return_path!(r; go_around_barriers = true)

    for side in (Nord, Ost, Sud, West)
        go_to_border_and_return_path!(r, side)
        putmarker!(r)
    end

    go_by_path!(r, path)
end


text:
export rule

@rule Bernoulli(:p, Marginalisation) (m_out::PointMass, ) = begin
    r = mean(m_out)
    Beta(one(r) + r, 2one(r) - r)
end

@rule Bernoulli(:p, Marginalisation) (q_out::Bernoulli, ) = begin
    r = first(probvec(q_out))
    return Beta(one(r) + r, 2one(r) - r)
end


text:
# ## InflationFixedExclusion.jl - Función de inflación de exclusión fija de gastos básicos sin utilizar fórmula de inflación IPC

"""
    InflationFixedExclusion{N} <: InflationFunction

Función de inflación para computar la inflación de exclusión fija de gastos
básicos para `N` bases del IPC. Esta versión no utiliza la fórmula del IPC, como
la función de inflación `InflationFixedExclusionCPI`.

## Utilización

    function (inflfn::InflationFixedExclusion)(base::VarCPIBase{T}, i::Int) where T

Define cómo opera InflationFixedExclusion sobre un objeto de tipo `VarCPIBase`,
con `i` listas de exclusión de gastos básicos (posiciones en el IPC) para las
bases 2000 y 2010.

    function (inflfn::InflationFixedExclusion)(cs::CountryStructure, ::CPIVarInterm)

Define cómo opera InflationFixedExclusion sobre un CountryStructure

## Ejemplos

Instanciamos la función y le damos dos listas de exclusión, una por cada base en
el CountryStructure. Se exclyen los gastos básicos con números de columnas
`[25,30,54,88]` del primer `VarCPIBase` y los gastos básicos `[65,95,85]` del
segundo `VarCPIBase`.
```julia-repl
julia> InfExc = InflationFixedExclusion([25,30,54,88], [65,95,85])
(::InflationFixedExclusion{2}) (generic function with 6 methods)
```

"""
Base.@kwdef struct InflationFixedExclusion{N} <: InflationFunction
    # Tupla con vectores de gastos básicos a excluir en cada base (tantos vectores como bases)
    v_exc::NTuple{N,Vector{Int}}
end

# Ampliación para que reciba solo los vectores, no necesariamente como una tupla.
InflationFixedExclusion(v_exc...) = InflationFixedExclusion(v_exc)

# Extender el método de nombre y de tag
measure_name(inflfn::InflationFixedExclusion) = "Exclusión fija de gastos básicos " * string(map(length, inflfn.v_exc))

# Método para obtener parámetros
params(inflfn::InflationFixedExclusion) = inflfn.v_exc

# Cómputo del resumen intermensual utilizando la lista de exclusión i
function (inflfn::InflationFixedExclusion)(base::VarCPIBase{T}, i::Int) where T
    # Elección del vector de exclusión a utlizar dependiendo de que base se está tomando
    # 1 -> base 2000, 2-> base 2010 (se puede generalizar para más bases)
    exc = inflfn.v_exc[i]
    # creación de una copia de la lista original de pesos desde base.w
    w_exc = copy(base.w)
    # Asignación de peso cero a los gastos básicos de la lista de exclusión (exc = inflfn.v_exc[i])
    w_exc[exc] .= 0
    # Renormalización de pesos
    w_exc = w_exc / sum(w_exc)
    # Obtener variación intermensual
    base.v * w_exc
end

# Redefinición del método para obtener variaciones intermensuales del CountryStructure
function (inflfn::InflationFixedExclusion)(cs::CountryStructure, ::CPIVarInterm)
    # Acá se llama a inflfn(base, i), en donde base es de tipo VarCPIBase e i es la posición del vector de exclusión.
    # Esta es la función que debe definirse para cualquier medida de inflación.
    l = length(cs.base)
    vm = mapfoldl(i -> inflfn(cs.base[i],i), vcat, 1:l)
    vm
end


text:
# download input files from: https://gist.github.com/isaacsas/7648fe065f1884307a2906a0a48ea12b
# then add ReactionNetworkImporter package (unregistered)
#]
# add https://github.com/isaacsas/ReactionNetworkImporter.jl.git

using DiffEqBase, DiffEqBiological, Plots, OrdinaryDiffEq, Sundials, DataFrames, CSVFiles, LinearAlgebra
using ReactionNetworkImporter
using TimerOutputs

# parameters
doplot = true
networkname = "testbcrbng"
tf = 10000.

# BNG simulation data
datadir  = joinpath(@__DIR__,"../data/bcr")
fname    = joinpath(datadir, "bcr.net")
cdatfile = joinpath(datadir, "bcr.cdat")
gdatfile = joinpath(datadir, "bcr.gdat")
print("getting cdat file...")
#cdatdf = CSV.File(cdatfile, delim=" ", ignorerepeated=true) |> DataFrame
cdatdf = DataFrame(load(File(format"CSV", cdatfile), header_exists=true, spacedelim=true) )
println("done")
print("getting gdat file...")
#gdatdf = CSV.File(gdatfile, delim=" ", ignorerepeated=true) |> DataFrame
gdatdf = DataFrame(load(File(format"CSV", gdatfile), header_exists=true, spacedelim=true) )
println("done")

const to = TimerOutput()
reset_timer!(to)

@timeit to "bionetgen" prnbng = loadrxnetwork(BNGNetwork(),string(networkname,"bng"), fname);
rnbng = prnbng.rn; u₀ = prnbng.u₀; p = prnbng.p;

u0 = u₀   #convert.(Float64, prn.u₀)

# get the bng reaction network
@timeit to "baddodes" addodes!(rnbng; build_jac=false, build_symfuncs=false)
@timeit to "bODEProb" boprob = ODEProblem(rnbng, u0, (0.,tf), p)
show(to)

# BNG simulation data testing
asykgroups = prnbng.groupstoids[:Activated_Syk]
asyksyms = findall(x -> x ∈ asykgroups, rnbng.syms_to_ints)
asynbng = zeros(length(cdatdf[:time]))
for sym in asyksyms
    global asynbng
    asynbng += cdatdf[sym]
end

# note solvers run _much_ faster the second time
reset_timer!(to); @timeit to "BNG_CVODE_BDF" begin bsol = solve(boprob, CVODE_BDF(),dense=false, saveat=1., abstol=1e-8, reltol=1e-8); end; show(to)

basyk = sum(bsol[asykgroups,:], dims=1)
#vars = sum(bsol[prnbng.groupstoids[:Ig_alpha_P],:],dims=1)

if doplot
    plotlyjs()
    plot(gdatdf[:time][2:end], gdatdf[:Activated_Syk][2:end], xscale=:log10, label=:AsykGroup)
    plot!(cdatdf[:time][2:end], asynbng[2:end], xscale=:log10, label=:AsykSum)
    plot!(bsol.t[2:end], basyk'[2:end], label=:AsykDEBio, xscale=:log10)
    # plot!(bsol.t[2:end], vars'[2:end], label=:Ig_alpha_P, xscale=:log10)
end

# test the error
norm(gdatdf[:Activated_Syk] - asynbng, Inf)
norm(gdatdf[:Activated_Syk] - basyk', Inf)
norm(asynbng - basyk', Inf)

max_issues_repo_name: isaacsas/ParseRxNetwork


text:
# # A sphere in a point charge field.

# In general, the LaplaceBIE allows finding the field on the body's surface if the potential in the absence of an object(s) is known. (Though one can not add objects one by one to get the solution directly. One, however, should be able to make an iterative scheme in the lines of multipole method.) To demonstrate that, let's consider a field on the dielectric sphere (radius r1) due to a point charge at distance ζ away. The solution inside the sphere is given as series [Straton page 221]:
# ```math
# \psi = \sum_{n=0}^{\infty} \frac{2n + 1}{\epsilon n + n + 1} \frac{r_1^n}{\zeta^{n+1}} L_n(\cos \theta)
# ```
# where $L_n$ is Legendre polynomial and $\theta$ angle from line connecting point charge and spehre.

using Jacobi
function ψt(cosθ,ϵ,ζ,r1)
    s = 0
    for n in 0:25
        s += (2*n + 1)/(n*ϵ + n + 1)*r1^n/ζ^(n+1)*legendre(cosθ,n)
    end
    return s
end

# The normal derivative can also be given for the sphere which coincides with derivative in the radial direction:
# ```math
# \frac{\partial \psi}{\partial n} = \sum_{n=0}^{\infty} \frac{n(2n + 1)}{\epsilon n + n + 1} \frac{r_1^{n-1}}{\zeta^{n+1}} L_n(\cos \theta)
# ```
function ∇ψn(cosθ,ϵ,ζ,r1)
    s = 0
    for n in 1:25
        s += n*(2*n + 1)/(n*ϵ + n + 1)*r1^(n-1)/ζ^(n+1) * legendre(cosθ,n)
    end
    return s
end

# First, we can define the exterior potential in the absence of objects:
using LinearAlgebra

ζ = 1.2
r1 = 1
y = [0,0,ζ]

ψ0(x) = 1/norm(x.-y)
∇ψ0(x) = -(x.-y)/norm(x.-y)^3

# Let's now place a sphere in the field with relative permittivity $\epsilon=10$.

ϵ = 10

include("sphere.jl")
msh = unitsphere(2)
vertices, faces = msh.vertices, msh.faces
n = vertices

# With LaplaceBIE to calculate the surface field, we execute these three lines:

using LaplaceBIE

ψ = surfacepotential(vertices,n,faces,ϵ,ψ0);
P∇ψ = tangentderivatives(vertices,n,faces,ψ);
n∇ψ = normalderivatives(vertices,n,faces,P∇ψ,ϵ,∇ψ0);

# Now to compare with analytics, we use azimuthal symmetry, which in numerics was chosen as the x-axis:

using Winston

cosθ = range(-1,1,length=100)
cosθs = [x[3]/1 for x in vertices]

# Potential on the surface vs analytics:

scatter(cosθs,ψ)
oplot(cosθ,ψt.(cosθ,ϵ,ζ,r1) )
savefig("potential.svg")

# ![](potential.svg)

# Normal derivatives on the surface vs analytics:

scatter(cosθs,n∇ψ)
oplot(cosθ,∇ψn.(cosθ,ϵ,ζ,r1) )
savefig("nderivatives.svg")

# ![](nderivatives.svg)


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

using Markdown
using InteractiveUtils

# ╔═╡ 81915028-18e7-11eb-0aed-57a75839d4a7
begin
	using CSV
	using DataFrames
	using Plots
	using StatsPlots
	using Pipe: @pipe
	using GLM
	using StatsBase
	using Statistics
	using Dates
end

# ╔═╡ 02eb9d36-18e8-11eb-24a9-db6f58080189
md"""
# Seoul Bike Sharing Demand

This notebook demonstrates how to use GLM to analyze the Seoul Bike Sharing Demand
data set.

## Data definitions

The following definitions are captured from the [UCI Datasets Archive page](https://archive.ics.uci.edu/ml/datasets/Seoul+Bike+Sharing+Demand).

- Date: year-month-day
- Rented Bike count: Count of bikes rented at each hour
- Hour: Hour of he day
- Temperature: Temperature in Celsius
- Humidity: %
- Windspeed: m/s
- Visibility: 10m
- Dew point temperature: Celsius
- Solar radiation: MJ/m2
- Rainfall: mm
- Snowfall: cm
- Seasons: Winter, Spring, Summer, Autumn
- Holiday: Holiday/No holiday
- Functional Day: NoFunc(Non Functional Hours), Fun(Functional hours)

As the dataset is targeted for bike sharing demand, the key metrics should be the Rented Bike Count field. The rest of the columns are variables that may be correlated to bike sharing demand.

Taking a quick peek at the variables, it makes sense that they may directly or indirectly affect people's decision of renting a bike. For example, when it's very hot or very cold, then I may take on some other transportation. Likewise, humidity, solar radiation, visibility, rainfall, and snowfall probably play a part as well.

Most variables are continuous variables. The last three columns seasons, holiday, and functional day are discrete variables. As for the hour of day variable, we may want to consider that as discrete because increasing the value (later in the day) does not always affect the response variable (rented bike count) in the same direction.
"""

# ╔═╡ eedb4bf8-19a9-11eb-106b-69b047d4118e
md"""
## Reading data
"""

# ╔═╡ dfa5fbf0-18e7-11eb-391b-1d878b9c96f0
file = "https://archive.ics.uci.edu/ml/machine-learning-databases/00560/SeoulBikeData.csv"

# ╔═╡ aed964bc-18e7-11eb-11f1-a3d9957a45f5
df = DataFrame(CSV.File(download(file); dateformats = Dict(:Date => "dd/mm/yyyy")));

# ╔═╡ 4d4a9158-18ec-11eb-2fb0-c178637f903a
md"""
## Data wrangling
"""

# ╔═╡ 09577b1e-19aa-11eb-3e6c-175705b4b951
md"""
Let's quickly examine the data frame and its columns.
"""

# ╔═╡ c5b4e26a-18e7-11eb-2bef-2d937faa252d
describe(df, :eltype, :min, :mean, :max, :nmissing)

# ╔═╡ 6ddeb944-18ec-11eb-05c7-b51e7b01578e
md"It looks pretty clean in general but let's rename some columns so they can be referenced more easily in the code below."

# ╔═╡ ebc9c09a-18e9-11eb-1f56-ef92258af11f
rename!(df,
	2 => :RentedBikeCount,
	4 => :Temperature, 5 => :Humidity, 6 => :WindSpeed, 7 => :Visibility,
	8 => :DewPointTemperature, 9 => :SolarRadiation, 10 => :Rainfall, 11 => :Snowfall,
	14 => :FunctioningDay
);

# ╔═╡ 16c31048-19b6-11eb-2ad8-1f984acfb06f
md"""
## Correlation analysis
"""

# ╔═╡ 9564ab04-19b4-11eb-3c4c-815fe8a6031a
names(df)

# ╔═╡ 3084d7a6-19b5-11eb-1fa5-9b9f7cd482f4
let
	cols = 3:11
	colnames = names(df)[cols]
	Σ = cor(Matrix(df[:, cols]))
	heatmap(colnames, colnames, Σ; xrotation = 45, seriescolor = :PiYG_4)
end

# ╔═╡ 2ef9d8a6-19b6-11eb-0cc5-e1831ce6c40a
md"""
It appears that temperature and dew point temperature columns are highly correlated. For that reason, we may want to exclude dew point temperature in our model below.
"""

# ╔═╡ 7deee71e-18ec-11eb-113d-fb0cdc6e07ef
md"""
## Seasonal analysis
Since this is a time series, let's check if there's a trend.
"""

# ╔═╡ 32cd8006-18ec-11eb-221a-d38281aadffb
plot(df.Date, df.RentedBikeCount;
	legend = :none,
	ylabel = "Rented Bike Count")

# ╔═╡ a9f3aeda-18ec-11eb-2cb9-9f12fd001d88
md"""
As we only have 1 year of data, we cannot really see any seasonal trends. However, my intuition is that people don't tend to ride bikes during winter due to the low temperature. When spring comes around March/April, the demand picked up.

Now, December 2018 still look like a higher demand when compared with December 2017. A possibility is that biking may have gotten popular? Not sure...
"""

# ╔═╡ 69199458-18ef-11eb-1a46-7d71379a4c04
md"""
## Time of day

My initial guess is that people probably do not rent bikes early morning or late night. Let's verify that.
"""

# ╔═╡ 5de23772-18ed-11eb-1105-85661397bcd3
let
	data = @pipe df |>
		select(_, :Date, :Hour, :RentedBikeCount) |>
		unstack(_, :Hour, :RentedBikeCount) |>
		sort(_, :Date)

	matrix = Matrix(data[:, 2:end])'

	heatmap(data[:, 1], 0:23, matrix;
		title = "Rented Bike Count by Date / Hour",
		ylabel = "Hour",
		xrotation = 45.0,
		xticks = Date(2017,12,1):Month(1):Date(2019,1,1),
		yticks = 0:4:24)
end

# ╔═╡ 8789d48e-18ef-11eb-2b87-2337472759c7
md"""
**The two horizontal stripes show that 8 AM and 6 PM are popular hours for bike rentals.** That makes sense because people might want to do some exercise right before work or after work. Or, they may use it for transportation.

**Korean seems to stay up late.** Bike rental demand continues after 6 PM until almost midnight. By 3 AM, it's total silence.

Now, let's take a look at the same data from a different angle using a boxplot.
"""

# ╔═╡ 5ef37b5a-18eb-11eb-3f75-cfb5947950fb
begin
	hline([1800], color = :green)
	@df df StatsPlots.boxplot!(:Hour, :RentedBikeCount;
		legend = :none,
		color = palette(:default)[1],  # reset to first color due to `hline` above
		xlabel = "Hour",
		ylabel = "Rented Bike Count")
end

# ╔═╡ 3dc1127e-19a4-11eb-2c09-4f2624f479ac
md"""
The boxplot shows the quartiles for each hour across all dates in the year. It confirms our findings that 8 AM and 6 AM are popular times.
"""

# ╔═╡ 716cb212-19b1-11eb-190a-7702595f6bab
md"""
## Categorical Variables

Some variables should be converted to categorical such that GLM.jl can encode dummy variables as such. Instead of mutating the existing `Hour` column, I will create a new `Hour2` column and make it categorical.
"""

# ╔═╡ 79056460-19b6-11eb-1c5c-bf9224a0b18e
df.Hour2 = df.Hour

# ╔═╡ 9506548a-19b6-11eb-02dc-ad2832999a5b
categorical!(df, :Hour2);

# ╔═╡ 7553988e-19a5-11eb-3549-bfd4ea98b2fd
md"""
## Linear Regression

Using linear regression, we can fit the data to a linear equation. Here, the response variable is `RentedBikeCount`. We can choose any of the other fields as explanatory variables. Let's start with something simple - an ordinary linear model.
"""

# ╔═╡ c18fa04e-19a5-11eb-2f5c-0bdba322a26f
hour_model = @formula(RentedBikeCount ~ 1 + Hour2)

# ╔═╡ 558893b4-19a6-11eb-3e84-ef27b2749b80
md"""
To fit the model, we can use the `lm` function.
"""

# ╔═╡ 66b6f680-19a6-11eb-2d2f-afe8d206ee5f
ols = lm(hour_model, df)

# ╔═╡ a46a544a-19a6-11eb-3de3-87dcf7bf431f
md"""
We can now use the fitted model to predict bike demand.
"""

# ╔═╡ bfe39e20-19a6-11eb-1b9c-b95be62dadb5
let
	y = df.RentedBikeCount
	yhat = round.(Int, predict(ols))
	(r_squared = r2(ols), rmsd = rmsd(y, yhat))
end

# ╔═╡ 079fe824-19a7-11eb-3f1d-b3aa9c9b0f33
md"""
The $R^2$ is very low, meaning that the current model cannot make very accurate prediction. That's understandable because we have only used a single explanatory variable.

The root mean squared deviation (`rmsd`) value shows how much the predicted values deviates from the actual values.

Now, let's design a more complex model but we will continue to use an ordinary linear model.
"""

# ╔═╡ 41fb7ff4-19a7-11eb-0945-017db54aa411
let
	model = @formula(RentedBikeCount ~
						1 + Hour2 + Temperature + Humidity + WindSpeed +
						Visibility + SolarRadiation + Rainfall + Snowfall +
						Seasons + Holiday + FunctioningDay)
	fitted = lm(model, df)
	y = df.RentedBikeCount
	yhat = round.(Int, predict(fitted))
	(r_squared = r2(fitted), rmsd = rmsd(y, yhat))
end

# ╔═╡ e3615da0-19a7-11eb-3dc7-d734a6da4420
md"""
That's great result. The $R^2$ has jumped to 0.66 now! The RMSD is also reduced quite significantly from 543 to 375.
"""

# ╔═╡ ab20e2f6-19a9-11eb-09cd-eff9e75608ee
md"""
## Generalized Linear Model (GLM)

I wonder if we can do better. The GLM.jl package supports Generalized Linear Model (GLM) which is more flexible that linear regression. We will demonstrate how it works below.

Given that the response variable is a count, the general wisdom (not mine) is to design the model with Poisson distribution.
"""

# ╔═╡ 1accea7a-19a8-11eb-2e37-b943977ad801
let
	model = @formula(RentedBikeCount ~
						1 + Hour2 + Temperature + Humidity + WindSpeed +
						Visibility + SolarRadiation + Rainfall + Snowfall +
						Seasons + Holiday + FunctioningDay)
	fitted = glm(model, df, Poisson(), LogLink())
	y = df.RentedBikeCount
	yhat = round.(Int, predict(fitted))
	(rmsd = rmsd(y, yhat), )
end

# ╔═╡ 61eb2598-19a8-11eb-055d-b166da58fa50
md"""
There is no `r2` function defined for GLM models, so we just show RMSD here. As you can see, RMSD is further reduced using this model. That's an improvement.

What if we build an even more complex model? So far, all variables are independent. If we introduce interaction terms (multiple variables interacting with each other) then we may create a more powerful predictor.

The question is how to choose the right variables for the interaction terms. My gut feeling is that it would be appropriate to choose variables that are "orthogonal" to each other. For example, humidity and rainfall should be highly correlated and the interaction between them would be somewhat uninteresting. Hence, I have chosen to mix hour of day, temperature, and humidity in the following experiment.
"""

# ╔═╡ 66967c54-19b3-11eb-1b13-856ef40132fa
let
	model = @formula(RentedBikeCount ~
						1 + Hour2 + Temperature + Humidity + WindSpeed +
						Visibility + SolarRadiation + Rainfall + Snowfall +
						Seasons + Holiday + FunctioningDay +
						Hour2 * Temperature * Humidity
	)
	fitted = glm(model, df, Poisson(), LogLink())
	y = df.RentedBikeCount
	yhat = round.(Int, predict(fitted))
	(rmsd = rmsd(y, yhat), )
end

# ╔═╡ 1de7c9b2-19b4-11eb-2085-7b7383899819
md"""
That's great! The RMSD is now further reduced although not by a whole lot.
"""

# ╔═╡ 34be34aa-19b4-11eb-1a17-e1253c48b050
md"""
## Todo's

So far, I have been fitting the model with the complete data set. In order to test the predictive power of the model, I should test it against unseen data. Of course, I don't have any more data than what I have downloaded. What I should do is to split the data set and do cross validation.
"""

# ╔═╡ 544e7ac0-19b2-11eb-3931-8f7f10c23212
md"""
## Resources

I don't know much about GLM before working on this. I found the following resources useful as I learn about the subject:

- Foundations of Linear and Generalized Linear Models by Alan Agresti (Wiley 2015)
- [Introduction to Generalized Linear Models](https://online.stat.psu.edu/stat504/node/216/)
- [MIT 18.650 Statistics for Applications, Phillipe Rigolle, Lecture 21-22](https://www.youtube.com/watch?v=X-ix97pw0xY)
"""

# ╔═╡ 7c1251d4-19a9-11eb-0450-8dcae5fdc135
md"""
Thanks you for reading. I hope you enjoy this notebook.

_Tom Kwong_,
_October 2020_
"""

# ╔═╡ Cell order:
# ╟─02eb9d36-18e8-11eb-24a9-db6f58080189
# ╟─eedb4bf8-19a9-11eb-106b-69b047d4118e
# ╠═81915028-18e7-11eb-0aed-57a75839d4a7
# ╠═dfa5fbf0-18e7-11eb-391b-1d878b9c96f0
# ╠═aed964bc-18e7-11eb-11f1-a3d9957a45f5
# ╟─4d4a9158-18ec-11eb-2fb0-c178637f903a
# ╟─09577b1e-19aa-11eb-3e6c-175705b4b951
# ╠═c5b4e26a-18e7-11eb-2bef-2d937faa252d
# ╟─6ddeb944-18ec-11eb-05c7-b51e7b01578e
# ╠═ebc9c09a-18e9-11eb-1f56-ef92258af11f
# ╟─16c31048-19b6-11eb-2ad8-1f984acfb06f
# ╠═9564ab04-19b4-11eb-3c4c-815fe8a6031a
# ╠═3084d7a6-19b5-11eb-1fa5-9b9f7cd482f4
# ╟─2ef9d8a6-19b6-11eb-0cc5-e1831ce6c40a
# ╟─7deee71e-18ec-11eb-113d-fb0cdc6e07ef
# ╠═32cd8006-18ec-11eb-221a-d38281aadffb
# ╟─a9f3aeda-18ec-11eb-2cb9-9f12fd001d88
# ╟─69199458-18ef-11eb-1a46-7d71379a4c04
# ╠═5de23772-18ed-11eb-1105-85661397bcd3
# ╟─8789d48e-18ef-11eb-2b87-2337472759c7
# ╠═5ef37b5a-18eb-11eb-3f75-cfb5947950fb
# ╟─3dc1127e-19a4-11eb-2c09-4f2624f479ac
# ╟─716cb212-19b1-11eb-190a-7702595f6bab
# ╠═79056460-19b6-11eb-1c5c-bf9224a0b18e
# ╠═9506548a-19b6-11eb-02dc-ad2832999a5b
# ╟─7553988e-19a5-11eb-3549-bfd4ea98b2fd
# ╠═c18fa04e-19a5-11eb-2f5c-0bdba322a26f
# ╟─558893b4-19a6-11eb-3e84-ef27b2749b80
# ╠═66b6f680-19a6-11eb-2d2f-afe8d206ee5f
# ╟─a46a544a-19a6-11eb-3de3-87dcf7bf431f
# ╠═bfe39e20-19a6-11eb-1b9c-b95be62dadb5
# ╟─079fe824-19a7-11eb-3f1d-b3aa9c9b0f33
# ╠═41fb7ff4-19a7-11eb-0945-017db54aa411
# ╟─e3615da0-19a7-11eb-3dc7-d734a6da4420
# ╟─ab20e2f6-19a9-11eb-09cd-eff9e75608ee
# ╠═1accea7a-19a8-11eb-2e37-b943977ad801
# ╟─61eb2598-19a8-11eb-055d-b166da58fa50
# ╠═66967c54-19b3-11eb-1b13-856ef40132fa
# ╟─1de7c9b2-19b4-11eb-2085-7b7383899819
# ╟─34be34aa-19b4-11eb-1a17-e1253c48b050
# ╟─544e7ac0-19b2-11eb-3931-8f7f10c23212
# ╟─7c1251d4-19a9-11eb-0450-8dcae5fdc135

max_stars_repo_name: Humans-of-Julia/Challenges

max_issues_repo_name: Humans-of-Julia/WeeklyContest

max_forks_repo_name: Humans-of-Julia/WeeklyContest


text:
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MathOptInterface.Test
const MOIU = MathOptInterface.Utilities
const MOIB = MathOptInterface.Bridges

include("utilities.jl")

include("simple_model.jl")

mock = MOIU.MockOptimizer(SimpleModel{Float64}())
config = MOIT.TestConfig()
config_with_basis = MOIT.TestConfig(basis = true)

@testset "Scalar functionize" begin
    MOI.empty!(mock)
    bridged_mock = MOIB.ScalarFunctionize{Float64}(mock)
    x = MOI.add_variable(bridged_mock)
    y = MOI.add_variable(bridged_mock)
    sx = MOI.SingleVariable(x)
    sy = MOI.SingleVariable(y)
    ci = MOI.add_constraint(bridged_mock, sx, MOI.GreaterThan(0.0))
    @test MOI.get(bridged_mock, MOI.ConstraintFunction(), ci) ≈ sx
    MOI.set(bridged_mock, MOI.ConstraintFunction(), ci, sy)
    @test MOI.get(bridged_mock, MOI.ConstraintFunction(), ci) ≈ sy
    @test MOI.get(bridged_mock, MOI.ConstraintSet(), ci) == MOI.GreaterThan(0.0)
    MOI.set(bridged_mock, MOI.ConstraintSet(), ci, MOI.GreaterThan(1.0))
    @test MOI.get(bridged_mock, MOI.ConstraintSet(), ci) == MOI.GreaterThan(1.0)
    test_delete_bridge(bridged_mock, ci, 2, ((MOI.ScalarAffineFunction{Float64}, MOI.Interval{Float64}, 0),))

    MOIT.basic_constraint_tests(bridged_mock, config,
                                include=[(F, S) for
                                F in [MOI.SingleVariable],
                                S in [MOI.GreaterThan{Float64},
                                    MOI.LessThan{Float64}]
                                    ])

    for T in [Int, Float64], S in [MOI.GreaterThan{T}, MOI.LessThan{T}]
        @test MOIB.added_constraint_types(MOIB.ScalarFunctionizeBridge{T, S}) ==
            [(MOI.ScalarAffineFunction{T}, S)]
    end

    @testset "linear2" begin
        MOI.empty!(bridged_mock)
        MOIU.set_mock_optimize!(mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1, 0],
                 (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64})    => [-1],
                 (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [0, 1],
            con_basis =
                [(MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64})    => [MOI.NONBASIC],
                 (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [MOI.BASIC, MOI.NONBASIC]]))
        MOIT.linear2test(bridged_mock, config_with_basis)

        for (i, ci) in enumerate(MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{MOI.SingleVariable,
                                        MOI.GreaterThan{Float64}}()))
            test_delete_bridge(bridged_mock, ci, 2,
                               ((MOI.ScalarAffineFunction{Float64},
                                 MOI.GreaterThan{Float64}, 0),),
                               num_bridged = 3 - i)
        end
    end
end

@testset "Vector functionize" begin
    MOI.empty!(mock)
    bridged_mock = MOIB.VectorFunctionize{Float64}(mock)
    x = MOI.add_variable(bridged_mock)
    y = MOI.add_variable(bridged_mock)
    z = MOI.add_variable(bridged_mock)
    v1 = MOI.VectorOfVariables([x,y,z])
    v2 = MOI.VectorOfVariables([x,y,y])
    ci = MOI.add_constraint(bridged_mock, v1, MOI.PowerCone(0.1))
    @test MOI.get(bridged_mock, MOI.ConstraintFunction(), ci) ≈ v1
    MOI.set(bridged_mock, MOI.ConstraintFunction(), ci, v2)
    @test MOI.get(bridged_mock, MOI.ConstraintFunction(), ci) ≈ v2
    @test MOI.get(bridged_mock, MOI.ConstraintSet(), ci) == MOI.PowerCone(0.1)
    MOI.set(bridged_mock, MOI.ConstraintSet(), ci, MOI.PowerCone(0.2))
    @test MOI.get(bridged_mock, MOI.ConstraintSet(), ci) == MOI.PowerCone(0.2)
    test_delete_bridge(bridged_mock, ci, 3, ((MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives, 0),))

    MOIT.basic_constraint_tests(bridged_mock, config,
                                include=[(F, S) for
                                F in [MOI.VectorOfVariables],
                                S in [MOI.Nonnegatives,
                                    MOI.Nonpositives]
                                    ])

    for T in [Int, Float64], S in [MOI.Nonnegatives, MOI.Nonpositives]
        @test MOIB.added_constraint_types(MOIB.VectorFunctionizeBridge{T, S}) ==
            [(MOI.VectorAffineFunction{T}, S)]
    end

    @testset "lin1v" begin
        mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock, [1.0, 0.0, 2.0],
            (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) => [[0, 2, 0]],
            (MOI.VectorAffineFunction{Float64}, MOI.Zeros)        => [[-3, -1]])
        MOIT.lin1vtest(bridged_mock, config)
        ci = first(MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{MOI.VectorOfVariables,
                                        MOI.Nonnegatives}()))
        test_delete_bridge(bridged_mock, ci, 3,
                           ((MOI.VectorAffineFunction{Float64},
                             MOI.Nonnegatives, 0),))
    end
end


text:
function Base.display(p::CmdDependency)
	print("CmdDependency\n  exec             →")
	display(p.exec)
	print("  test_args        →")
	display(p.test_args)
	println("  validate_success → $(p.validate_success)\n  validate_stdout  → $(p.validate_stdout)\n  validate_stderr  → $(p.validate_stderr)\n  exit_when_fail   → $(p.exit_when_fail)")
end

function Base.print(io::IO, p::CmdDependency)
	print(io, p.exec)
end

function Base.show(io::IO, p::CmdDependency)
	show(io, p.exec)
end

@eval function Base.display(p::CmdProgram)
    fs = $(fieldnames(CmdProgram))
    fs_string = $(map(string, fieldnames(CmdProgram)))
    max_byte = $(maximum(length, map(string, fieldnames(CmdProgram))))
    println("$CmdProgram:")
    for (i,f) in enumerate(fs)
		if f == :inputs
			print("  ", f, " " ^ (max_byte - length(fs_string[i])), " → ")
			display_xxputs(max_byte, p.inputs, p.input_types, p.default_inputs)
		elseif f == :outputs
			print("  ", f, " " ^ (max_byte - length(fs_string[i])), " → ")
			display_xxputs(max_byte, p.outputs, p.output_types, p.default_outputs)
		elseif f in Symbol[:input_types, :default_inputs, :output_types, :default_outputs]
		else
			print("  ", f, " " ^ (max_byte - length(fs_string[i])), " → ")
	        println(getfield(p, f))
		end
    end
end

@eval function Base.display(p::JuliaProgram)
    fs = $(fieldnames(JuliaProgram))
    fs_string = $(map(string, fieldnames(JuliaProgram)))
    max_byte = $(maximum(length, map(string, fieldnames(JuliaProgram))))
    println("$JuliaProgram:")
    for (i,f) in enumerate(fs)
		if f == :inputs
			print("  ", f, " " ^ (max_byte - length(fs_string[i])), " → ")
			display_xxputs(max_byte, p.inputs, p.input_types, p.default_inputs)
		elseif f == :outputs
			print("  ", f, " " ^ (max_byte - length(fs_string[i])), " → ")
			display_xxputs(max_byte, p.outputs, p.output_types, p.default_outputs)
		elseif f in Symbol[:input_types, :default_inputs, :output_types, :default_outputs]
		else
			print("  ", f, " " ^ (max_byte - length(fs_string[i])), " → ")
	        println(getfield(p, f))
		end
    end
end

function display_xxputs(max_bype::Int, xxputs::Vector{String}, xxput_types::Vector{DataType}, default_xxputs::Vector)
	n = length(xxputs)
	if n == 0
		println("<empty>")
		return
	end
	max_bype_xxputs = maximum(length, xxputs)
	max_bype_xxput_types = maximum(length, map(string, xxput_types))
	for i in 1:n
		name = xxputs[i]
		type = string(xxput_types[i])
		default = default_xxputs[i]
		println(i == 1 ? "" : " " ^ (max_bype + 5),
			"\"$(name)\" ", " " ^ (max_bype_xxputs - length(name)),
			":: $(type) ", " " ^ (max_bype_xxput_types - length(type)),
			isnothing(default) ? "(required)" : "(default: $default)"
		)
	end
	return
end


text:
using Test, SpatioTemporalSystems, DrWatson
# Test:

sts = STS{:bk}(N = 100, T = 1000, Δt = 0.1)
d = makesim(sts)

U = d["U"]

# using PyPlot
# figure()
# imshow(U[16])

# %% Henon map
sts = STS{:henon1d}(ic = 4, N = 200, T = 400, S = 1000)
d = makesim(sts)
u = d["U"]
#
# using PyPlot
# figure()
# U = hcat(u...)
# imshow(U)

# %% Kuramoto
L = 200  # physical size of system
N = 512  # number of gridpoints
Δt = 1/4
T = 100.0
x = L*(0:N-1)/N # position vector
u0 = cos.(x) + 0.1*sin.(x/8) + 0.01*cos.((2*pi/L)*x);

sts = STS{:ksiva}(N = N, ic = u0, T = T, Δt = Δt, p = (L = L,))
savename(sts)

sim = makesim(sts)
u = sim["u"]
# @unpack u, t = sim

# using PyPlot
# ax = gca()
# ax.clear()
# U = hcat(u...)
# ax.imshow(U, extent = (0, T, 0, L), cmap = "inferno")
# ax.set_aspect("auto")
# ax.set_ylabel("x")
# ax.set_xlabel("t")
# gcf().tight_layout()


text:
abstract type AbstractScheduler end

struct StepDecay <: AbstractScheduler
    xmax
    xmin
    Δ
    T
end

(d::StepDecay)(t) = max(d.xmin, d.xmax - div(t, d.T) * d.Δ)

d = StepDecay(1.0, 0.1, 0.1, 5)
plot(1:100, d.(1:100))

struct ExponentialDecay <: AbstractScheduler
    xmax
    xmin
    ρ
end

(d::ExponentialDecay)(t) = max(d.xmin, d.xmax * d.ρ^t)

d = ExponentialDecay(1.0, 0.1, 0.99)
plot(1:250, d.(1:250))


text:
using Revise
using AdFem
using PyCall
using LinearAlgebra
using PyPlot
using SparseArrays
using MAT
np = pyimport("numpy")

# mode = "data" generate data
# mode != "data" train
# mode = "train"
mode = "data"

## alpha-scheme
β = 1/4; γ = 1/2
a = b = 0.1

m = 20
n = 10
h = 0.01
NT = 500
it0 = NT÷2
Δt = 20/NT
ηmax = 1
ηmin = 0.5
bdedge = []
for j = 1:n
  push!(bdedge, [(j-1)*(m+1)+m+1 j*(m+1)+m+1])
end
bdedge = vcat(bdedge...)

bdnode = Int64[]
for j = 1:n+1
  push!(bdnode, (j-1)*(m+1)+1)
end

# λ = Variable(1.0)
# μ = Variable(1.0)
# invη = Variable(1.0)

function eta_model(idx)
  if idx == 1
    out = ηmin * ones(4, m, n)
    out[:, :, 1:div(n,3)] .= ηmax
    # out[:, 1:div(m,3), :] .= 12.0
    # out[:, :, div(n,3)*2:end] .= 8.0
    out[:]
  elseif idx == 2
    out = ones(4, m, n)
    for i = 1:m
        for j = 1:n
            x = i*h; y = j*h
            out[1,i,j] = (x^2 + y^2)
            out[2,i,j] = (x^2 + y^2)
            out[3,i,j] = (x^2 + y^2)
            out[4,i,j] = (x^2 + y^2)
        end
    end
    out[:]
  end
end

function visualize_inv_eta(X, k)
    x = LinRange(0.5h,m*h, m)
    y = LinRange(0.5h,n*h, n)
    V = zeros(m, n)
    for i = 1:m
        for j = 1:n
            elem = (j-1)*m + i
            V[i, j] = mean(X[4(elem-1)+1:4elem])
        end
    end
    close("all")
    pcolormesh(x, y, V', vmin=ηmin-(ηmax-ηmin)/4, vmax=ηmax+(ηmax-ηmin)/4)
    colorbar(shrink=0.5)
    xlabel("x")
    ylabel("y")
    # title("Iteration = $k")
    title("True model")
    axis("scaled")
    gca().invert_yaxis()
    savefig("iter$k.png")
end

λ = constant(2.0)
μ = constant(0.5)
if mode=="data"
  global invη = constant(eta_model(1))
else
    # global invη_ = Variable(10.0*ones(4*m*n))

    # invη1d_ = Variable(10*ones(2m))
    # invη_ = repeat(invη1d_, 1, 2)[:]
    # global invη_ = repeat(invη_, n)

    invη_ = Variable((ηmin + ηmax)/2*ones(n))
    # invη1d_ = repeat(invη1d_, 1, 2)[:]
    invη_ = tf.reshape(repeat(invη_, 1, 4m), (-1,))
    # invη_ = repeat(invη_, m)
    # global invη_ = tf.transpose(invη2d_)[:]
    global invη_ = invη_

    # invη1d_ = Variable(10*ones(1,2m))
    # invη_ = repeat(invη1d_, 2, 1)[:]
    # global invη_ = repeat(invη_, n, 1)

    # invη1d_ = Variable(10*ones(1,2n))
    # invη_ = repeat(invη1d_, 2, 1)[:]
    # global invη_ = repeat(invη_, m, 1)

    # global invη = repeat(invη_, 1, 4)[:]
    # global invη = [invη_[1] * ones(2*m*n); invη_[2] * ones(2*m*n)]

    # global invη_ = Variable(8.0)*ones(4*m*n)
    global invη = invη_
end

fn_G = invη->begin
  G = tensor([1/Δt+μ*invη -μ/3*invη 0.0
    -μ/3*invη 1/Δt+μ*invη-μ/3*invη 0.0
    0.0 0.0 1/Δt+μ*invη])
  invG = inv(G)
end
invG = map(fn_G, invη)
S = tensor([2μ/Δt+λ/Δt λ/Δt 0.0
    λ/Δt 2μ/Δt+λ/Δt 0.0
    0.0 0.0 μ/Δt])

# error()
# invG = Variable([1.0 0.0 0.0
#                 0.0 1.0 0.0
#                 0.0 0.0 1.0])
H = invG*S

M = compute_fem_mass_matrix1(m, n, h)
Zero = spzeros((m+1)*(n+1), (m+1)*(n+1))
M = SparseTensor([M Zero;Zero M])

K = compute_fem_stiffness_matrix(H, m, n, h)
C = a*M + b*K # damping matrix

L = M + γ*Δt*C + β*Δt^2*K
L, Lbd = fem_impose_Dirichlet_boundary_condition(L, bdnode, m, n, h)
# error()

a = TensorArray(NT+1); a = write(a, 1, zeros(2(m+1)*(n+1))|>constant)
v = TensorArray(NT+1); v = write(v, 1, zeros(2(m+1)*(n+1))|>constant)
d = TensorArray(NT+1); d = write(d, 1, zeros(2(m+1)*(n+1))|>constant)
U = TensorArray(NT+1); U = write(U, 1, zeros(2(m+1)*(n+1))|>constant)
Sigma = TensorArray(NT+1); Sigma = write(Sigma, 1, zeros(4*m*n, 3)|>constant)
Varepsilon = TensorArray(NT+1); Varepsilon = write(Varepsilon, 1,zeros(4*m*n, 3)|>constant)

Forces = zeros(NT, 2(m+1)*(n+1))
for i = 1:NT
  T = eval_f_on_boundary_edge((x,y)->0.1, bdedge, m, n, h)
  T = [-T zeros(length(T))]
#   T = [T T]
  rhs = compute_fem_traction_term(T, bdedge, m, n, h)

#   if i*Δt>0.5
#     rhs = zero(rhs)
#   end
  Forces[i, :] = rhs
end
Forces = constant(Forces)

function condition(i, tas...)
  i <= NT
end

function body(i, tas...)
  a_, v_, d_, U_, Sigma_, Varepsilon_ = tas
  a = read(a_, i)
  v = read(v_, i)
  d = read(d_, i)
  U = read(U_, i)
  Sigma = read(Sigma_, i)
  Varepsilon = read(Varepsilon_, i)

  # Sigma * (invG/Δt) Sigma 800x3, invG/Δt: 800x3x3
  res = squeeze(tf.matmul(tf.reshape(Sigma, (size(Sigma,1), 1, 3)),(invG/Δt)))
  F = compute_strain_energy_term(res, m, n, h) - K * U
  rhs = Forces[i] - F

  td = d + Δt*v + Δt^2/2*(1-2β)*a
  tv = v + (1-γ)*Δt*a
  rhs = rhs - C*tv - K*td

  # rhs[[bdnode; bdnode.+(m+1)*(n+1)]] .= 0.0

  rhs = scatter_update(rhs, constant([bdnode; bdnode.+(m+1)*(n+1)]), constant(zeros(2*length(bdnode))))

  ## alpha-scheme
  a = L\rhs # bottleneck
  d = td + β*Δt^2*a
  v = tv + γ*Δt*a
  U_new = d
  Varepsilon_new = eval_strain_on_gauss_pts(U_new, m, n, h)

  res2 = squeeze(tf.matmul(tf.reshape(Varepsilon_new-Varepsilon, (size(Sigma,1), 1, 3)),
                        tf.matmul(invG, S)))
  Sigma_new = res +  res2

  i+1, write(a_, i+1, a), write(v_, i+1, v), write(d_, i+1, d), write(U_, i+1, U_new),
        write(Sigma_, i+1, Sigma_new), write(Varepsilon_, i+1, Varepsilon_new)
end

i = constant(1, dtype=Int32)
_, _, _, _, u, sigma, varepsilon = while_loop(condition, body,
                  [i, a, v, d, U, Sigma, Varepsilon])

U = stack(u)
Sigma = stack(sigma)
Varepsilon = stack(varepsilon)

if mode!="data"
    data = matread("U.mat")
  global Uval,Sigmaval, Varepsilonval = data["U"], data["Sigma"], data["Varepsilon"]
  U.set_shape((NT+1, size(U, 2)))
  idx0 = 1:4m*n
  Sigma = map(x->x[idx0,:], Sigma)

  idx = 1:m+1
  # idx = [idx; idx .+ (m+1)*(n+1)]
  # global loss = sum((U[:,idx] - Uval[:,idx])^2)  + sum((Sigma - Sigmaval[:, idx0, :])^2) #+ sum((Varepsilon - Varepsilonval)^2)
  # global loss = sum((U - Uval)^2)  + sum((Sigma - Sigmaval[:, idx0, :])^2) #+ sum((Varepsilon - Varepsilonval)^2)
  # global loss = sum((U[:,idx] - Uval[:,idx])^2) #+ sum((Varepsilon - Varepsilonval)^2)
  # global loss = sum((U - Uval)^2)
  # global loss = sum((U[it0:end, :] - Uval[it0:end, :])^2)
  global loss = sum((U[it0:end, idx] - Uval[it0:end, idx])^2)
  # global loss = sum((Sigma - Sigmaval[:, idx0, :])^2)
end

# opt = AdamOptimizer(0.1).minimize(loss)
sess = Session(); init(sess)
cb = (v, i, l)->begin
  println("[$i] loss = $l")
  inv_eta = v[1]
  visualize_inv_eta(inv_eta, i)
end
if mode=="data"
    Uval,Sigmaval, Varepsilonval = run(sess, [U, Sigma, Varepsilon])
  matwrite("U.mat", Dict("U"=>Uval, "Sigma"=>Sigmaval, "Varepsilon"=>Varepsilonval))
  # visualize_displacement(U, m, n, h; name = "_eta$η", xlim_=[-0.01,0.5], ylim_=[-0.05,0.22])
    # # visualize_displacement(U, m, n, h;  name = "_viscoelasticity")
    # # visualize_stress(H, U, m, n, h;  name = "_viscoelasticity")

    close("all")

    figure(figsize=(15,5))
    subplot(1,3,1)
    idx = div(n,2)*(m+1) + m+1
    # idx = (m+1)÷2
    # plot((it0-1:NT)*Δt, Uval[it0:end, idx])
    plot((0:NT)*Δt, Uval[:, idx])
    xlabel("time")
    ylabel("\$u_x\$")
    plt.gca().ticklabel_format(style="sci", scilimits=(0,0), axis="y")

    ax = plt.gca().inset_axes([0.3, 0.2, 0.6, 0.7])
    ax.plot((it0-1:NT)*Δt, Uval[it0:end, idx])
    ax.ticklabel_format(style="sci", scilimits=(0,0), axis="y")

    subplot(1,3,2)
    idx = 4*(div(n,2)*m + m)
    # idx = (m+1)*2
    # plot((it0-1:NT)*Δt, Sigmaval[it0:end,idx,1])
    plot((0:NT)*Δt, Sigmaval[:,idx,1])
    xlabel("time")
    ylabel("\$\\sigma_{xx}\$")
    plt.gca().ticklabel_format(style="sci", scilimits=(0,0), axis="y")

    ax = plt.gca().inset_axes([0.4, 0.1, 0.5, 0.5])
    ax.plot((it0-1:NT)*Δt, Sigmaval[it0:end,idx,1])
    ax.ticklabel_format(style="sci", scilimits=(0,0), axis="y")

    subplot(1,3,3)
    idx = 4*(div(n,2)*m + m)
    # idx = (m+1)*2
    # plot((it0-1:NT)*Δt, Varepsilonval[it0:end,idx,1])
    plot((0:NT)*Δt, Varepsilonval[:,idx,1])
    xlabel("time")
    ylabel("\$\\varepsilon_{xx}\$")
    plt.gca().ticklabel_format(style="sci", scilimits=(0,0), axis="y")

    ax = plt.gca().inset_axes([0.3, 0.2, 0.6, 0.7])
    ax.plot((it0-1:NT)*Δt, Varepsilonval[it0:end,idx,1])
    ax.ticklabel_format(style="sci", scilimits=(0,0), axis="y")

    savefig("visco_eta.png")

    cb([run(sess, invη)], "true", 0)
  error()
end

@info run(sess, loss)
v_ = []
i_ = []
l_ = []

loss_ = BFGS!(sess, loss*1e10, vars=[invη], callback=cb, var_to_bounds=Dict(invη_=>(5.0,15.0)))
# matwrite("R.mat", Dict("V"=>v_, "L"=>l_))

# for i = 1:1000
#   _, l, invη_ = run(sess, [opt, loss, invη])
#   @show i, l #, invη_
#   mod(i,5)==1 && cb([invη_], i, l)
# end

# ηs = 11:0.1:13
# losses = []
# for η in ηs
#    push!(losses,run(sess, loss, invη=>η))
# end
# plot(ηs, losses)

# Uval, Sigmaval, Varepsilonval = run(sess, [U, Sigma, Varepsilon])
# Uval[idx]