text:
module DICOMTools

	export get_tag
	export get_transformation_matrix
	export dcm2array
	export load_dcms, load_dcms!
	export merge_multiseries
	export DICOMDict
	export DICOMSeries

	using DICOM
	using LinearAlgebra

	# Sorry for the bad commenting/documenting ... when I find the time I'll add it

	# Define a types to store DICOMs found in a directory structure
	struct DICOMSeries
		paths::Vector{String} # Individual dicom files
		dcms::Vector{DICOM.DICOMData} # The read in dicom structures
		pos::Vector{Int64} # Position in the files where reading stopped (see enum WhatToRead)
	end
	const DICOMDict = Dict{String, DICOMSeries} # SeriesInstanceUID => DICOMSeries

	# This enum can be used to define up to which group to read
	@enum WhatToRead read_all read_uid read_voxelsize skip_pixels
	# Convert above enum to a DICOM group
	function get_max_group(what::WhatToRead)::UInt16
		# Get DICOM.dcm_parse's max_group argument to stop reading early
		max_group = (what == read_all) ? 0xffff :
					(what == read_uid) ? 0x0020 :
					(what == read_voxelsize) ? 0x0028 : 0x7FDF # last option must be skip_pixels
		return max_group
	end

	# Load Dicoms recursively
	function load_dcms(directory::String, what::WhatToRead=read_all)::DICOMDict

		max_group = get_max_group(what)

		dcmdict = DICOMDict()
		for (root, _, files) in walkdir(directory)
			for filename in files
				fullpath = joinpath(root, filename)
				DICOM.isdicom(fullpath) || continue
				local dcm, pos
				open(fullpath, "r") do file
					dcm = dcm_parse(file; max_group)
					if eof(file)
						pos = -1
					else
						pos = position(file) - sizeof(UInt16) # Group tag was read in, revert position
					end
				end
				uid = dcm[tag"SeriesInstanceUID"]
				if haskey(dcmdict, uid)
					push!(dcmdict[uid].paths, fullpath)
					push!(dcmdict[uid].dcms, dcm)
					push!(dcmdict[uid].pos, pos)
				else
					dcmdict[uid] = DICOMSeries([fullpath], [dcm], [pos])
				end
			end
		end
		return dcmdict
	end

	# Read tags retrospectively
	function load_dcms!(dcmdict::DICOMDict, what::WhatToRead=read_all)::DICOMDict
		max_group = get_max_group(what)
		for series in values(dcmdict)
			load_dcms!(series, max_group)
		end
		return dcmdict
	end
	function load_dcms!(series::DICOMSeries, what::WhatToRead=read_all)::DICOMSeries
		max_group = get_max_group(what)
		load_dcms!(series::DICOMSeries, max_group)
		return series
	end
	function load_dcms!(series::DICOMSeries, max_group::UInt16)::DICOMSeries
		for (path, dcm, pos) in zip(series.paths, series.dcms, series.pos)
			pos == -1 && continue
			open(path, "r") do file
				seek(file, pos)
				DICOM.read_body!(file, dcm; max_group)
			end
		end
		return series
	end

	function filter(f::Function, dcmdict::DICOMDict)
		Base.filter(series -> f(series.second.dcms), dcmdict)
	end
	function filter!(f::Function, dcmdict::DICOMDict)
		Base.filter!(series -> f(series.second.dcms), dcmdict)
	end

	function get_tag(
		multiseries::AbstractVector{<: DICOMSeries},
		tag::Tuple{UInt16,UInt16},
		outtype::Val{T} = Val(typeof(multiseries[1].dcms[1][tag]))
	) where T
		T[ series.dcms[1][tag] for series in multiseries ]
	end
	function get_tags(
		multiseries::AbstractVector{<: DICOMSeries},
		tags::NTuple{N, Tuple{UInt16,UInt16}},
		outtypes::Val{T} = Val(ntuple(i -> typeof(multiseries[1].dcms[1][tags[i]]), N))
	) where {N, T}
		T[ series.dcms[1][tag] for tag in tags for series in multiseries ]
	end

	function get_transformation_matrix(series::DICOMSeries)::Matrix{Float64}
		# This breaks if InstanceNumber tag does not equal slice number
		# Gaps between slices are handled correclty if volume is rectangular

		dcms = series.dcms

		# Get origin of first and second slice
		# Translation of first slice is origin of the volume
		# Translation of second minus translation of first is vector pointing in slice direction.

		# Get "step"-vectors in row and column direction
		local rows_vector::Vector{Float64}, columns_vector::Vector{Float64}
		let
			orientation = dcms[1][(0x0020, 0x0037)] # [3 row cosines and 3 column cosines]
			# Note: Direction cosines of row point along x-direction
			pixelsize = dcms[1][(0x0028, 0x0030)] # [row pixel spacing and column pixel spacing]
			rows_vector = pixelsize[1] .* orientation[1:3]
			columns_vector = pixelsize[2] .* orientation[4:6]
		end

		local translation::Vector{Float64}, slices_vector::Vector{Float64}
		let first_slice = 0, second_slice = 0
			for (i, dcm) in enumerate(dcms)
				instance_number = dcm[(0x0020, 0x0013)]
				if instance_number == 1
					first_slice = i
				elseif instance_number == 2
					second_slice = i
				end
				first_slice != 0 && second_slice != 0 && break
			end
			translation = dcms[first_slice][(0x0020, 0x0032)] # [x, y, z]
			if second_slice > 0
				slices_vector = dcms[second_slice][(0x0020, 0x0032)] - translation
			else
				# Need to compute cross product, there is no second slice
				slices_vector = rows_vector × columns_vector
			end
			# Do not normalise, so that slice thickness is included
		end

		transformation = Matrix{Float64}(undef, 4, 4)
		transformation[1:3, 1] = columns_vector # This order is correct,
		# because the first column is multiplied with the row index
		transformation[1:3, 2] = rows_vector
		transformation[1:3, 3] = slices_vector
		transformation[1:3, 4] = translation
		transformation[4, 1:3] .= 0
		transformation[4, 4] = 1
		return transformation
	end

	function dcm2array(
		series::DICOMSeries,
		outtype::Val{T} = Val(Float64)
	)::Array{T, 3} where T
		# Note: This breaks if InstanceNumber tag does not equal slice number

		dcms = series.dcms

		# Get size
		reference = dcms[1]
		array_size = (
			get(reference, tag"Rows", 1),
			get(reference, tag"Columns", 1),
			maximum( get(dcm, tag"InstanceNumber", 1) for dcm in dcms )
		)

		# Fill array
		voxels = Array{T, 3}(undef, array_size...)
		for dcm in dcms
			data = dcm[tag"PixelData"]
			voxels[:, :, dcm[tag"InstanceNumber"]] = data
		end

		return voxels
	end

	"""

		Get an array extracted from multiple DICOM series.
		For example useful for MRI: fitting relaxation curves to series
		acquired with different inversion/echo times.

	"""
	function merge_multiseries(
		multiseries::AbstractVector{<: DICOMSeries},
		data_type::Val{T} = Val(eltype(multiseries[1].dcms[1]))
	)::Array where T
		# TODO: Would expect to return Array{..., 4} because dcm2array gives 3D volume.
		# However what happens with timeseries in 3D volumes?

		# Get arrays from dcms
		tojoin = [ dcm2array(dcms) for dcms in multiseries ]

		# Join arrays
		# Rearrange all volumes into a single array, first axis iterates over parameters,
		# this makes access to signal(parameter) easier.
		# Create empty array from the size of the first element in tojoin
		arrays = Array{
			T,
			1 + ndims(tojoin[1]) # parameter (1) and spatial dimensions
		}(
			undef,
			length(tojoin),
			size(tojoin[1])... # Use first element as reference
		)
		array_size = tuple((size(arrays, d) for d = 2:4)...)
		for (i, array) in enumerate(tojoin)
			@assert size(array) == array_size # Does size match?
			arrays[i, :, :, :] = array
		end
		return arrays
	end

end


text:
module print_information_example

using GcodeParser;

function run()

    # create any data object
    # it will be passed as a second parameter to your callbacks
    # here simple dictionary is used to store information during the print
    myPrinter = Dict{String,Any}();
    myPrinter["positioning"] = "absolute";
    myPrinter["x"] = 0.;
    myPrinter["y"] = 0.;
    myPrinter["z"] = 0.;
    myPrinter["e"] = 0.;
    myPrinter["filamentUsage"] = 0.; # store total filament usage (printed length of filament)
    myPrinter["distanceMoved"] = 0.; # store total distance moved

    # Setup a dictionary of callbacks for specified commands
    callbacks = Dict{String,Function}();
    callbacks["G0"] = movement; # just move the printhead
    callbacks["G1"] = extrude;  # move the printhead as well as extrude material

    # watch out for relative and absolute positioning
    callbacks["G90"] = (cmds, dataobject)->dataobject["positioning"] = "absolute";
    callbacks["G91"] = (cmds, dataobject)->dataobject["positioning"] = "relative";

    # parse g-code file and simulate print using our own callbacks and data object
    parseFile("examples/3dprint/gcodes/AI3M_test.gcode", callbacks, myPrinter);

    # Show printer data after print with some interesting stats
    @show myPrinter;

    return;
end

"""
    movement(cmds, dataobject)

Example movement callback for `G0` and `G1` which calculates the total distance moved in all axes.

It is calculated by watching the `X`, `Y` and `Z` axes movement.
"""
function movement(cmds, dataobject)

    dx = 0.;
    dy = 0.;
    dz = 0.;

    x = findfirst((x->lowercase(x.first) == "x"), cmds);
    if x !== nothing
        val = parse(Float64, cmds[x].second);

        if dataobject["positioning"] === "absolute"
            dx = val - dataobject["x"];
            dataobject["x"] = val;
        else
            dx = val;
        end
    end

    y = findfirst((x->lowercase(x.first) == "y"), cmds);
    if y !== nothing
        val = parse(Float64, cmds[y].second);

        if dataobject["positioning"] === "absolute"
            dy = val - dataobject["y"];
            dataobject["y"] = val;
        else
            dy = val;
        end
    end

    z = findfirst((x->lowercase(x.first) == "z"), cmds);
    if z !== nothing
        val = parse(Float64, cmds[z].second);

        if dataobject["positioning"] === "absolute"
            dz = val - dataobject["z"];
            dataobject["z"] = val;
        else
            dz = val;
        end
    end

    dataobject["distanceMoved"] += sqrt(dx * dx + dy * dy + dz * dz);
end

"""
    extrude(cmds, dataobject)

Example extrusion callback for `G1` which calculates total length of filament extruded.

The extruded filament length is obtained by watching the `E` axis movement in the g-code file.
"""
function extrude(cmds, dataobject)
    movement(cmds, dataobject);

    # calculate used filament length
    e = findfirst((x->lowercase(x.first) == "e"), cmds);
    if e !== nothing
        # Current E axis value
        e = parse(Float64, cmds[e].second);

        # Printed length of a current move
        if dataobject["positioning"] === "absolute"
            de = e - dataobject["e"];
            dataobject["e"] = e;
        else
            de = e;
        end

        # Used filament
        dataobject["filamentUsage"] += de;
        # println(dataobject["filamentUsage"]);5
    end
end

end # module

ext: jl


text:
# ---
# title: 679. 24 Game
# id: problem679
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Depth-first Search
# link: <https://leetcode.com/problems/24-game/description/>
# hidden: true
# ---
#
# You have 4 cards each containing a number from 1 to 9. You need to judge
# whether they could operated through `*`, `/`, `+`, `-`, `(`, `)` to get the
# value of 24.
#
# **Example 1:**
#
#
#
#     Input: [4, 1, 8, 7]
#     Output: True
#     Explanation: (8-4) * (7-1) = 24
#
#
# **Example 2:**
#
#
#
#     Input: [1, 2, 1, 2]
#     Output: False
#
#
# **Note:**
#
#   1. The division operator `/` represents real division, not integer division. For example, 4 / (1 - 2/3) = 12.
#   2. Every operation done is between two numbers. In particular, we cannot use `-` as a unary operator. For example, with `[1, 1, 1, 1]` as input, the expression `-1 - 1 - 1 - 1` is not allowed.
#   3. You cannot concatenate numbers together. For example, if the input is `[1, 2, 1, 2]`, we cannot write this as 12 + 12.
#
#
## @lc code=start
using LeetCode

