text:
"""
    Gain(k=1; name)

Output the product of a gain value with the input signal.

# Parameters:
- `k`: Scalar gain
"""
function Gain(k=1; name)
    @named siso = SISO()
    @unpack u, y = siso
    pars = @parameters k=k
    eqs = [
        y ~ k * u
    ]
    extend(ODESystem(eqs, t, [], pars; name=name), siso)
end

"""
    MatrixGain(K::AbstractArray; name)

Output the product of a gain matrix with the input signal vector.

# Parameters:
- `K`: Matrix gain
"""
function MatrixGain(K::AbstractArray; name)
    nout, nin = size(K)
    @named input = RealInput(;nin=nin)
    @named output = RealOutput(;nout=nout)
    eqs = [
        output.u[i] ~ sum(K[i,j] * input.u[j] for j in 1:nin) for i in 1:nout # FIXME: if array equations work
    ]
    compose(ODESystem(eqs, t, [], []; name=name), [input, output])
end

"""
    Sum(n::Int; name)

Output the sum of the elements of the input port vector.

# Parameters:
- `n`: Input port dimension
"""
function Sum(n::Int; name)
    @named input = RealInput(;nin=n)
    @named output = RealOutput()
    eqs = [
        output.u ~ sum(input.u)
    ]
    compose(ODESystem(eqs, t, [], []; name=name), [input, output])
end

"""
    Feedback(;name)

Output difference between reference input (input1) and feedback input (input2).
"""
function Feedback(;name)
    @named input1 = RealInput()
    @named input2 = RealInput()
    @named output = RealOutput()
    eqs= [
        output.u ~ input1.u - input2.u
    ]
    return compose(ODESystem(eqs, t, [], []; name=name), input1, input2, output)
end

"""
    Add(;name, k1=1, k2=1)

Output the sum of the two scalar inputs.

# Parameters:
- `k1`: Gain for first input
- `k2`: Gain for second input
"""
function Add(;name, k1=1, k2=1)
    @named input1 = RealInput()
    @named input2 = RealInput()
    @named output = RealOutput()
    pars = @parameters begin
        k1=k1
        k2=k2
    end
    eqs= [
        output.u ~ k1 * input1.u + k2 * input2.u
    ]
    return compose(ODESystem(eqs, t, [], pars; name=name), input1, input2, output)
end

"""
    Add(;name, k1=1, k2=1,k3=1)

Output the sum of the three scalar inputs.

# Parameters:
- `k1`: Gain for first input
- `k2`: Gain for second input
- `k3`: Gain for third input
"""
function Add3(;name, k1=1, k2=1, k3=1)
    @named input1 = RealInput()
    @named input2 = RealInput()
    @named input3 = RealInput()
    @named output = RealOutput()
    pars = @parameters begin
        k1=k1
        k2=k2
        k3=k3
    end
    eqs= [
        output.u ~ k1 * input1.u + k2 * input2.u + k3 * input3.u
    ]
    return compose(ODESystem(eqs, t, [], pars; name=name), input1, input2, input3, output)
end

"""
    Product(;name)

Output product of the two inputs.
"""
function Product(;name)
    @named input1 = RealInput()
    @named input2 = RealInput()
    @named output = RealOutput()
    eqs= [
        output.u ~ input1.u * input2.u
    ]
    return compose(ODESystem(eqs, t, [], []; name=name), input1, input2, output)
end

"""
    Division(;name)

Output first input divided by second input.
"""
function Division(;name)
    @named input1 = RealInput()
    @named input2 = RealInput()
    @named output = RealOutput()
    eqs= [
        output.u ~ input1.u / input2.u
    ]
    return compose(ODESystem(eqs, t, [], []; name=name), input1, input2, output)
end

"""
    StaticNonLinearity(func ;name)

Applies the given function to the input.

If the given function is not composed of simple core methods (e.g. sin, abs, ...), it has to be registered via `@register_symbolic func(u)`
"""
function StaticNonLinearity(func; name)
    @named siso = SISO()
    @unpack u, y = siso
    eqs = [y ~ func(u)]
    extend(ODESystem(eqs, t, [], []; name=name), siso)
end

"""
    Abs(;name)

Output the absolute value of the input.
"""
Abs(;name) = StaticNonLinearity(abs; name)

"""
    Sign(;name)

Output the sign of the input
"""
Sign(;name) = StaticNonLinearity(sign; name)

"""
    Sqrt(;name)

Output the square root of the input (input >= 0 required).
"""
Sqrt(;name) = StaticNonLinearity(sqrt; name)

"""
    Sin(;name)

Output the sine of the input.
"""
Sin(;name) = StaticNonLinearity(sin; name)

"""
    Cos(;name)

Output the cosine of the input.
"""
Cos(;name) = StaticNonLinearity(cos; name)

"""
    Tan(;name)

Output the tangent of the input.
"""
Tan(;name) = StaticNonLinearity(tan; name)

"""
    Asin(;name)

Output the arc sine of the input.
"""
Asin(;name) = StaticNonLinearity(asin; name)

"""
    Acos(;name)

Output the arc cosine of the input.
"""
Acos(;name) = StaticNonLinearity(acos; name)

"""
    Atan(;name)

Output the arc tangent of the input.
"""
Atan(;name) = StaticNonLinearity(atan; name)

"""
    Atan2(;name)

Output the arc tangent of the input.
"""
function Atan2(;name)
    @named input1 = RealInput()
    @named input2 = RealInput()
    @named output = RealOutput()
    eqs = [
        output.u ~ atan(input1.u, input2.u)
    ]
    compose(ODESystem(eqs, t, [], []; name=name), [input1, input2, output])
end

"""
    Sinh(;name)

Output the hyperbolic sine of the input.
"""
Sinh(;name) = StaticNonLinearity(sinh; name)

"""
    Cosh(;name)

Output the hyperbolic cosine of the input.
"""
Cosh(;name) = StaticNonLinearity(cosh; name)

"""
    Tanh(;name)

Output the hyperbolic tangent of the input.
"""
Tanh(;name) = StaticNonLinearity(tanh; name)

"""
    Exp(;name)

Output the exponential (base e) of the input.
"""
Exp(;name) = StaticNonLinearity(exp; name)

"""
    Log(;name)

Output the natural (base e) logarithm of the input.
"""
Log(;name) = StaticNonLinearity(log; name)

"""
    Log10(;name)

Output the base 10 logarithm of the input.
"""
Log10(;name) = StaticNonLinearity(log10; name)


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions related to searching text.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

const _search_highlighting = Dict{Bool, Decoration}(
    false => Decoration(foreground = "30", background = "47"),
    true  => Decoration(foreground = "30", background = "43")
)

"""
    reset_highlighting()

Reset the search highlighting to the default one.

"""
function reset_highlighting()
    _search_highlighting[false] = Decoration(
        foreground = "30",
        background = "47"
    )

    _search_highlighting[true] = Decoration(
        foreground = "30",
        background = "43"
    )

    return nothing
end

"""
    _change_active_match!(pagerd::Pager, forward::Bool = true)

Change the active matches in `pagerd`. If `forward` is `true`, then the search
is performed forward. Otherwise, it is performed backwards.

"""
function _change_active_match!(pagerd::Pager, forward::Bool = true)
    search_matches         = pagerd.search_matches
    active_search_match_id = pagerd.active_search_match_id

    num_matches = length(search_matches)

    if num_matches == 0
        active_number_match = 0
    else
        # Activate the next match according to the user preference.
        if forward
            active_search_match_id += 1
        else
            active_search_match_id -= 1
        end

        if active_search_match_id > num_matches
            active_search_match_id = 1
        elseif active_search_match_id < 1
            active_search_match_id = num_matches
        end
    end

    pagerd.active_search_match_id = active_search_match_id

    return nothing
end

"""
    _find_matches!(pagerd::Pager, regex::Regex)

Find all matches of `regex` in the text of the pager `pager`. The vector with
the matches will be written to `pagerd`.

"""
function _find_matches!(pagerd::Pager, regex::Regex)
    lines          = pagerd.lines
    num_lines      = pagerd.num_lines
    search_matches = pagerd.search_matches

    # Reset the previous search.
    empty!(search_matches)
    pagerd.active_search_match_id = 0

    # Regex to remove the ANSI escape sequence.
    regex_ansi = r"\x1B(?:[@-Z\\-_]|\[[0-?]*[ -/]*[@-~])"

    # For each line, find matches based on regex.
    for i = 1:num_lines
        # We need to filter the escape sequences from the line before searching.
        # TODO: Should we maintain a version of the input without the escape to
        # improve performance?

        line = String(*(split(lines[i], regex_ansi)...))

        matches_i = eachmatch(regex, line)

        for m in matches_i
            # `m.offset` contains the byte in which the match starts. However,
            # we need to obtain the character. Hence, it is necessary to compute
            # the text width from the beginning to the offset.
            push!(search_matches,
                  (i, textwidth(line[1:m.offset]), textwidth(m.match))
            )
        end
    end

    return nothing
end

"""
    _move_view_to_match!(pagerd::Pager)

Move the view of the pager `pagerd` to ensure that the current highlighted match
is inside it.

"""
function _move_view_to_match!(pagerd::Pager)
    # Unpack.
    start_row              = pagerd.start_row
    start_col              = pagerd.start_col
    active_search_match_id = pagerd.active_search_match_id
    search_matches         = pagerd.search_matches
    freeze_columns         = pagerd.freeze_columns
    freeze_rows            = pagerd.freeze_rows
    title_rows             = pagerd.title_rows

    rows, cols = _get_pager_display_size(pagerd)

    # Compute the last row and columns that is displayed.
    end_row = (start_row - 1) + (rows - freeze_rows)
    end_col = start_col + (cols - freeze_columns)

    # Get the active match.
    hl_i = active_search_match_id
    hl_i == 0 && return nothing

    # Get the position of the highlight.
    m = search_matches[hl_i]
    hl_line = m[1]
    hl_col_beg = m[2]
    hl_col_end = hl_col_beg + m[3] - 1

    # Check if the highlight row is visible.
    if (hl_line < start_row)
        start_row = max(hl_line, freeze_rows + 1)
    elseif hl_line > end_row
        start_row = (hl_line + 1) - (rows - freeze_rows)
    end

    # If the highlight is outsidde the title rows, then we can move the view to
    # display it.
    if title_rows < hl_line
        # Check if the highlight column is visible.
        if hl_col_beg < start_col
            start_col = hl_col_beg
        elseif hl_col_end > end_col
            start_col = (hl_col_end + 1) - (cols - freeze_columns)
        end
    end

    pagerd.start_row = start_row
    pagerd.start_col = start_col

    return nothing
end

"""
    _quit_search!(pagerd::Pager)

Quit search mode of pager `pagerd`.

"""
function _quit_search!(pagerd::Pager)
    empty!(pagerd.search_matches)
    pagerd.active_search_match_id = 0
    return nothing
end

ext: jl


text:
# ---
# title: 778. Swim in Rising Water
# id: problem778
# author: Indigo
# date: 2021-01-28
# difficulty: Hard
# categories: Binary Search, Heap, Depth-first Search, Union Find
# link: <https://leetcode.com/problems/swim-in-rising-water/description/>
# hidden: true
# ---
#
# On an N x N `grid`, each square `grid[i][j]` represents the elevation at that
# point `(i,j)`.
#
# Now rain starts to fall. At time `t`, the depth of the water everywhere is
# `t`. You can swim from a square to another 4-directionally adjacent square if
# and only if the elevation of both squares individually are at most `t`. You
# can swim infinite distance in zero time. Of course, you must stay within the
# boundaries of the grid during your swim.
#
# You start at the top left square `(0, 0)`. What is the least time until you
# can reach the bottom right square `(N-1, N-1)`?
#
# **Example 1:**
#
#
#
#     Input: [[0,2],[1,3]]
#     Output: 3
#     Explanation:
#     At time 0, you are in grid location (0, 0).
#     You cannot go anywhere else because 4-directionally adjacent neighbors have a higher elevation than t = 0.
#
#     You cannot reach point (1, 1) until time 3.
#     When the depth of water is 3, we can swim anywhere inside the grid.
#
#
# **Example 2:**
#
#
#
#     Input: [[0,1,2,3,4],[24,23,22,21,5],[12,13,14,15,16],[11,17,18,19,20],[10,9,8,7,6]]
#     Output: 16
#     Explanation:
#     **0  1  2  3  4**
#     24 23 22 21  **5**
#     **12 13 14 15 16**
#     **11** 17 18 19 20
#     **10  9  8  7  6**
#
#     The final route is marked in bold.
#     We need to wait until time 16 so that (0, 0) and (4, 4) are connected.
#
#
# **Note:**
#
#   1. `2 <= N <= 50`.
#   2. grid[i][j] is a permutation of [0, ..., N*N - 1].
#
#
## @lc code=start
using LeetCode

function swim_in_water(grid::Vector{Vector{Int}})
    function neighbors(pos, len)
        res = Tuple{Int, Int}[]
        (pos[1] > 1) && (push!(res, (pos[1] - 1, pos[2])))
        (pos[1] < len) && (push!(res, (pos[1] + 1, pos[2])))
        (pos[2] > 1) && (push!(res, (pos[1], pos[2] - 1)))
        (pos[2] < len) && (push!(res, (pos[1], pos[2] + 1)))
        res
    end
    len = length(grid)
    visited = fill(false, len, len)
    pq = PriorityQueue{Tuple{Int, Int}, Int}()
    pq[(1, 1)] = grid[1][1]
    while !isempty(pq)
        pos, cost = dequeue_pair!(pq)
        visited[pos[1], pos[2]] = true
        if pos == (len, len)
            return cost
        end
        for neighbor in neighbors(pos, len)
            if !visited[neighbor[1], neighbor[2]]
                pq[neighbor] = min(max(grid[neighbor[1]][neighbor[2]], cost), get(pq, neighbor, typemax(Int)))
            end
        end
    end
    ## return 0
end

## @lc code=end


text:
function graphloadertest()

	# create a graph
	# with no edges or vertices loaded
	g = Graph()
	# create a graph loader for the Graph
	gl = TestGraphLoader(g)

	# use another test method to obtain graph data for the test
	sourcegraph = buildsocialtestgraph()
	# copy this data to the test graph loader
	# this is to act in place of a persistent store from which to load
	gl.sourceedges = sourcegraph.edges
	gl.sourcevertices = sourcegraph.vertices

	# count the countries
	@test query(gl, (vertices, v->v.typelabel == "Country"), count) == 30

	@test query(gl, (vertices, v->v.typelabel == "Country" && get(v,"Name") == "Germany"), count) == 1

	# start with the vertices for Germany and France
	# follow the edges of type "LivesIn" to find all the people within the graph who live in Germany or France
	#     only edges that have a ForYears property of >= 5 (indicating the person has lived in the associated country for 5 or more years) will be followed
	# the tail of each edge is a person, filter the people selected to those with an age property >= 30
	part1 = query(gl,
				(vertices, v->v.typelabel == "Country" && (get(v,"Name") == "Germany" || get(v,"Name") == "France")),
			 	(incoming, e->e.typelabel == "LivesIn" && get(e,"ForYears",0) >= 5),
			  	(tail, v->get(v,"Age",0) >= 30))
	# get a count of these results
	@test query(part1, count) == 10

	# start with the results of the first query
	# follow edges of type "Follows" (in an incoming direction)
	# the tail of each edge is a person who follows a person selected in the first query
	# select a distinct set of vertices (as there may have been multiple edges to each person..as a person could follow more than one person selected by the first query)
	part2 = query(part1,
				(incoming, e->e.typelabel == "Follows"),
			  	tail,
				distinct
			 )
	# get a count of the results
	testcount = query(part2, count)
	@test testcount == 78
end

@time graphloadertest()


text:
using Test

include("playMadoka.jl")

println()
q1=[[:x は :y が 好き],[:y は 上条 が 好き],[:x は 魔法少女 です]]
println("Q:$(stringT(q1))")
a1=solveQ([:x,:y,:z],q1,魔法会議)
println("A: $a1")

println()
q2=[[:x は :y が 好き],[:y は 上条 が 好き],[:x は 魔法少女 です]]
println("Q:$(stringT(q2))")
a2=solveQ([:x,:y,:z],q2,魔法会議)
println("A: $a2")

println()
q3=[[:x は :y が 好き],[:y は :z が 好き],[:z は :x が 好き]]
println("Q:$(stringT(q3))")
a3=solveQ([:x,:y,:z],q3,魔法会議)
println("A: $a3")

println()
q4=[[:x は :y が 好き]]
println("Q:$(stringT(q4))")
a4=solveQ([:x,:y],q4,魔法会議)
println("A: $a4")

println()
q5=[[:x は :y が 好き],[:y は :x が 好き]]
println("Q:$(stringT(q5))")
a5=solveQ([:x,:y],q5,魔法会議)
println("A: $a5") #; printQM([:x,:y],q5,a5)

println()
q6=[[:y が 魔女 に なる],[:x は :y です],[:z は :y を 作る]]
#q6=[[:x は :y です],[:z は :y を 作る],[:y が 魔女 に なる]]
println("Q:$(stringT(q6))")
a6=solveQ([:x,:y,:z],q6,魔法会議)
println("A: $a6") #; printQM([:x,:y,:z],q6,a6)

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
module Types

using MySQL

const help_msg = """
+-------------------------------------------------------------+
| /register <name> <password> <password> - Create an account  |
| /login <name> <password> - Connect with an existing account |
| /who - Show your name, id                                   |
| /create <channel_name> [password]- Create a channel         |
| /join <channel_name> [password]- Join a channel             |
| /leave - Leave a channel                                    |
| /help - Display this message                                |
+-------------------------------------------------------------+

"""

const welcome_msg = """
Use /help if you need more informations

"""

mutable struct Client
    id::String
    name::String
    ip_addr::String
    current_channel_id::String
    conn::IO
    is_auth::Bool
end

function Client(ip_addr::String, conn::IO)
    Client("", "", ip_addr, "", conn, false)
end

mutable struct Channel
    id::String
    name::String
    description::String
    is_protected::Int
    password::String
    owner_id::String
end

mutable struct Storage
    listener::Any
    sql_conn::MySQL.Connection
    # ip => Client
    active_clients::Dict{String, Any}
    # id => Channel
    active_channels::Dict{String, Any}
end

function Storage(listener::Any, sql_conn::MySQL.Connection)
    Storage(listener, sql_conn, Dict(), Dict())
end

abstract type AbstractCommand{A <: Number, B <: Bool} end

function addCommandType(name::String, args_amount::String, auth::String)
    expr = """
    struct $name{A, B} <: AbstractCommand{A, B}
        args_amount::A
        auth::B
        function $name()
            if ($args_amount < 0)
                error(\"The arguments amount must be >= 0\")
            end
            new{Int, Bool}($args_amount, $auth)
        end
    end"""
    expr = Meta.parse(expr)
    eval(expr)
end

addCommandType("Register", "0", "false")
addCommandType("Login", "2", "false")
addCommandType("Help", "0", "false")
addCommandType("Who", "0", "true")
addCommandType("Create", "1", "true")
addCommandType("Join", "1", "true")
addCommandType("Leave", "0", "true")
addCommandType("Unknown", "0", "false")

const table = Dict(
    "register" => Register,
    "login" => Login,
    "help" => Help,
    "who" => Who,
    "create" => Create,
    "join" => Join,
    "leave" => Leave
)

function getCommandInfos(name::String)
    if (any([name == key for (key, _) in table]) == false)
        return (Unknown())
    end
    table[name]()
end

export Client, Channel, Storage, getCommandInfos

end # Data


text:
function givenscheck(T,n,tol)
    e = eps(real(T))
    A = rand(T,n,n)
    Ab = big.(A)
    G,_ = givens(A[4,5],A[5,5],4,5)
    Gb,_ = givens(Ab[4,5],Ab[5,5],4,5)
    B = copy(A)
    Bb = copy(Ab)
    lmul!(G,B)
    lmul!(Gb,Bb)
    @test norm(B-Bb) < tol * n * e
    rmul!(A,G)
    rmul!(Ab,Gb)
    @test norm(A-Ab) < tol * n * e
end

@testset "givens $T" for T in (Double64, Double32, Complex{Double64})
    n = 67
    givenscheck(T,n,4)
end

max_stars_repo_name: RalphAS/DoubleBLAS


text:
module Simulate_Pendulum

using  Modia3D
import Modia3D.ModiaMath

# Solids
groundMaterial = Modia3D.Material(color="DarkGreen", transparency=0.5)
# massProperties = Modia3D.MassProperties()
vmat1 = Modia3D.Material(color="LightBlue" , transparency=0.5)
vmat2 = Modia3D.Material(color="Red")

@assembly Pendulum(;Lx = 1.0, Ly=0.2*Lx, Lz=0.2*Lx, m=1.0) begin
   world  = Modia3D.Object3D(Modia3D.CoordinateSystem(0.5*Lx))

   # Pendulum
   body   = Modia3D.Object3D(Modia3D.Solid(Modia3D.SolidBeam(Lx,Ly,Lz), Modia3D.MassProperties(m=m), vmat1))
   frame1 = Modia3D.Object3D(body; r=[-Lx/2, 0.0, 0.0])
   cyl    = Modia3D.Object3D(frame1,Modia3D.Cylinder(Ly/2,1.2*Ly; material=vmat2))

   # Connect pendulum to world
   rev    = Modia3D.Revolute(world, frame1)
end

const Lx = 1.6
const m  = 0.5
# Modia3D.visualizeAssembly!( Pendulum(Lx=Lx, m=m) )

pendulum = Pendulum(Lx=Lx, m=m, sceneOptions=Modia3D.SceneOptions(visualizeFrames=true, defaultFrameLength=0.3))
model    = Modia3D.SimulationModel( pendulum )
result   = ModiaMath.simulate!(model, stopTime=4.5, tolerance=1e-6,log=true)

ModiaMath.plot(result, ["rev.phi", "rev.w", "rev.a"])

println("... success of Simulate_Pendulum.jl!")
end


text:
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

function compute(n::Int, digits::AbstractString)::Int
    digits = [parse(Int, digit) for digit ∈ replace(digits, "\n" => "")]
    max_product = 0
    for i ∈ 1:length(digits) - n + 1
        if 0 ∈ digits[i:i+n-1] continue end
        cur_product = prod(digits[i:i+n-1])
        if cur_product > max_product max_product = cur_product end
    end
    return max_product
end

number =
"""73167176531330624919225119674426574742355349194934
96983520312774506326239578318016984801869478851843
85861560789112949495459501737958331952853208805511
12540698747158523863050715693290963295227443043557
66896648950445244523161731856403098711121722383113
62229893423380308135336276614282806444486645238749
30358907296290491560440772390713810515859307960866
70172427121883998797908792274921901699720888093776
65727333001053367881220235421809751254540594752243
52584907711670556013604839586446706324415722155397
53697817977846174064955149290862569321978468622482
83972241375657056057490261407972968652414535100474
82166370484403199890008895243450658541227588666881
16427171479924442928230863465674813919123162824586
17866458359124566529476545682848912883142607690042
24219022671055626321111109370544217506941658960408
07198403850962455444362981230987879927244284909188
84580156166097919133875499200524063689912560717606
05886116467109405077541002256983155200055935729725
71636269561882670428252483600823257530420752963450"""

compute(4, number)

compute(13, number)

@benchmark compute(13, number)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
# cudnnSoftmaxAlgorithm_t
const CUDNN_SOFTMAX_FAST = Cint(0)
const CUDNN_SOFTMAX_ACCURATE = Cint(1)
const CUDNN_SOFTMAX_LOG = Cint(2)

# cudnnSoftmaxMode_t
const CUDNN_SOFTMAX_MODE_INSTANCE = Cint(0)
const CUDNN_SOFTMAX_MODE_CHANNEL = Cint(1)

function softmax(x::CuArray{T}, algo=CUDNN_SOFTMAX_ACCURATE) where T
    h = gethandle()
    mode = CUDNN_SOFTMAX_MODE_CHANNEL
    xdesc = TensorDesc(x, 4)
    y = similar(x)
    @cudnn(:cudnnSoftmaxForward,
        (Cptr,Cint,Cint,
        Cptr,Cptr,Cptr,
        Cptr,Cptr,Cptr),
        h, algo, mode, T[1], xdesc, x, T[0], xdesc, y)
    y
end

function ∇softmax!(y::CuArray{T}, dy, dx, algo=CUDNN_SOFTMAX_ACCURATE) where T
    h = gethandle()
    ydesc = TensorDesc(y, 4)
    mode = CUDNN_SOFTMAX_MODE_CHANNEL
    @cudnn(:cudnnSoftmaxBackward,
        (Cptr,Cint,Cint,
        Cptr,Cptr,Cptr,Cptr,Cptr,
        Cptr,Cptr,Cptr),
        h, algo, mode, T[1], ydesc, y, ydesc, dy, T[1], ydesc, dx)
end


text:
export pack_input

# one il + protein + water
function pack_input(data, pdb_dir::String, nil, nwater, sides)

  protein = data.protein
  cation  = data.cation
  anion   = data.anion

  lx = round(Int64,sides[1] / 2)
  ly = round(Int64,sides[2] / 2)
  lz = round(Int64,sides[3] / 2)

  io = open("box.inp","w")
  println(io,"tolerance 2.0")
  println(io,"output system.pdb")
  println(io,"add_box_sides 1.0")
  println(io,"filetype pdb")
  println(io,"seed -1")
  println(io,"                  ")
  println(io,"structure $(protein)")
  println(io," number 1")
  println(io," center")
  println(io," fixed 0. 0. 0. 0. 0. 0.")
  println(io,"end structure")
  println(io,"                  ")
  println(io,"structure $pdb_dir/WATER.pdb")
  println(io," number $nwater")
  println(io," inside box -$(lx). -$(ly). -$(lz). $(lx). $(ly). $(lz).")
  println(io,"end structure")
  println(io,"                  ")
  println(io,"structure $pdb_dir/$(cation)_VSIL.pdb")
  println(io," number $nil")
  println(io," inside box -$(lx). -$(ly). -$(lz). $(lx). $(ly). $(lz).")
  println(io,"end structure")
  println(io,"                  ")
  println(io,"structure $pdb_dir/$(anion)_VSIL.pdb")
  println(io," number $nil")
  println(io," inside box -$(lx). -$(ly). -$(lz). $(lx). $(ly). $(lz).")
  println(io,"end structure")

  close(io)