function judgePoint24(nums::Vector)::Bool
    n = length(nums)
    if n == 1
        return nums[1] == 24
    else
        for i in 1:n
            for j in 1:n
                if j != i
                    rest = [nums[k] for k in 1:n if k != i && k != j]
                    if judgePoint24([rest..., nums[i] + nums[j]]) ||
                       judgePoint24([rest..., nums[i] - nums[j]]) ||
                       judgePoint24([rest..., nums[j] - nums[i]]) ||
                       judgePoint24([rest..., nums[i] * nums[j]]) ||
                       (!iszero(nums[j]) && judgePoint24([rest..., nums[i]//nums[j]])) ||
                       (!iszero(nums[i]) && judgePoint24([rest..., nums[j]//nums[i]]))
                        return true
                    end
                end
            end
        end
        false
    end
end
## @lc code=end


text:
# stoichiometric matrix: fluxes (in order) are v1-v6 (v6 is reverse of v5), b1-b20
S = [-1 0 0 1 2 -2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; #citrulline
     -1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;  #aspartate
     1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; #arginosuccinate
     0 1 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; #fumarate
     0 1 -1 0 -2 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; #arginine
     0 0 1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;  #ornithine
     0 0 1 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; #urea
     0 0 0 -1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; #carbamoyl phosphate
     -1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; #ATP
     1 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0; #AMP
     1 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0; #diphosphate
     0 0 0 1 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0; #phosphate
     0 0 0 0 -3 3 0 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0 0 0; #NADPH
     0 0 0 0 -3 3 0 0 0 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0; #H+
     0 0 0 0 -4 4 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0; #O2
     0 0 0 0 2 -2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 0 0 0 0; #NO2
     0 0 0 0 3 -3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 0 0; #NADP+
     0 0 -1 0 4 -4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1]; #H2O

#atom matrix: C H N O P S
A = [6 13 3 3 0 0; #citrulline
    4 7 1 4 0 0;  #aspartate
    10 18 4 6 0 0; #arginosuccinate
    4 4 0 4 0 0; #fumarate
    6 14 4 2 0 0; #arginine
    5 12 2 2 0 0; #ornithine
    1 4 2 1 0 0; #urea
    1 4 1 5 1 0; #carbamoyl phosphate
    10 16 5 13 3 0; #ATP
    10 14 5 7 1 0;  #AMP
    0 4 0 7 2 0;  #diphosphate
    0 3 0 4 1 0;  #phosphate
    21 30 7 17 3 0; #NADPH
    0 1 0 0 0 0; #H+
    0 0 0 2 0 0; #O2
    0 0 1 1 0 0; #NO
    21 29 7 17 3 0; #NADP+
    0 2 0 1 0 0] #H2O

e = transpose(A)*S

#look at only the balances for v1-v6, species are not balanced over boundaries
print(e[:,1:6])


text:
Constant(x)=1
ConstantLearned(train_data, train_labels) = Constant

ConstantLearnerVerbose(train_data, train_labels) = Constant , rand(), ones(size(train_data)[2])


text:
using Documenter, PowerSystems
using InfrastructureSystems
const PSYPATH = dirname(pathof(PowerSystems))

include(joinpath(@__DIR__, "src", "generate_validation_table.jl"))

makedocs(
    modules = [PowerSystems],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true",
    assets = [
        "assets/logo.png"
    ]),
    sitename = "PowerSystems.jl",
    authors = "Jose Daniel Lara, Daniel Thom and Clayton Barrows",
    pages = Any[ # Compat: `Any` for 0.4 compat
        "Introduction" => "index.md",
        "User Guide" => Any["man/data.md"],
        "Developer" => Any["Tests"=>"developer/tests.md",
                           "Logging"=>"developer/logging.md",
                           "Style Guide" => "developer/style.md"],
        "API" => Any[
                        "PowerSystems" => "api/PowerSystems.md"
        ]
    ]
)

deploydocs(
    repo = "github.com/NREL/PowerSystems.jl.git",
    branch = "gh-pages",
    target = "build",
    deps = Deps.pip("pygments", "mkdocs", "python-markdown-math"),
    make = nothing,
)


text:
module ForecastPlots

using Distributions, LinearAlgebra, Statistics
using Plots, Dates

include("acf.jl")
include("ccf.jl")
include("pacf.jl")
export acf, ccf, pacf

include("candle.jl")
export candle

include("dplot.jl")
export dplot

include("fplot.jl")
export fplot

include("splot.jl")
export splot

"""
Package: ForecastPlots

Collection of plot functionalities for time series analysis. The available plots are:

- `acf`:     Auto-Correlation plot
- `candle`:  Candelstick plot for stock prices
- `ccf`:     Cross-Correlation plot
- `dplot`:   Decomposition plot for Data, Trend, Seasonality and Remainder.
- `fplot`:   Multivariate forecasting plots with with prediction intervals.
- `pacf`:    Partial Auto-Correlation plot
- `splot`:   Seasonal plot, similar to `monthplot` in R.
"""
ForecastPlots

end


text:
# Particularly useful or beautiful grammar of graphics invocations.

"""
    plot(fs::Vector{T}, lower::Number, upper::Number, elements::ElementOrFunction...;
         mapping...) where T <: Base.Callable
"""
function plot(fs::Vector{T}, lower::Number, upper::Number, elements::ElementOrFunction...; mapping...) where T <: Base.Callable
    if isempty(elements)
        elements = ElementOrFunction[]
    elseif isa(elements, Tuple)
        elements = ElementOrFunction[elements...]
    end

    element_types = Set(map(typeof, elements))

    if !in(Guide.xlabel, element_types)
        push!(elements, Guide.xlabel("x"))
    end

    if !in(Guide.ylabel, element_types)
        push!(elements, Guide.ylabel("f(x)"))
    end

    if upper < lower
        push!(elements, Coord.cartesian(xflip=true))
    end

    mappingdict = Dict{Symbol, Any}(:y => fs, :xmin => [lower], :xmax => [upper])
    for (k, v) in mapping
        mappingdict[k] = v
    end

    plot(Stat.func, Geom.line, elements...; mappingdict...)
end

"""
    plot(f::Function, lower::Number, upper::Number, elements::ElementOrFunction...;
         mapping...)

Plot the function or expression `f`, which takes a single
argument or operates on a single variable, respectively, between the `lower`
and `upper` bounds.  See [`Stat.func`](@ref) and [`Geom.line`](@ref).
"""
plot(f::Function, lower::Number, upper::Number, elements::ElementOrFunction...; mapping...) =
        plot(Function[f], lower, upper, elements...; mapping...)

"""
    plot(f::Function, xmin::Number, xmax::Number, ymin::Number, ymax::Number,
         elements::ElementOrFunction...; mapping...)

Plot the contours of the 2D function or expression in `f`.
See [`Stat.func`](@ref) and [`Geom.contour`](@ref).
"""
function plot(f::Function, xmin::Number, xmax::Number, ymin::Number, ymax::Number,
              elements::ElementOrFunction...; mapping...)
    default_elements = ElementOrFunction[]
    element_types = Set(map(typeof, elements))

    if !in(Guide.xlabel, element_types)
        push!(default_elements, Guide.xlabel("x"))
    end

    if !in(Guide.ylabel, element_types)
        push!(default_elements, Guide.ylabel("y"))
    end

    if !in(Guide.ColorKey, element_types) && !in(Guide.ManualColorKey, element_types)
        push!(default_elements, Guide.colorkey(title="f(x,y)"))
    end

    push!(default_elements, Coord.cartesian(xflip=xmin > xmax, yflip=ymin > ymax))

    plot(layer(f, xmin, xmax, ymin, ymax,
               elements...; mapping...), default_elements...)
end

"""
    layer(fs::Vector{T}, lower::Number, upper::Number,
          elements::ElementOrFunction...) where T <: Base.Callable -> [Layers]

Create a layer from a list of functions or expressions in `fs`.
"""
layer(fs::Vector{T}, lower::Number, upper::Number, elements::ElementOrFunction...) where T <: Base.Callable =
    layer(y=fs, xmin=[lower], xmax=[upper], Stat.func, Geom.line, elements...)

"""
    layer(f::Function, lower::Number, upper::Number,
          elements::ElementOrFunction...) -> [Layers]

Create a layer from the function or expression `f`, which takes a single
argument or operates on a single variable, respectively, between the `lower`
and `upper` bounds.  See [`Stat.func`](@ref) and [`Geom.line`](@ref).
"""
layer(f::Function, lower::Number, upper::Number, elements::ElementOrFunction...) =
        layer(Function[f], lower, upper, elements...)

"""
    layer(f::Function, xmin::Number, xmax::Number, ymin::Number, ymax::Number,
          elements::ElementOrFunction...; mapping...) -> [Layers]

Create a layer of the contours of the 2D function or expression in `f`.
See [`Stat.func`](@ref) and [`Geom.contour`](@ref).
"""
function layer(f::Function, xmin::Number, xmax::Number, ymin::Number, ymax::Number,
               elements::ElementOrFunction...; mapping...)
    if isempty(elements)
        elements = ElementOrFunction[]
    elseif isa(elements, Tuple)
        elements = ElementOrFunction[elements...]
    end

    mappingdict = Dict{Symbol, Any}(:z => f,
                                    :xmin => [xmin], :xmax => [xmax],
                                    :ymin => [ymin], :ymax => [ymax])
    for (k, v) in mapping
        mappingdict[k] = v
    end

    layer(Geom.contour, elements...; mappingdict...)
end

### why `spy` and not `plot`?
"""
    spy(M::AbstractMatrix, elements::ElementOrFunction...; mapping...) -> Plot

Plots a heatmap of `M`, with M[1,1] in the upper left.  `NaN` values are
left blank, and an error is thrown if all elements of `M` are `NaN`.  See
[`Geom.rectbin`](@ref) and [`Coord.cartesian(fixed=true)...)`](@ref
Gadfly.Coord.cartesian).
"""
function spy(M::AbstractMatrix, elements::ElementOrFunction...; mapping...)
    is, js, values = _findnz(x->!isnan(x), M)
    n,m = size(M)
    plot(x=js, y=is, color=values,
        Coord.cartesian(yflip=true, fixed=true, xmin=0.5, xmax=m+.5, ymin=0.5, ymax=n+.5),
        Scale.color_continuous,
        Geom.rectbin,
        Scale.x_continuous,
        Scale.y_continuous,
        elements...; mapping...)
end

# Finds the subscripts and values where the predicate returns true.
#
# It takes a predicate (`testf`), a matrix (`A`) and returns a tuple `(is, js, zs)`,
# where `is`,'js' and 'zs' are arrays of subscripts and values where the predicate
# returned true.
#
# This function is used by spy()
# Hopefully at some point something like this will be in the standard library
# and then this can be removed (https://github.com/JuliaLang/julia/pull/9340).
#
function _findnz(testf::Function, A::AbstractMatrix{T}) where T
    N = Base.count(testf, A)
    is = zeros(Int, N)
    js = zeros(Int, N)
    zs = Array{T}(undef, N)
    if N == 0
        return (is, js, zs)
    end
    count = 1
    for j=1:size(A,2), i=1:size(A,1)
        Aij = A[i,j]
        if testf(Aij)
            is[count] = i
            js[count] = j
            zs[count] = Aij
            count += 1
        end
    end
    return (is, js, zs)
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export progress

"""
    progress(;kwargs...)

A Progress component.
Give user feedback for status of the task. For more information, see: https://mantine.dev/core/progress/
Keyword arguments:
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks
- `className` (String; optional): Often used with CSS to style elements with common properties
- `color` (a value equal to: "dark", "gray", "red", "pink", "grape", "violet", "indigo", "blue", "cyan", "teal", "green", "lime", "yellow", "orange"; optional): Progress color from theme
- `radius` (a value equal to: "xs", "sm", "md", "lg", "xl" | Real; optional): Predefined progress radius from theme.radius or number for height in px
- `sections` (optional): Replaces value if present, renders multiple sections instead of single one. sections has the following type: Array of lists containing elements 'value', 'color'.
Those elements have the following types:
  - `value` (Real; optional)
  - `color` (a value equal to: "dark", "gray", "red", "pink", "grape", "violet", "indigo", "blue", "cyan", "teal", "green", "lime", "yellow", "orange"; optional)s
- `size` (a value equal to: "xs", "sm", "md", "lg", "xl"; optional): Predefined progress height or number for height in px
- `striped` (Bool; optional): Adds stripes
- `value` (Real; optional): Current value for controlled slider
"""
function progress(; kwargs...)
        available_props = Symbol[:id, :className, :color, :radius, :sections, :size, :striped, :value]
        wild_props = Symbol[]
        return Component("progress", "Progress", "dash_mantine_components", available_props, wild_props; kwargs...)
end


text:
const input = Int[0,1,4,13,15,12,16]

function init(input::Vector{Int})
    state = Dict{Int,Vector{Int}}()
    for (i, val) in enumerate(input)
        state[val] = [i]
    end
    return state
end

function turn!(state::Dict{Int,Vector{Int}}, value::Int, now::Int)
    new_val = 0
    already_said = keys(state)
    if value ∈ already_said
        prev_turns = state[value]
        if length(prev_turns) >= 2
            new_val = prev_turns[end] - prev_turns[end - 1]
        end
    end
    if new_val ∉ already_said
        state[new_val] = [now]
    else
        push!(state[new_val], now)
    end
    return state, new_val
end

function play(input::Vector{Int}, until::Int)
    state = init(input)
    val = last(input)
    start = length(input) + 1
    for now in start:until
        state, val = turn!(state, val, now)
    end
    return val
end

println("Part 1: $(play(input, 2020))")
println("Part 2: $(play(input, 30000000))")


text:
#########
# use KissMCMC
#########
"""
Numerical parameters for MCMC

```
struct MCMCNum<:APara @deftype Int
    alg=[:met, :emcee, :emceep][2] # which MCMC stepper to use
    niter = 10^5 # number of MC steps
    nthin = 10     # only store theta every so often
    nburnin = niter÷3 # how much burn-in
    nchains = max(min(niter÷1000,100),2) # number of chains used in the emcee sampler
    sigma_ppdf::Vector{F}
    sample_ppdf::Function = theta_vec -> [randn()*sigma_ppdf[i]+theta_vec[i] for i=1:length(theta_vec)]

    # emcee
    drop_low_accept_ratio::Bool = true
    ball_radius::Any = 0.1 # can be a number or vector

    # output
    plotyes::Bool = false
    verbose::Bool = false
end
```
"""
@with_kw struct MCMCNum<:APara @deftype Int
    alg::Symbol=[:metro, :emcee, :emceep][2] # which MCMC stepper to use
    niter = 10^5 # number of MC steps
    nthin = 10     # only store theta every so often
    nburnin = niter÷2 # how much burn-in
    @assert nburnin<niter
    nchains = max(min(niter÷3000,100),11) # number of chains used in the emcee sampler
    sigma_ppdf::Vector{F} = F[] # needs to be length(theta0.th0)
    sample_ppdf::Function = theta_vec -> [randn()*sigma_ppdf[i]+theta_vec[i] for i=1:length(theta_vec)]

    # emcee
    drop_low_accept_ratio::Bool = true
    ball_radius::Any = 0.1 # can be a number or vector
    a_scale::Float64 = 2.0 # jump size, probably needn't be adjusted

    # output
    plotyes::Bool = false
    verbose::Bool = false
    use_progress_meter::Bool = true
end

"""
Results of an MCMC run.

Parameterised on T to allow nothing, when only the prior was sampled.

```
@with_kw struct MCMCSol{T} @deftype T
    thetas::Matrix{Float64}
    reshape_revert::Any
    theta_expect::Vector{Float64}
    theta_mode::Vector{Float64}
    thetas_std::Vector{Float64}
    thetas_minmax::Matrix{Float64}
    blobs::Any
    logdensity_vals::Vector{Float64}

    vols::Union{Void,Vector{Float64}}
    vols_above::Union{Void,Vector{Float64}}
    mc_hs1d
    mc_ivs1d
    mh2d
    sh2d
    miv2d
    siv2d
    min_h2d
    max_h2d
    min_iv2d
    max_iv2d

    # solver performance
    solverstats::Dict{Symbol,Any} # wall_time, eff_sample_size, τ, τ_convergence, accept_ratio, # of samples in thetas

    # inputs
    theta0::Theta0{Float64}
    pmcmc::MCMCNum
    logdensity::Function # density which was sampled

    # Holds various things, such as:
    # misc=Dict(:logposterior=>logposterior,
    #       :logprior=>logprior,
    #       :logprior_debug=>logprior_debug,
    #       :loglikelihood=>loglikelihood,
    #       :fwdm_fn=>fwdm_fn,
    #       :fit_target=>fit_target))
    misc::Dict{Symbol,Any}
end
```
"""
@with_kw struct MCMCSol{T} @deftype T
    thetas::Matrix{Float64}
    reshape_revert::Any
    theta_expect::Vector{Float64}
    theta_mode::Vector{Float64}
    thetas_std::Vector{Float64}
    thetas_minmax::Matrix{Float64}
    blobs::Any
    logdensity_vals::Vector{Float64}

    vols::Union{Void,Vector{Float64}}
    vols_above::Union{Void,Vector{Float64}}
    mc_hs1d
    mc_ivs1d
    mh2d
    sh2d
    miv2d
    siv2d
    min_h2d
    max_h2d
    min_iv2d
    max_iv2d

    # solver performance
    solverstats::Dict{Symbol,Any}  # wall_time, eff_sample_size, τ, τ_convergence, accept_ratio, # of samples in thetas

    # inputs
    theta0::Theta0{Float64}
    pmcmc::MCMCNum
    logdensity::Function # density which was sampled

    # Holds various things, such as:
    # misc=Dict(:logposterior=>logposterior,
    #       :logprior=>logprior,
    #       :logprior_debug=>logprior_debug,
    #       :loglikelihood=>loglikelihood,
    #       :fwdm_fn=>fwdm_fn,
    #       :fit_target=>fit_target))
    misc::Dict{Symbol,Any}
end

function bite(logdensity, theta0::Theta0, pmcmc::MCMCNum;
              resample_logdensity=nothing,
              n_resamples=300,
              misc=Dict{Symbol,Any}())

    sol = mcmc(logdensity, theta0, pmcmc; misc=misc)
    if resample_logdensity!=nothing
        sol = BM.retrace_samples(resample_logdensity, sol, pmcmc, nsamples=n_resamples, misc=sol.misc)
    end

    return sol
end

"""
    mcmc(logposterior, theta0::Theta0, pmcmc::MCMCNum; verbose=true, use_progress_meter=true)

Runs the MCMC.

Input:
- logdensity: log-density to sample
- theta0: one value for metropolis, tuple for emcee (theta0, ball_radius)
KW:
- misc : if passed in it will be stored in the sol.misc dictionary

Output:
- thetas, theta_expect, theta_mode, accept_ratio, blobs, logdensity_value: results
- mc_hs1d, mc_ivs1d : all 1D results
- mh2d, sh2d, miv2d, siv2d : 2D mean and standard dev
- eff_sample_size, τ, τ_convergence: from KissMCMC.eff_samples
- reshape_revert: to undo the chain squash
"""
function mcmc(logdensity, theta0::Theta0, pmcmc::MCMCNum;
              verbose=pmcmc.verbose, use_progress_meter=pmcmc.use_progress_meter,
              misc=Dict{Symbol,Any}())
    tic = time()
    @unpack alg,a_scale = pmcmc
    gl = theta0.gb.gl
    p = pmcmc
    blob_reduce! = if KissMCMC.hasblob(logdensity, theta0.th0)
        _blob_reduce!
    else
        KissMCMC.default_blob_reduce!
    end

    if alg==:emcee
        res_ec = emcee(logdensity, (theta0.th0, p.ball_radius), niter=p.niter,
                       nburnin=p.nburnin, nchains=p.nchains,
                       blob_reduce! = blob_reduce!,
                       use_progress_meter=use_progress_meter,
                       a_scale=a_scale)

        # get metrics
        eff_sample_size, τ, τ_convergence, eff_sample_sizes = eff_samples(res_ec[1])

        thetas, accept_ratio, blobs, logdensity_vals, reshape_revert = squash_chains(res_ec... ,
                                                                                drop_low_accept_ratio=p.drop_low_accept_ratio,
                                                                                blob_reduce! = blob_reduce!,
                                                                                verbose=verbose)
    elseif alg==:metro
        thetas, accept_ratio, blobs, logdensity_vals = metropolis(logdensity, p.sample_ppdf, theta0.th0,
                                                                  niter=p.niter, nthin=p.nthin, nburnin=p.nburnin,
                                                                  blob_reduce! = blob_reduce!)
        eff_sample_size, τ, τ_convergence, nthin, reshape_revert = nothing, nothing, nothing, nothing, nothing
    else
        error("alg unknown.  Note parallel algs may not work currently due to closures.")
    end
    mc_vol, mc_vol_above, mc_hs1d, mc_ivs1d,
    mh2d, sh2d, miv2d, siv2d,
    min_h2d, max_h2d, min_iv2d, max_iv2d #=,hquant, ivquant=# = unpack_blob(blobs, size(gl.dem.v))

    theta_expect = squeeze(mean(thetas,2),2)
    theta_mode = thetas[:,findmax(logdensity_vals)[2]]

    solverstats = Dict{Symbol,Any}()
    nsamples = size(thetas, 2)
    @pack solverstats = eff_sample_size, τ, τ_convergence, eff_sample_sizes, accept_ratio, nsamples

    solverstats[:wall_time] = time() - tic

    return MCMCSol(thetas, reshape_revert,
                   theta_expect, theta_mode, std(thetas,2)[:], hcat([[e...] for e in extrema(thetas,2)]...)',
                   blobs, logdensity_vals,
                   mc_vol, mc_vol_above,
                   mc_hs1d, mc_ivs1d,
                   mh2d, sh2d, miv2d, siv2d, min_h2d, max_h2d, min_iv2d, max_iv2d,
                   solverstats,
                   theta0, pmcmc, logdensity, misc)
end

"""
    retrace_samples(logdensity, sol_in, pmcmc::MCMCNum;
                         nsamples = 200,
                         verbose=true, use_progress_meter=true)

When just sampling the prior, this function will re-sample a certain number of
samples with running the forward model (thus logdensity will tipically be logposterior).
This thus allows to do the predictions
of ice thickness whilst only sampling the priors (which is much cheaper).

Note:
- that of the logdensity, only the blob is used.
- this means also that the statistics, e.g. theta_expect, is calculated from sol_in.theta
  and not from just the thetas used in re-sampling.  This is maybe slightly inconsistent.
"""
function retrace_samples(logdensity, sol_in, pmcmc::MCMCNum;
                         nsamples=200,
                         verbose=pmcmc.verbose, use_progress_meter=pmcmc.use_progress_meter,
                         misc=Dict{Symbol,Any}())

    @unpack theta0 = sol_in
    thetas_in = sol_in.thetas
    p = pmcmc
    demsize = size(theta0.gb.gl.dem.v)

    out = _retrace_samples(logdensity, thetas_in, pmcmc, demsize,
                           nsamples, verbose, use_progress_meter,
                           sol_in.logdensity_vals, sol_in.solverstats)

    return MCMCSol(out..., theta0, pmcmc, logdensity, misc)
end
# Split out logic such that it can be used in retrace_samples_condensed_MCMCSol as well:
function _retrace_samples(logdensity, thetas_orig, pmcmc::MCMCNum, demsize,
                          nsamples, verbose, use_progress_meter,
                          logdensity_vals,  # we want the density of the original pdf not the resampled one
                          solverstats)
    tic = time()

    # Those we want of the original thetas as that is more accurate:
    theta_expect = squeeze(mean(thetas_orig,2),2)
    theta_mode = thetas_orig[:,findmax(logdensity_vals)[2]]

    # Now thin the samples:
    inds = 1:(max(1, size(thetas_orig,2)÷nsamples)):size(thetas_orig,2)
    thetas_in = thetas_orig[:,inds]

    @assert KissMCMC.hasblob(logdensity, thetas_in[:,1])

    thetas, accept_ratio, blobs, logdensity_vals_ = KissMCMC.retrace_samples(logdensity, thetas_in,
                                                                            blob_reduce! = _blob_reduce!,
                                                                            use_progress_meter=use_progress_meter)

    @assert thetas_in == thetas # otherwise something went wrong

    (mc_vol, mc_vol_above, mc_hs1d, mc_ivs1d,
     mh2d, sh2d, miv2d, siv2d,
     min_h2d, max_h2d, min_iv2d, max_iv2d #=,hquant, ivquant=#) = unpack_blob(blobs, demsize)

    solverstats[:wall_time_resample] = time() - tic

    # return as the first arguments to MCMCSol
    return (thetas, nothing,
            theta_expect, theta_mode, std(thetas_orig,2)[:], hcat([[e...] for e in extrema(thetas_orig,2)]...)',
            blobs,
            logdensity_vals[inds], # logdensity_vals needs to be original density
            mc_vol, mc_vol_above,
            mc_hs1d, mc_ivs1d,
            mh2d, sh2d, miv2d, siv2d, min_h2d, max_h2d, min_iv2d, max_iv2d,
            solverstats)
end

############
# Blobs
###########
# see https://github.com/mauro3/KissMCMC.jl/blob/master/examples/notebooks/bayesian-ex-predictions-v2.ipynb

# Two methods needed to intitialize the storage-blob:
function _blob_reduce!(new_blob, niter::Int, nburnin::Int, nthin::Int)
    vol, vol_above, hs1d, ivs1d, hs2d, ivs2d = new_blob
    n = length(hs2d)
    # store 0D results
    vol_out = Array{eltype(hs1d)}((niter-nburnin)÷nthin)
    vol_above_out = Array{eltype(hs1d)}((niter-nburnin)÷nthin)
    # store 1D results
    h1d_out = Array{eltype(hs1d)}(length(hs1d), (niter-nburnin)÷nthin)
    iv1d_out = Array{eltype(ivs1d)}(length(ivs1d), (niter-nburnin)÷nthin)
    # store online-stats: mean, variance, extremas and quantiles
    h2d_out = Series(Group([Mean() for i=1:n]),
                     Group([Variance() for i=1:n]),
                     Group([Extrema() for i=1:n]),
                     #Group([Quantile([0.25,0.5,0.75]) for i=1:n]) # uses too much RAM
                     )
    iv2d_out = Series(Group([Mean() for i=1:n]),
                      Group([Variance() for i=1:n]),
                      Group([Extrema() for i=1:n]),
                      #Group([Quantile([0.25,0.5,0.75]) for i=1:n])
                      )
    return vol_out, vol_above_out, h1d_out, iv1d_out, h2d_out, iv2d_out
end
_blob_reduce!(new_blob, niter::Int, nburnin::Int, nthin::Int, nchains::Int) =
             [_blob_reduce!(new_blob,niter,nburnin,nthin) for i=1:nchains]

# Two methods to update it
function _blob_reduce!(stored_blob, new_blob, ni::Int)
    vol_out, vol_above_out, h1d_out, iv1d_out, h2d_out, iv2d_out = stored_blob
    vol, vol_above, hs1d, ivs1d, hs2d, ivs2d = new_blob
    vol_out[ni] = vol
    vol_above_out[ni] = vol_above
    h1d_out[:,ni] = hs1d
    iv1d_out[:,ni] = ivs1d
    n = length(hs2d)
    fit!(h2d_out, reshape(hs2d, n))
    fit!(iv2d_out, reshape(ivs2d, n))
    nothing
end
_blob_reduce!(stored_blob, new_blob, ni::Int, nc::Int) = _blob_reduce!(stored_blob[nc], new_blob, ni);

# One method to merge several chains
function _blob_reduce!(stored_blobs, chains2keep)
    # chain 1 has one stored_blobs[1], etc.
    vol_out, vol_above_out, h1d_out, iv1d_out, h2d_out, iv2d_out = stored_blobs[chains2keep[1]]
    n = size(h1d_out,1)
    h1d_out = reshape(h1d_out, length(h1d_out))
    iv1d_out = reshape(iv1d_out, length(iv1d_out))
    for i in chains2keep[2:end]
        append!(vol_out, stored_blobs[i][1])
        append!(vol_above_out, stored_blobs[i][2])
        append!(h1d_out, stored_blobs[i][3])
        append!(iv1d_out, stored_blobs[i][4])
        merge!(h2d_out, stored_blobs[i][5])
        merge!(iv2d_out, stored_blobs[i][6])
    end
    h1d_out = reshape(h1d_out, n, length(h1d_out)÷n)
    iv1d_out = reshape(iv1d_out, n, length(iv1d_out)÷n)
    return vol_out, vol_above_out, h1d_out, iv1d_out, h2d_out, iv2d_out
end

"""
Returns
- all samples of hs1d and iv1d
- mean, std, min, and max of the 2D fields
- h2d_stats, iv2d_stats: OnlineStats Series

hs1d, iv1d, mh2d, sh2d, miv2d, siv2d,
min_h2d, max_h2d, min_iv2d, max_iv2d,
# hquant, ivquant,
h2d_stats, iv2d_stats  = unpack_blob(blob, sz)
"""
function unpack_blob(blob, sz)
    vol, vol_above, h1d, iv1d, h2d, iv2d = blob
    hmean, hvar, hext#=, hquant=# = [value.(u) for u in value(h2d)]
    ivmean, ivvar, ivext#=, ivquant=# = [value.(u) for u in value(iv2d)]
    hmin = reshape([h[1] for h in hext], sz)
    hmax = reshape([h[2] for h in hext], sz)
    ivmin = reshape([iv[1] for iv in ivext], sz)
    ivmax = reshape([iv[2] for iv in ivext], sz)
    return tuple(vol, vol_above, h1d, iv1d,
                 map(x->reshape(x,sz), (hmean, sqrt.(hvar), ivmean, sqrt.(ivvar)))...,
                 hmin,hmax,ivmin,ivmax, #hquant, ivquant,
                 h2d, iv2d)
end
# nice:
unpack_blob(blob::Void, sz) = (nothing, nothing,
                               nothing,nothing,nothing,nothing,nothing,nothing,
                               nothing,nothing,nothing,nothing,nothing,nothing)

##########
# Parameters
##########


text:
# image plots
function imagesc(x::Coord,y::Coord,Z::Coord;
                 handle::Handle    = gnuplot_state.current,
                 clim::Vector{Int} = [0, 255],
                 args...)

    # process arguments
    PA = PlotArgs()
    for k in keys(args)
        # substitute synonyms
        key = k
        for s in syn
            key ∈ s && (key = s[1]; break)
        end
        # store arguments
        for f in fieldnames(PlotArgs)
            f == key && (setfield!(PA, f, string(args[k])); break)
        end
    end

    # validation
    valid_range(PA.xrange)
    valid_range(PA.yrange)
    length(clim) == 2 || throw(DomainError(length(clim), "clim must be a 2-element vector."))
    length(x) == Base.size(Z)[2] || throw(DimensionMismatch("invalid coordinates."))
    length(y) == Base.size(Z)[1] || throw(DimensionMismatch("invalid coordinates."))
    2 <= ndims(Z) <= 3 || throw(DimensionMismatch(ndims(Z), "Z must have two or three dimensions"))

    handle = figure(handle, redraw = false)
    clearfigure(handle)

    # create figure configuration
    ndims(Z) == 2 ? plotstyle = "image" : plotstyle = "rgbimage"
    if ndims(Z) == 3
        Z[:] = Z.-clim[1]
        Z[Z.<0] .= 0.0
        Z[:] = Z.*255.0/(clim[2]-clim[1])
        Z[Z.>255] .= 255
    end

    term = config[:term][:terminal]
    PA.font == "" && (PA.font = TerminalDefaults[term][:font])
    PA.size == "" && (PA.size = TerminalDefaults[term][:size])
    PA.background == "" && (PA.background = "white")

    tc = TermConf(font       = PA.font,
                  size       = PA.size,
                  background = "white")

    ac = AxesConf(title  = PA.title,
                  xlabel = PA.xlabel,
                  ylabel = PA.ylabel,
                  xrange = PA.xrange,
                  yrange = PA.yrange)

    cc = CurveConf(plotstyle = plotstyle)

    c = Curve(x    = x,
              y    = y,
              Z    = Z,
              conf = cc)

    push_figure!(handle,tc,ac,c,PA.gpcom)
    return gnuplot_state.figs[findfigure(handle)]
end
imagesc(Z::Coord;args...) = imagesc(1:size(Z)[2],1:size(Z)[1],Z;args...)


text:
visualize(args...; kwargs...) = nothing
visualize!(args...; kwargs...) = nothing
visualize_sequence(args...; kwargs...) = nothing

is_kwarg(arg_or_kwarg::Symbol) = false
is_kwarg(arg_or_kwarg::Expr) = (arg_or_kwarg.head == :parameters)

function has_kwargs(args_kwargs::Vector)
  isempty(args_kwargs) && return false
  return is_kwarg(first(args_kwargs))
end

function get_kwargs(args_kwargs::Vector)
  @assert has_kwargs(args_kwargs)
  return first(args_kwargs)
end

function get_kwarg(kwargs::Expr, key::Symbol)
  n = findfirst(kw -> kw.args[1] == :sequence, kwargs.args)
  if !isnothing(n)
    @assert kwargs.args[n].head == :kw
    return esc(kwargs.args[n].args[2])
  end
  return nothing
end

function args_kwargs(ex::Vector)
  kwargs = has_kwargs(ex) ? get_kwargs(ex) : :()
  args = has_kwargs(ex) ? ex[2:end] : ex
  return args, kwargs
end

function function_args_kwargs(ex::Symbol)
  func = :identity
  args = [ex]
  kwargs = :()
  iscollection = true
  return func, args, kwargs, iscollection
end

function function_args_kwargs(ex::Expr)
  if ex.head == :call
    func = first(ex.args)
    args, kwargs = args_kwargs(ex.args[2:end])
    iscollection = true
  elseif ex.head == :ref
    #func, args, kwargs, iscollection = function_args_kwargs(Symbol(ex.args))
    func = :identity
    args = [ex]
    kwargs = :()
    iscollection = false
  else
    dump(ex)
    error("Visualizing expression $ex not supported right now.")
  end
  return func, args, kwargs, iscollection
end

expr_to_string(s::Symbol) = String(s)
expr_to_string(ex::Expr) = String(repr(ex))[3:(end - 1)]

# Take the symbols of the arguments and output
# the labels if there are multiple inputs or
# the prefix for the labels if there is only
# one input.
function vertex_labels_kwargs(args, iscollection)
  if iscollection && isone(length(args))
    vertex_labels_kw = :vertex_labels_prefix
    vertex_labels_arg = string(only(args))
  else
    vertex_labels_kw = :vertex_labels
    vertex_labels_arg = string.(args)
  end
  return vertex_labels_kw, vertex_labels_arg
end

function func_args_sequence_kwargs(ex, vis_kwargs...)
  func, args, kwargs, iscollection = function_args_kwargs(ex)
  sequence = get_kwarg(kwargs, :sequence)
  vertex_labels_kw, vertex_labels_arg = vertex_labels_kwargs(args, iscollection)
  # Merge labels kwarg with kwargs
  vis_kwargs_dict = Dict([
    vis_kwarg.args[1] => vis_kwarg.args[2] for vis_kwarg in vis_kwargs
  ])
  vertex_labels_kwarg_dict = Dict(vertex_labels_kw => vertex_labels_arg)
  merged_kwargs_dict = merge(vertex_labels_kwarg_dict, vis_kwargs_dict)
  merged_kwargs_expr = [:($k = $v) for (k, v) in pairs(merged_kwargs_dict)]
  return func, esc.(args), sequence, esc.(merged_kwargs_expr)
end

function visualize_expr(vis_func, ex::Union{Symbol,Expr}, vis_kwargs::Expr...)
  func, args, sequence, kwargs = func_args_sequence_kwargs(ex, vis_kwargs...)
  e = quote
    $(vis_func)($(func), ($(args...),), $(sequence); $(kwargs...))
  end
  return e
end

function visualize_expr!(fig, vis_func!, ex::Union{Symbol,Expr}, vis_kwargs::Expr...)
  func, args, sequence, kwargs = func_args_sequence_kwargs(ex, vis_kwargs...)
  e = quote
    $(vis_func!)($(esc(fig)), $(func), ($(args...),), $(sequence); $(kwargs...))
  end
  return e
end

"""
    @visualize

Visualize a contraction of ITensors, returning the result of the contraction.

The contraction should be written in terms of a series of ITensors contracted with `*`.

# Examples
```julia
using ITensors
using ITensorVisualization

i = Index(2, "index_i")
j = Index(10, "index_j")
k = Index(40, "index_k")
l = Index(40, "index_l")
m = Index(40, "index_m")
A = randomITensor(i, j, k)
B = randomITensor(i, j, l, m)
C = randomITensor(k, l)

# Contract the tensors over the common indices
# and visualize the results
ABC = @visualize A * B * C

AB = @visualize A * B
# Use readline() to pause between plots
readline()
ABC = @visualize AB * C vertex_labels = ["A*B", "C"]
readline()

# Save the results to figures for viewing later
AB = @visualize fig1 A * B
ABC = @visualize fig2 AB * C vertex_labels = ["A*B", "C"]

display(fig1)
readline()
display(fig2)
readline()
```

# Keyword arguments:
- `vertex_labels`: Custom tensor labels to display on the vertices of the digram. If not specified, they are determined automatically from the input to the macro.
- `edge_labels=IndexLabels()`: A list of the edge labels or an `AbstractEdgeLabels` object specifying how they should be made.
- `arrow_show`: Whether or not to show arrows on the edges.
"""
macro visualize(fig::Symbol, ex::Symbol, kwargs::Expr...)
  e = quote
    $(esc(fig)) = $(visualize_expr(visualize, ex, kwargs...))
    $(esc(ex))
  end
  return e
end

macro visualize!(fig, ex::Symbol, kwargs::Expr...)
  e = quote
    $(visualize_expr!(fig, visualize!, ex, kwargs...))
    $(esc(ex))
  end
  return e
end

macro visualize(ex::Symbol)
  e = quote
    display($(visualize_expr(visualize, ex)))
    $(esc(ex))
  end
  return e
end

macro visualize(ex_or_fig::Symbol, ex_or_kwarg::Expr, last_kwargs::Expr...)
  if ex_or_kwarg.head == :(=)
    # The second input is a keyword argument which means that the
    # first input is the collection to visualize (no figure output binding specified)
    ex = ex_or_fig
    kwargs = (ex_or_kwarg, last_kwargs...)
    e = quote
      display($(visualize_expr(visualize, ex, kwargs...)))
      $(esc(ex))
    end
  else
    # The second input is not a keyword argument which means that the
    # first input is the binding for the figure output, the second is the expression
    # to visualize
    fig = ex_or_fig
    ex = ex_or_kwarg
    kwargs = last_kwargs
    e = quote
      $(esc(fig)) = $(visualize_expr(visualize, ex, kwargs...))
      $(esc(ex))
    end
  end
  return e
end

macro visualize!(fig, ex::Expr, kwargs::Expr...)
  e = quote
    $(visualize_expr!(fig, visualize!, ex, kwargs...))
    $(esc(ex))
  end
  return e
end

macro visualize(ex::Expr, kwargs::Expr...)
  e = quote
    display($(visualize_expr(visualize, ex, kwargs...)))
    $(esc(ex))
  end
  return e
end

macro visualize_noeval(ex::Symbol, kwargs::Expr...)
  e = quote
    $(visualize_expr(visualize, ex, kwargs...))
  end
  return e
end

macro visualize_noeval(ex::Expr, kwargs::Expr...)
  e = quote
    $(visualize_expr(visualize, ex, kwargs...))
  end
  return e
end

macro visualize_noeval!(fig, ex::Symbol, kwargs::Expr...)
  e = quote
    $(visualize_expr!(fig, visualize!, ex, kwargs...))
  end
  return e
end

macro visualize_noeval!(fig, ex::Expr, kwargs::Expr...)
  e = quote
    $(visualize_expr!(fig, visualize!, ex, kwargs...))
  end
  return e
end

macro visualize_sequence(fig::Symbol, ex::Expr, kwargs::Expr...)
  e = quote
    $(esc(fig)) = $(visualize_expr(visualize_sequence, ex, kwargs...))
    $(esc(ex))
  end
  return e
end

macro visualize_sequence(ex::Expr, kwargs::Expr...)
  e = quote
    display($(visualize_expr(visualize_sequence, ex, kwargs...)))
    $(esc(ex))
  end
  return e
end

macro visualize_sequence_noeval(ex::Expr, kwargs::Expr...)
  e = quote
    $(visualize_expr(visualize_sequence, ex, kwargs...))
  end
  return e
end


text:
import Base: show, push!

"""
    CLIException

Command line interface exceptions, contains a member
`cmdstack` for command backtrace.
"""
abstract type CLIException <: Exception end
push!(e::CLIException, cmdname::String) = push!(e.cmdstack, cmdname)
cmdtrace(e::CLIException) = join(e.cmdstack[end:-1:1], " ")

function helpguide(io::IO, msg::CLIException)
    cmd = cmdtrace(msg)
    println(io, "Try '$(cmd) --help/-h' for more information")
end

function errorinfo(io::IO, e::Exception)
    println(io, e)
    println(io, "Try 'option --help/-h' for more information")
end

function errorinfo(io::IO, msg::CLIException)
    println(io, msg)
    helpguide(io, msg)
end

errorinfo(msg::CLIException) = errorinfo(STDERR, msg)

struct OptionError <: CLIException
    cmdstack::Vector{String}
    option::String
end

OptionError(option::String) = OptionError([], option)
OptionError(init::String, option::String) = OptionError([init], option)

function show(io::IO, msg::OptionError)
    print(io, "unrecognized option `$(msg.option)`")
end

"""
Input command error.
"""
struct CmdError <: CLIException
    cmdstack::Vector{String}
    cmdname::String
end

CmdError(cmdname::String) = CmdError([], cmdname)
CmdError(init::String, cmdname::String) = CmdError([init], cmdname)

function show(io::IO, msg::CmdError)
    print(io, "unrecognized command `$(msg.cmdname)`")
end

"""
Invalid CLI argument (type)
"""
struct CLIArgError{SIG, GOT} <: CLIException
    cmdstack::Vector{String}
end

CLIArgError(::Type{S}, ::Type{T}) where {S, T} = CLIArgError{S, T}([])
CLIArgError(init::String, ::Type{S}, ::Type{T}) where {S, T} = CLIArgError{S, T}([init])

function show(io::IO, msg::CLIArgError{S, T}) where {S, T}
    print(io, "Unknown type: $(T) (expect $(S))")
end

"""
General CLI error message
"""
struct CLIError <: CLIException
    cmdstack::Vector{String}
    msg::String
end

CLIError(msg::String) = CLIError([], msg)

function show(io::IO, msg::CLIError)
    print(msg.msg)
end


text:
using Simplify: PatternRule, EvalRule, OrderRule
using Simplify: diff
using SpecialSets

@syms f g
@syms a b c d e
@vars x y z

@testset "Rule" begin

    @testset "PatternRule" begin
        @test normalize(@term(f(a, 0)), PatternRule(@term(f(x, 0)), @term(x))) == @term(a)
        @test normalize(@term(a + 0), PatternRule(@term(x + 0), @term(x))) == @term(a)
        @test normalize(@term(b + 1), PatternRule(@term(x + 0), @term(x))) == @term(b + 1)
        @test normalize(@term(b), PatternRule(@term(x + 0), @term(x))) == @term(b)
        @test normalize(@term(f(a, b)), Rules(@term(f(x, y)) => @term(g(x)))) == @term(g(a))

        @testset "Predicates" begin
            trs = Rules(PatternRule(
                @term(x / x),
                @term(one(x)),
                [σ -> isvalid(Image(σ[x], Nonzero))]
            ))
            odd = Symbolic(:odd)

            with_context([get_context(); Image(x, TypeSet(Number)); Image(odd, Odd)]) do
                @test normalize(@term(3 / 3)        , trs) == @term(one(3))
                @test normalize(@term(2 / 3)        , trs) == @term(2 / 3)
                @test normalize(@term(x / x)        , trs) == @term(x / x)
                @test normalize(@term((2^x) / (2^x)), trs) == @term(one(2^x))
                @test normalize(@term(odd / odd)    , trs) == @term(one(odd))
                @test_skip normalize(@term((odd + 2) / (odd + 2)), trs) == @term(one(odd + 2))
            end

            @vars h
            associative_rs = @term RULES [
                (h(h(x, y), z) => h(x, y, z)) where {σ -> isvalid(Associative(σ[h]))}
                (h(x, h(y, z)) => h(x, y, z)) where {σ -> isvalid(Associative(σ[h]))}
            ]
            @test normalize(@term(a * (b * c)), associative_rs) == @term(a * b * c)
            @test_skip normalize(@term(((a + b) + c) * (d * e)), associative_rs) == @term((a + b + c) * d * e)
            with_context(Context(Associative(f))) do
                @test normalize(@term(a * (b * c))  , associative_rs) == @term(a * (b * c))
                @test normalize(@term(f(f(a, b), c)), associative_rs) == @term(f(a, b, c))
                @test_skip normalize(@term(f(f(1, f(2, 3), 4), f(5, f(6, 7)))), associative_rs) ==
                    @term(f(1, 2, 3, 4, 5, 6, 7))
            end
        end
    end
    @testset "EvalRule" begin
        @test normalize(@term(f(2, 3)), EvalRule(f, *)) == @term(6)
        @test normalize(@term(f(a, 3)), EvalRule(f, *)) == @term(f(a, 3))
        @test normalize(@term(f(2, b)), EvalRule(f, *)) == @term(f(2, b))
        @test normalize(@term("a" * "b"), EvalRule(*)) == @term("ab")
        @test normalize(@term(a + 2 * 3), EvalRule(*)) == @term(a + 2 * 3)
        @test normalize(@term(a + 2 * 3), Rules(EvalRule(*))) == @term(a + 6)
        @test normalize(@term(2 * 3 + 4 * 5), Rules(EvalRule(*))) == @term(6 + 20)
        @test normalize(@term(2 * 3 + 4 * 5), Rules(EvalRule(+), EvalRule(*))) == @term(26)

        with_context(Context(Associative(f))) do
            rule = EvalRule(f, +)
            @test normalize(@term(f(a, 1, 2, b, 3, c)), rule) == @term(f(a, 3, b, 3, c))
            @test normalize(@term(f(1, 2, 3, 4, 5)), rule) == @term(15)
            @test normalize(@term(f(1, 2, x, 3, 4, 5)), rule) == @term(f(3, x, 12))
            @test normalize(@term(f(1, 2, x, y, 3, 4, 5)), rule) == @term(f(3, x, y, 12))
        end

        with_context(Context(Associative(f), Commutative(f))) do
            rule = EvalRule(f, +)
            @test normalize(@term(f(a, 1, 2, b, 3, c)), rule) == @term(f(a, b, c, 6))
            @test normalize(@term(f(1, 2, 3, 4, 5)), rule) == @term(15)
            @test normalize(@term(f(1, 2, x, 3, 4, 5)), rule) == @term(f(x, 15))
            @test normalize(@term(f(1, 2, x, y, 3, 4, 5)), rule) == @term(f(x, y, 15))
        end
    end
    @testset "OrderRule" begin
        with_context(Context(Commutative(f))) do
            rule = OrderRule(x -> sprint(show, x))
            @test normalize(@term(f(a, b)), rule) == @term(f(a, b))
            @test normalize(@term(f(b, a)), rule) == @term(f(a, b))
            @test normalize(@term(f(2, 1)), rule) == @term(f(1, 2))
            @test normalize(@term(f(1, 1.0)), rule) == @term(f(1, 1.0))
            @test normalize(@term(f(1.0, 1)), rule) == @term(f(1, 1.0))
        end
    end
end

@testset "accuracy of standard rules" begin
    @testset "$set" for set ∈ [:BASIC, :ABSOLUTE_VALUE, :BOOLEAN, :CALCULUS, :LOGARITHM, :TRIGONOMETRY, :TYPES]
        CASES = Any[-10:10; -5:0.1:5; false:true]

        set ∈ [:CALCULUS, :LOGARITHM] && continue

        for rule ∈ rules(set)
            rule isa PatternRule || continue
            l, r = rule.left, rule.right

            @testset "$rule" begin
                for case ∈ CASES
                    vars = Simplify.vars(l)

                    σ = Dict(var => case for var ∈ vars)
                    Simplify._preds_match(rule.ps, σ) || continue

                    lres = replace(l, σ) |> get |> eval
                    rres = replace(r, σ) |> get |> eval

                    success = isapprox(lres, rres, atol = 1e-9)
                    success || @error "Case" case
                    @test success
                end
            end
        end
    end
end

@testset "normalize" begin

    @testset "BASIC" begin
        @syms a b

        with_context([get_context(); Image.([a, b], Ref(TypeSet(Number)))]) do
            @test normalize(@term(37)) == @term(37)
            @test normalize(@term(a)) == @term(a)
            @test normalize(@term(a + 0)) == @term(a)
            @test normalize(@term(b + 0 + 0)) == @term(b)
            @test normalize(@term(0 + 0 + b)) == @term(b)
            @test normalize(@term(b * (1 + 2 - 3))) == @term(0)
            @test normalize(@term(0 + b + 0)) == @term(b)
            @test normalize(@term(a * (2 * b))) == @term(a * 2 * b)
            @test normalize(@term(a * (b * c))) == @term(a * b * c)
            @test normalize(@term((a * b) * c)) == @term(a * b * c)
            @test_broken normalize(@term(a * ((b * c) * d))) == @term(a * b * c * d)
            @test_broken normalize(@term(a * ((b * c) * d) * e)) == @term(a * b * c * d * e)
            @test normalize(@term(a + ((2 + b) + 3))) == @term(a + b + 5)
            @test normalize(@term(a + 2)) == normalize(@term(2 + a))
        end
    end

    @testset "ABSOLUTE_VALUE" begin
        @syms a b

        with_context([get_context(); Image(a, TypeSet(Int)); Image(b, Nonzero)]) do
            @test normalize(@term(abs(a))) == @term(abs(a))
            @test normalize(@term(abs(-a))) == @term(abs(a))
            @test normalize(@term(abs(0))) == @term(0)
            @test normalize(@term(abs(-3))) == @term(3)
            @test normalize(@term(abs(3))) == @term(3)
            @test normalize(@term(abs(2a))) == @term(2abs(a))
            @test normalize(@term(abs(-(5a)))) == @term(5abs(a))
            @test normalize(@term(abs(a * b))) == @term(abs(a) * abs(b))
            @test normalize(@term(abs(a / b))) == @term(abs(a) * inv(abs(b)))
            @test normalize(@term(abs(a / 1))) == @term(abs(a))
            @test normalize(@term(abs(abs(a)))) == @term(abs(a))
            @test normalize(@term(abs(a^2))) == @term(a^2)
        end

        @syms d1 d2
        with_context([get_context(); Image(d1, Set([1, 2])); Image(d2, Set([-1, 1]))]) do
            @test normalize(@term(abs(d1))) == @term(d1)
            @test normalize(@term(abs(d2))) == @term(abs(d2))
        end
    end

    @testset "BOOLEAN" begin
        @syms x y

        with_context([get_context(); Image.([x, y], Ref(TypeSet(Bool)))]) do
            @test normalize(@term(x & true)) == @term(x)
            @test normalize(@term(x ⊻ false)) == @term(x)
            @test normalize(@term(x | (x & y))) == @term(x)
            @test normalize(@term(y | !y)) == @term(true)
            @test normalize(@term(!y | y)) == @term(true)
            @test normalize(@term(y & y)) == @term(y)
            @test normalize(@term(y ⊻ y)) == @term(false)
            @test normalize(@term(!(!x))) == @term(x)
            @test normalize(@term(!(!x & !x))) == @term(x)
            @test normalize(@term(!(!x & !x) & !x)) == @term(false)
            @test normalize(@term(!(!x & !x) | x)) == @term(x)
            @test normalize(@term(!x & x | (y & (y | true)))) == @term(y)
        end
    end

    @testset "CALCULUS" begin
        @syms x y z

        ctx = [
            Image.([x, y, z], Ref(TypeSet(Int)))
            Signature(f, [TypeSet(Number)], TypeSet(Number))
        ]

        with_context([get_context(); ctx]) do
            @test normalize(@term diff(2, x)) == @term(0)
            @test normalize(@term diff(x * y, x)) == @term(x*diff(y, x) + y)
            @test normalize(@term diff(sin(2x + 3y), x)) == @term(cos(2x + 3y) * (3diff(y, x) + 2))
            @test normalize(@term diff(x * y + sin(x^z), x)) ==
                  normalize(@term(y + x*diff(y, x) + cos(x^z)*(z*x^(z-1) + (x^z * log(x) * diff(z, x)))))
            @test normalize(@term diff(2x + tan(x), x)) == @term(tan(x)^2 + 3)
            @test normalize(@term diff(f(x) + 3x, x)) == @term(diff(f(x), x) + 3)
        end

        @syms w
        with_context([get_context(); Image(w, Nonzero ∩ TypeSet(Float64))]) do
            @test normalize(@term diff(log(w), w)) == @term(inv(w))
        end
    end

    @testset "LOGARITHM" begin
        @syms b x y n

        ctx = Context(
            Associative(+),
            Commutative(+),
            Associative(*),
            Commutative(*),
            Image(b, TypeSet(Real)),
            Image(x, TypeSet(Real)),
            Image(y, TypeSet(Real)),
            Image(n, GreaterThan(3)),
        )

        with_context([get_context(); ctx]) do
            @test normalize(@term(log(b, x * y))) == @term(log(b, x) + log(b, y))
            @test normalize(@term(log(n, 1))) == @term(0)
            @test normalize(@term(log(n, n ^ x))) == @term(x)
            @test normalize(@term(log(n, n * y))) == @term(log(n, y) + 1)
            @test normalize(@term(log(n, y * n))) == @term(log(n, y) + 1)
            @test normalize(@term(log(b, 1/n))) == @term(-log(b, n))
            @test_skip normalize(@term(b ^ log(b, x*b))) == @term(x * b)
            @test_skip normalize(@term(b ^ (log(b, x) + log(b, b)))) == @term(x * b)
        end
    end

    @testset "TRIGONOMETRY" begin
        @syms α β θ

        @test normalize(@term(sin(0) * tan(π / 4))) == @term(0)
        @test normalize(@term(sin(π/3)cos(0) + cos(π/3)sin(0))) == @term(√3 * inv(2))
        @test normalize(@term(one(θ) + tan(θ) ^ 2)) == @term(sec(θ) ^ 2)
        @test normalize(@term(tan(π / 6))) == @term(√3 * inv(3))
        @test_broken normalize(@term(1 / (sin(-3) / cos(-3)))) == @term(-cot(3))
        @test normalize(@term(2 * cos((α + β) / 2) * cos((α - β) / 2))) == @term(cos(α) + cos(β))
        @test normalize(@term((tan(α) - tan(β)) * inv(1 + tan(α) * tan(β)))) == @term(tan(-β + α))
        @test normalize(@term(csc(π/2 - θ))) == @term(sec(θ))

        with_context(Context(Image(a, TypeSet(Int)))) do
            @test_broken normalize(@term sin(a)^2 + cos(a)^2 + 1) == 2
        end
    end

    @testset "custom" begin
        @syms a
        @vars x y

        @test normalize(@term(f(a, a)), @term RULES [
            f(x, x) => x
        ]) == @term(a)

        @test normalize(@term(f(a, b)), @term RULES [
            f(x, x) => x
        ]) == @term(f(a, b))

        @test normalize(@term(f(f(a), a)), @term RULES [
            f(x, x) => x
            f(x)     => x
        ]) == @term(a)

        @test normalize(@term(f(f(a), g(b))), @term RULES [
            f(f(x), y) => y
            g(x) => x
        ]) == @term(b)

        @test normalize(@term(f(a, a)), @term RULES [
            f(b, b) => b
        ]) == @term(f(a, a))
    end
end


text:
#
# Implementation of scope operations for the toy backend.
#

compile(binding::AbstractBinding, base::Scope, args::Vector{AbstractSyntax}) =
    compile(binding, base, args...)

immutable SimpleBinding <: AbstractBinding
    query::Query
end

compile(binding::SimpleBinding, base::Scope) = binding.query

immutable ParamBinding <: AbstractBinding
    tag::Symbol
    output::Output
end

compile(binding::ParamBinding, base::Scope) =
    let scope = nest(base, domain(binding.output)),
        pipe = ParamPipe(domain(base), binding.tag, binding.output)
        Query(scope, pipe)
    end

immutable ClassBinding <: AbstractBinding
    db::ToyDatabase
    classname::Symbol
end

compile(binding::ClassBinding, base::Scope) =
    let scope = nest(base, Entity{binding.classname}),
        class = binding.db.schema.name2class[binding.classname],
        set = binding.db.instance.sets[binding.classname],
        pipe = SetPipe(binding.classname, set, ismonic=true, iscovering=true)
        Query(scope, pipe)
    end

immutable ArrowBinding <: AbstractBinding
    db::ToyDatabase
    classname::Symbol
    arrowname::Symbol
end

compile(binding::ArrowBinding, base::Scope) =
    let class = binding.db.schema.name2class[binding.classname],
        arrow = class.name2arrow[binding.arrowname],
        scope = nest(base, odomain(arrow)),
        scope = settag(scope, binding.arrowname),
        tag = Symbol(binding.classname, "/", binding.arrowname),
        map = binding.db.instance.maps[binding.classname, binding.arrowname],
        pipe = EntityMapPipe(tag, Entity{binding.classname}, map, omode(arrow))
        if arrow.select != nothing
            out = mkselect(scope, arrow.select)
            out = Query(out, scope=settag(out.scope, binding.arrowname))
            scope = addlocal(scope, :__out, out)
        end
        Query(scope, pipe)
    end

immutable EntityIdBinding <: AbstractBinding
    classname::Symbol
end

compile(binding::EntityIdBinding, base::Scope) =
    let scope = settag(nest(base, Int), :id),
        pipe = ItemPipe(Entity{binding.classname}, :id, runique=true)
        Query(scope, pipe)
    end

immutable GlobalBinding <: AbstractBinding
    fn::Symbol
end

compile(binding::GlobalBinding, base::Scope, args::Vector{AbstractSyntax}) =
    compile(Fn{binding.fn}, base, args...)

scope(db::ToyDatabase) =
    begin
        globals = SymbolTable()
        for m in methods(RBT.compile)
            if length(m.sig.parameters) < 1
                continue
            end
            p = m.sig.parameters[2]
            ps =
                isa(p, Union) ? collect(Type, p.types) :
                p <: Type ? Type[p.parameters[1]] : Type[]
            for p in ps
                if p <: Type
                    p = p.parameters[1]
                end
                if p <: Fn
                    name = fnname(p)
                    if isa(name, Symbol)
                        globals = assoc(globals, (name, -1), GlobalBinding(name))
                    end
                end
            end
        end
        locals = SymbolTable()
        for c in db.schema.classes
            locals = assoc(locals, (c.name, 0), ClassBinding(db, c.name))
        end
        return Scope(db, Unit, globals=globals, locals=locals)
    end

nest(db::ToyDatabase, base::Scope, domain::Type) =
    begin
        tag = NullableSymbol()
        locals = SymbolTable()
        if isentity(domain)
            class = db.schema.name2class[classname(domain)]
            tag = NullableSymbol(class.name)
            locals = assoc(locals, (:id, 0), EntityIdBinding(class.name))
            for arrow in class.arrows
                locals = assoc(locals, (arrow.name, 0), ArrowBinding(db, class.name, arrow.name))
            end
        elseif domain == Unit
            for c in db.schema.classes
                locals = assoc(locals, (c.name, 0), ClassBinding(db, c.name))
            end
        end
        scope = Scope(db, domain, globals=base.globals, locals=locals, tag=tag)
        if isentity(domain)
            scope = addlocal(scope, (:__id, 0), locals[:id, 0])
            if class.select != nothing
                out = mkselect(scope, class.select)
                out = Query(out, scope=settag(out.scope, classname(domain)))
                scope = addlocal(scope, (:__out, 0), out)
            end
            if isa(class.select, Tuple)
                items = tuple([SimpleBinding(mkselect(scope, name, false)) for name in class.select]...)
                scope = setitems(scope, items)
            end
        end
        return scope
    end

mkselect(base::Scope, name::Symbol, finish::Bool=true) =
    let binding = get(lookup(base, (name, 0)))
        finish ? select(compile(binding, base)) : compile(binding, base)
    end

mkselect(base::Scope, names::Tuple, finish::Bool=true) =
    let fields = Query[mkselect(base, name, finish) for name in names]
        pipe = TuplePipe([field.pipe for field in fields])
        scope = nest(base, odomain(pipe))
        scope = setitems(
            scope,
            tuple([SimpleBinding(Query(
                fields[k],
                scope=settag(scope, names[k]),
                pipe=ItemPipe(pipe, k))) for k = eachindex(names)]...))
        Query(scope, pipe)
    end

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*-----------------------------------------------------------------------
#* Author: Tiago Os򱨯rio
#* $Id$
#*---------------------------------------------------------------------
type Lag_1
	Lag_1()=begin
		[
			:(diff(output)=0/"s"),
		],
		[
			"",
		],
		new(
			positive(Dict{Symbol,Any}(
				:Brief=>"model gain"
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"model time constant",
				:Unit=>"s"
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"input signal"
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"output signal"
			)),
			[
				:(gain*input= tau*diff(output) + output),
			],
			[
				"Calculate output",
			],
			[:gain,:tau,],
			[:input,:output,]
		)
	end
	gain::positive
	tau::positive
	input::DanaReal
	output::DanaReal
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	initials::Array{Expr,1}
	initialNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Lag_1
function setEquationFlow(in::Lag_1)
	addEquation(1)
end
function initial(in::Lag_1)
	addEquation(1)
end
function atributes(in::Lag_1,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Icon]="icon/Lag_1"
	fields[:Brief]="Model Lag."
	fields[:Info]="== Inputs ==
* One input signal.

== Outputs ==
* One output signal.
"
	drive!(fields,_)
	return fields
end
Lag_1(_::Dict{Symbol,Any})=begin
	newModel=Lag_1()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
int main() {
    p();
    printInt(1);
    return 0;
}

void p() {}

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
using LinearAlgebra

import Libdl
const depsfile = joinpath(@__DIR__, "deps.jl")

settings = joinpath(first(DEPOT_PATH), "prefs", "FFTW")
mkpath(dirname(settings))
if haskey(ENV, "JULIA_FFTW_PROVIDER")
    provider = ENV["JULIA_FFTW_PROVIDER"]
    open(f -> println(f, provider), settings, "w")
elseif isfile(settings)
    provider = readchomp(settings)
else
    provider = "FFTW"
    open(f -> println(f, provider), settings, "w")
end

if provider == "MKL"
    const mkllib = Sys.iswindows() ? "mkl_rt" : "libmkl_rt"
    # If BLAS was compiled with MKL and the user wants MKL-based FFTs, we'll oblige.
    open(depsfile, "w") do io
        println(io, """
            using IntelOpenMP_jll, MKL_jll
            check_deps() = nothing
            const libfftw3 = MKL_jll.libmkl_rt_path
            const libfftw3f = libfftw3
        """)
    end
elseif provider == "FFTW"
    open(depsfile, "w") do io
        println(io, """
            using FFTW_jll
            check_deps() = nothing
        """)
    end
else
    error("Unrecognized JULIA_FFTW_PROVIDER \"$provider\".\n",
          "To fix this, set ENV[\"JULIA_FFTW_PROVIDER\"] to \"FFTW\" or \"MKL\"\n",
          "and rerun Pkg.build(\"FFTW\").")
end


text:
@testset "W2V network builder" begin
    @test_broken 1 == 2
end


text:
primitive type MyInt 32 end

struct Person
    id::Int
    name::String
end

@testset "ArrowTypes" begin

@test ArrowTypes.ArrowKind(MyInt) == ArrowTypes.PrimitiveKind()
@test ArrowTypes.ArrowKind(Person) == ArrowTypes.StructKind()
@test ArrowTypes.ArrowKind(Person(0, "bob")) == ArrowTypes.StructKind()

@test ArrowTypes.ArrowType(Int) == Int
@test ArrowTypes.ArrowType(Union{Int, Missing}) == Union{Int, Missing}
@test ArrowTypes.ArrowType(Missing) == Missing

@test ArrowTypes.toarrow(1) === 1

@test ArrowTypes.arrowname(Int) == Symbol()
@test !ArrowTypes.hasarrowname(Int)

@test ArrowTypes.arrowmetadata(Int) == ""

@test ArrowTypes.JuliaType(1) === nothing
@test ArrowTypes.JuliaType(1, Int) === nothing
@test ArrowTypes.JuliaType(1, Int, nothing) === nothing

@test ArrowTypes.fromarrow(Int, 1) === 1
@test ArrowTypes.fromarrow(Person, 1, "bob") == Person(1, "bob")
@test ArrowTypes.fromarrow(Union{Int, Missing}, missing) === missing
@test ArrowTypes.fromarrow(Union{Int, Missing}, 1) === 1
@test ArrowTypes.fromarrow(Union{Float64, Missing}, 1) === 1.0

@test ArrowTypes.ArrowKind(Missing) == ArrowTypes.NullKind()
@test ArrowTypes.ArrowKind(Nothing) == ArrowTypes.NullKind()
@test ArrowTypes.ArrowType(Nothing) == Missing
@test ArrowTypes.toarrow(nothing) === missing
@test ArrowTypes.arrowname(Nothing) == ArrowTypes.NOTHING
@test ArrowTypes.JuliaType(Val(ArrowTypes.NOTHING)) == Nothing
@test ArrowTypes.fromarrow(Nothing, missing) === nothing

@test ArrowTypes.ArrowKind(Int) == ArrowTypes.PrimitiveKind()
@test ArrowTypes.ArrowKind(Float64) == ArrowTypes.PrimitiveKind()

@test ArrowTypes.ArrowType(Char) == UInt32
@test ArrowTypes.toarrow('1') == UInt32('1')
@test ArrowTypes.arrowname(Char) == ArrowTypes.CHAR
@test ArrowTypes.JuliaType(Val(ArrowTypes.CHAR)) == Char
@test ArrowTypes.fromarrow(Char, UInt32('1')) == '1'

@test ArrowTypes.ArrowKind(Bool) == ArrowTypes.BoolKind()

@test ArrowTypes.ListKind() == ArrowTypes.ListKind{false}()
@test !ArrowTypes.isstringtype(ArrowTypes.ListKind())
@test !ArrowTypes.isstringtype(typeof(ArrowTypes.ListKind()))
@test ArrowTypes.ArrowKind(String) == ArrowTypes.ListKind{true}()

hey = collect(b"hey")
@test ArrowTypes.fromarrow(String, pointer(hey), 3) == "hey"

@test ArrowTypes.ArrowType(Symbol) == String
@test ArrowTypes.toarrow(:hey) == "hey"
@test ArrowTypes.arrowname(Symbol) == ArrowTypes.SYMBOL
@test ArrowTypes.JuliaType(Val(ArrowTypes.SYMBOL)) == Symbol
@test ArrowTypes.fromarrow(Symbol, pointer(hey), 3) == :hey

@test ArrowTypes.ArrowKind(Vector{Int}) == ArrowTypes.ListKind()
@test ArrowTypes.ArrowKind(Set{Int}) == ArrowTypes.ListKind()
@test ArrowTypes.ArrowType(Set{Int}) == Vector{Int}
@test typeof(ArrowTypes.toarrow(Set([1,2,3]))) <: Vector{Int}
@test ArrowTypes.arrowname(Set{Int}) == ArrowTypes.SET
@test ArrowTypes.JuliaType(Val(ArrowTypes.SET), Vector{Int}) == Set{Int}
@test ArrowTypes.fromarrow(Set{Int}, [1,2,3]) == Set([1,2,3])

K = ArrowTypes.ArrowKind(NTuple{3, UInt8})
@test ArrowTypes.gettype(K) == UInt8
@test ArrowTypes.getsize(K) == 3
@test K == ArrowTypes.FixedSizeListKind{3, UInt8}()

u = UUID(rand(UInt128))
ubytes = ArrowTypes._cast(NTuple{16, UInt8}, u.value)
@test ArrowTypes.ArrowKind(u) == ArrowTypes.FixedSizeListKind{16, UInt8}()
@test ArrowTypes.ArrowType(UUID) == NTuple{16, UInt8}
@test ArrowTypes.toarrow(u) == ubytes
@test ArrowTypes.arrowname(UUID) == ArrowTypes.UUIDSYMBOL
@test ArrowTypes.JuliaType(Val(ArrowTypes.UUIDSYMBOL)) == UUID
@test ArrowTypes.fromarrow(UUID, ubytes) == u

nt = (id=1, name="bob")
@test ArrowTypes.ArrowKind(NamedTuple) == ArrowTypes.StructKind()
@test ArrowTypes.fromarrow(typeof(nt), nt) === nt
@test ArrowTypes.fromarrow(Person, nt) == Person(1, "bob")
@test ArrowTypes.ArrowKind(Tuple) == ArrowTypes.StructKind()
@test ArrowTypes.ArrowKind(Tuple{}) == ArrowTypes.StructKind()
@test ArrowTypes.arrowname(Tuple{Int, String}) == ArrowTypes.TUPLE
@test ArrowTypes.arrowname(Tuple{}) == ArrowTypes.TUPLE
@test ArrowTypes.JuliaType(Val(ArrowTypes.TUPLE), NamedTuple{(Symbol("1"), Symbol("2")), Tuple{Int, String}}) == Tuple{Int, String}
@test ArrowTypes.fromarrow(Tuple{Int, String}, nt) == (1, "bob")
@test ArrowTypes.fromarrow(Union{Missing, typeof(nt)}, nt) == nt

v = v"1"
v_nt = (major=1, minor=0, patch=0, prerelease=(), build=())
@test ArrowTypes.ArrowKind(VersionNumber) == ArrowTypes.StructKind()
@test ArrowTypes.arrowname(VersionNumber) == ArrowTypes.VERSION_NUMBER
@test ArrowTypes.JuliaType(Val(ArrowTypes.VERSION_NUMBER)) == VersionNumber
@test ArrowTypes.fromarrow(typeof(v), v_nt) == v

@test ArrowTypes.ArrowKind(Dict{String, Int}) == ArrowTypes.MapKind()
@test ArrowTypes.ArrowKind(Union{String, Int}) == ArrowTypes.UnionKind()

@test ArrowTypes.default(Int) == Int(0)
@test ArrowTypes.default(Symbol) == Symbol()
@test ArrowTypes.default(Char) == '\0'
@test ArrowTypes.default(String) == ""
@test ArrowTypes.default(Union{Int, Missing}) == Int(0)

@test ArrowTypes.promoteunion(Int, Float64) == Float64
@test ArrowTypes.promoteunion(Int, String) == Union{Int, String}

@test ArrowTypes.concrete_or_concreteunion(Int)
@test !ArrowTypes.concrete_or_concreteunion(Union{Real, String})
@test !ArrowTypes.concrete_or_concreteunion(Any)

@test ArrowTypes.ToArrow([1,2,3]) == [1,2,3]
@test ArrowTypes.ToArrow([:hey, :ho]) == ["hey", "ho"]
@test ArrowTypes.ToArrow(Any[1, 3.14]) == [1.0, 3.14]
@test ArrowTypes.ToArrow(Any[1, 3.14, "hey"]) == [1.0, 3.14, "hey"]

end


text:
using OnlineAI, FactCheck

n = 10_000
x = rand(n, 10)
y = rand(n, 1)

nin = ncols(x)
nout = ncols(y)

updater = AdadeltaModel(1e-6, 0.95, 1e-6)
dropout = Dropout()
mloss = L2CostModel()
params = NetParams(updater, dropout, mloss)


text:
using Test

@testset "Tests" begin
    tests = [
        "moment_matching.jl",
    ]

    for t in tests
        include(t)
    end
end


text:
using IRTools: varargs!, inlineable!, pis!, slots!
using IRTools.Inner: argnames!, update!

ignore_sig(T) = all(T -> T <: Type, T.parameters)

function edge!(m::IRTools.Meta, edge::Core.MethodInstance)
  m.code.edges === nothing && (m.code.edges = Core.MethodInstance[])
  push!(m.code.edges, edge)
  return
end

@generated function _pullback(ctx::AContext, f, args...)
  T = Tuple{f,args...}
  ignore_sig(T) && return :(f(args...), Pullback{$T}(()))

  iskw = is_kwfunc(f, args...)
  # if it is_kw then `args[1]` are the keyword args, `args[2]` is actual function
  base_T = iskw ? Tuple{args[2:end]...} : T
  hascr, cr_edge = has_chain_rrule(base_T)
  chain_rrule_f = iskw ? :chain_rrule_kw : :chain_rrule
  hascr && return :($chain_rrule_f(f, args...))

  g = try _lookup_grad(T) catch e e end
  !(g isa Tuple) && return :(f(args...), Pullback{$T}((f,)))
  meta, forw, _ = g
  argnames!(meta, Symbol("#self#"), :ctx, :f, :args)
  forw = varargs!(meta, forw, 3)
  # IRTools.verify(forw)
  forw = slots!(pis!(inlineable!(forw)))
  @static if VERSION >= v"1.3" # no edges pre-1.3
    # be ready to swap to using chainrule if one is declared
    cr_edge != nothing && edge!(meta, cr_edge)
  end
  return update!(meta.code, forw)
end

@generated function (j::Pullback{T})(Δ) where T
  ignore_sig(T) && return :nothing
  g = try _lookup_grad(T)
  catch e
    rethrow(CompileError(T,e))
  end
  if g === nothing
    Δ == Nothing && return :nothing
    return :(error("Non-differentiable function $(repr(j.t[1]))"))
  end
  meta, _, back = g
  argnames!(meta, Symbol("#self#"), :Δ)
  # IRTools.verify(back)
  back = slots!(inlineable!(back))
  return update!(meta.code, back)
end


text:
## Multi-element airfoil using constant-strength doublet-source panel method
using AeroMDAO
using CoordinateTransformations, Rotations

## Geometry variables
alpha_u   = [0.2, 0.3, 0.2, 0.15, 0.2]
alpha_l   = [-0.2, -0.1, -0.1, -0.001]
dzs       = (0., 0.)

## Airfoils
airfoil_1 = kulfan_CST(alpha_u, alpha_l, dzs, (0., 0.), 60);      # Kulfan CST airfoil
airfoil_2 = affine(scale(naca4((0,0,1,2)), 0.2); vector = [1.0, -0.1], angle = deg2rad(5.))  # NACA 4-digit airfoil

## Case setup
panels  = ComponentVector(
                            kulfan = make_panels(airfoil_1, 60),
                            naca   = make_panels(airfoil_2, 60)
                          )
uniform = Uniform2D(1., 0.)

## Obviously wrong because Kutta condition is hand-coded -_-
system  = @time solve_system(
                     panels, uniform, 15, 1e5
                    );

##
@time ues = surface_velocities(system);
@time cl  = lift_coefficient(system)
@time cls, cms, cps = surface_coefficients(system)

##
using CairoMakie

af_1 = coordinates(airfoil_1)
af_2 = coordinates(airfoil_2)

f = Figure()
Axis(f[1,1])
lines!(af_1[:,1], af_1[:,2])
lines!(af_2[:,1], af_2[:,2])
f


text:
module Test16

import Automa
import Automa.RegExp: @re_str
using Test

@testset "Test16" begin
    @test_throws ArgumentError Automa.CodeGenContext(generator=:table,  loopunroll=1)
    @test_throws ArgumentError Automa.CodeGenContext(generator=:inline, loopunroll=1)
    @test_throws ArgumentError Automa.CodeGenContext(generator=:goto,   loopunroll=-1)

    re = re"A+(B+C)*(D|E)+"
    machine = Automa.compile(re)
    for checkbounds in (true, false), loopunroll in 0:4
        ctx = Automa.CodeGenContext(generator=:goto, checkbounds=checkbounds, loopunroll=loopunroll)
        init_code = Automa.generate_init_code(ctx, machine)
        exec_code = Automa.generate_exec_code(ctx, machine)
        validate = @eval function (data)
            $(init_code)
            p_end = p_eof = lastindex(data)
            $(exec_code)
            return cs == 0
        end
        @test validate(b"ABCD")
        @test validate(b"AABCD")
        @test validate(b"AAABBCD")
        @test validate(b"AAAABCD")
        @test validate(b"AAAABBBBBCD")
        @test validate(b"AAAAAAAAAAAAAABBBBBCBCBBBBBBBBCDE")
        @test validate(b"AAAAAAAAAAAAAABBBBBCBCBBBBBBBBCDEDDDDDDEDEDDDEEEEED")
    end
end

end


text:
#=
    The File contains parser tests
=#

#------------------------------------------------------------------------------

@testset "parser - JSON to DOK_Sparsik" begin

    set_fn = "matrices_1.json"
    matrices = [from_COO(arr..., Nemo.QQ) for arr in load_COO_set(set_fn)]

    @test length(matrices) == 1
    @test matrices[1] == from_dense([2 0; -1 0;], Nemo.QQ)

    # and the second..
    set_fn = "matrices_2.json"
    matrices = [from_COO(arr..., Nemo.QQ) for arr in load_COO_set(set_fn)]

    @test length(matrices) == 3
    @test matrices[1] == from_dense([3 2 0; 0 1 0; 0 0 0;], Nemo.QQ)
    @test matrices[2] == from_dense([1 0 0;], Nemo.QQ)
    @test iszero(matrices[3])

end

@testset "parser - MTX to DOK_Sparsik" begin

    set_group = "test"
    set_name  = "test_set.mtx"

    fn, dim, sz, matrix = load_MTX(set_group, set_name)

    # hmmm
    matrix = matrix[2:end]

    @test fn == set_group*'/'*set_name && dim == 3 && sz == 1
    @test from_COO(dim, dim, matrix, Nemo.QQ) == from_dense([1 3 0; 3//2 1 0; 0 0 0;], Nemo.QQ)

end

#=
@testset "parser - *.ode to polynomials" begin

    # path to CLUE parser
    true_parser_path = "python_parser/true_parser.py"
    # path to testsets *.ode files
    testsets_path = "python_parser/testsets"

    true_parser_path, testsets_path = map(
        name -> replace("$(normpath(joinpath(@__FILE__, "..", "..")))test/$name",
                        "\\" => "/"),
        (true_parser_path, testsets_path)
    )

    # for test_filename in readdir(testsets_path)
    # for test_filename in ["e2.ode", "simple.ode", "test1.ode", "test2.ode"]
    for test_filename in ["where_is_the_bug.ode"]
        test_path = testsets_path * "/" * test_filename
        true_parser_cmd = `python $true_parser_path $test_path`

        # Now that we have output polys from CLUE parser..
        true_output = @capture_out run(true_parser_cmd)
        # ..we beautify them a bit and then compare with ours

        true_output = collect(
            filter(!isempty, map(
                        strip ∘ String,
                        split(replace(true_output, "+ -" => "- "), r"\r?\n")
            )
        ))

        test_output = collect(sort(map(
                        strip ∘ repr,
                        values(load_ODEs(test_path))
        )))

        # Browse all equations and find similar ones in true_output <--> test_output
        # with greedy matching
        # If there is a bijection established, we consider outputs equal
        test_to_true = Dict{Int, Int}()
        for (i, true_eq) in enumerate(true_output)
            for (j, test_eq) in enumerate(test_output)
                ### Controversial, but ok for now
                if Set(true_eq) == Set(test_eq) && !haskey(test_to_true, j)
                    test_to_true[j] = i
                    break
                end
            end
        end

        @test length(test_to_true) == length(test_output) == length(true_output)

    end

end
=#

@info "OK"


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*--------------------------------------------------------------------
#* Author: Gerson Balbueno Bicca
#* $Id: HeatExchangerDetailed.mso 197 2007-03-08 14:31:57Z bicca $
#*------------------------------------------------------------------
type ShellandTubes_NTU
	ShellandTubes_NTU()=begin
		new(
			ShellandTubesBasic(),
			NTU_Basic(Dict{Symbol,Any}(
				:Brief=>"NTU Method"
			)),
			[
				:(Method.NTU*Method.Cmin = _base_1.Details.Ud*_base_1.Pi*_base_1.Tubes.TubeOD*_base_1.Tubes.NumberOfTubes*_base_1.Tubes.TubeLength),
				:(Method.Cmin = min([Method.Ch,Method.Cc])),
				:(Method.Cmax = max([Method.Ch,Method.Cc])),
				:(Method.Cr = Method.Cmin/Method.Cmax),
				:(_base_1.Details.Q = Method.Eft*Method.Cmin*(_base_1.InletShell.T-_base_1.InletTube.T)),
				:(Method.Ch = _base_1.InletShell.F*_base_1.Shell.Properties.Average.Cp),
				:(Method.Cc = _base_1.InletTube.F*_base_1.Tubes.Properties.Average.Cp),
				:(_base_1.Details.Q = Method.Eft*Method.Cmin*(_base_1.InletTube.T-_base_1.InletShell.T)),
				:(Method.Cc = _base_1.InletShell.F*_base_1.Shell.Properties.Average.Cp),
				:(Method.Ch = _base_1.InletTube.F*_base_1.Tubes.Properties.Average.Cp),
				:(Method.Eft1 = 2*(1+Method.Cr+sqrt(1+Method.Cr^2)*((1+exp(-Method.NTU*sqrt(1+Method.Cr^2)))/(1-exp(-Method.NTU*sqrt(1+Method.Cr^2)))) )^(-1)),
				:(Method.Eft = ( ((1-Method.Eft1*Method.Cr)/(1-Method.Eft1))^2 -1 )*( ((1-Method.Eft1*Method.Cr)/(1-Method.Eft1))^2 - Method.Cr )^(-1)),
				:(Method.Eft = 2*(1+Method.Cr+sqrt(1+Method.Cr^2)*((1+exp(-Method.NTU*sqrt(1+Method.Cr^2)))/(1-exp(-Method.NTU*sqrt(1+Method.Cr^2)))) )^(-1)),
				:(Method.Eft1 = 1),
			],
			[
				"Number of Units Transference","Minimum Heat Capacity","Maximum Heat Capacity","Thermal Capacity Ratio","Duty","Hot Stream Heat Capacity","Cold Stream Heat Capacity","Duty","Hot Stream Heat Capacity","Cold Stream Heat Capacity","Effectiveness Correction for 2 pass shell side","TEMA F Shell Effectiveness","TEMA E Shell Effectiveness","Variable not in use when 1 Pass Shell Side",
			],
			[:Method,]
		)
	end
	_base_1::ShellandTubesBasic
	Method::NTU_Basic
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export ShellandTubes_NTU
function setEquationFlow(in::ShellandTubes_NTU)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	let switch=HotSide
		if HotSide==InletTube.T > InletShell.T
			set(switch,"tubes")
		end
		if HotSide==InletTube.T < InletShell.T
			set(switch,"shell")
		end
		if switch=="shell"
			addEquation(5)
			addEquation(6)
			addEquation(7)
		elseif switch=="tubes"
			addEquation(8)
			addEquation(9)
			addEquation(10)
		end
	end
	let switch=ShellType
		if switch=="Fshell"
			addEquation(11)
			addEquation(12)
		elseif switch=="Eshell"
			addEquation(13)
			#	Method.Eft = 1;
			addEquation(14)
		end
	end
end
function atributes(in::ShellandTubes_NTU,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/STHE"
	fields[:Brief]="Shell and Tubes Heat Exchangers"
	fields[:Info]="to be documented"
	drive!(fields,_)
	return fields
end
ShellandTubes_NTU(_::Dict{Symbol,Any})=begin
	newModel=ShellandTubes_NTU()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
#
# Plot the results of a SIRD simulation
#

using Plots

function plotsird( input :: SIRDInput, traj :: SIRDTraj, filename :: String;
                   size=(400,400) )

  ENV["GKSwstype"]="nul"

  time = traj.time
  U = traj.U
  S = traj.S
  D = traj.D
  I = traj.I
  nsteps = traj.nsteps

  xmin = time[1]
  xmax = time[nsteps]

  plot(size=(400,400),framestyle=:box)
  plot!(time,U,label="",linewidth=2,color=:blue)
  plot!(time,S,label="",linewidth=2,color=:red)
  plot!(time,I,label="",linewidth=2,color=:darkgreen)
  plot!(time,D,label="",linewidth=2,color=:black)
  plot!(xlabel="time")
  plot!(ylabel="Fraction of the population")
  plot!(xlim=[0,xmax],ylim=[0,1])

  plot!(rectangle(-0.01*xmax,0.4*xmax,0.76,1.02), opacity=0.9,label="",color=:white)

  fontsize=8
  x = 0.01*(xmax-xmin) + xmin
  d = 0.05 ; y = [ 0.99 + d ]
  annotate!(x,yd!(y,d),text("Final data:",:left,fontsize,:serif,:black))
  annotate!(x,yd!(y,d),text("Susceptible: $(@sprintf("%2.1f",U[nsteps]*100))%",:left,fontsize,:serif,:blue))
  annotate!(x,yd!(y,d),text("Sick: $(@sprintf("%2.1f",S[nsteps]*100))%",:left,fontsize,:serif,:red))
  annotate!(x,yd!(y,d),text("Immune: $(@sprintf("%2.1f",I[nsteps]*100))%",:left,fontsize,:serif,:darkgreen))
  annotate!(x,yd!(y,d),text("Dead: $(@sprintf("%2.1f",D[nsteps]*100))%",:left,fontsize,:serif,:black))

  plot!(rectangle(0.44*xmax, 1.01*xmax, 0.915, 1.02), opacity=0.9,label="",color=:white)

  x = 0.98*(xmax-xmin)+xmin
  d = 0.05 ; y = [ 0.99 + d ]
  s = @sprintf("kc = %7.5f ; kd = %7.5f",input.kc,input.kd)
  annotate!(x,yd!(y,d),text(s,:right,fontsize,:serif,:black))
  s = @sprintf("ki = %7.5f ; kiu = %7.5f",input.ki,input.kiu)
  annotate!(x,yd!(y,d),text(s,:right,fontsize,:serif,:black))

  savefig(filename)

end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#    Tests of errors in text backend.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Errors in keywords" begin
    @test_throws Exception pretty_table([1 2], header_crayon = [crayon"white"])
    @test_throws Exception pretty_table([1 2], subheader_crayon = [crayon"white"])
end


text:
using aoc2019.computer: load_program
using aoc2019.thrusters: output, feedback_loop
using Combinatorics

input = joinpath(@__DIR__, "input")
p = load_program(input)

perms = collect(permutations([4,3,2,1,0]))

p1 = maximum(asyncmap(perm -> output(copy(p), perm, 0), perms))
@assert p1 == 338603

perms = collect(permutations([9,8,7,6,5]))

p2 = maximum(asyncmap(perm -> feedback_loop(copy(p), perm, 0), perms))
@assert p2 == 63103596

print("-----------------------------------------------------------------------\n")
print("amplification circuit -- part one\n    max output : $p1\n")
print("amplification circuit -- part two\n    max output : $p2\n")
print("-----------------------------------------------------------------------\n")


text:
function part1(filename::AbstractString)
    sum(length(line) - length(unescape_string(line[begin+1:end-1])) for line in eachline(filename))
end

function part2(filename::AbstractString)
    sum(length(escape_string(line)) - length(line) + 2 for line in eachline(filename))
end

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code