end

# one il + water
function pack_input_sol(data, pdb_dir::String, nil, nwater, sides)

  cation  = data.cation
  anion   = data.anion

  lx = round(Int64,sides[1] / 2)
  ly = round(Int64,sides[2] / 2)
  lz = round(Int64,sides[3] / 2)

  io = open("box.inp","w")
  println(io,"tolerance 2.0")
  println(io,"output system.pdb")
  println(io,"add_box_sides 1.0")
  println(io,"filetype pdb")
  println(io,"seed -1")
  println(io,"                  ")
  println(io,"structure $pdb_dir/WATER.pdb")
  println(io," number $nwater")
  println(io," inside box -$(lx). -$(ly). -$(lz). $(lx). $(ly). $(lz).")
  println(io,"end structure")
  println(io,"                  ")
  println(io,"structure $pdb_dir/$(cation)_VSIL.pdb")
  println(io," number $nil")
  println(io," inside box -$(lx). -$(ly). -$(lz). $(lx). $(ly). $(lz).")
  println(io,"end structure")
  println(io,"                  ")
  println(io,"structure $pdb_dir/$(anion)_VSIL.pdb")
  println(io," number $nil")
  println(io," inside box -$(lx). -$(ly). -$(lz). $(lx). $(ly). $(lz).")
  println(io,"end structure")

  close(io)

end

# two ils + protein + water
function pack_input(data, pdb_dir::String, ncat, nan, nwater, sides)

  protein = data.protein
  cation  = data.cation
  anion   = data.anion1
  anion2  = data.anion2

  lx = round(Int64,sides[1] / 2)
  ly = round(Int64,sides[2] / 2)
  lz = round(Int64,sides[3] / 2)

  io = open("box.inp","w")
  println(io,"tolerance 2.0")
  println(io,"output system.pdb")
  println(io,"add_box_sides 1.0")
  println(io,"filetype pdb")
  println(io,"seed -1")
  println(io,"                  ")
  println(io,"structure $(protein)")
  println(io," number 1")
  println(io," center")
  println(io," fixed 0. 0. 0. 0. 0. 0.")
  println(io,"end structure")
  println(io,"                  ")
  println(io,"structure $pdb_dir/WATER.pdb")
  println(io," number $nwater")
  println(io," inside box -$(lx). -$(ly). -$(lz). $(lx). $(ly). $(lz).")
  println(io,"end structure")
  println(io,"                  ")
  println(io,"structure $pdb_dir/$(cation)_VSIL.pdb")
  println(io," number $ncat")
  println(io," inside box -$(lx). -$(ly). -$(lz). $(lx). $(ly). $(lz).")
  println(io,"end structure")
  println(io,"                  ")
  println(io,"structure $pdb_dir/$(anion)_VSIL.pdb")
  println(io," number $nan")
  println(io," inside box -$(lx). -$(ly). -$(lz). $(lx). $(ly). $(lz).")
  println(io,"end structure")
  println(io,"                  ")
  println(io,"structure $pdb_dir/$(anion2)_VSIL.pdb")
  println(io," number $nan")
  println(io," inside box -$(lx). -$(ly). -$(lz). $(lx). $(ly). $(lz).")
  println(io,"end structure")

  close(io)

end

# two ils + water
function pack_input_elec(data, pdb_dir::String, ncat, nan, nwater, sides)

  cation  = data.cation
  anion   = data.anion1
  anion2  = data.anion2

  lx = round(Int64,sides[1] / 2)
  ly = round(Int64,sides[2] / 2)
  lz = round(Int64,sides[3] / 2)

  io = open("box.inp","w")
  println(io,"tolerance 2.0")
  println(io,"output system.pdb")
  println(io,"add_box_sides 1.0")
  println(io,"filetype pdb")
  println(io,"seed -1")
  println(io,"                  ")
  println(io,"structure $pdb_dir/WATER.pdb")
  println(io," number $nwater")
  println(io," inside box -$(lx). -$(ly). -$(lz). $(lx). $(ly). $(lz).")
  println(io,"end structure")
  println(io,"                  ")
  println(io,"structure $pdb_dir/$(cation)_VSIL.pdb")
  println(io," number $ncat")
  println(io," inside box -$(lx). -$(ly). -$(lz). $(lx). $(ly). $(lz).")
  println(io,"end structure")
  println(io,"                  ")
  println(io,"structure $pdb_dir/$(anion)_VSIL.pdb")
  println(io," number $nan")
  println(io," inside box -$(lx). -$(ly). -$(lz). $(lx). $(ly). $(lz).")
  println(io,"end structure")
  println(io,"                  ")
  println(io,"structure $pdb_dir/$(anion2)_VSIL.pdb")
  println(io," number $nan")
  println(io," inside box -$(lx). -$(ly). -$(lz). $(lx). $(ly). $(lz).")
  println(io,"end structure")

  close(io)

end


text:
module TestDay08

import AOC2021.Day08.part1, AOC2021.Day08.part2
using Test

test_input_raw = """
be cfbegad cbdgef fgaecd cgeb fdcge agebfd fecdb fabcd edb | fdgacbe cefdb cefbgd gcbe
edbfga begcd cbg gc gcadebf fbgde acbgfd abcde gfcbed gfec | fcgedb cgb dgebacf gc
fgaebd cg bdaec gdafb agbcfd gdcbef bgcad gfac gcb cdgabef | cg cg fdcagb cbg
fbegcd cbd adcefb dageb afcb bc aefdc ecdab fgdeca fcdbega | efabcd cedba gadfec cb
aecbfdg fbg gf bafeg dbefa fcge gcbea fcaegb dgceab fcbdga | gecf egdcabf bgf bfgea
fgeab ca afcebg bdacfeg cfaedg gcfdb baec bfadeg bafgc acf | gebdcfa ecba ca fadegcb
dbcfg fgd bdegcaf fgec aegbdf ecdfab fbedc dacgb gdcebf gf | cefg dcbef fcge gbcadfe
bdfegc cbegaf gecbf dfcage bdacg ed bedf ced adcbefg gebcd | ed bcgafe cdgba cbgef
egadfb cdbfeg cegd fecab cgb gbdefca cg fgcdab egfdb bfceg | gbdfcae bgc cg cgb
gcafb gcf dcaebfg ecagb gf abcdeg gaef cafbge fdbac fegbdc | fgae cfgab fg bagce"""

test_input = [string(x) for x in split(test_input_raw, "\n")]

@testset "Day 08" begin
    @testset "part 1" begin
        @test part1(test_input) == 26
    end
    @testset "part 2" begin
        @test part2(test_input) == 61229
    end
end

end # module


text:
import YAML

const WRONG_FORMAT_CONFIG_FILE =
    joinpath(dirname(pathof(PowerSystems)), "descriptors", "config.yml")
include(joinpath(BASE_DIR, "test", "data_5bus_pu.jl"))

@testset "Test reading in config data" begin
    data = IS.read_validation_descriptor(PSY.POWER_SYSTEM_STRUCT_DESCRIPTOR_FILE)
    @test data isa Vector
    @test !isempty(data)
    function find_struct()
        for item in data
            if item["struct_name"] == "ThermalStandard"
                return true
            end
        end
        return false
    end
    @test find_struct()
    @test_throws(ErrorException, IS.read_validation_descriptor("badfile.toml"))
end

@testset "Test adding custom validation YAML file to System" begin
    nodes = nodes5()
    sys_no_config = System(
        nodes,
        thermal_generators5(nodes),
        loads5(nodes),
        nothing,
        nothing,
        100.0,
        nothing,
        nothing;
        runchecks = true,
    )
    @test !isempty(sys_no_config.data.validation_descriptors)

    nodes = nodes5()
    sys_no_runchecks = System(
        nodes,
        thermal_generators5(nodes),
        loads5(nodes),
        nothing,
        nothing,
        100.0,
        nothing,
        nothing;
        runchecks = false,
    )
    @test isempty(sys_no_runchecks.data.validation_descriptors)
end

@testset "Test extracting struct info from validation_descriptor vector" begin
    data = [
        Dict(
            "fields" => Dict{Any, Any}[
                Dict(
                    "name" => "curtailpenalty",
                    "valid_range" => Dict{Any, Any}("max" => nothing, "min" => 0.0),
                ),
                Dict(
                    "name" => "variablecost",
                    "valid_range" => Dict{Any, Any}("max" => nothing, "min" => 0.0),
                ),
                Dict("name" => "internal"),
            ],
            "struct_name" => "EconHydro",
        ),
        Dict(
            "fields" => Dict{Any, Any}[
                Dict(
                    "name" => "curtailpenalty",
                    "valid_range" => Dict{Any, Any}("max" => nothing, "min" => 0.0),
                ),
                Dict(
                    "name" => "variablecost",
                    "valid_range" => Dict{Any, Any}("max" => nothing, "min" => 0.0),
                ),
                Dict("name" => "internal"),
            ],
            "struct_name" => "EconLoad",
        ),
    ]
    struct_name = "EconHydro"
    descriptor = IS.get_config_descriptor(data, struct_name)
    @test descriptor isa Dict{String, Any}
    @test haskey(descriptor, "struct_name")
    @test haskey(descriptor, "fields")
    @test descriptor["struct_name"] == struct_name
end

@testset "Test extracting field info from struct descriptor dictionary" begin
    config = Dict{Any, Any}(
        "fields" => Dict{Any, Any}[
            Dict("name" => "name", "data_type" => "String"),
            Dict("name" => "available", "data_type" => "Bool"),
            Dict("name" => "bus", "data_type" => "Bus"),
            Dict("name" => "tech", "data_type" => "Union{Nothing, TechThermal}"),
            Dict("name" => "econ", "data_type" => "Union{Nothing, EconThermal}"),
            Dict("name" => "internal", "data_type" => "IS.InfrastructureSystemsInternal"),
        ],
        "struct_name" => "ThermalStandard",
    )
    field_name = "econ"
    field_descriptor = IS.get_field_descriptor(config, field_name)
    @test field_descriptor isa Dict{Any, Any}
    @test haskey(field_descriptor, "name")
    @test field_descriptor["name"] == field_name
end

@testset "Test retrieving validation action" begin
    warn_descriptor = Dict{Any, Any}(
        "name" => "ramplimits",
        "valid_range" => Dict{Any, Any}("max" => 5, "min" => 0),
        "validation_action" => "warn",
    )
    error_descriptor = Dict{Any, Any}(
        "name" => "ramplimits",
        "valid_range" => Dict{Any, Any}("max" => 5, "min" => 0),
        "validation_action" => "error",
    )
    typo_descriptor = Dict{Any, Any}(
        "name" => "ramplimits",
        "valid_range" => Dict{Any, Any}("max" => 5, "min" => 0),
        "validation_action" => "asdfasdfsd",
    )
    @test IS.get_validation_action(warn_descriptor) == IS.validation_warning
    @test IS.get_validation_action(error_descriptor) == IS.validation_error
    @test_throws(ErrorException, IS.get_validation_action(typo_descriptor))
end

@testset "Test field validation" begin
    #test recursive call of validate_fields and a regular valid range
    nodes = nodes5()
    bad_therm_gen_rating = thermal_generators5(nodes)
    bad_therm_gen_rating[1].rating = -10
    @test_logs(
        (:error, r"Invalid range"),
        @test_throws(
            PSY.InvalidRange,
            System(
                nodes,
                bad_therm_gen_rating,
                loads5(nodes),
                nothing,
                nothing,
                100.0,
                nothing,
                nothing;
                runchecks = true,
            )
        )
    )

    #test custom range (activepower and activepowerlimits)
    bad_therm_gen_act_power = thermal_generators5(nodes)
    bad_therm_gen_act_power[1].activepower = 10
    nodes = nodes5()
    @test_logs (:warn, r"Invalid range") System(
        nodes,
        bad_therm_gen_act_power,
        loads5(nodes),
        nothing,
        nothing,
        100.0,
        nothing,
        nothing;
        runchecks = true,
    )

    #test validating named tuple
    nodes = nodes5()
    bad_therm_gen_ramp_lim = thermal_generators5(nodes)
    bad_therm_gen_ramp_lim[1].ramplimits = (up = -10, down = -3)
    @test_logs(
        (:error, r"Invalid range"),
        match_mode = :any,
        @test_throws(
            PSY.InvalidRange,
            System(
                nodes,
                bad_therm_gen_ramp_lim,
                loads5(nodes),
                nothing,
                nothing,
                100.0,
                nothing,
                nothing;
                runchecks = true,
            )
        )
    )
end

@testset "Test field validation" begin
    nodes = nodes5()
    sys = System(
        nodes,
        thermal_generators5(nodes),
        loads5(nodes),
        nothing,
        nothing,
        100.0,
        nothing,
        nothing;
        runchecks = true,
    )

    add_component!(
        sys,
        Bus(11, "11", BusTypes.PQ, 1, 1, (min = 0.9, max = 1.1), 123, nothing, nothing),
    )
    B = collect(get_components(Bus, sys))
    a = Arc(B[1], B[6])
    badline = Line(
        "badline",
        true,
        0.01,
        0.01,
        a,
        0.002,
        0.014,
        (from = 0.015, to = 0.015),
        5.0,
        (min = -1, max = 1),
    )
    @test_logs(
        (:error, r"cannot create Line"),
        match_mode = :any,
        @test_throws(PSY.InvalidValue, add_component!(sys, badline))
    )
end

@testset "Test field validation after deserialization" begin
    nodes = nodes5()
    sys = System(
        nodes,
        thermal_generators5(nodes),
        loads5(nodes),
        nothing,
        nothing,
        100.0,
        nothing,
        nothing,
    )

    add_component!(
        sys,
        Bus(11, "11", BusTypes.PQ, 1, 1, (min = 0.9, max = 1.1), 123, nothing, nothing),
    )
    path = joinpath(mktempdir(), "test_validation.json")
    try
        PSY.to_json(sys, path)
    catch
        rm(path)
        rethrow()
    end

    try
        sys2 = PSY.System(path)

        B = collect(get_components(Bus, sys2))
        a = Arc(B[1], B[6])
        badline = Line(
            "badline",
            true,
            0.01,
            0.01,
            a,
            0.002,
            0.014,
            (from = 0.015, to = 0.015),
            5.0,
            (min = -1, max = 1),
        )
        @test_logs(
            (:error, r"cannot create Line"),
            match_mode = :any,
            @test_throws(PSY.InvalidValue, add_component!(sys2, badline))
        )
    finally
        rm(path)
    end
end


text:
# Per Lyndon. Needs adaptation, and/or needs redefinition of map to use functions rather
# than AbstractOp.
#function rrule(map, f, xs)
#    # Rather than 3 maps really want 1 multimap
#    ys_and_pullbacks = map(x->rrule(f, x), xs) #Take this to ys = map(f, x)
#    ys = map(first, ys_and_pullbacks)
#    pullbacks = map(last, ys_and_pullbacks)
#    function map_pullback(dys)
#        _call(f, x) = f(x)
#        dfs_and_dxs = map(_call, pullbacks, dys)
#        # but in your case you know it will be NoTangent() so can  skip
#        df = sum(first, dfs_and_dxs)
#        dxs = map(last, dfs_and_dxs)
#        return NoTangent(), df, dxs
#    end
#    return ys, map_pullback
#end
macro scalarmaprule(func, derivative)
    return ChainRulesCore.@strip_linenos quote
        function ChainRulesCore.frule(
            (_, _, $(esc(:ΔA))),
            ::typeof(Base.map),
            ::typeof($(func)),
            $(esc(:A))::GBArray
        )
            $(esc(:Ω)) = map($(esc(func)), $(esc(:A)))
            return $(esc(:Ω)), $(esc(derivative)) .* unthunk($(esc(:ΔA)))
        end
        function ChainRulesCore.rrule(
            ::typeof(Base.map),
            ::typeof($(func)),
            $(esc(:A))::GBArray
        )
            $(esc(:Ω)) = map($(esc(func)), $(esc(:A)))
            function mapback($(esc(:ΔA)))
                NoTangent(), NoTangent(), $(esc(derivative)) .* $(esc(:ΔA))
            end
            return $(esc(:Ω)), mapback
        end
    end
end

function ChainRulesCore.frule(
    (_,_,ΔA),
    ::typeof(map),
    ::typeof(sqrt),
    A::Array
)
    Ω = map(sqrt, A)
    return Ω, inv.(2 .* Ω)
end

#Trig
@scalarmaprule UnaryOps.SIN cos.(A)
@scalarmaprule UnaryOps.COS -sin.(A)
@scalarmaprule UnaryOps.TAN @. 1 + (Ω ^ 2)

#Hyperbolic Trig
@scalarmaprule UnaryOps.SINH cosh.(A)
@scalarmaprule UnaryOps.COSH sinh.(A)
@scalarmaprule UnaryOps.TANH @. 1 - (Ω ^ 2)

@scalarmaprule UnaryOps.MINV -(Ω .^ 2)
@scalarmaprule UnaryOps.EXP Ω

@scalarmaprule UnaryOps.ABS sign.(A)
#Anything that uses MINV fails the isapprox tests :().
# Since in the immortal words of Miha - "FiniteDiff is smarter than you", these shouldn't be enabled.
#@scalarmaprule UnaryOps.ASIN @. inv(sqrt.(1 - A ^ 2))
#@scalarmaprule UnaryOps.ACOS @. inv(sqrt.(1 - A ^ 2))
#@scalarmaprule UnaryOps.ATAN @. inv(1 + A ^ 2)
#@scalarmaprule UnaryOps.SQRT inv.(2 .* Ω)

function frule(
    (_, _, ΔA),
    ::typeof(map),
    ::typeof(UnaryOps.IDENTITY),
    A::GBArray
)
    return (A, ΔA)
end
function rrule(::typeof(map), ::typeof(UnaryOps.IDENTITY), A::GBArray)
    return A, (ΔΩ) -> (NoTangent(), NoTangent(), ΔΩ)
end


text:
# TODO: the input + idx combination should probably be replaced by a stream

"""
    TokenBuffer("foo bar")

Turns a string into a readable stream, used for building tokenisers. Utility
parsers such as `spaces` and `number` read characters from the stream and
into an array of tokens.

Parsers return `true` or `false` to indicate whether they matched anything
in the input stream. They can therefore be combined easily, e.g.

    spacesornumber(ts) = spaces(ts) || number(ts)

either skips whitespace or parses a number token, if possible.

The simplest possible tokeniser accepts any `character` with no token breaks:

    function tokenise(input)
        ts = TokenBuffer(input)
        while !isdone(ts)
            character(ts)
        end
        return ts.tokens
    end

    tokenise("foo bar baz") # ["foo bar baz"]

The second simplest splits only on spaces:

    function tokenise(input)
        ts = TokenBuffer(input)
        while !isdone(ts)
            spaces(ts) || character(ts)
        end
        return ts.tokens
    end

    tokenise("foo bar baz") # ["foo", "bar", "baz"]

See `nltk_word_tokenize` for a more advanced example.
"""
mutable struct TokenBuffer
    input::Vector{Char}
    buffer::Vector{Char}
    tokens::Vector{String}
    idx::Int
end

TokenBuffer(input) = TokenBuffer(input, [], [], 1)

TokenBuffer(input::AbstractString) = TokenBuffer(collect(input))

Base.getindex(ts::TokenBuffer, i = ts.idx) = ts.input[i]

function isdone(ts::TokenBuffer)
    done = ts.idx > length(ts.input)
    done && flush!(ts)
    return done
end

"""
    flush!(::TokenBuffer, tokens...)

TokenBuffer builds the current token as characters are read from the input. When
the end of the current token is detected, call `flush!` to finish it and append
it to the token stream. Optionally, give additional tokens to be added to the
stream after the current one.
"""
function flush!(ts::TokenBuffer, s...)
    if !isempty(ts.buffer)
        push!(ts.tokens, String(ts.buffer))
        empty!(ts.buffer)
    end
    push!(ts.tokens, s...)
    return
end

"""
    lookahead(::TokenBuffer, s; boundary = false)

Peek at the input to see if `s` is coming up next. `boundary` specifies whether
a word boundary should follow `s`.

    julia> lookahead(TokenBuffer("foo bar"), "foo")
    true

    julia> lookahead(TokenBuffer("foo bar"), "bar")
    false

    julia> lookahead(TokenBuffer("foo bar"), "foo", boundary = true)
    true

    julia> lookahead(TokenBuffer("foobar"), "foo", boundary = true)
    false
"""
function lookahead(ts::TokenBuffer, s; boundary = false)
    ts.idx + length(s) - 1 > length(ts.input) && return false
    for j = 1:length(s)
        ts.input[ts.idx-1+j] == s[j] || return false
    end
    if boundary
        next = ts.idx + length(s)
        next > length(ts.input) && return true
        (isletter(ts[next]) || ts[next] == '-') && return false
    end
    return true
end

"""
    character(::TokenBuffer)

Push the next character in the input into the buffer's current token.
"""
function character(ts)
    push!(ts.buffer, ts[])
    ts.idx += 1
    return true
end

"""
    spaces(::TokenBuffer)

If there is whitespace in the input, skip it, and flush the current token.
"""
function spaces(ts)
    isspace(ts[]) || return false
    flush!(ts)
    ts.idx += 1
    return true
end

"""
    atoms(::TokenBuffer, ["--", "...", ...])

Matches a set of atomic tokens, such as `...`, which should always be treated
as a single token, regardless of word boundaries.
"""
function atoms(ts, as)
    for a in as
        lookahead(ts, a) || continue
        flush!(ts, String(a))
        ts.idx += length(a)
        return true
    end
    (ispunct(ts[]) && ts[] != '.' && ts[] != '-') || return false
    flush!(ts, string(ts[]))
    ts.idx += 1
    return true
end

"""
    suffixes(::TokenBuffer, ["'ll", "'re", ...])

Matches tokens with suffixes, such as `you're`, that should be treated as
separate tokens.
"""
function suffixes(ts, ss)
    isempty(ts.buffer) && return false
    for s in ss
        lookahead(ts, s, boundary=true) || continue
        flush!(ts, String(s))
        ts.idx += length(s)
        return true
    end
    return false
end

"""
    replaces(::TokenBuffer, ["cannot"=>("can", "not"), ("freeeee"=>("free",), ...])

Matches tokens, and flushs their replacement to the stream.
The replacements should be a tuple of strings (potentially a 1-tuple), as this
can be used for splitting tokens.
For example, `cannot` would be split into `can` and `not`.
`freeeee` would be replaced with `free`
"""
function replaces(ts, ss)
    for (pat, subs) in ss
        lookahead(ts, pat, boundary=true) || continue
        flush!(ts, subs...)
        ts.idx += length(pat)
        return true
    end
    return false
end

"""
    openquote(::TokenBuffer)

Matches " used as an opening quote, and tokenises it as ``.
"""
function openquote(ts)
    ts[] == '"' || return false
    flush!(ts, "``")
    ts.idx += 1
    return true
end

"""
    closingquote(::TokenBuffer)

Matches " used as a closing quote, and tokenises it as ''.
"""
function closingquote(ts)
    ts[] == '"' || return false
    flush!(ts, "''")
    ts.idx += 1
    return true
end

"""
    number(::TokenBuffer)

Matches numbers such as `10,000.5`, preserving formatting.
"""
function number(ts, sep = (':', ',', '\'', '.'); check_sign = false)
    i = ts.idx
    if check_sign && ts[] ∈ ['+', '-'] && ( i == 1 || isspace(ts[i-1]))
        i += 1
    end

    i <= length(ts.input) && isdigit(ts[i]) || return false
    while i <= length(ts.input) && (isdigit(ts[i]) ||
                (ts[i] in sep && i < length(ts.input) && isdigit(ts[i+1])))
        i += 1
    end
    flush!(ts, String(ts[ts.idx:i-1]))
    ts.idx = i
    return true
end


text:
module debug

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
push!(LOAD_PATH, "./src/")   # TODO: Refactor module into folder? Also path is
push!(LOAD_PATH, "./src/CommonUtils/")
using Documenter
using CommonUtils
using Basis1D,Basis2DTri,Basis2DQuad,Basis3DHex
using UniformQuadMesh,UniformTriMesh,UniformHexMesh

makedocs(modules=[CommonUtils,Basis1D,Basis2DTri,UniformTriMesh], sitename="Documentation for JuliaDG")

max_stars_repo_name: jlchan/ESDG-CNS

max_issues_repo_name: jlchan/ESDG-CNS

max_forks_repo_name: jlchan/ESDG-CNS


text:
using GitHub, Base64

GitHub.Content(::Nothing) = nothing

(@isdefined AUTH) || (const AUTH = GitHub.authenticate(ENV["GITHUB_AUTH"]))
(@isdefined REPO) || (const REPO = Repo("crstnbr/dqmc-benchmarks"))

function push(file::AbstractString;
            filename=file,
            message="push $filename"
            )
    content = read(file, String)
    d = Dict(
        "content" => base64encode(content), # maybe escape_string
        "message" => string(message, " (api)")
    )
    create_file(REPO, filename; auth=AUTH, params=d)
    nothing
end

function rm(filename::AbstractString;
            message="rm $filename"
            )
    d = Dict(
        "message" => string(message, " (api)"),
        "sha" => GitHub.file(REPO, filename; auth=AUTH).sha
    )
    delete_file(REPO, filename; auth=AUTH, params=d)
    nothing
end

function rm(filenames::AbstractVector{T}) where T<:AbstractString
    for f in filenames
        @info "Removing $f"
        rm(f)
    end
    nothing
end

function download(filename)
    f = GitHub.file(REPO, filename; auth=AUTH)
    c = replace(f.content, "\n"=>"")
    s = String(base64decode(c))
    write(filename, s)
    nothing
end

getfilelist() = GitHub.name.(directory(REPO, "."; auth=AUTH)[1])
function getresultfiles(; branch="")
    filter(x->occursin("result", x)&&occursin(branch, x)&&endswith(x, "json"), getfilelist())
end

clear_repo() = rm(filter(f->occursin("result", f), getfilelist()))

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
# ---
# title: noiseonf
# id: demo_noiseonf
# cover: assets/noiseonf.png
# author: Peter Kovesi
# date: 2018-10-26
# ---

using Images
using ImageContrastAdjustment
using ImagePhaseCongruency

## Noise images with amplitude decay exponents of 1.5 and 2.5
img1 = noiseonf(512, 1.5)
img2 = noiseonf(512, 2.5)

mosaic(
    adjust_histogram(Gray.(img1), LinearStretching()),
    adjust_histogram(img2, LinearStretching());
    nrow=1
)

# save cover image #src
isdir("assets") || mkdir("assets") #src
save(joinpath("assets", "noiseonf.png"), adjust_histogram(Gray.(img1), LinearStretching())) #src

max_stars_repo_name: peterkovesi/ImagePhaseCongruency

max_issues_repo_name: peterkovesi/ImagePhaseCongruency

max_forks_repo_name: peterkovesi/ImagePhaseCongruency


text:
using Revise
using IsopycnalSurfaces, Test

@testset "IsopycnalSurfaces.jl" begin
    # test input from a 3D array (failed Dec 09 2021)
    @testset "3d_array" begin

        # currently empty
        pz = collect(0.:500.:4000.) # pressure levels
        nz = length(pz)

        θz = collect(range(20,stop=10,length=nz))
        Sz = collect(range(36,stop=35,length=nz))

        ztest = sort(rand(2:8,2))
        σ₁true = sigma1column(θz[ztest],Sz[ztest],pz[ztest])

        # put θ,S,p into 3D array
        nx = 10; ny = 10;
        θ = Array{Float64,3}(undef,nx,ny,nz)
        S = Array{Float64,3}(undef,nx,ny,nz)
        p = Array{Float64,3}(undef,nx,ny,nz)

        [θ[i,j,k] = θz[k] for i = 1:nx for j = 1:ny for k = 1:nz ]
        [S[i,j,k] = Sz[k] for i = 1:nx for j = 1:ny for k = 1:nz ]
        [p[i,j,k] = pz[k] for i = 1:nx for j = 1:ny for k = 1:nz ]

        σ₁grid = collect(range(minimum(σ₁true),stop=maximum(σ₁true),length=20))

        vars = Dict("θ" => θ, "Sp" => S)

        @testset "3d_array_spline" begin
            splorder = 3
            varsσ = vars2sigma1(vars,pz,σ₁grid,splorder)
            xx = rand(1:nx); yy = rand(1:ny)
            @test isapprox(varsσ["p"][xx,yy,begin],pz[ztest[begin]])
            @test isapprox(varsσ["p"][xx,yy,end],pz[ztest[end]])
        end

        @testset "3d_array_linear" begin
            splorder = 3
            linearinterp = true
            varsσ = vars2sigma1(vars,pz,σ₁grid,splorder,linearinterp)
            xx = rand(1:nx); yy = rand(1:ny)
            @test isapprox(varsσ["p"][xx,yy,begin],pz[ztest[begin]])
            @test isapprox(varsσ["p"][xx,yy,end],pz[ztest[end]])
        end
    end
end


text:
function HMMgamma(alpha,phghm)
    #HMMGAMMA HMM Posterior smoothing using the Rauch-Tung-Striebel correction method
    # gamma=HMMbackward(alpha,phghm)
    #
    # Inputs:
    # alpha : alpha forward messages (see HMMforward.m)
    # phghm : transition distribution in a matrix
    #
    # Outputs: gamma(i,t) is p(h(t)=i|v(1:T))
    # See also HMMbackward.m, HMMviterbi.m, demoHMMinference.m

    T=size(alpha,2); H=size(phghm, 1);
    # gamma recursion
    gamma=zeros(size(alpha))
    gamma[:,T]=alpha[:,T];
    for t=T-1:-1:1
        phghp=condp(phghm'.*repmat(alpha[:,t],1,H));
        gamma[:,t]=condp(phghp*gamma[:,t+1]);
    end

    if 1==0 # gamma recursion: More human readable
        gamma[:, T]=alpha[:, T]./sum(alpha[:, T])
        for t = T-1:-1:1
            phghp=phghm'.*repmat(alpha[:,t],1,H)
            phghp=phghp./repmat(sum(phghp, 1),H,1)
            gamma[:,t]=phghp*gamma[:,t+1]
        end
    end
    return gamma
end


text:
open("input_file","r") do f
    for line in eachline(f)
      println("read line: ", line)
    end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# # Support Vector Machine
#
# In this notebook we show how you can use KernelFunctions.jl to generate
# kernel matrices for classification with a support vector machine, as
# implemented by [LIBSVM](https://github.com/JuliaML/LIBSVM.jl).

using Distributions
using KernelFunctions
using LIBSVM
using LinearAlgebra
using Plots
using Random

## Set seed
Random.seed!(1234);

# ## Generate half-moon dataset

# Number of samples per class:
n1 = n2 = 50;

# We generate data based on SciKit-Learn's sklearn.datasets.make_moons function:

angle1 = range(0, π; length=n1)
angle2 = range(0, π; length=n2)
X1 = [cos.(angle1) sin.(angle1)] .+ 0.1 .* randn.()
X2 = [1 .- cos.(angle2) 1 .- sin.(angle2) .- 0.5] .+ 0.1 .* randn.()
X = [X1; X2]
x_train = RowVecs(X)
y_train = vcat(fill(-1, n1), fill(1, n2));

# ## Training
#
# We create a kernel function:
k = SqExponentialKernel() ∘ ScaleTransform(1.5)

# LIBSVM can make use of a pre-computed kernel matrix.
# KernelFunctions.jl can be used to produce that using `kernelmatrix`:
model = svmtrain(kernelmatrix(k, x_train), y_train; kernel=LIBSVM.Kernel.Precomputed)

# ## Prediction
#
# For evaluation, we create a 100×100 2D grid based on the extent of the training data:
test_range = range(floor(Int, minimum(X)), ceil(Int, maximum(X)); length=100)
x_test = ColVecs(mapreduce(collect, hcat, Iterators.product(test_range, test_range)));

# Again, we pass the result of KernelFunctions.jl's `kernelmatrix` to LIBSVM:
y_pred, _ = svmpredict(model, kernelmatrix(k, x_train, x_test));

# We can see that the kernelized, non-linear classification successfully separates the two classes in the training data:
plot(; lim=extrema(test_range), aspect_ratio=1)
contourf!(
    test_range,
    test_range,
    y_pred;
    levels=1,
    color=cgrad(:redsblues),
    alpha=0.7,
    colorbar_title="prediction",
)
scatter!(X1[:, 1], X1[:, 2]; color=:red, label="training data: class –1")
scatter!(X2[:, 1], X2[:, 2]; color=:blue, label="training data: class 1")


text:
__precompile__(false)

module PixhawkData

using Cxx
using Libdl

const depfile = joinpath(@__DIR__, "..", "deps", "deps.jl")
if isfile(depfile)
    include(depfile)
else
    error("libPixhawkData is not properly installed. Please run Pkg.build(\"PixhawkData\") to rectify this error.")
end

const path_to_header = joinpath(@__DIR__, "..", "deps", "usr", "include")

function __init__()
    addHeaderDir(path_to_header, kind=C_System)
    Libdl.dlopen(libpixhawkdata, Libdl.RTLD_GLOBAL)

    cxx"#include <iostream>"
    cxx"#include <vector>"
    cxx"#include <tuple>"
end

export
    top,
    commands,
    parse_commandline,
    quit_handler,
    main

include("mavlink_control.jl")

export
    Serial_Port,
    Serial_Port1,
    initialize_defaults,
    read_message,
    write_message,
    open_serial,
    close_serial,
    start,
    stop,
    handle_quit,
    _open_port,
    _setup_port,
    _read_port,
    _write_port

include("serial_port.jl")

export
    get_time_usec,
    set_position,
    set_velocity,
    set_acceleration,
    set_yaw,
    set_yaw_rate
    start_autopilot_interface_read_thread,
    start_autopilot_interface_write_thread,
    Time_Stamps,
    Mavlink_Messages,
    get_time_usec,
    set_position,
    set_velocity,
    set_acceleration,
    set_yaw,
    set_yaw_rate,
    Autopilot_Interface,
    Autopilot_Interface1,
    update_setpoint,
    read_messages,
    write_messages,
    write_setpoint,
    enable_offboard_control,
    disable_offboard_control,
    toggle_offboard_control,
    start,
    stop,
    start_read_thread,
    start_write_thread,
    handle_quit,
    read_thread,
    write_thread,
    start_autopilot_interface_read_thread,
    start_autopilot_interface_write_thread

include("autopilot_interface.jl")

end


text:
push!(LOAD_PATH, joinpath(@__DIR__, ".."))

using PkgBenchmark
using Oceananigans
using Benchmarks

baseline = BenchmarkConfig(id="master")
script = joinpath(@__DIR__, "benchmarkable_incompressible_model.jl")
resultfile = joinpath(@__DIR__, "regression_benchmarks.json")

judge(Oceananigans, baseline, script=script, resultfile=resultfile, verbose=true)


text:
# Univariate

const VectorOfUnivariate = Distributions.Product

function arraydist(dists::AbstractVector{<:UnivariateDistribution})
    return product_distribution(dists)
end
function arraydist(dists::AbstractVector{<:Normal})
    m = mapvcat(mean, dists)
    s = mapvcat(std, dists)
    return TuringMvNormal(m, s)
end

function Distributions.logpdf(dist::VectorOfUnivariate, x::AbstractVector{<:Real})
    return sum(map((d, x) -> logpdf(d, x), dist.v, x))
end
function Distributions.logpdf(dist::VectorOfUnivariate, x::AbstractMatrix{<:Real})
    # eachcol breaks Zygote, so we need an adjoint
    return mapvcat(dist.v, eachcol(x)) do dist, c
        sum(map(c) do x
            logpdf(dist, x)
        end)
    end
end
@adjoint function Distributions.logpdf(dist::VectorOfUnivariate, x::AbstractMatrix{<:Real})
    # Any other more efficient implementation breaks Zygote
    f(dist, x) = [sum(logpdf.(dist.v, view(x, :, i))) for i in 1:size(x, 2)]
    return pullback(f, dist, x)
end

struct MatrixOfUnivariate{
    S <: ValueSupport,
    Tdist <: UnivariateDistribution{S},
    Tdists <: AbstractMatrix{Tdist},
} <: MatrixDistribution{S}
    dists::Tdists
end
Base.size(dist::MatrixOfUnivariate) = size(dist.dists)
function arraydist(dists::AbstractMatrix{<:UnivariateDistribution})
    return MatrixOfUnivariate(dists)
end
function Distributions.logpdf(dist::MatrixOfUnivariate, x::AbstractMatrix{<:Real})
    # Broadcasting here breaks Tracker for some reason
    # A Zygote adjoint is defined for mapvcat to use broadcasting
    return sum(map(dist.dists, x) do dist, x
        logpdf(dist, x)
    end)
end
function Distributions.logpdf(dist::MatrixOfUnivariate, x::AbstractArray{<:AbstractMatrix{<:Real}})
    return mapvcat(x -> logpdf(dist, x), x)
end
function Distributions.logpdf(dist::MatrixOfUnivariate, x::AbstractArray{<:Matrix{<:Real}})
    return mapvcat(x -> logpdf(dist, x), x)
end
function Distributions.rand(rng::Random.AbstractRNG, dist::MatrixOfUnivariate)
    return rand.(Ref(rng), dist.dists)
end

# Multivariate

struct VectorOfMultivariate{
    S <: ValueSupport,
    Tdist <: MultivariateDistribution{S},
    Tdists <: AbstractVector{Tdist},
} <: MatrixDistribution{S}
    dists::Tdists
end
Base.size(dist::VectorOfMultivariate) = (length(dist.dists[1]), length(dist))
Base.length(dist::VectorOfMultivariate) = length(dist.dists)
function arraydist(dists::AbstractVector{<:MultivariateDistribution})
    return VectorOfMultivariate(dists)
end
function Distributions.logpdf(dist::VectorOfMultivariate, x::AbstractMatrix{<:Real})
    # eachcol breaks Zygote, so we define an adjoint
    return sum(logpdf.(dist.dists, eachcol(x)))
end
function Distributions.logpdf(dist::VectorOfMultivariate, x::AbstractArray{<:AbstractMatrix{<:Real}})
    return mapvcat(x -> logpdf(dist, x), x)
end
function Distributions.logpdf(dist::VectorOfMultivariate, x::AbstractArray{<:Matrix{<:Real}})
    return mapvcat(x -> logpdf(dist, x), x)
end
@adjoint function Distributions.logpdf(dist::VectorOfMultivariate, x::AbstractMatrix{<:Real})
    f(dist, x) = sum(mapvcat(i -> logpdf(dist.dists[i], view(x, :, i)), 1:size(x, 2)))
    return pullback(f, dist, x)
end
function Distributions.rand(rng::Random.AbstractRNG, dist::VectorOfMultivariate)
    init = reshape(rand(rng, dist.dists[1]), :, 1)
    return mapreduce(i -> rand(rng, dist.dists[i]), hcat, 2:length(dist); init = init)
end


text:
function aggregate(obj::ThirdMoment, arg0::SecondMoment)
    return jcall(obj, "aggregate", void, (SecondMoment,), arg0)
end

function clear(obj::ThirdMoment)
    return jcall(obj, "clear", void, ())
end

function copy(obj::ThirdMoment)
    return jcall(obj, "copy", ThirdMoment, ())
end

function get_result(obj::ThirdMoment)
    return jcall(obj, "getResult", jdouble, ())
end

function increment(obj::ThirdMoment, arg0::jdouble)
    return jcall(obj, "increment", void, (jdouble,), arg0)
end


text:
using TraceCalls
using TraceCalls: tree_size
using NBInclude
using Base.Test

@testset "TraceCalls" begin

@testset "Basics" begin
    include("basics.jl")
end

################################################################################
# Testing with popular packages
root_of(mod::Module) =
    joinpath(dirname(TraceCalls.module_definition_file(string(mod))), "..")

@testset "JuMP" begin
    using JuMP
    example = joinpath(root_of(JuMP), "examples", "basic.jl")
    @test ctree_size(@trace JuMP include(example)) > 10
end

# Disabled because I think it brought in Cairo, which caused homebrew issues compiling
# GCC...? Not sure - July'18
# @testset "Gadfly" begin
#     using Gadfly
#     # Most of the complexity in Gadfly comes from Base.show(::Gadfly.PlotObject, ...), and
#     # that's not tested here.
#     @test ctree_size(@trace Gadfly plot(x=rand(10), y=rand(10))) > 5
# end

# Knet should be tested again either once they've updated to 0.6, or once
# https://github.com/JuliaLang/julia/issues/22729 is fixed.
# using KNet
# @trace Knet include(joinpath(root_of(Knet), "examples", "linreg.jl"))

@testset "DataStructures" begin
    using DataStructures: binary_minheap, OrderedDict
    @test ctree_size(@trace DataStructures binary_minheap([2,11, 14])) > 3
    # This looks innocent, but it helped us find a problem with typed slurped args
    @test ctree_size(@trace DataStructures.OrderedDict OrderedDict(1=>2, 3=>4)) >= 1
end

@testset "PyCall" begin
    using PyCall
    @test ctree_size(@trace PyCall pyimport(:math)[:pi]) > 5
end

@testset "Optim" begin
    using Optim   # From https://github.com/JuliaNLSolvers/Optim.jl
    rosenbrock(x) =  (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
    @test ctree_size(@trace Optim optimize(rosenbrock, zeros(2), BFGS())) > 50
end

# August'18: this example no longer works (even without TraceCalls usage), and
# I don't have time to fix it ATM.
# @testset "DiffEq" begin
#     # We can use DifferentialEquations.jl, but that's a huge download and install
#     # Example from http://docs.juliadiffeq.org/stable/tutorials/ode_example.html
#     using OrdinaryDiffEq, DiffEqBase
#     # interp_summary is to fix some issue - lost the link
#     DiffEqBase.interp_summary(::OrdinaryDiffEq.InterpolationData) = "42"
#     f(t,u) = 1.01*u
#     u0=1/2
#     tspan = (0.0,1.0)
#     prob = ODEProblem(f,u0,tspan)
#     @test ctree_size(@trace OrdinaryDiffEq OrdinaryDiffEq.solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)) > 20
# end

@testset "LightGraphs" begin
    using LightGraphs
    A = [
        0 1 1
        1 0 1
        1 1 0
    ]
    graph = Graph(A)
    @test ctree_size(@trace LightGraphs adjacency_matrix(graph)) > 5
end

################################################################################

@testset "README" begin
    @nbinclude("../README.ipynb")
end

end


text:
##
using Pkg
Pkg.activate(".")
Pkg.instantiate()
include("src/EasyStocks.jl")

## Setup
using CairoMakie # For plotting

## Solve
h1 = 1.0
h2 = 2.0
y2 = 3.0
mp1 = ModPar(q=0.003,ψ=0.0,xstar=.0,η=0.3)
mp1.γ = 2 + 1/(1-mp1.η)
np1 = NumPar(mp1,nx=101,nsav=101,nq=31,nα=101,nh=1,xmax=10,ygrd=[0.,y2])
np1.hgrd[:] .= h1
MS1 = ModelSolution(mp1,np1)
@time SolveModel!(MS1)

# With housing or jump?
# Jump:
mp2 = ModPar(q=mp1.q,ψ=0.015,xstar=mp1.xstar,η =mp1.η,γ = mp1.γ)
np2 = NumPar(mp1,nx=np1.nx,nsav=np1.nsav,nq=np1.nq,nα=np1.nα,nh=1,xmax=np1.xmax,ygrd=np1.ygrd)
# Housing:
mp2 = ModPar(q=mp1.q,ψ=mp1.ψ,xstar=mp1.xstar,η =mp1.η,γ = mp1.γ)
np2 = NumPar(mp2,nx=np1.nx,nsav=np1.nsav,nq=np1.nq,nα=np1.nα,nh=2,xmax=np1.xmax,ygrd=np1.ygrd)

np2.hgrd[:] .= h2
MS2 = ModelSolution(mp2,np2)
@time SolveModel!(MS2)

# mp2.ms=0.3
np2.hgrd[:] .= [h1, h2]
MS3 = ModelSolution(mp2,np2)
@time SolveModel!(MS3)
#
fig = Figure()
g2 = fig[1,1] = GridLayout()
g1 = fig[2,1] = GridLayout()
ax1  = Axis(g1[1,1],xlabel = "", ylabel = "")
ax2  = Axis(g1[1,2],xlabel = "", ylabel = "")
ax3  = Axis(g1[1,3],xlabel = "", ylabel = "")

v2 = Axis(g2[1,1], xlabel ="", ylabel = "Utils")
h2 = Axis(g2[1,2], xlabel ="", ylabel = "Housing")
s3 = Axis(g2[1,3])

ih = 1

lines!(s3,MS1.s[:,ih] + MS1.b[:,ih],linestyle=:dash,color=:gray,)
lines!(s3,MS2.s[:,ih] + MS2.b[:,ih],color=:orange,)

ia = 2
lines!(h2,MS1.np.xgrd,MS1.h[:,ih,ia],linestyle=:dash,color=:gray,)
lines!(h2,MS2.np.xgrd,MS2.h[:,ih,ia],color=:orange,)
lines!(h2,MS2.np.xgrd,MS3.h[:,ih,ia],color=:black,)
ylims!(h2,np2.hgrd[1]-0.5,np2.hgrd[end]+0.5)

ia = 2
ih = 1
lines!(v2,MS2.np.xgrd,MS2.V[:,1,ia],color=:orange,)
lines!(v2,MS3.np.xgrd,MS3.V[:,1,ia],color=:black,)
lines!(v2,MS1.np.xgrd,MS1.V[:,1,ia],linestyle=:dash,color=:gray,)
ylims!(v2,max(minimum(MS3.V[:,1,ia]),-5),maximum(MS3.V[:,1,ia]))

# rowsize!(g2,2,Relative(1.33))

ih = 1
lines!(ax1,MS1.np.xgrd, MS1.α[:,ih])
ylims!(ax1,-0.05,1.05)
# text!(ax1, "Not willing\nto pay cost \n ↓", position=(0.0,0.5),align=(:left,:center))

ih = 1
lines!(ax2,MS2.np.xgrd, MS1.α[:,ih], linestyle=:dash,color=:gray,)
lines!(ax2,MS2.np.xgrd, MS2.α[:,ih], color=:orange)
ylims!(ax2,-0.05,1.05)
# text!(ax2, "↑\n Gambling\n for housing", position=(3.2,0.2),align=(:right,:center))
# text!(ax2, "Stay safe \nfor housing\n↓", position=(3.8,0.2),align=(:left,:center))
# text!(ax2, "Gradually re-enter\n↓", position=(4.7,0.8),align=(:left,:center))

ih = 2
lines!(ax3,MS2.np.xgrd, MS1.α[:,1],color=:gray)
lines!(ax3,MS2.np.xgrd, MS2.α[:,1],color=:orange)
lines!(ax3,MS2.np.xgrd, MS3.α[:,1],color=:black)
ylims!(ax3,-0.05,1.05)

fig
# ## Run-code

# fig = Figure()
# ax1 = Axis(fig[1,1])
# MS = MS1
# hitp = LinearInterpolation(MS.np.xgrd,MS2.h[:,2],extrapolation_bc=Flat())

# probh2 = fill(-1.0,MS.np.nx)
# for ix in eachindex(MS.np.xgrd)
#   xns = MS.s[ix,1]*(1.0.+MS.np.rsgrd) .+ MS.b[ix,1] *(1.0 + MS.mp.r)
#   probh2[ix] = sum(MS.np.πrs.*(hitp.(xns).==MS2.np.hgrd[2]))
# end

# scatter!(ax1,probh2)

# MS = MS2
# hitp = LinearInterpolation(MS2.np.xgrd,MS.h[:,2],extrapolation_bc=Flat())

# probh2 = fill(-1.0,MS.np.nx)
# for ix in eachindex(MS.np.xgrd)
#   xns = MS.s[ix,1]*(1.0.+MS.np.rsgrd) .+ MS.b[ix,1] *(1.0 + MS.mp.r)
#   probh2[ix] = sum(MS.np.πrs.*(hitp.(xns).==MS2.np.hgrd[2]))
# end

# scatter!(ax1,probh2.+1.0)
# fig


text:
#=
Circuit builders for MPSCircuit.jl.
=#
export DCbuilder, MPSbuilder

"""
    DCbuilder(nBit::Int64, depth::Int64)
Structure of elements may needed for a Quantum differentiable circuit.
\n
\nFields:
\n`block::ChainBlock`:    The block for 1 depth.
\n`circuit::ChainBlock`:  The standard structure of a differentiable circuit without the modifications on the head and the tail.
\n`fullbody::ChainBlock`: The modified version of a differentiable circuit with structures of a head and a tail (same number of blocks).
\n`head::ChainBlock`:     The "head" of the differentiable circuit(if directly input bits).
\n`tail::ChainBlock`:     The "tail" of the differentiable circuit(2 layers of rotaional gates).
"""
struct DCbuilder
    block::ChainBlock    # The block for 1 depth.
    circuit::ChainBlock  # The standard structure of a differentiable circuit without the modifications on the head and the tail..
    fullbody::ChainBlock # The modified version of a differentiable circuit with structures of a head and a tail (same number of blocks).
    head::ChainBlock     # The "head" of the differentiable circuit(if directly input bits).
    tail::ChainBlock     # The "tail" of the differentiable circuit(2 layers of rotaional gates).

    function DCbuilder(nBit::Int64, depth::Int64) #DON'T USE repeat to avoid POINTER side-effect!
        block = chain(nBit, vcat([chain(nBit,[put(nBit, i=>Rz(0)) for i=nBit:-1:1])],
                                 [chain(nBit,[put(nBit, i=>Rx(0)) for i=nBit:-1:1])],
                                 [chain(nBit,[put(nBit, i=>Rz(0)) for i=nBit:-1:1])],
                                 [chain(nBit, vcat([control(nBit, i, (i-1)=>X) for i=nBit:-1:2],
                                                   [control(nBit, 1, nBit=>X)]))] ))
        tail = chain(nBit, vcat([chain(nBit,[put(nBit, i=>Rz(0)) for i=nBit:-1:1])],
                                [chain(nBit,[put(nBit, i=>Rx(0)) for i=nBit:-1:1])]))
        head = deepcopy(block[2:end])
        blocks = ChainBlock[]
        for i=1:depth push!(blocks, deepcopy(block)) end
        Cblock = chain(nBit, blocks)
        body = chain(nBit, vcat([head], blocks[2:end], [tail]))
        new(block, Cblock, body, head, tail)
    end
end

"""
    MPSc{circuit::ChainBlock, cExtend::ChainBlock}
\n
Fields:
\n`circuit::ChainBlock`: QMPS circuit.
\n`cExtend::ChainBlock`: QMPS-extended circuit.
"""
struct MPSc
    circuit::ChainBlock # MPS-qubit-reusable circuit.
    cExtend::ChainBlock # MPS-extended circuit.
end

"""
    MPSbuilder(nBitA::Int64, vBit::Int64, rBit::Int64, blockT::Union{String, Tuple{String, Int64}})
    ->
    MPSc{circuit::ChainBlock, cExtend::ChainBlock}
Function for creating different types of MPS circuits.
\n`blockT::Union{String, Tuple{String, Int64}}`
\n1) `blockT = ("DC", depth)`
\n   "DC" stands for "differentiable circuit".
\n2) `blockT = "CS"`
\n   "CS" stands for "cluster state".
\n`dBlocksPar::Array{Float64,1}`
\n   The array of the parameters (for differentiable blocks) in the circuit.
"""
function MPSbuilder(nBitA::Int64, vBit::Int64, rBit::Int64, blockT::Union{String, Tuple{String, Int64}})
    par2nd = MPSpar(nBitA, vBit, rBit)
    nBlock = par2nd.nBlock
    nBit = par2nd.nBit
    if typeof(blockT) == Tuple{String, Int64}
        if blockT[1] == "DC" # differentiable circuit.
            depth = blockT[2]
            swapV1 = chain(nBit, [put(nBit, (inBit,inBit+1)=>SWAP) for irBit=rBit:-1:1 for inBit=irBit  :(vBit+irBit-1)])
            swap1V = chain(nBit, [put(nBit, (inBit+1,inBit)=>SWAP) for irBit=1:   rBit for inBit=(vBit+irBit-1):-1:irBit])
            MeasureBlock = Measure(nBit, locs=(vBit+1):nBit, resetto=0)
            cBlockHead = chain(nBit, chain(nBit, DCbuilder(nBit, depth).circuit, swapV1) |> markDiff, MeasureBlock)
            cBlocks = [cBlockHead]
            cEBlocks = [subroutine( nBitA, cBlockHead[1], (nBitA-rBit-vBit+1):nBitA )]
            for i=2:nBlock-1
                cBlockHead = chain(nBit, DCbuilder(nBit, depth).circuit, swapV1) |> markDiff
                cBlock =  chain(nBit, chain(nBit, vcat([swap1V], cBlockHead.blocks)), MeasureBlock)
                push!(cBlocks, cBlock)
                # push!(cEBlocks, put(nBitA, Tuple((nBitA-i*rBit-vBit+1):(nBitA-(i-1)*rBit),)=>cBlockHead)) #Use `subroutine` instead of `put` for CuYao compatibility and efficiency.
                push!(cEBlocks, subroutine( nBitA, cBlockHead, (nBitA-i*rBit-vBit+1):(nBitA-(i-1)*rBit) ))
            end
            cBlockHead = chain(nBit, DCbuilder(nBit, depth).circuit, swapV1) |> markDiff
            cBlock =  chain(nBit, vcat([swap1V], cBlockHead.blocks))
            push!(cBlocks, cBlock)
            push!(cEBlocks, subroutine( nBitA, cBlockHead, (nBitA-nBlock*rBit-vBit+1):(nBitA-(nBlock-1)*rBit) ))
            circuit = chain(nBit, cBlocks)
            cExtend = chain(nBitA, cEBlocks)
        end
    end
    if typeof(blockT) == String
        if blockT == "CS" # MPS blocks for 1D cluster state.
            if vBit !=1 || rBit !=1
                error("The nBit=vBit+rBit of cluster state MPS blocks should be 2.\n")
            end
            cBlocks = ChainBlock[]
            MeasureBlock = Measure(nBit, locs=(vBit+1):nBit, resetto=0)
            push!(cBlocks, chain(nBit, chain(nBit, repeat(nBit,H,(2,1)), control(nBit, 1, nBit=>Z)), MeasureBlock))
            for i =2:nBlock-1
                push!(cBlocks, chain(nBit, chain(nBit,put(nBit, (2,1)=>SWAP), put(nBit, 1=>H), control(nBit, 1, nBit=>Z)), MeasureBlock))
            end
            push!(cBlocks, chain(nBit,put(nBit, (2,1)=>SWAP), put(nBit, 1=>H), control(nBit, 1, nBit=>Z)))
            circuit = chain(nBit,cBlocks)
            cExtend = chain(nBitA, vcat([repeat(nBitA,H,Tuple(nBitA:-1:1,))], [control(nBitA, i-1, i=>Z) for i=nBitA:-1:2]))
        end
    end

    MPSc(circuit, cExtend)
end


text:
struct MeshAccretionGeometry{T} <: AbstractAccretionGeometry{T}
    mesh::Vector{Tuple{SVector{3,T},SVector{3,T},SVector{3,T}}}
    x_extent::Tuple{T,T}
    y_extent::Tuple{T,T}
    z_extent::Tuple{T,T}
end

function MeshAccretionGeometry(mesh)
    static_mesh = map(mesh) do triangle
        Tuple(SVector(p[1], p[2], p[3]) for p in triangle)
    end
    MeshAccretionGeometry(static_mesh, bounding_box(mesh)...)
end

# naive implementation
function bounding_box(mesh::Union{GeometryBasics.Mesh{3,T}}) where {T}
    bounding_box(T, mesh)
end

function bounding_box(mesh::Vector{Tuple{SVector{3,T},SVector{3,T},SVector{3,T}}}) where {T}
    bounding_box(T, mesh)
end

function bounding_box(T, mesh)
    xmin = typemax(T)
    xmax = -typemax(T)
    ymin = typemax(T)
    ymax = -typemax(T)
    zmin = typemax(T)
    zmax = -typemax(T)
    for t in mesh
        for p in t
            let x = p[1], y = p[2], z = p[3]
                (x > xmax) && (xmax = x)
                (x < xmin) && (xmin = x)
                (y > ymax) && (ymax = y)
                (y < ymin) && (ymin = y)
                (z > zmax) && (zmax = z)
                (z < zmin) && (zmin = z)
            end
        end
    end
    (xmin, xmax), (ymin, ymax), (zmin, zmax)
end

function in_nearby_region(m::MeshAccretionGeometry{T}, line_element) where {T}
    p = line_element[2]
    @inbounds m.x_extent[1] < p[1] < m.x_extent[2] &&
              m.y_extent[1] < p[2] < m.y_extent[2] &&
              m.z_extent[1] < p[3] < m.z_extent[2]
end

function has_intersect(m::MeshAccretionGeometry{T}, line_element) where {T}
    for triangle in m.mesh
        dist_sq = sum(@.((triangle[1] - line_element[2])^2))
        # assume line element and mesh triangle are small; check if we're within a
        # certain distance before running jsr
        if dist_sq < 3.0 && jsr_algorithm(triangle..., line_element...)
            return true
        end
    end
    false
end

function collision_callback(m::MeshAccretionGeometry{T}) where {T}
    (u, λ, integrator) -> intersects_geometry(m, cartesian_line_element(u, integrator))
end

export MeshAccretionGeometry


text:
#Put your code here Sangchu.

using CSV, DataFrames

#Import monotonic data from CSV file.
data = CSV.File("/Users/crismoen/.julia/dev/StructuresKit/src/assets/Connections/Zhang_2020_steel_steel_screw_monotonic.csv")

#Interpolation

import XLSX

# Read fitst lab data Excel file. Please change the path if necessary
el = XLSX.readxlsx("D:/JHU Spring 2021/Fastener Testing Data_Zhidong/Trail1/Central_Summary_20200108_PhaseI.xlsx")

# Relate the experiment serial number with its filename
filenames1 = el["Sheet1!G2:G94"]
fileorders1 = el["Sheet1!A2:A94"]

# Read txt. files by order
for i in 1:size(filenames1,1)
    # Creat two temporary recording current data, one for displacement, the other for force
    disp = []
    force = []

    # Concatenate the file name with orders
    file2read = filenames1[i]
    # Please change the path if necessary
    filepath = "D:/JHU Spring 2021/Fastener Testing Data_Zhidong/Fastener Testing Data_Zhidong/"
    filename = filepath * file2read

    # Read the txt. lab data file
    f = open(filename, "r")
    lines = readlines(f, keep = true)

    # Only pick the displacement column and force column
    for line in lines

        spliteddata = split(line)
        choosedisp = spliteddata[[8],:]
        chooseforce = spliteddata[[10],:]

        # Record them
        append!(disp,choosedisp)
        append!(force,chooseforce)

    end

    close(f)

    # load the displacement data and force data in to DataFrame
    using DataFrames

    df = DataFrame(displacement = disp, forcemag = force)

    using CSV

    outputorder = string(i)
    outputtype = ".csv"
    outputfilename = outputorder * outputtype

    CSV.write(outputfilename, df, header = false)

end

##################################################
#     Same process for the second Excel file     #
##################################################

# Please change the path if necessary
el = XLSX.readxlsx("D:/JHU Spring 2021/Fastener Testing Data_Zhidong/Trail1/Central_Summary_20200108_Phase2&3.xlsx")

filenames2 = el["Sheet1!G2:G64"]

fileorders2 = el["Sheet1!A2:A64"]

for i in 1:size(filenames2,1)

    disp = []
    force = []
    file2read = filenames1[i]
    # Please change the path if necessary
    filepath = "D:/JHU Spring 2021/Fastener Testing Data_Zhidong/Fastener Testing Data_Zhidong/"

    filename = filepath * file2read

    f = open(filename, "r")

    lines = readlines(f, keep = true)

    for line in lines

        spliteddata = split(line)
        choosedisp = spliteddata[[8],:]
        chooseforce = spliteddata[[10],:]
        append!(disp,choosedisp)
        append!(force,chooseforce)
    end

    close(f)

    using DataFrames
    df = DataFrame(displacement = disp, forcemag = force)

    for k in 1:size(disp,1)
        disp[k] = parse(Float64,disp[k])
    end

    for j in 1:size(force,1)
        force[j] = parse(Float64,force[j])
    end

    Maxforce = findmax(force)
    println(Maxforce)
    println(disp[Maxforce[2]])
    println(i)

    using CSV

    outputorder = string(i+88)
    outputtype = ".csv"
    outputfilename = outputorder * outputtype
    CSV.write(outputfilename, df, header = false)

end

#############################
#  Plot the Monotonic data  #
#############################

using DataFrames
using CSV
using Plots

# Load the CSV file. Please change the path if necessary
el = DataFrame(CSV.File("D:/JHU Spring 2021/Fastener Testing Data_Zhidong/Trail1/Mon.csv"))

# Convert the already read CSV file into array
df_matrix = Matrix(el)

# Visualize the chosen node
for i in 1:size(df_matrix,1)

    # Load the 4 node's displacement and force
    F = df_matrix[i,6:9]
    D = df_matrix[i,2:5]
    #
    name = df_matrix[i,1]

    df = DataFrame(force = F, displacement = D)
    plot(D,F)
    Plots.display(plot(D,F,color ="red",title = name,xlabel = "disp",ylabel = "force"))
    sleep(1)
    println(i)
end


text:
"""
    function reduce_image_2D( image::Array{<:Real},
                              x_pixels::Int64, y_pixels::Int64)

Unflattens an image array to a 2D array of pixels.
"""
@inline @fastmath function reduce_image_2D( image::Array{<:Real},
                                            x_pixels::Int64, y_pixels::Int64)

    im_plot = zeros(y_pixels, x_pixels)
    k = 1
    @inbounds for i = 1:y_pixels, j = 1:x_pixels
        im_plot[j,i] = image[k,1]
        if image[k,2] > 0.0
            im_plot[j, i] /= image[k,2]
        end
        k += 1
    end
    return im_plot
end

"""
    function reduce_image_3D( image::Array{<:Real}, w_image::Array{<:Real},
                                            x_pixels::Int64, y_pixels::Int64, z_pixels::Int64)

Unflattens an image array to a 3D array of pixels.
"""
@inline @fastmath function reduce_image_3D( image::Array{<:Real},
                                            x_pixels::Int64, y_pixels::Int64, z_pixels::Int64)

    im_plot = zeros(z_pixels, y_pixels, x_pixels)
    m = 1
    @inbounds for i = 1:z_pixels, j = 1:y_pixels, k = 1:x_pixels

        im_plot[k,j,i] = image[m,1]

        if image[m,1] > 0.0
            im_plot[k, j, i] /= image[m,2]
        end
        m += 1
    end
    return im_plot
end


text: # This is the content


text:
module Export

export generate_html, default_index

import Pluto
import Pluto: ServerSession, generate_html
using HTTP
using Base64
using SHA
import Pkg

## CACHE

export try_fromcache, try_tocache

cache_filename(cache_dir::String, current_hash::String) = joinpath(cache_dir, HTTP.URIs.escapeuri(string(try_get_exact_pluto_version(), current_hash)) * ".plutostate")

function try_fromcache(cache_dir::String, current_hash::String)
    p = cache_filename(cache_dir, current_hash)
    if isfile(p)
        try
            open(Pluto.unpack, p, "r")
        catch e
            @warn "Failed to load statefile from cache" current_hash exception=(e,catch_backtrace())
        end
    end
end
try_fromcache(cache_dir::Nothing, current_hash) = nothing

function try_tocache(cache_dir::String, current_hash::String, state)
    mkpath(cache_dir)
    try
        open(cache_filename(cache_dir, current_hash), "w") do io
            Pluto.pack(io, state)
        end
    catch e
        @warn "Failed to write to cache file" current_hash exception=(e,catch_backtrace())
    end
end
try_tocache(cache_dir::Nothing, current_hash, state) = nothing

## FINDING THE PLUTO VERSION

function try_get_exact_pluto_version()
    try
        deps = Pkg.API.dependencies()

        p_index = findfirst(p -> p.name == "Pluto", deps)
        p = deps[p_index]

        if p.is_tracking_registry
            p.version
        elseif p.is_tracking_path
            error("Do not add the Pluto dependency as a local path, but by specifying its VERSION or an exact COMMIT SHA.")
        else
            # ugh
            is_probably_a_commit_thing = all(in(('0':'9') ∪ ('a':'f')), p.git_revision)
            if !is_probably_a_commit_thing
                error("Do not add the Pluto dependency by specifying its BRANCH, but by specifying its VERSION or an exact COMMIT SHA.")
            end

            p.git_revision
        end
    catch e
        @error "Failed to get exact Pluto version from dependency. Your website is not guaranteed to work forever." exception=(e, catch_backtrace())
        Pluto.PLUTO_VERSION
    end
end

function default_index(paths)
    """
    <!DOCTYPE html>

        <meta charset="utf-8">
        <meta name="viewport" content="width=device-width, initial-scale=1">

        <link rel="stylesheet" href="index.css">

Notebooks

        $(join(
            if link === nothing
                """
$(name) (Loading...)
"""
            else
                """
$(name)
"""
            end
            for (name,link) in paths
        ))

    """
end

end


text:
using GeoPhyInv
using Statistics
#md using Plots; gr();

# ### Medium
medium=Medium(:acou_homo1); # load a simple homogeneous acoustic medium from the gallery
update!(medium, [:vp,:rho], randn_perc=5); # add some random noise to the medium

# ### AGeom
ageom=AGeom(medium.mgrid,:xwell, SSrcs(2)); # load a simple acquisition using `mgrid` of the medium

# ### SrcWav
tgrid = range(0.0,stop=2.0,length=2000); # generate a time grid
wav = ricker(10.0, tgrid, tpeak=0.25,); # ricker wavelet
srcwav = SrcWav(tgrid, ageom, [:P]);
update!(srcwav, [:P], wav);

# ### Plotting #1
#md p1=heatmap(medium, :vp)
#md scatter!(ageom, SSrcs())
#md scatter!(ageom, Recs())
#md plot(p1)

# ### SeisForwExpt
# Now we have all the required variables to create `SeisForwExpt` object and
# prepare the forward modeling.
# While creating, we switched the `snaps_flag` on, and instruct recording field at
# `tsnaps`.
# Once the `Expt` object is created, do the modeling "without approximately any
# memory allocations" using `update!`

pa=SeisForwExpt(Fdtd(),medium=medium, ageom=ageom, srcwav=srcwav,
	snaps_flag=true, tsnaps=[0.3, 0.4, 0.5], tgrid=tgrid, verbose=true);

# ### Modeling
@time update!(pa);

# ### Extracting snaps
snaps=pa[:snaps,1]; # extracting snaps of the first supersource
snaps=pa[:snaps,2]; # second supersource

# ### Plotting #2

#md p=[]
#md for ii in 1:3
#md 	push!(p,heatmap(medium.mgrid[2], medium.mgrid[1], snaps[:,:,ii], xlabel="x", ylabel="z"))
#md end
#md plot(p..., layout=(1,3), size=(1100,250))


text:
function triangulate_with_cell_data(dataset::AbstractVTKUnstructuredData, OutputT = NTuple{3,Int}, reduction = :max)
    filter_cells!(dataset, [POINT_CELLS; LINE_CELLS])

    point_coords = dataset.point_coords
    point_data = dataset.point_data

    #Sorted inds and cell type are keys
    #Cell connectivity, cell count, and cell data are the values

    cell_register = Dict{NTuple{3, Int}, Tuple{NTuple{3, Int}, _Counter, typeof(dataset.cell_data)}}()
    for i in 1:length(dataset.cell_connectivity)
        _cells = triangulate_cell(dataset.cell_connectivity[i], dataset.cell_types[i])
        for j in 1:length(_cells)
            _key = sort(SVector{3,Int}(_cells[j])).data
            if haskey(cell_register, _key)
                cell_register[_key][2].a += 1
                for m in keys(dataset.cell_data)
                    _var_dim = var_dim(dataset, m, "Cell")
                    if _var_dim == 1
                        if reduction == :mean
                            cell_register[_key][3][m] += [dataset.cell_data[m][i]]
                        else
                            cell_register[_key][3][m] = max(cell_register[_key][3][m], [dataset.cell_data[m][i]])
                        end
                    else
                        if reduction == :mean
                            cell_register[_key][3][m] += dataset.cell_data[m][:,i]
                        else
                            cell_register[_key][3][m] = max.(cell_register[_key][3][m], dataset.cell_data[m][:,i])
                        end
                    end
                end
            else
                cell_register[_key] = (_cells[j], _Counter(1), empty(dataset.cell_data))
                for m in keys(dataset.cell_data)
                    _var_dim = var_dim(dataset, m, "Cell")
                    if _var_dim == 1
                        cell_register[_key][3][m] = [dataset.cell_data[m][i]]
                    else
                        cell_register[_key][3][m] = dataset.cell_data[m][:,i]
                    end
                end
            end
        end
    end

    ncells = length(cell_register)
    cell_data = empty(dataset.cell_data)
    _cell_connectivity = OutputT[]

    for m in keys(dataset.cell_data)
        _var_dim = var_dim(dataset, m, "Cell")
        if _var_dim == 1
            cell_data[m] = zeros(length(cell_register))
        else
            cell_data[m] = zeros(_var_dim, length(cell_register))
        end
    end

    for (i, kv) in enumerate(cell_register)
        k, v = kv
        for m in keys(dataset.cell_data)
            if reduction == :mean
                v[3][m] = v[3][m] ./ v[2].a
            end
            _var_dim = var_dim(dataset, m, "Cell")
            if _var_dim == 1
                cell_data[m][i] = v[3][m][1]
            else
                cell_data[m][:,i] = v[3][m]
            end
        end
        cc = OutputT <: Tuple ? v[1] : OutputT(v[1]...)
        push!(_cell_connectivity, cc)
    end

    _cell_types = [5 for i in 1:ncells]

    return VTKPolyData(point_coords, _cell_types, _cell_connectivity, point_data, cell_data)
end

function triangulate_no_cell_data(dataset::AbstractVTKUnstructuredData, OutputT = NTuple{3, Int})
    filter_cells!(dataset, [POINT_CELLS; LINE_CELLS])

    point_coords = dataset.point_coords
    point_data = dataset.point_data

    #Sorted inds and cell type are keys
    #Cell connectivity is the value
    cell_register = Dict{NTuple{3, Int}, NTuple{3, Int}}()
    for i in 1:length(dataset.cell_connectivity)
        _cells = triangulate_cell(dataset.cell_connectivity[i], dataset.cell_types[i])
        for j in 1:length(_cells)
            _key = sort(SVector{3,Int}(_cells[j])).data
            if !haskey(cell_register, _key)
                cell_register[_key] = _cells[j]
            end
        end
    end

    ncells = length(cell_register)
    cell_data = empty(dataset.cell_data)
    _cell_connectivity = OutputT[]

    for (i, kv) in enumerate(cell_register)
        k, v = kv
        cc = OutputT <: Tuple ? v : OutputT(v...)
        push!(_cell_connectivity, cc)
    end
    _cell_types = [5 for i in 1:ncells]

    return VTKPolyData(point_coords, _cell_types, _cell_connectivity, point_data, cell_data)
end

function remove_unused_vertices(dataset::AbstractVTKUnstructuredData)
    if length(dataset.cell_connectivity) == 0
        N = size(dataset.point_coords, 1)
        return typeof(dataset)(zeros(dim, 0), Int[], copy(dataset.cell_connectivity), empty(dataset.point_data), empty(dataset.cell_data))
    end
    I = length(dataset.cell_connectivity)
    J = length(dataset.cell_connectivity[1].data)
    all_node_inds = unique!(vec([Int(dataset.cell_connectivity[i].data[j]) for j in 1:J, i in 1:I]))
    nnodes = length(all_node_inds)
    ind_map = Dict(all_node_inds[i] => i for i in 1:nnodes)
    diff = size(dataset.point_coords, 2) - nnodes
    cell_connectivity = map(dataset.cell_connectivity) do cc
        getindex.(Ref(ind_map), Int.(cc))
    end
    point_coords = dataset.point_coords[:, all_node_inds]
    point_data = Dict(
        k => (
            temp = dataset.point_data[k];
            temp isa Matrix ? temp[:, all_node_inds] : temp[all_node_inds]
        ) for k in keys(dataset.point_data)
    )
    return typeof(dataset)(point_coords, copy(dataset.cell_types), cell_connectivity, point_data, deepcopy(dataset.cell_data))
end

function duplicate_vertices(dataset::AbstractVTKUnstructuredData)
    @assert all(isequal(5), dataset.cell_types)
    N = sum(length(dataset.cell_connectivity[i]) for i in 1:num_of_cells(dataset))
    point_coords = zeros(size(dataset.point_coords, 1), N)
    cell_connectivity = deepcopy(dataset.cell_connectivity)
    point_data = empty(dataset.point_data)
    for m in keys(dataset.point_data)
        _var_dim = var_dim(dataset, m, "Point")
        if _var_dim == 1
            point_data[m] = zeros(N)
        else
            point_data[m] = zeros(_var_dim, N)
        end
    end

    point_counter = 0
    for i in 1:num_of_cells(dataset)
        for j in dataset.cell_connectivity[i]
            point_counter += 1
            @views point_coords[:,point_counter] += dataset.point_coords[:,j]
            for m in keys(dataset.point_data)
                _var_dim = var_dim(dataset, m, "Point")
                if _var_dim == 1
                    point_data[m][point_counter] += dataset.point_data[m][j]
                else
                    @views point_data[m][:,point_counter] += dataset.point_data[m][:,j]
                end
            end
        end
        cell_connectivity[i] = (point_counter-2, point_counter-1, point_counter)
    end

    for m in keys(dataset.point_data)
        _var_dim = var_dim(dataset, m, "Point")
        if _var_dim == 1
            point_data[m] = point_data[m] ./ max.(point_counter, 1)
        else
            point_data[m] = point_data[m] ./ max.(point_counter', 1)
        end
    end
    cell_data = deepcopy(dataset.cell_data)
    cell_types = copy(dataset.cell_types)

    return VTKPolyData(point_coords, cell_types, cell_connectivity, point_data, cell_data)
end

function triangulate(dataset::AbstractVTKUnstructuredData, triangulate_cell_data=false, OutputT = NTuple{3,Int})
    if triangulate_cell_data
        return triangulate_with_cell_data(dataset, OutputT)
    else
        return triangulate_no_cell_data(dataset, OutputT)
    end
end


text:
for cd in (
    CategoricalDiscretizer(Dict(:A=>1, :B=>2, :C=>3)),
    CategoricalDiscretizer([:A, :B, :C], Int),
    CategoricalDiscretizer([:A, :B, :C])
    )

    @test encode(cd, :B) == 2
    @test encode(cd, "B") == 2
    @test encode(cd, [:C, :B, :A, :B, :C]) == [3,2,1,2,3]
    @test encode(cd, [:A :B; :C :A]) == [1 2; 3 1]
    @test_throws KeyError encode(cd, :D)
    @test_throws KeyError encode(cd, [:A, :D])

    @test decode(cd, 2) == :B
    @test decode(cd, 2 % UInt8 ) == :B
    @test decode(cd, [3,2,1,2,3]) == [:C, :B, :A, :B, :C]
    @test decode(cd, [1 2; 3 1]) == [:A :B; :C :A]

    @test nlabels(cd) == 3
    @test encoded_type(cd) == Int
    @test decoded_type(cd) == Symbol

    @test supports_encoding(cd, :A)
    @test supports_encoding(cd, :B)
    @test supports_encoding(cd, :C)
    @test !supports_encoding(cd, :D)

    @test supports_decoding(cd, 1)
    @test supports_decoding(cd, 2)
    @test supports_decoding(cd, 3)
    @test !supports_decoding(cd, 0)
    @test !supports_decoding(cd, 4)
end