text:
using NotInferenceDontLookHere

struct Argument
    n::Int
end

#=
struct PhiNode
    edges::Vector{Int}
    values::Vector{Any}
end
=#
Core.PhiNode() = Core.PhiNode(Any[], Any[])

#=
struct PiNode
    val::Any
    typ::Type
end
=#

#Maybe:
#struct LabelNode
#    ninstrs::Int
#end
#const Stmt = Union{Call, PhiNode, GotoIfNot, GlobalRef, ReturnNode}

struct GotoIfNot{T}
    cond::T
    dest::Int
end

struct ReturnNode{T}
    val::T
    ReturnNode{T}(val::T) where {T} = new{T}(val)
    ReturnNode{T}() where {T} = new{T}()
end

using Base.Meta
function normalize(expr)
    if isexpr(expr, :gotoifnot)
        return GotoIfNot{Any}(expr.args...)
    elseif isexpr(expr, :return)
        return ReturnNode{Any}(expr.args...)
    elseif isa(expr, LabelNode)
        return nothing
    end
    expr
end

function block_for_inst(index, inst)
    searchsortedfirst(index, inst, lt=(<=))
end
block_for_inst(cfg::CFG, inst) = block_for_inst(cfg.index, inst)

function compute_basic_blocks(stmts::Vector{Any})
    jump_dests = Set{Int}(1)
    terminators = Vector{Int}()
    # First go through and compute jump destinations
    for (idx, stmt) in pairs(stmts)
        # Terminators
        if isa(stmt, Union{GotoIfNot, GotoNode, ReturnNode})
            push!(terminators, idx)
            isa(stmt, ReturnNode) && continue
            if isa(stmt, GotoIfNot)
                push!(jump_dests, idx+1)
                push!(jump_dests, stmt.dest)
            else
                push!(jump_dests, stmt.label)
            end
        end
    end
    bb_starts = sort(collect(jump_dests))
    for i = length(stmts):-1:1
        if stmts[i] != nothing
            push!(bb_starts, i+1)
            break
        end
    end
    # Compute ranges
    basic_block_index = Int[]
    blocks = map(zip(bb_starts, Iterators.drop(bb_starts, 1))) do (first, last)
        push!(basic_block_index, first)
        BasicBlock(StmtRange(first, last-1))
    end
    popfirst!(basic_block_index)
    # Compute successors/predecessors
    for (num, b) in pairs(blocks)
        terminator = stmts[last(b.stmts)]
        # Conditional Branch
        if isa(terminator, GotoIfNot)
            block′ = block_for_inst(basic_block_index, terminator.dest)
            push!(blocks[block′].preds, num)
            push!(b.succs, block′)
        end
        if isa(terminator, GotoNode)
            block′ = block_for_inst(basic_block_index, terminator.label)
            push!(blocks[block′].preds, num)
            push!(b.succs, block′)
        elseif !isa(terminator, ReturnNode)
            if num + 1 <= length(blocks)
                push!(blocks[num+1].preds, num)
                push!(b.succs, num+1)
            end
        end
    end
    CFG(blocks, basic_block_index)
end

# Run iterated dominance frontier
function idf(cfg, defuse, domtree, slot)
    # This should be a priority queue, but TODO - sorted array for now
    pq = [(n, domtree.nodes[n].level) for n in defuse[slot].defs]
    sort!(pq, by=x->x[2])
    phiblocks = Int[]
    processed = Set{Int}()
    while !isempty(pq)
        node, level = pop!(pq)
        worklist = Int[]
        visited = Set{Int}()
        push!(worklist, node)
        while !isempty(worklist)
            active = pop!(worklist)
            for succ in cfg.blocks[active].succs
                succ_level = domtree.nodes[succ].level
                succ_level > level && continue
                succ in processed && continue
                push!(processed, succ)
                # <- TODO: Use liveness here
                push!(phiblocks, succ)
                if !(succ in defuse[slot].defs)
                    push!(pq, (succ, succ_level))
                    sort!(pq, by=x->x[2])
                end
            end

            for child in domtree.nodes[active].children
                child in visited && continue
                push!(visited, child)
                push!(worklist, child)
            end
        end
    end
    phiblocks
end

function first_insert_for_bb(code, cfg, block)
    for idx in cfg.blocks[block].stmts
        stmt = code[idx]
        if !isa(stmt, LabelNode) && !isa(stmt, PhiNode)
            return idx
        end
    end
end

function is_call(stmt, sym)
    isexpr(stmt, :call) || return false
    isa(stmt.args[1], GlobalRef) || return false
    (stmt.args[1].name === sym) || return false
    return true
end

function stmt_effect_free(stmt, src, mod)
    isa(stmt, Union{PiNode, PhiNode}) && return true
    isa(stmt, Union{ReturnNode, PhiNode, GotoNode, GotoIfNot}) && return false
    NI.statement_effect_free(stmt, src, mod)
end

function get_def(ir::IRCode, val::SSAValue)
    ir.stmts[val.id]
end

function annotate_pred!(ir, cfg, domtree, root, val, replacement)
    for bb in dominated(domtree, root)
        for idx in cfg.blocks[bb].stmts
            ir.stmts[idx] = ssamap(ir.stmts[idx]) do use
                use.id == val.id ? replacement : use
            end
        end
    end
end

function insert_node!(ir::IRCode, pos, typ, val)
    push!(ir.new_nodes, (pos, typ, val))
    return SSAValue(length(ir.stmts) + length(ir.new_nodes))
end

function predicate_insertion_pass!(ir::IRCode, domtree)
    cfg = ir.cfg
    for (idx, stmt) in pairs(ir.stmts)
        !isa(stmt, GotoIfNot) && continue
        !isa(stmt.cond, SSAValue) && continue
        def = get_def(ir, stmt.cond)
        is_call(def, :(===)) || continue
        (val, cmp) = def.args[2:3]
        if !isa(val, SSAValue)
            (cmp, val) = (val, cmp)
        end
        isa(val, SSAValue) || continue
        isa(cmp, SSAValue) && continue
        if isa(cmp, GlobalRef)
            cmp = getfield(cmp.mod, cmp.name)
        end
        isdefined(typeof(cmp), :instance) || continue
        true_type = typeof(cmp)
        isa(ir.types[val.id], Type) || continue
        false_type = Core.Compiler.typesubtract(ir.types[val.id], true_type)
        true_block = block_for_inst(cfg, idx) + 1
        false_block = stmt.dest
        true_val = insert_node!(ir, first_insert_for_bb(ir.stmts, cfg, true_block), true_type, PiNode(val, true_type))
        false_val = insert_node!(ir, first_insert_for_bb(ir.stmts, cfg, false_block), false_type, PiNode(val, false_type))
        annotate_pred!(ir, cfg, domtree, true_block, val, true_val)
        annotate_pred!(ir, cfg, domtree, false_block, val, false_val)
    end
    ir
end

function run_passes(ci::CodeInfo, mod::Module, nargs::Int)
    ci.code = map(normalize, ci.code)
    ci.code = strip_trailing_junk(ci.code)
    cfg = compute_basic_blocks(ci.code)
    defuse_insts = scan_slot_def_use(nargs, ci)
    domtree = construct_domtree(cfg)
    # @show ci.code
    ir = construct_ssa!(ci, mod, cfg, domtree, defuse_insts, nargs)
    #@show ("pre_compact", ir)
    ir = compact!(ir)
    # @show ("pre_verify", ir)
    verify_ir(ir)
    ir = predicate_insertion_pass!(ir, domtree)
    ir = compact!(ir)
    #@show ("pre_getfield_elim", ir)
    ir = getfield_elim_pass!(ir)
    ir = compact!(ir)
    # @show ("pre_lift", ir)
    ir = type_lift_pass!(ir)
    ir = compact!(ir)
    ir
end

function construct_ssa(f, args, mod::Module=Core.Main)
    ci = NI.code_typed(f, args)[1].first
    ci.code = map(normalize, ci.code)
    ci.code = strip_trailing_junk(ci.code)
    cfg = compute_basic_blocks(ci.code)
    defuse_insts = scan_slot_def_use(nargs, ci)
    domtree = construct_domtree(cfg)
    ir = construct_ssa!(ci, mod, cfg, domtree, defuse_insts, nargs)
    ir
end

function run_passes(f, args, mod::Module=Core.Main)
    ci = NI.code_typed(f, args)[1].first
    run_passes(ci, mod,  length(args.parameters))
end

function run_passes_ci(f, args)
    ci = NI.code_typed(f, args)[1].first
    ir = run_passes(ci, length(args.parameters))
    replace_code!(ci, ir, length(args.parameters))
    ci
end


text:
module Backward

export Var, grad!

mutable struct Var{T<:AbstractFloat} <: AbstractFloat
    value::T
    grad::Union{Nothing, T}
    children::Vector{Tuple{Var{T}, T}}
end

Var(x::T) where T<:AbstractFloat = Var(x, nothing, Tuple{Var{T}, T}[])

function grad!(v::Var)
    if isnothing(v.grad)
        v.grad = sum(partial * grad!(u) for (u, partial) in v.children)
    end
    v.grad
end

import Base: +, -, *, /, ^
function +(u::Var{T}, v::Var{T}) where T<:AbstractFloat
    z = Var(u.value + v.value)
    push!(u.children, (z, one(T)))
    push!(v.children, (z, one(T)))
    z
end
function -(u::Var{T}, v::Var{T}) where T<:AbstractFloat
    z = Var(u.value - v.value)
    push!(u.children, (z, one(T)))
    push!(v.children, (z, -one(T)))
    z
end
function *(u::Var{T}, v::Var{T}) where T<:AbstractFloat
    z = Var(u.value * v.value)
    push!(u.children, (z, v.value))
    push!(v.children, (z, u.value))
    z
end
function /(u::Var{T}, v::Var{T}) where T<:AbstractFloat
    z = Var(u.value / v.value)
    # dz = d(u / v) = (du * v - u * dv) / v^2
    push!(u.children, (z, one(T) / v.value))
    push!(v.children, (z, -u.value / (v.value ^ 2)))
    z
end

function ^(u::Var, k::Int)
    z = Var(u.value ^ k)
    push!(u.children, (z, k * (u.value ^ (k - 1))))
    z
end

function Base.sin(u::Var)
    z = Var(sin(u.value))
    push!(u.children, (z, cos(u.value)))
    z
end
function Base.cos(u::Var)
    z = Var(cos(u.value))
    push!(u.children, (z, -sin(u.value)))
    z
end

function Base.exp(u::Var)
    exp_u = exp(u.value)
    z = Var(exp_u)
    push!(u.children, (z, exp_u))
    z
end
function Base.log(u::Var)
    z = Var(log(u.value))
    push!(u.children, (z, inv(u.value)))
    z
end

end

ext: jl

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


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

using Markdown
using InteractiveUtils

# ╔═╡ 902eeb2e-71d2-11eb-3322-53451cb23bf2
# dependencies
begin
	using PlutoUI
	using Logging
	Logging.global_logger(Logging.SimpleLogger(stdout, LogLevel(-5000)  ) )
	Logging.disable_logging(LogLevel(-3001))
end

# ╔═╡ d0fa06da-71d1-11eb-2d00-97ea31bbb5ed
md"""# Constraint Satisfaction Problems (CSP)"""

# ╔═╡ 7321b9ee-71d2-11eb-1338-b35f914af6bf
md"""
## Application 1
How many solutions are there for the map-coloring problem in Figure 6.1 using RGB? How many solutions if four colors are allowed? Two colors? What about $n$ colors?

$(PlutoUI.LocalResource("./img/map.png", :width => 500, :align => "middle"))

"""

# ╔═╡ 13e93050-71d3-11eb-14a9-c145d4e13c8e
md"""
## Application 2
Give precise formulations for each of the following as constraint satisfaction problems:
* Rectilinear floor-planning:

    The problem of allocating given rooms within a given predefined contour shape, while satisfying given adjacency and size constraints. This work is significant with respect to complex building structures, such as hospitals, schools, offices, etc. You need to find non-overlapping places in a large rectangle for a number of smaller rectangles.

$(PlutoUI.LocalResource("./img/floorplan.gif", "width" => 500))

* Hamiltonian tour: given a network of cities connected by roads, choose an order to visit all cities in a country without repeating any. You can consider this as a relaxed version of the TSP, as the "cost" of the entire trip is not taken into account.
"""

# ╔═╡ 297786ea-71d5-11eb-3401-05d1c135d4c6
md"""
## Application 3
Consider the exam scheduling problem. We want to schedule 5 exams in the period 1-19 June. Each exam has a minimum number of days required to be able to study for it. The faculty members have some other things to do as well, so their availability is limited. The table below gives you an overview. The assignment of a course to a given date $D$ is considered as a block that starts at $D-\text{preparation days}$. The actual exam period starts on the first, and preparation days should take place in the exam period as well (e.g. SE426 can take place NET the 5th). Exams are never scheduled on the same day.

| Course  | available dates staff  |  required preparation days |
|---|---|---|
| DS425	| $1\rightarrow19$ | 3 |
|SE422	| $2\rightarrow6$ |	2 |
|SE423	| $8\rightarrow19$ |	4 |
|SE426	| $1\rightarrow16$ |	4 |
|TN421	| $1\rightarrow19$ |	1|

1. Sketch the constraint graph for this problem.
2. Use the backtracking algorithm to solve this problem. Use the minimum remaining values heuristic for variable selection, the least constraining value for value ordering and forward checking for inference. In the case of ties, use the value/variable that occurs first using numerical/alphabetical order.

"""

# ╔═╡ e3db6424-71d6-11eb-0b85-db4e9e6f4bcf
md"""
## Application 4  (Limited coding)
Consider the following logic puzzle: In five houses, each with a different color, live five persons of different nationalities, each of whom prefers a different brand of candy, a different drink, and a different pet. Given the following facts, the questions to answer are “Where does the zebra live, and in which house do they drink water?”
* The Englishman lives in the red house.
* The Spaniard owns the dog.
* The Norwegian lives in the first house on the left.
* The green house is immediately to the right of the ivory house.
* The man who eats Hershey bars lives in the house next to the man with the fox.
* Kit Kats are eaten in the yellow house.
* The Norwegian lives next to the blue house.
* The Smarties eater owns snails.
* The Snickers eater drinks orange juice.
* The Ukrainian drinks tea.
* The Japanese eats Milky Ways.
* Kit Kats are eaten in a house next to the house where the horse is kept. Coffee is drunk in the green house.
* Milk is drunk in the middle house.

Questions:
1. What could be possible variables, domains and constraints?
2. Try to understand the implementation that is given below and come up with a solution. Hint: rubber ducking it might be a good way of doing this.

"""

# ╔═╡ 6f572702-71d9-11eb-2871-d707f3ee1df7
begin

	# to be able to add new methods to the existing ones in Julia Base
	import Base.getindex, Base.get, Base.haskey, Base.in, Base.copy, Base.deepcopy

	abstract type AbstractProblem end;
	abstract type AbstractCSP <: AbstractProblem end;

	# ___________________________________________________________________ #
	# Everything related to dictionaries for domain and neighbor tracking #
	# ___________________________________________________________________ #
	struct ConstantFunctionDict{V}
		value::V
		function ConstantFunctionDict{V}(val::V) where V
			return new(val);
		end
	end

	ConstantFunctionDict(val) = ConstantFunctionDict{typeof(val)}(val);

	copy(cfd::ConstantFunctionDict) = ConstantFunctionDict{typeof(cfd.value)}(cfd.value);

	deepcopy(cfd::ConstantFunctionDict) = ConstantFunctionDict{typeof(cfd.value)}(deepcopy(cfd.value));

	mutable struct CSPDict
		dict::Union{Nothing, Dict, ConstantFunctionDict}

		function CSPDict(dictionary::Union{Nothing, Dict, ConstantFunctionDict})
			return new(dictionary);
		end
	end

	function getindex(dict::CSPDict, key)
		if (typeof(dict.dict) <: ConstantFunctionDict)
			return dict.dict.value;
		else
			return getindex(dict.dict, key);
		end
	end

	function getkey(dict::CSPDict, key, default)
		if (typeof(dict.dict) <: ConstantFunctionDict)
			return dict.dict.value;
		else
			return getkey(dict.dict, key, default);
		end
	end

	function get(dict::CSPDict, key, default)
		if (typeof(dict.dict) <: ConstantFunctionDict)
			return dict.dict.value;
		else
			return get(dict.dict, key, default);
		end
	end

	function haskey(dict::CSPDict, key)
		if (typeof(dict.dict) <: ConstantFunctionDict)
			return true;
		else
			return haskey(dict.dict, key);
		end
	end

	function in(pair::Pair, dict::CSPDict)
		if (typeof(dict.dict) <: ConstantFunctionDict)
			if (getindex(pair, 2) == dict.dict.value)
				return true;
			else
				return false;
			end
		else
			return in(pair, dict.dict);
		end
	end

	# ___________________________________________________________________ #
	# Everything related to assigning variables #
	# ___________________________________________________________________ #

	"""
		assign(problem, key, val, assignment)
	Overwrite (if an value exists already) assignment[key] with 'val'.
	"""
	function assign(problem::T, key, val, assignment::Dict) where {T <: AbstractCSP}
		assignment[key] = val;
		problem.nassigns = problem.nassigns + 1;
		nothing;
	end

	"""
		unassign(problem, key, val, assignment)
	Delete the existing (key, val) pair from 'assignment'.
	"""
	function unassign(problem::T, key, assignment::Dict) where {T <: AbstractCSP}
		if (haskey(assignment, key))
			delete!(assignment, key);
		end
		nothing;
	end

	function nconflicts(problem::T, key, val, assignment::Dict) where {T <: AbstractCSP}
		return count(
					(function(second_key)
						return (haskey(assignment, second_key) &&
							!(problem.constraints(key, val, second_key, assignment[second_key])));
					end),
					problem.neighbors[key]);
	end

	function support_pruning(problem::T) where {T <: AbstractCSP}
		if (problem.current_domains === nothing)
			problem.current_domains = Dict(collect(Pair(key, collect(problem.domains[key])) for key in problem.vars));
		end
		nothing;
	end

	function suppose(problem::T, key, val) where {T <: AbstractCSP}
		support_pruning(problem);
		local removals::AbstractVector = collect(Pair(key, a) for a in problem.current_domains[key]
												if (a != val));
		problem.current_domains[key] = [val];
		return removals;
	end

	function parse_neighbors(neighbors::String; vars::AbstractVector=[])
		local new_dict = Dict();
		for var in vars
			new_dict[var] = [];
		end
		local specs::AbstractVector = collect(map(String, split(spec, [':'])) for spec in split(neighbors, [';']));
		for (A, A_n) in specs
			A = String(strip(A));
			if (!haskey(new_dict, A))
				new_dict[A] = [];
			end
			for B in map(String, split(A_n))
				push!(new_dict[A], B);
				if (!haskey(new_dict, B))
					new_dict[B] = [];
				end
				push!(new_dict[B], A);
			end
		end
		return new_dict;
	end

	"""
		CSP

	CSP is a Constraint Satisfaction Problem implementation of AbstractProblem and AbstractCSP. This problem contains an unused initial state field to accommodate the requirements of some search algorithms.

	Note: you will need to define the constraint function for each specific problem! cf.  zebra_constraint, different_values_constraint etc.

	"""
	mutable struct CSP <: AbstractCSP
		vars::AbstractVector
		domains::CSPDict
		neighbors::CSPDict
		constraints::Function
		initial::Tuple
		current_domains::Union{Nothing, Dict}
		nassigns::Int64

		function CSP(vars::AbstractVector, domains::CSPDict, neighbors::CSPDict, constraints::Function;
					initial::Tuple=(), current_domains::Union{Nothing, Dict}=nothing, nassigns::Int64=0)
			return new(vars, domains, neighbors, constraints, initial, current_domains, nassigns)
		end
	end

	# Initialising the problem's constants
	struct ZebraInitialState
		colors::Array{String, 1}
		pets::Array{String, 1}
		drinks::Array{String, 1}
		countries::Array{String, 1}
		candies::Array{String, 1}

		function ZebraInitialState()
			local colors = map(String, split("Red Yellow Blue Green Ivory"));
			local pets = map(String, split("Dog Fox Snails Horse Zebra"));
			local drinks = map(String, split("OJ Tea Coffee Milk Water"));
			local countries = map(String, split("Englishman Spaniard Norwegian Ukranian Japanese"));
			local candies = map(String, split("KitKat Hershey Smarties Snickers Milkyway"));
			return new(colors, pets, drinks, countries, candies);
		end
	end

	zebra_constants = ZebraInitialState();

	"""
		ZebraCSP

	ZebraCSP is a Zebra Constraint Satisfaction Problem implementation of AbstractProblem and AbstractCSP.
	"""
	mutable struct ZebraCSP <: AbstractCSP
		vars::AbstractVector
		domains::CSPDict
		neighbors::CSPDict
		constraints::Function
		initial::Tuple
		current_domains::Union{Nothing, Dict}
		nassigns::Int64

		function ZebraCSP(;initial::Tuple=(), current_domains::Union{Nothing, Dict}=nothing, nassigns::Int64=0)
			# 25 variables associated with the problem
			local vars = vcat(zebra_constants.colors,
							zebra_constants.pets,
							zebra_constants.drinks,
							zebra_constants.countries,
							zebra_constants.candies);
			local domains = Dict();
			# Set all possible domain values equal to 1,2, ... , 5
			for var in vars
				domains[var] = collect(1:5);
			end
			# impose constraints
			domains["Norwegian"] = [1];
			domains["Milk"] = [3];
			neighbors = parse_neighbors("Englishman: Red;
					Spaniard: Dog; KitKat: Yellow; Hershey: Fox;
					Norwegian: Blue; Smarties: Snails; Snickers: OJ;
					Ukranian: Tea; Japanese: Milkyway; KitKat: Horse;
					Coffee: Green; Green: Ivory", vars=vars);
			for category in [zebra_constants.colors,
							zebra_constants.pets,
							zebra_constants.drinks,
							zebra_constants.countries,
							zebra_constants.candies]
				for A in category
					for B in category
						if (A != B)
							if (!(B in neighbors[A]))
								push!(neighbors[A], B);
							end
							if (!(A in neighbors[B]))
								push!(neighbors[B], A);
							end
						end
					end
				end
			end
			# return correctyl defined problem
			return new(vars, CSPDict(domains), CSPDict(neighbors), zebra_constraint, initial, current_domains, nassigns);
		end
	end

	# deal with the constraints
	function zebra_constraint(A::String, a, B::String, b; recursed::Bool=false)
		local same::Bool = (a == b); # identical numerical value
		local next_to::Bool = (abs(a - b) == 1); # neigbors
		if (A == "Englishman" && B == "Red")
			return same;
		elseif (A == "Spaniard" && B == "Dog")
			return same;
		elseif (A == "Hershey" && B == "Fox")
			return next_to;
		elseif (A == "Norwegian" && B == "Blue")
			return next_to;
		elseif (A == "KitKat" && B == "Yellow")
			return same;
		elseif (A == "Smarties" && B == "Snails")
			return same;
		elseif (A == "Snickers" && B == "OJ")
			return same;
		elseif (A == "Ukranian" && B == "Tea")
			return same;
		elseif (A == "Japanese" && B == "Milkyway")
			return same;
		elseif (A == "KitKat" && B == "Horse")
			return next_to;
		elseif (A == "Coffee" && B == "Green")
			return same;
		elseif (A == "Green" && B == "Ivory")
			return ((a - 1) == b);
		elseif (!recursed)
			return zebra_constraint(B, b, A, a, recursed=true);
		elseif ((A in zebra_constants.colors && B in zebra_constants.colors) ||
				(A in zebra_constants.pets && B in zebra_constants.pets) ||
				(A in zebra_constants.drinks && B in zebra_constants.drinks) ||
				(A in zebra_constants.countries && B in zebra_constants.countries) ||
				(A in zebra_constants.candies && B in zebra_constants.candies))
			return !same;
		else
			error("ZebraConstraintError: This constraint could not be evaluated on the given arguments!");
		end
	end

	function solve_zebra(problem::ZebraCSP, algorithm::Function; kwargs...)
		local answer = algorithm(problem; kwargs...);
		for house in collect(1:5)
			print("House ", house);
			for (key, val) in collect(answer)
				if (val == house)
					print(" ", key);
				end
			end
			println();
		end
		return answer["Zebra"], answer["Water"], problem.nassigns, answer;
	end
end;

# ╔═╡ adc27b3a-71da-11eb-325e-011cf335877e
begin
	# ___________________________________________________________________ #
	#               Actual backtracking search algorithm                  #
	# ___________________________________________________________________ #

	"""
		backtracking_search(problem)
	Search the given problem by using the backtracking search algorithm (Fig 6.5) and return the solution
	if any are found.

	kwargs
	------
		- select_unassigned_variable : Function
			by default returns a list of unasigned variables
		- order_domain_values : Function
			by default returns a list of possible domain value (no specific order)
		- inference : Function
			defaults to something that does nothing. Can be something else, e.g. forward_checking

	"""
	function backtracking_search(problem::T;
								select_unassigned_variable::Function=first_unassigned_variable,
								order_domain_values::Function=unordered_domain_values,
								inference::Function=no_inference) where {T <: AbstractCSP}
		local result = backtrack(problem, Dict(),
								select_unassigned_variable=select_unassigned_variable,
										order_domain_values=order_domain_values,
										inference=inference);
		if (!(typeof(result) <: Nothing || goal_test(problem, result)))
			error("BacktrackingSearchError: Unexpected result!")
		end
		return result;
	end

	function backtrack(problem::T, assignment::Dict;
						select_unassigned_variable::Function=first_unassigned_variable,
						order_domain_values::Function=unordered_domain_values,
						inference::Function=no_inference) where {T <: AbstractCSP}
		if (length(assignment) == length(problem.vars))
			return assignment;
		end
		local var = select_unassigned_variable(problem, assignment);
		for value in order_domain_values(problem, var, assignment)
			if (nconflicts(problem, var, value, assignment) == 0)
				assign(problem, var, value, assignment);
				removals = suppose(problem, var, value);
				if (inference(problem, var, value, assignment, removals))
					result = backtrack(problem, assignment,
										select_unassigned_variable=select_unassigned_variable,
										order_domain_values=order_domain_values,
										inference=inference);
					if (!(typeof(result) <: Nothing))
						return result;
					end
				end
				restore(problem, removals);
			end
		end
		unassign(problem, var, assignment);
		return nothing;
	end

	function restore(problem::T, removals::AbstractVector) where {T <: AbstractCSP}
		for (key, val) in removals
			push!(problem.current_domains[key], val);
		end
		nothing;
	end

	function goal_test(problem::T, state::Tuple) where {T <: AbstractCSP}
		let
			local assignment = Dict(state);
			return (length(assignment) == length(problem.vars) &&
					all((function(key)
								return nconflicts(problem, key, assignment[key], assignment) == 0;
							end)
							,
							problem.vars));
		end
	end

	function goal_test(problem::T, state::Dict) where {T <: AbstractCSP}
		let
			local assignment = deepcopy(state);
			return (length(assignment) == length(problem.vars) &&
					all((function(key)
								return nconflicts(problem, key, assignment[key], assignment) == 0;
							end),
							problem.vars));
		end
	end

	function first_unassigned_variable(problem::T, assignment::Dict) where {T <: AbstractCSP}
		return getindex(problem.vars, findfirst((function(var)
							return !haskey(assignment, var);
						end),
						problem.vars));
	end

	function unordered_domain_values(problem::T, var, assignment::Dict) where {T <: AbstractCSP}
		return choices(problem, var);
	end

	function choices(problem::T, key) where {T <: AbstractCSP}
		if (!(problem.current_domains === nothing))
			return problem.current_domains[key];
		else
			return problem.domains[key];
		end
	end

	function no_inference(problem::T, var, value, assignment::Dict, removals::Union{Nothing, AbstractVector}) where {T <: AbstractCSP}
		return true;
	end

	# ___________________________________________________________________ #
	#  inference option (forward checking, not used standardbacktracking) #
	# ___________________________________________________________________ #

	function forward_checking(problem::T, var, value, assignment::Dict, removals::Union{Nothing, AbstractVector}) where {T <: AbstractCSP}
		for B in problem.neighbors[var]
			if (!haskey(assignment, B))
				for b in copy(problem.current_domains[B])
					if (!problem.constraints(var, value, B, b))
						prune(problem, B, b, removals);
					end
				end
				if (length(problem.current_domains[B]) == 0)
					return false;
				end
			end
		end
		return true;
	end

	function prune(problem::T, key, value, removals) where {T <: AbstractCSP}
		local not_removed::Bool = true;
		for (i, element) in enumerate(problem.current_domains[key])
			if (element == value)
				deleteat!(problem.current_domains[key], i);
				not_removed = false;
				break;
			end
		end
		if (not_removed)
			error("Could not find ", value, " in ", problem.current_domains[key], " for key '", key, "' to be removed!");
		end
		if (!(typeof(removals) <: Nothing))
			push!(removals, Pair(key, value));
		end
		nothing;
	end
end;

# ╔═╡ be0f324e-71da-11eb-12da-43a7c3b865ec
begin
	prob = ZebraCSP()
	# Solve zebra problem with backtracking search
	@time solve_zebra(prob, backtracking_search);
end

# ╔═╡ 8a65d932-71df-11eb-0a50-a30132237e47
md"""
## Application 5 - another look at map coloring
"""

# ╔═╡ 9b4a7ab6-71df-11eb-2324-b115dd57c628
begin
	# function to identify if two values are the same (in case of node coloring)
	function different_values_constraint(A::T1, a::T2, B::T1, b::T2) where {T1, T2}
		return (a != b);
	end

	function MapColoringCSP(colors::AbstractVector, neighbors::String)
		local parsed_neighbors = parse_neighbors(neighbors);
		return CSP(collect(keys(parsed_neighbors)), CSPDict(ConstantFunctionDict(colors)), CSPDict(parsed_neighbors), different_values_constraint);
	end

	function MapColoringCSP(colors::AbstractVector, neighbors::Dict)
		return CSP(collect(keys(neighbors)), CSPDict(ConstantFunctionDict(colors)), CSPDict(neighbors), different_values_constraint);
	end

end;

# ╔═╡ b9bafd2c-71df-11eb-05ee-55793f3a78fa
begin
	@info "Comparing performance for map coloring"

	australia_csp = MapColoringCSP(["R", "G", "B"], "SA: WA NT Q NSW V; NT: WA Q; NSW: Q V; T: ");
	@time backtracking_search(australia_csp)
	@info *("Backtracking search without inference, $(australia_csp.nassigns) variables assigned\n", ["\t$(key): $(val)\n" for (key, val) in australia_csp.current_domains]...)

	australia_csp = MapColoringCSP(["R", "G", "B"], "SA: WA NT Q NSW V; NT: WA Q; NSW: Q V; T: ");
	@time backtracking_search(australia_csp; inference=forward_checking)
	@info *("Backtracking search with inference (forward checking), $(australia_csp.nassigns) variables assigned\n", ["\t$(key): $(val)\n" for (key, val) in australia_csp.current_domains]...)
end

# ╔═╡ 4f044672-7fe1-11eb-2a21-dbd705daf5e8
australia_csp.vars

# ╔═╡ Cell order:
# ╠═902eeb2e-71d2-11eb-3322-53451cb23bf2
# ╟─d0fa06da-71d1-11eb-2d00-97ea31bbb5ed
# ╠═7321b9ee-71d2-11eb-1338-b35f914af6bf
# ╟─13e93050-71d3-11eb-14a9-c145d4e13c8e
# ╟─297786ea-71d5-11eb-3401-05d1c135d4c6
# ╟─e3db6424-71d6-11eb-0b85-db4e9e6f4bcf
# ╠═6f572702-71d9-11eb-2871-d707f3ee1df7
# ╠═adc27b3a-71da-11eb-325e-011cf335877e
# ╠═be0f324e-71da-11eb-12da-43a7c3b865ec
# ╠═8a65d932-71df-11eb-0a50-a30132237e47
# ╠═9b4a7ab6-71df-11eb-2324-b115dd57c628
# ╠═b9bafd2c-71df-11eb-05ee-55793f3a78fa
# ╠═4f044672-7fe1-11eb-2a21-dbd705daf5e8

max_stars_repo_path: PS/PS04 - CSP.jl

max_issues_repo_path: PS/PS04 - CSP.jl

max_forks_repo_path: PS/PS04 - CSP.jl


text:
type NeuralLayer
    w::Matrix{Float64}   # weights
    b::Vector{Float64}   # biases
    a_func::Function     # activation function
    a_derv::Function     # activation funciton derivative

    # The following must be tracked for back propagation
    hx::Vector{Float64}  # input values
    pa::Vector{Float64}  # pre activation values
    pr::Vector{Float64}  # predictions (activation values)
    # Gradients
    wgr::Matrix{Float64} # weight gradient
    bgr::Vector{Float64} # bias gradient
end

# sigmoidal activation function
function sigm(a::Vector{Float64})
    1. ./ (1. + exp(-a))
end

# sigmoid derivative
function sigm_derv(a::Vector{Float64})
    s_a = sigm(a)
    s_a .* (1. - s_a)
end

function NeuralLayer(in_dim::Integer,out_dim::Integer)
    # Glorot & Bengio, 2010
    b = sqrt(6) / sqrt(in_dim + out_dim)
    NeuralLayer(2b * rand(out_dim,in_dim) - b,
                zeros(out_dim),
                sigm,
                sigm_derv,

                zeros(in_dim),
                zeros(out_dim),
                zeros(out_dim),

                zeros(out_dim,in_dim),
                zeros(out_dim)
    )
end

type ArtificialNeuralNetwork
    layers::Vector{NeuralLayer}
    hidden_sizes::Vector{Int64} # Number of nodes in each hidden level
    classes::Vector{Int64}
end

function softmax(ann_output::Vector{Float64})
    # Takes the output of a neural network and produces a valid
    # probability distribution
    ann_output = exp(ann_output)
    ann_output / sum(ann_output)
end

function ArtificialNeuralNetwork(n_hidden_units::Integer)
    ann = ArtificialNeuralNetwork(Array(NeuralLayer,0),
                                 [n_hidden_units],
                                 Array(Int64,0))
end

function forward_propagate(nl::NeuralLayer,x::Vector{Float64})
    nl.hx = x
    nl.pa = nl.b + nl.w * x
    nl.pr = nl.a_func(nl.pa)
end

function back_propagate(nl::NeuralLayer,output_gradient::Vector{Float64})
    nl.wgr = output_gradient * nl.hx' # compute weight gradient
    nl.bgr = output_gradient # compute bias gradient
    nl.w' * output_gradient # return gradient of level below
end

function init!(ann::ArtificialNeuralNetwork,
                     x::Matrix{Float64},
                     y::Vector{Int64})
    layers = Array(NeuralLayer,length(ann.hidden_sizes) + 1)
    ann.classes = unique(y)
    sort!(ann.classes)
    input_dim = size(x)[2]
    for i = 1:length(ann.hidden_sizes)
        out_dim = ann.hidden_sizes[i]
        layers[i] = NeuralLayer(input_dim,out_dim)
        input_dim = out_dim
    end
    layers[length(layers)] = NeuralLayer(input_dim,length(ann.classes))
    layers[length(layers)].a_func = softmax
    ann.layers = layers
    ann
end

function fit!(ann::ArtificialNeuralNetwork,
              x::Matrix{Float64},
              y::Vector{Int64};
              epochs::Int64 = 5,
              alpha::Float64 = 0.1,
              lambda::Float64 = 0.1)
    init!(ann,x,y)
    n_obs, n_feats = size(x)
    layers = ann.layers
    n_layers = length(layers)
    for _ = 1:epochs
        for i = 1:n_obs
            y_hat = zeros(length(ann.classes))
            y_hat[findfirst(ann.classes,y[i])] = 1.

            y_pred = predict(ann,x[i,:][:])
            output_gradient = -(y_hat - y_pred)
            for j = n_layers:-1:2
                # This returns the gradient of the j-1 layer
                next_layer_gr = back_propagate(layers[j],output_gradient)
                next_layer = layers[j-1]
                output_gradient = next_layer_gr .* next_layer.a_derv(next_layer.pa)
            end
            back_propagate(layers[1],output_gradient)

            # Compute delta and step
            for j = 1:n_layers
                nl = layers[j]
                # Computer L2 weight penatly
                weight_delta = (-nl.wgr) - (lambda * (2 * nl.w))
                nl.w = alpha * weight_delta + nl.w
                nl.b = alpha * (-nl.bgr) + nl.b
            end
        end
    end
end

# Predict class probabilities for a given observation
function predict(ann::ArtificialNeuralNetwork,x::Vector{Float64})
    for i in 1:length(ann.layers)
        x = forward_propagate(ann.layers[i],x)
    end
    x
end

# Handle a Matrix input
function predict(ann::ArtificialNeuralNetwork,x::Matrix{Float64})
    n_obs,n_feats = size(x)
    y_proba = zeros((n_obs,length(ann.classes)))
    for i = 1:n_obs
        y_proba[i,:] = predict(ann,x[i,:][:])
    end
    y_proba
end


text:
module Quadratic

using STMO: Tracker, notrack, PathTrack, trace

export fquad, plot_quadratic, solve_quadratic, quadratic_ls, gradient_descent

using LinearAlgebra, Plots

"""Compute a quadratic function."""
fquad(x::Vector, P::Matrix, q::Vector, r::Real=0.0) = 0.5x' * P * x + q' * x + r

"""Plot a 1-D quadratic function."""
function plot_quadratic(p::Real, q::Real, r::Real, (xl, xu), stepsize=0.1; kwargs...)
    return plot(x -> p*x^2 + q*x+r, xl:stepsize:xu, xlabel="\$x\$"; kwargs...)
end

"""Plot a 2-D quadratic function."""
function plot_quadratic(P::AbstractMatrix, q::AbstractVector, r::Real,
                            (x1l, x1u), (x2l, x2u), stepsize=0.1; kwargs...)
    fun = (x1, x2) -> [x1,x2] |> x -> 0.5x' * P * x + q' * x + r
    return contour(x1l:stepsize:x1u, x2l:stepsize:x2u, fun, xlabel="\$x_1\$",
                                            ylabel="\$x_2\$"; kwargs...)
end

"""
    solve_quadratic(p::Real, q::Real[, r::Real=0.0])

Finds the minimizer of an 1-D quadratic system. If the quadratic term `p` is
negative, the function raises an error.

Inputs:
    - `p`, `q`, `r`: the terms of the 1-D quadratic system

Output:
    - xstar: the minimizer, a number
"""
function solve_quadratic(p::Real, q::Real, r::Real=0.0)
    @assert p > 0.0
    return - q / p
end

"""
    solve_quadratic(P::AbstractMatrix, q::AbstractVector, r::Real=0)

Finds the minimizer of an N-D quadratic system.
P is assumed to be a symmetric positive-definite matrix.

Inputs:
    - P, q, r: the terms of the nD quadratic system

Output:
    - xstar: the minimizer, an (n x 1) vector
"""
function solve_quadratic(P::AbstractMatrix, q::AbstractVector, r::Real=0.0)
    return - P \ q
end

"""
    quadratic_ls(P::AbstractMatrix, q::AbstractVector, Δx::AbstractVector,
                                    x::AbstractVector)

Find the exact step size that minimized a quadratic system in
a given point x for a given search direction Dx

Inputs:
    - P, q: the terms of the nD quadratic system
    - x: starting point
    - Δx: search direction

Output:
    - t: optimal step size
"""
function quadratic_ls(P::AbstractMatrix, q::AbstractVector, Δx::AbstractVector,
                                    x::AbstractVector)
    t = solve_quadratic(Δx' * P * Δx, Δx' * P * x + Δx ⋅ q)
    return t
end

"""
    gradient_descent(P::AbstractArray, q::AbstractVector,
            x₀::AbstractVector; β::Real=0.0, ϵ::Real=1e-6,
            tracker::Tracker=notrack)

Computes the minimizes of a quadratic system using gradient descent. Optionally
provide momentum.

Inputs:
    - P, q: the terms of the nD quadratic system
    - x₀: starting point
    - ϵ: convergence parameter
    - β: momentum parameter
    - tracker: object of the type `Tracker` to save the steps

Outputs:
    - xstar: the found minimum
"""
function gradient_descent(P::AbstractArray, q::AbstractVector,
            x₀::AbstractVector; β::Real=0.0, ϵ::Real=1e-6,
            tracker::Tracker=notrack)
    @assert 0 ≤ β < 1
    x = x₀  # initial value
    Δx = zero(x)  # pre-allocate a vector for the gradient
    while true
        Δx .= (1.0 - β) * (- P * x .- q) .+ β *  Δx
        if norm(Δx) < ϵ
            break
        end
        # determine stepsize using exact line search
        t = quadratic_ls(P, q, Δx, x)
        # perform step
        x .+= t * Δx
        trace(tracker, x)  # saves the steps
    end
    return x
end

# SIGNAL RECOVERY
# ---------------

module SignalRecovery
using LinearAlgebra: I

signalfun(x, n) = 3sin(x * 2 * pi / n) +
            2cos(x * 4 * pi / n) +
            sin(x * 4 * pi / n) + 0.8 * cos(x * 12 * pi / n)

"""
    generate_noisy_measurements(m, n; σ=1.0)

Generate noisy measurements according to some function f

Inputs:
    - m : number of observations
    - n : dim of x
    - σ : normally distributed noise (default = 1)

Output:
    - y : vector of noisy measurements
    - ind : vector of indices
"""
function generate_noisy_measurements(m, n; σ=1.0)
    ind = rand(1:n, m)
    y = signalfun.(ind, n) .+ randn(m) * σ
    return y, ind
end

"""
    make_connection_matrix(n, γ=100)

Generates the kernel matrix and the inverse kernel matrix

Uses γ as a characteristic length scale of the radial basis kernel,
assumes periodic boundaries.

Inputs:
    - n : number of points
    - γ : length scale

Output:
    - K
    - Kinv
"""
function make_connection_matrix(n, γ=100)
    # cyclic distance
    d(i, j) = min((i-j)^2, ((i-n)-j)^2, ((j-n)-i)^2)
    # distance matrix
    D = [d(i, j) for i in 1:n, j in 1:n]
    # kernel matrix
    K = exp.(-D / γ^2) + 1e-2I
    return K, inv(K)
end

"""
    make_bookkeeping(I, n)

Constructs the bookkeeping matrix R

Inputs:
    - ind : the indices of the measurements
    - n : dimensionality of signal vector

Output:
    - R : m x n bookkeeping matrix
"""
function make_bookkeeping(ind, n)
    m = length(ind)
    R = zeros(Int, m, n)
    for (i, j) in enumerate(ind)
        R[i,j] = 1
    end
    return R
end

export generate_noisy_measurements, make_connection_matrix, make_bookkeeping, signalfun

end # module SignalRecovery

export SignalRecovery

end  # module Quadratic


text:
"""
    Model(mesh, ρ, D, T₂, κ, γ)

Discretized biological model with initial spin densities `ρ`, diffusion tensors `D`,
T₂-relaxation times `T₂`, wall permeabilities `κ`, and gyromacnetic ratio `γ` (defaults to
the one of water protons).
"""
Base.@kwdef struct Model{T}
    mesh::FEMesh{T}
    ρ::Vector{Complex{T}}
    D::Vector{SMatrix{3,3,T,9}}
    T₂::Vector{T}
    κ::Vector{T}
    γ::T
end


text:
const gbm_study_uids = [
    "1.3.6.1.4.1.14519.5.2.1.4591.4001.278082550121070125285213632206",
    "1.3.6.1.4.1.14519.5.2.1.4591.4001.335353575986269052491315637674",
    "1.3.6.1.4.1.14519.5.2.1.4591.4001.365805576275232517344053939830",
    "1.3.6.1.4.1.14519.5.2.1.4591.4001.100057969162276274933613772317",
    "1.3.6.1.4.1.14519.5.2.1.4591.4001.269887096484012292940330991126",
    "1.3.6.1.4.1.14519.5.2.1.4591.4001.961791689281776173751323306588",
    "1.3.6.1.4.1.14519.5.2.1.4591.4001.763554173270318063812534542847",
    "1.3.6.1.4.1.14519.5.2.1.4591.4001.304604545029494418165835320551"
    ]

const alt_gbm_name = Dict(
    "TCGA-06-0185-1" => "1.3.6.1.4.1.14519.5.2.1.4591.4001.278082550121070125285213632206",
    "TCGA-06-0185-2" => "1.3.6.1.4.1.14519.5.2.1.4591.4001.335353575986269052491315637674",
    "TCGA-06-0185-3" => "1.3.6.1.4.1.14519.5.2.1.4591.4001.365805576275232517344053939830",
    "TCGA-06-0881-1" => "1.3.6.1.4.1.14519.5.2.1.4591.4001.100057969162276274933613772317",
    "TCGA-06-0881-2" => "1.3.6.1.4.1.14519.5.2.1.4591.4001.269887096484012292940330991126",
    "TCGA-06-1802-1" => "1.3.6.1.4.1.14519.5.2.1.4591.4001.961791689281776173751323306588",
    "TCGA-06-2570-1" => "1.3.6.1.4.1.14519.5.2.1.4591.4001.763554173270318063812534542847",
    "TCGA-06-5417-1" => "1.3.6.1.4.1.14519.5.2.1.4591.4001.304604545029494418165835320551"
    )

function download_invivo_studies(studies::AbstractVector{String}=gbm_study_uids; destination::AbstractString, overwrite = false)
    make_folder(destination)

    number_of_studies = length(studies)
    vfa_folders, dce_folders = [fill("", number_of_studies) for _=1:2]
    for (index, study) in enumerate(studies)
        (vfa_folders[index], dce_folders[index]) = download_vfa_and_dce(study, destination; overwrite=overwrite)
    end
    return (vfa_folders = vfa_folders, dce_folders = dce_folders)
end

download_invivo_studies(study::AbstractString; kwargs...) = download_invivo_studies([study]; kwargs...)

function download_vfa_and_dce(study_id, destination; overwrite = false)
    vfa_folder = joinpath(destination, study_id, "vfa")
    dce_folder = joinpath(destination, study_id, "dce")

    if isdir(vfa_folder) && isdir(dce_folder)
        if !overwrite
            return (vfa_folder, dce_folder)
        end
    end

    gbm_series = series(study = study_id)
    vfa_series = find_vfa_series(gbm_series)
    dce_series = find_dce_series(gbm_series)

    if !isdir(vfa_folder) || overwrite == true
        download_series(vfa_series, destination = vfa_folder)
    end

    if !isdir(dce_folder) || overwrite == true
        download_series(dce_series, destination = dce_folder)
    end
    return (vfa_folder, dce_folder)
end

find_vfa_series(series_dataframe) = find_in_description("MAP", series_dataframe)
find_dce_series(series_dataframe) = find_in_description("DYN", series_dataframe)

function find_in_description(word_to_find::AbstractString, series_dataframe)
    descriptions = series_dataframe.SeriesDescription
    found_indices = findall(occursin.(word_to_find, descriptions))
    if length(found_indices) < 1
        errant_id = series_dataframe.PatientID[1]
        errant_study = series_dataframe.StudyInstanceUID[1]
        error("No single series in $errant_id found with $word_to_find in their description.
              The full study UID is $errant_study\n")
    elseif length(found_indices) > 1
        errant_id = series_dataframe.PatientID[1]
        @warn "Multiple series in $errant_id found containing $word_to_find in their description.
              An arbitrary series will be used among the found series"
        found_index = found_indices[end]
    else
        found_index = found_indices[1]
    end
    found_series = series_dataframe.SeriesInstanceUID[found_index]
    return found_series
end

function download_series(series_id::AbstractString; destination)
    make_folder(destination; remove_existing = true)
    zip_file = joinpath(destination, "downloaded.zip")
    images(series = series_id, file = zip_file)
    unzip_command = `unzip -o $zip_file -d $destination`
    run(unzip_command)
    rm(zip_file)
    return destination
end

function download_invivo_masks(; destination, overwrite = false)
    if isdir(destination) && overwrite == false
        return destination
    end
    make_folder(destination; remove_existing = true)
    zip_file = joinpath(destination, "invivo_masks.zip")
    download("https://osf.io/uxe3p/download", zip_file)
    unzip_cmd = `unzip -o $zip_file -d $destination`
    run(unzip_cmd)
    rm(zip_file)
    return destination
end

function download_invivo_preprocessed(; destination, overwrite = false)
    if isdir(destination) && overwrite == false
        return destination
    end
    make_folder(destination; remove_existing = true)
    zip_file = joinpath(destination, "preprocessed.zip")
    download("https://osf.io/d3ext/download", zip_file)
    unzip_cmd = `unzip -o $zip_file -d $destination`
    run(unzip_cmd)
    rm(zip_file)
    return destination
end


text:
using LPC
using Base.Test # can I use FactCheck instead ?

a, err, relf = lpc_burg(rand(10), 3)

@test false # TODO


text:
module Geom

using Color
using Compat
using Compose
using DataArrays
using DataStructures
using Gadfly

import Compose.combine # Prevent DataFrame.combine from taking over.
import Gadfly: render, layers, element_aesthetics, inherit, escape_id,
               default_statistic, default_scales, element_coordinate_type,
               ScaleElement, svg_color_class_from_label, isconcrete,
               concretize
import Iterators
import Iterators: cycle, product, distinct, takestrict, chain, repeated

# Geometry that renders nothing.
immutable Nil <: Gadfly.GeometryElement
end

const nil = Nil

function render(geom::Nil, theme::Gadfly.Theme, aes::Gadfly.Aesthetics,
                data::Gadfly.Data, scales::Dict{Symbol, ScaleElement},
                subplot_layer_aess::Vector{Gadfly.Aesthetics})
end

# Subplot geometries require some more arguments to render. A simpler render
# function is defined and passed through to here for non-subplot geometries.
function render(geom::Gadfly.GeometryElement, theme::Gadfly.Theme, aes::Gadfly.Aesthetics,
                subplot_layer_aess::Union(Nothing, Vector{Gadfly.Aesthetics}),
                subplot_layer_datas::Union(Nothing, Vector{Gadfly.Data}),
                scales::Dict{Symbol, ScaleElement})
    render(geom, theme, aes)
end

# Catchall
function default_statistic(::Gadfly.GeometryElement)
    return Gadfly.Stat.identity()
end

include("geom/bar.jl")
include("geom/boxplot.jl")
include("geom/errorbar.jl")
include("geom/hexbin.jl")
include("geom/hline.jl")
include("geom/label.jl")
include("geom/line.jl")
include("geom/point.jl")
include("geom/rectbin.jl")
include("geom/subplot.jl")
include("geom/vline.jl")
include("geom/ribbon.jl")
include("geom/violin.jl")
include("geom/polygon.jl")

end # module Geom


text:
	struct Scaling <: AbstractScaling
		W::Matrix{Float64}
		iW::Matrix{Float64}
		iWiW::Matrix{Float64}
		l::Vector{Float64}

		#per cone
		mu::Vector{Float64}
		wbs::Vector{Float64} # sqrt(s/z) for POC, wb for SOC
		zik::Vector{Float64} # storage for compute_scaling
		sik::Vector{Float64}
		function Scaling(p::Problem)
			return Scaling(zeros(p.k, p.k), zeros(p.k,p.k), zeros(p.k,p.k), zeros(p.k), p.cones)
		end
		function Scaling(W,iW,iWiW,l,cones)
			md = maximum(conedim.(cones))
			zik,sik = zeros(md),zeros(md)
			return new(W, iW, iWiW, l, zeros(length(cones)), zeros(cones[end].offs + conedim(cones[end])), zik, sik)
		end
	end

	function compute_scaling(cone::POC{dim}, cind, scaling::Scaling, s, z) where dim
		for i = 1:dim
			ii = cti(cone, i)
			scaling.W[ii,ii] = sqrt(s[ii] / z[ii])
			scaling.iW[ii,ii] = sqrt(z[ii] / s[ii])
			scaling.l[ii] = sqrt(s[ii] * z[ii])
			scaling.wbs[ii] = sqrt(s[ii] / z[ii])
		end
	end

	function compute_scaling(c::SOC{dim}, cind, scaling::Scaling, s, z) where dim
		sik,zik = scaling.sik,scaling.zik
		for ii = 1:dim
			i = cti(c, ii)
			sik[ii] = s[i]
			zik[ii] = z[i]
		end
		onrmz = zik[1]^2
		onrms = sik[1]^2
		wl = dim
		for i=2:wl
			onrmz -= zik[i] * zik[i]
			onrms -= sik[i] * sik[i]
		end
		nrmz = sqrt(onrmz)
		nrms = sqrt(onrms)
		rmul!(zik, 1/nrmz)
		rmul!(sik, 1/nrms)
		zbk = zik
		sbk = sik

		nsum = 0.0
		for i = 1:wl
			nsum += zbk[i] * sbk[i]
		end
		gamma = sqrt((1 + nsum)/2)

		wb = @view scaling.wbs[cti(c,1):cti(c,dim)]
		wb[1] = sbk[1] + zbk[1]
		for i = 2:wl
			wb[i] = sbk[i] - zbk[i]
		end
		rmul!(wb, 1/(2*gamma))

		bl = length(wb)-1
		denom = wb[1] + 1
		mu = sqrt(nrms/nrmz)
		scaling.mu[cind] = mu
		for i=1:bl, j=1:bl
			cellv = ((i==j) ? 1.0 : 0.0) + wb[i+1]*wb[j+1]/denom
			ii,ji = cti(c, i), cti(c, j)
			scaling.W[ii+1,ji+1] = cellv*mu
			scaling.iW[ii+1,ji+1] = cellv/mu
		end

		ic1 = cti(c,1)
		for i=1:wl
			scaling.W[ic1,cti(c, i)] = wb[i]*mu
		end
		scaling.iW[ic1,ic1] = wb[1]/mu

		for i=2:wl
			ii = cti(c, i)
			scaling.W[ii,ic1] = wb[i]*mu
			scaling.iW[ic1,ii] = -wb[i]/mu
			scaling.iW[ii,ic1] = -wb[i]/mu
		end

		ziv,siv = zbk[1], sbk[1]
		tmv1 = sqrt(nrms*nrmz)
		mult = tmv1/(ziv + siv + 2*gamma)
		rmul!(sbk, gamma + ziv)
		rmul!(zbk, gamma + siv)
		for i=2:wl
			scaling.l[cti(c, i)] = (sbk[i] + zbk[i])*mult
		end
		scaling.l[ic1] = gamma*tmv1
	end

	@unroll function compute_scaling(cones::Tuple{Vararg{Cone}}, scaling::Scaling, s, z)
		# assume that cones is in the order POC ... POC SOC ... SOC
		i=1
		@unroll for cone in cones
			compute_scaling(cone, i, scaling, s, z)
			i += 1
		end
		mul!(scaling.iWiW, scaling.iW, scaling.iW')
		return scaling
	end

	function scale!(cone::POC{dim}, cind, scl, s, op) where dim
		wb = scl.wbs
		for i=cti(cone,1):cti(cone, dim)
			op[i] = wb[i]*s[i]
		end
	end

	function iscale!(cone::POC{dim}, cind, scl, s, op) where dim
		wb = scl.wbs
		for i=cti(cone,1):cti(cone, dim)
			op[i] = 1/wb[i]*s[i]
		end
	end

	function scale!(cone::SOC{dim}, cind, scl, s, op) where dim
		wb = scl.wbs
		mu = scl.mu[cind]
		del = 0.0
		for i=cti(cone, 2):cti(cone, dim)
			del += wb[i]*s[i]
		end

		ind1 = cti(cone, 1)
		cst = (s[ind1] + del/(1+wb[ind1]))
		op[ind1] = mu * (wb[ind1]*s[ind1] + del)
		for i=cti(cone, 2):cti(cone, dim)
			op[i] = mu * (s[i] + cst*wb[i])
		end
	end

	function iscale!(cone::SOC{dim}, cind, scl, s, op) where dim
		wb = scl.wbs
		mu = scl.mu[cind]
		del = 0.0
		for i=cti(cone, 2):cti(cone, dim)
			del += wb[i]*s[i]
		end

		ind1 = cti(cone, 1)
		cst = (-s[ind1] + del/(1+wb[ind1]))
		op[ind1] = 1/mu * (wb[ind1]*s[ind1] - del)
		for i=cti(cone, 2):cti(cone, dim)
			op[i] = 1/mu * (s[i] + cst*wb[i])
		end
	end

	# computes W s
	@unroll function scale!(cones::Tuple{Vararg{Cone}}, scl::Union{SqrScaling, Scaling}, s, op)
		i=1
		@unroll for cone in cones
			scale!(cone, i, scl, s, op)
			i+=1
		end
	end
	# computes W^-1 s
	@unroll function iscale!(cones::Tuple{Vararg{Cone}}, scl::Union{SqrScaling, Scaling}, s, op)
		i=1
		@unroll for cone in cones
			iscale!(cone, i, scl, s, op)
			i+=1
		end
	end

	# update NT scalings

	function update_scaling!(cone::POC, scaling, s, z)
		mat,imat,sv = scaling
	end

	function update_scaling!(cone::SOC, scaling, s, z)
	end

	function update_scaling!(cones::Union{Vector{Cone},Vector{Cone{d}} where d}, scaling, s, z)

	end


text:
using AutoPreallocation
using BenchmarkTools
using Test

# functions to test on
f_ones() = ones(64)
f_matmul() = ones(32,64) * ones(64, 2)

f_matmul_noprealloc() = ones(32,64) * @no_prealloc(ones(64, 2))

@testset "ones example" begin
    @assert (@ballocated f_ones()) === 624
    val, record = record_allocations(f_ones)
    @test (@ballocated avoid_allocations($record, f_ones)) <= 64
end

@testset "matmul example" begin
    @assert (@ballocated f_matmul()) === 18_304
    val, record = record_allocations(f_matmul)
    # NOTE: (@Roger-luo) not sure why this is 256 on my machine
    @test (@ballocated avoid_allocations($record, f_matmul)) <= 352

    _, pf = preallocate(f_matmul)
    @test (@ballocated $pf()) <= 352
end

@testset "noprealloc example" begin
    @assert (@ballocated f_matmul_noprealloc()) === 18_304
    val, record = record_allocations(f_matmul_noprealloc)
    @test length(record.allocations) == 2
    @test record.initial_sizes == [(32, 64), (32, 2)]
    @test (@ballocated avoid_allocations($record, f_matmul_noprealloc)) <= 1520
end

@testset "check thread-safe" begin
    f(x, y) = x * y
    n = Threads.nthreads()
    results = Vector{Any}(undef, n)
    As = [rand(4, 4) for _ in 1:n]
    Bs = [rand(4, 4) for _ in 1:n]
    _, pf = preallocate(f, As[1], Bs[1])
    Threads.@threads for k in 1:n
        results[k] = pf(As[k], Bs[k])
    end

    # if it's not thread-safe, the result
    # will be modified by other multiplications
    # since they will share the same memory
    for k in 1:n
        @test results[k] ≈ f(As[k], Bs[k])
    end
end

@testset "resizing operations" begin
    function push_pop_test(a)
        x = [0, 0, 0]
        @test length(x) == 3
        push!(x, a)
        push!(x, 10a)
        @test length(x) == 5
        pop!(x)
        @test length(x) == 4
    end

    _, p_push_pop_test = preallocate(push_pop_test, 1)
    p_push_pop_test(2)  # check works on second call
    p_push_pop_test(3)  # check works on third call
    p_push_pop_test(4)  # check works on fourth call
end


text:
export initDropWav
export dropWav

"""
    y = dropWav(wav::Array, fs::Real=16000.0; ratio::Real=0.05)
droping frames to simulate network packet loss.
"""
function dropWav(wav::Array, fs::Real=16000.0; ratio::Real=0.05)
    ZERO = eltype(wav)(0.0)
    winlen = floor(Int, 0.016 * fs)                # 0.016毫秒一帧
    frames = div(length(wav)-winlen, winlen) + 1   # 帧数
    firstIds = (0:(frames-1)) .* winlen .+ 1       # 帧起始下标
    lasstIds = firstIds .+ (winlen - 1)            # 帧结束下标
    for t = 1:frames
        if rand()<ratio
            wav[firstIds[t]:lasstIds[t]] .= ZERO;
        end
    end
    return wav
end

"""
    initDropWav(fs::Real, ratioSpan::NTuple{2,Number}) -> dropwav(wav::Array)
init packet-loss-simulation function.
+ `fs` sampling rate e.g. 16000.0
+ `ratioSpan` e.g. (0.02, 0.09)
"""
function initDropWav(fs::Real, ratioSpan::NTuple{2,Number})
    MinR, MaxR = ratioSpan
    @assert MinR <= MaxR
    function dropwav(wav::Array)
        ratio = rand()*(MaxR - MinR) + MinR
        return dropWav(wav, fs; ratio=ratio)
    end
    return dropwav
end

max_stars_repo_name: sonosole/SpeechAugment

max_issues_repo_name: sonosole/SpeechAugment

max_forks_repo_name: sonosole/SpeechAugment


text:
"""
    nat_spec_filter(l, D; σ = 0.25 * maximum(D), method = :regular, thres = 0.2)

assemble the natural spectral graph filter centered at the l-th eigenvector via
the distance matrix `D`.

# Input Arguments
- `l::Int64`: index of the centered eigenvector
- `D::Matrix{Float64}`: non-trivial distance matrix of the eigenvectors
- `σ::Float64`: Gaussian window width parameter (default: `0.25 * maximum(D)`)
- `method::Symbol`: `:regular` or `:reduced` (default: `:regular`)
- `thres::Float64`: cutoff threshold ∈ (0, 1).

# Output Argument
- `𝛍::Vector{Float64}`: the natural spectral graph filter

"""
function nat_spec_filter(l, D; σ = 0.25 * maximum(D), method = :regular, thres = 0.2)
    d = D[:, l]
    𝛍 = exp.(-(d ./ σ).^2)
    𝛍 ./= sum(𝛍)
    if method == :reduced
        𝛍 .*= (𝛍 .> (thres * maximum(𝛍)))
    end
    return 𝛍
end

"""
    ngwf_all_vectors(D, 𝚽; σ = 0.2 * maximum(D))

assemble the whole NGWF dictionary.

# Input Arguments
- `D::Matrix{Float64}`: non-trivial distance matrix of the eigenvectors
- `𝚽::Matrix{Float64}`: graph Laplacian eigenvectors
- `σ::Float64`: Gaussian window width parameter (default: `0.25 * maximum(D)`)

# Output Argument
- `𝓤::Matrix{Float64}`: the NGWF dictionary

"""
function ngwf_all_vectors(D, 𝚽; σ = 0.2 * maximum(D))
    N = size(D, 1)
    𝓤 = zeros(N, 0)
    for l = 1:N
        𝛍 = nat_spec_filter(l, D; σ = σ)
        P = 𝚽 * Diagonal(𝛍) * 𝚽'
        𝓤 = hcat(𝓤, P)
    end
    return 𝓤
end

"""
    rngwf_all_vectors(D, 𝚽; σ = 0.2 * maximum(D), thres = 0.2)

assemble the reduced NGWF (rNGWF) dictionary.

# Input Arguments
- `D::Matrix{Float64}`: non-trivial distance matrix of the eigenvectors
- `𝚽::Matrix{Float64}`: graph Laplacian eigenvectors
- `σ::Float64`: Gaussian window width parameter (default: `0.25 * maximum(D)`)
- `thres::Float64`: cutoff threshold ∈ (0, 1).

# Output Argument
- `𝓤::Matrix{Float64}`: the rNGWF dictionary
- `dic_l2x::Dict`: a dictionary to store the filtered locations by QR at the l-th
    centered eigenvector

"""
function rngwf_all_vectors(D, 𝚽; σ = 0.2 * maximum(D), thres = 0.2)
    N = size(D, 1)
    𝓤 = zeros(N, 0)
    dic_l2x = Dict()
    for l = 1:N
        𝛍 = nat_spec_filter(l, D; σ = σ, method = :reduced, thres = thres)
        P = 𝚽 * Diagonal(𝛍) * 𝚽'
        dic_l2x[l] = qr(P, Val(true)).p[1:rank(P, rtol = 10^4 * eps())]
        𝓤 = hcat(𝓤, P[:, dic_l2x[l]])
    end
    return 𝓤, dic_l2x
end

function ngwf_vector(D, l, x, 𝚽; σ = 0.1 * maximum(D))
    N = size(𝚽, 1)
    P = 𝚽 * Diagonal(nat_spec_filter(l, D; σ = σ)) * 𝚽'
    ψ = P * spike(x, N)
    ψ ./= norm(ψ, 2)
    return ψ
end

"""
    frame_approx(f, U, V; num_kept = length(f))

approximate signal `f` by the frame `U`.

# Input Arguments
- `f::Vector{Float64}`: input graph signal
- `U::Matrix{Float64}`: a frame operator (matrix or dictionary)
- `V::Matrix{Float64}`: the dual frame operator of `U`
- `num_kept::Int64`: number of kept coefficients (NCR)

# Output Argument
- `rel_error::Vector{Float64}`: the relative errors
- `f_approx::Vector{Float64}`: the approximated signal

"""
function frame_approx(f, U, V; num_kept = length(f))
    g = U' * f
    ind = sortperm(g.^2; rev = true)[1:num_kept]
    f_approx = zeros(length(f))
    rel_error = [1.0]
    for γ in ind
        f_approx += g[γ] * V[:, γ]
        f_res = f - f_approx
        push!(rel_error, norm(f_res)/norm(f))
    end
    return rel_error, f_approx
end

"""
    rngwf_lx(dic_l2x)

find the sequential subindices of rNGWF vectors.

# Input Arguments
- `dic_l2x::Dict`: a dictionary to store the filtered locations by QR at the l-th
    centered eigenvector

# Output Argument
- `Γ::Vector{Tuple{Int64,Int64}}`: the sequential subindices of rNGWF vectors.

"""
function rngwf_lx(dic_l2x)
    N = length(dic_l2x)
    Γ = (Tuple{Int64,Int64})[]
    for l = 1:N
        for x in dic_l2x[l]
            push!(Γ, (l - 1, x))
        end
    end
    return Γ
end


text:
using ITensors
using ITensorVisualizationBase

s = siteinds("S=1/2", 5; conserve_qns=true)
ψ = randomMPS(s, n -> isodd(n) ? "↑" : "↓"; linkdims=2)
orthogonalize!(ψ, 2)
ψdag = prime(linkinds, dag(ψ))
tn = [ψ..., ψdag...]

edge_labels = (; plevs=true, qns=true)
@visualize fig tn edge_labels = edge_labels edge_textsize = 20

fig


text:
# First lets define a fast natural inner join, before worrying about other kinds
# of joins.

# Our initial implementation will use a simple hashing algorithm
@generated function Base.join{Names1, Names2}(t1::Table{Names1}, t2::Table{Names2})
    int_names = (intersect(Names1, Names2)...)
    int_indices = collect(columnindex(Names1, int_names))
    int_types = Expr(:curly, :Tuple, eltypes(t1).parameters[int_indices]...)

    quote
        # Make a hashmap of the overlapping part of t1, indicating the relevant
        # rows`
        hashmap = Dict{Row{$int_names, $int_types}, Vector{Int}}()
        for i1 = 1:nrow(t1)
            @inbounds int_row = t1[i1, Val{$int_names}]
            #=
            if haskey(hashmap, int_row) # Can we avoid querying twice? Yes
                push!(hashmap[int_row], i1)
            else
                hashmap[int_row] = [i1]
            end
            =#
            push!(get!(hashmap, int_row, Vector{Int}(0)), i1)
        end

        # Make an empty output table
        out = jointype(t1, t2)()

        # Now see if t2 overlaps with the hashmap and push onto the solution as
        # appropriate
        for i2 = 1:nrow(t2)
            @inbounds int_row = t2[i2, Val{$int_names}]
            if haskey(hashmap, int_row)
                for i1 in hashmap[int_row]
                    # Combine matching rows and push onto output table
                    push!(out, joinrow(t1[i1], t2[i2]))
                end
            end
        end
        return out
    end
end

@generated function jointype{Names1,Types1,Names2,Types2}(t1::Table{Names1,Types1}, t2::Table{Names2,Types2})
    names_ab = (intersect(Names1, Names2)...)
    indices_ab = collect(columnindex(Names1, names_ab))
    types_ab = storagetypes(t1).parameters[indices_ab]

    names_a = (setdiff(Names1, Names2)...)
    indices_a = collect(columnindex(Names1, names_a))
    types_a = storagetypes(t1).parameters[indices_a]

    names_b = (setdiff(Names2, Names1)...)
    indices_b = collect(columnindex(Names2, names_b))
    types_b = storagetypes(t2).parameters[indices_b]

    new_names = (names_ab..., names_a..., names_b...)
    new_types = Expr(:curly, :Tuple, types_ab..., types_a..., types_b...)

    return :(Table{$new_names, $new_types})
end

@generated function joinrow{Names1,Types1,Names2,Types2}(r1::Row{Names1,Types1}, r2::Row{Names2,Types2})
    names_ab = (intersect(Names1, Names2)...)
    indices_ab = collect(columnindex(Names1, names_ab))
    types_ab = eltypes(r1).parameters[indices_ab]

    names_a = (setdiff(Names1, Names2)...)
    indices_a = collect(columnindex(Names1, names_a))
    types_a = eltypes(r1).parameters[indices_a]

    names_b = (setdiff(Names2, Names1)...)
    indices_b = collect(columnindex(Names2, names_b))
    types_b = eltypes(r2).parameters[indices_b]

    new_names = (names_ab..., names_a..., names_b...)
    new_types = Expr(:curly, :Tuple, types_ab..., types_a..., types_b...)

    exprs = vcat([:(r1.data[$(indices_ab[j])]) for j = 1:length(indices_ab)],
                 [:(r1.data[$(indices_a[j])]) for j = 1:length(indices_a)],
                 [:(r2.data[$(indices_b[j])]) for j = 1:length(indices_b)])

    return Expr(:call, :(Row{$new_names, $new_types}), Expr(:tuple, exprs...))
end

#=
# Built-in join types
abstract AbstractJoin
immutable InnerJoin <: AbstractJoin; end
immutable LeftJoin <: AbstractJoin; end
immutable RightJoin <: AbstractJoin; end
immutable OuterJoin <: AbstractJoin; end

# A trait for the join is whether it can be tested element-wise or not.
immutable DispatchableJoin; end
isdispatchable(x) = DispatchableJoin()

"""
    testjoin(a, b, [jointype])

Determine if two pieces of raw data are considered equal for the purpose of a
join, and return true or false. By default, uses jointype `InnerJoin()`, but
also other popular joins are implemented, including `LeftJoin()`, `RightJoin()` and
`OuterJoin()`.

Users may define their own jointype as a subtype of `AbstractJoin` as well as the
corresponding `testjoin()` and `dojoin()` functions to implement
`join(table1, table2, jointype)`.
"""
testjoin(a, b, jointype::InnerJoin = InnerJoin()) = a == b

"""
    dojoin(a, b, [jointype])

Given two pieces of raw data that are previously tested as being joinable, and
return the correct data. By default, uses jointype `InnerJoin()`, but
also other popular joins are implemented, including `LeftJoin()`, `RightJoin()` and
`OuterJoin()`.

Users may define their own jointype as a subtype of `AbstractJoin` as well as the
corresponding `testjoin()` and `dojoin()` functions to implement
join(table1, table2, jointype).
"""
dojoin(a, b, jointype::InnerJoin = InnerJoin()) = a

"""
    join(a::Row, b::Row, [jointype])

Determine if two rows data are considered equal for the purpose of a
join, and returns a nullable of the joined rows. By default, this will be
performed individually on each element of the row.

Users may define their own jointype as a subtype of `AbstractJoin` and the
corresponding method

  `join{Names}(a::Row{Names},b::Row{Names}, jointype)`

for joining (the overlapping potion of Rows) with more complex behaviour.
"""
Base.join{Names}(a::Row{Names}, b::Row{Names}, jointype::AbstractJoin = InnerJoin()) = join(a, b, jointype, isdispatchable(jointype))

@generated function Base.join{Names, Types}(a::Row{Names, Types}, b::Row{Names}, jointype::AbstractJoin = InnerJoin(), ::DispatchableJoin = DispatchableJoin())
    expr1 = Vector{Expr}(length(Names))
    for i = 1:length(Names)
        expr1[i] = :(if !testjoin(a.data[$i], b.data[$i], jointype); return Nullable{Row{Names,Types}}(); end)
    end
    expr1 = Expr(:block, expr1...)

    expr2 = Vector{Expr}(length(Names))
    for i = 1:length(Names)
        expr2[i] = :( dojoin(a.data[$i], b.data[$i], jointype) )
    end
    expr2 = Expr(:tuple, expr2...)
    :($expr1; Nullable{Row{Names,Types}}(Row{Names,Types}($expr2)) )
end

@generated function Base.join(a::Row, b::Row, jointype::AbstractJoin = InnerJoin())
    idx_a = (setdiff(names(a), names(b))...)
    idx_ab = (intersect(names(a), names(b))...)
    idx_b = (setdiff(names(b), names(a))...)
    idx_out = (idx_ab..., idx_a..., idx_b...)

    t_a  = eltypes(a).parameters[columnindex(names(a), idx_a)]
    t_ab = eltypes(a).parameters[columnindex(names(a), idx_ab)]
    t_b  = eltypes(b).parameters[columnindex(names(b), idx_b)]

    types = Expr(:curly, :Tuple, t_ab..., t_a..., t_b...)

    quote
        tmp = join(a[Val{$idx_ab}], b[Val{$idx_ab}], jointype)
        if isnull(tmp)
            return Nullable{Row{$idx_out, types}}()
        else
            return Nullable{Row{$idx_out, types}}(hcat(get(tmp), a[Val{$idx_a}], b[Val{$idx_b}]))
        end
    end
end

"""
    join(a::Table, b::Table, [jointype])

Perform the relational join operation on two tables, generating a new `Table`.

A typical join operation searches for `Row`s which are equal in their overalapping
fields, and pushes a concatenation of the matches to the output table. However,
the join test may be something other than equality. Built-in joins include:

  `InnerJoin()` - equality tested with `==` [default]
  ...

Users may define their own jointype as a subtype of `TypedTables.AbstractJoin` and the
corresponding method for either (a) `TypedTables.testjoin()` and `TypedTables.dojoin()` for
individual data elements, or (b) `join()` on `Row`s for more complex behaviour (e.g.
multi-cell tests).
"""
@generated function Base.join(a::Table, b::Table, jointype::AbstractJoin = InnerJoin())
    NewNames = ...
    NewTypes = ...
    quote
        out = Table{$NewNames, $NewTypes}()
        for i = 1:length(a)
            for j = 1:length(b)
                row = join(a[i], b[j], jointype)
                if !isnull(row)
                    push!(out, get(row))
                end
            end
        end
        return out
    end
end
=#


text:
import Knet
using Knet.Ops20: softmax
using FileIO, ImageCore, ImageMagick, ImageTransformations, PyCall, Artifacts, Base.Threads

function torch_preprocess(file::String; atype=Knet.atype, resolution=224, format="whcn")
    transforms = pyimport("torchvision.transforms")
    torchvision_transform = transforms.Compose([
        transforms.Resize(resolution * 8 ÷ 7),
        transforms.CenterCrop(resolution),
        transforms.ToTensor(),
        transforms.Normalize(mean=[0.485, 0.456, 0.406], std=[0.229, 0.224, 0.225])
    ])
    PIL = pyimport("PIL")
    img = nothing
    @pywith PIL.Image.open(file) as f begin
        img = torchvision_transform(f.convert("RGB"))
    end
    img = img.numpy()
    img = reshape(img, (1, size(img)...)) # n,c,h,w=1,3,224,224
    if format != "nchw"
        perm = findfirst.(collect(format), "nchw")
        img = permutedims(img, (perm...,))
    end
    return atype(img)
end

function keras_preprocess_buggy(file::String; atype=Knet.atype, resolution=224, format="whcn")
    tf = pyimport("tensorflow")
    img = tf.io.read_file(file)
    img = tf.image.decode_jpeg(img, channels=3)
    img = tf.cast(img, tf.float32)
    img = tf.image.resize_with_crop_or_pad(img, resolution, resolution)
    img = tf.keras.applications.imagenet_utils.preprocess_input(img, mode="tf")
    img = img.numpy() # hwc
    img = reshape(img, (1, size(img)...)) # nhwc
    if format != "nhwc"
        perm = findfirst.(collect(format), "nhwc")
        img = permutedims(img, (perm...,))
    end
    return atype(img)
end

function imagenet_preprocess(file::String; o...)
    img = occursin(r"^http", file) ? mktemp() do fn,io
        load(download(file,fn))
    end : load(file)            # size=(h,w)
    imagenet_preprocess(img; o...)
end

function imagenet_preprocess(img::Matrix{<:Gray}; o...)
    imagenet_preprocess(RGB.(img); o...)
end

function imagenet_preprocess(img::Matrix{<:RGB}; normalize=nothing, resolution=224, format="whcn", atype=Knet.atype)
    if normalize === nothing; normalize = identity; end
    minsize = resolution * 8 ÷ 7                          # 224 => 256
    img = imresize(img, ratio=minsize/minimum(size(img))) # min(h,w)=256
    hcenter,vcenter = size(img) .>> 1                     # h÷2, w÷2
    half = resolution ÷ 2
    img = img[hcenter-half+1:hcenter+half, vcenter-half+1:vcenter+half] # hw=224,224
    img = Float32.(channelview(img))                                    # chw=3,224,224
    img = normalize(img)
    img = reshape(img, (1, size(img)...)) # nchw=1,3,224,224
    if format != "nchw"
        perm = findfirst.(collect(format), "nchw")
        img = permutedims(img, (perm...,))
    end
    return atype(img)
end

# fallback
imagenet_preprocess(x; atype=Knet.atype, o...) = atype(x)

# ImageNet meta-information
function imagenet_labels()
    global _imagenet_labels
    if !@isdefined(_imagenet_labels)
        _imagenet_labels = [ replace(x, r"\S+ ([^,]+).*"=>s"\1") for x in
                             readlines(joinpath(artifact"imagenet_labels","LOC_synset_mapping.txt")) ]
    end
    _imagenet_labels
end

function imagenet_synsets()
    global _imagenet_synsets
    if !@isdefined(_imagenet_synsets)
        _imagenet_synsets = [ split(s)[1] for s in
                              readlines(joinpath(artifact"imagenet_labels", "LOC_synset_mapping.txt")) ]
    end
    _imagenet_synsets
end

function imagenet_val()
    global _imagenet_val
    if !@isdefined(_imagenet_val)
        synset2index = Dict(s=>i for (i,s) in enumerate(imagenet_synsets()))
        _imagenet_val = Dict(x=>synset2index[y] for (x,y) in (z->split(z,[',',' '])[1:2]).(
            readlines(joinpath(artifact"imagenet_labels", "LOC_val_solution.csv"))[2:end]))
    end
    _imagenet_val
end

# Human readable predictions from tensors, images, files
function imagenet_predict(model, img="ILSVRC2012_val_00000001.JPEG"; preprocess=identity, apply_softmax=true)
    img = preprocess(img)
    cls = model(img)
    if cls isa PyObject; cls = cls.cpu().numpy(); end
    cls = convert(Array, vec(cls))
    cls = apply_softmax ? softmax(cls) : cls
    idx = sortperm(cls, rev=true)
    [ idx cls[idx] imagenet_labels()[idx] ]
end

# Recursive walk and predict each image in a directory
function imagenet_walkdir(
    model, dir;
    n=typemax(Int),
    b=32,
    preprocess=model[1],
    format="whcn",
    batchinput = identity,
    batchoutput = (x->convert(Array,x)),
    usethreads = false,
    o...
)
    files = []
    for (root, dirs, fs) in walkdir(dir)
        for f in fs
            push!(files, joinpath(root, f))
            length(files) > n && break
        end
        length(files) > n && break
    end
    n = min(n, length(files))
    images = Array{Any}(undef, b)
    preds = []
    for i in Knet.progress(1:b:n)
        j = min(n, i+b-1)
        if usethreads
            @threads for k in 0:j-i
                images[1+k] = preprocess(files[i+k])
            end
        else
            for k in 0:j-i
                images[1+k] = preprocess(files[i+k])
            end
        end
        batch = cat(images[1:j-i+1]...; dims=findfirst('n', format))
        p = (batch |> batchinput |> model |> batchoutput)
        append!(preds, vec((i->i[1]).(argmax(p; dims=1))))
    end
    [ preds imagenet_labels()[preds] files[1:n] ]
end

# Top-1 error in validation set
function imagenet_top1(model, valdir; o...)
    pred = imagenet_walkdir(model, valdir; o...)
    gold = imagenet_val()
    error, count = 0, size(pred, 1)
    for i in 1:count
        image = match(r"ILSVRC2012_val_\d+", pred[i,3]).match
        if pred[i,1] != gold[image]
            error += 1
        end
    end
    error = error / count
    (accuracy = 1-error, error = error)
end

# torch model top-1 example:
# r18a = models.resnet18(pretrained=true).eval()
# imagenet_top1(r18a, val; n=1000, preprocess=torch_preprocess, format="nchw", batchinput=torch.tensor, batchoutput=(x->permutedims(x.detach().cpu().numpy())), atype=Array{Float32}, usethreads=false)


text:
# FIXME this thing eats up macro calls!
"""
Remove LineNumberNode from quoted blocks of code
"""
rmlines(e::Expr) = Expr(e.head, map(rmlines, filter(x -> !(x isa LineNumberNode), e.args))...)
rmlines(a) = a

"""
HARD FIX of n-arity of operators in `Expr` trees
"""
function binarize!(e, op::Symbol)
    f(e) = if (isexpr(e, :call) && e.args[1] == op && length(e.args) > 3)
        foldl((x,y) -> Expr(:call, op, x, y), e.args[2:end])
    else e end

    df_walk!(f, e)
end

"""
Binarize n-ary operators (`+` and `*`) and call [`rmlines`](@ref)
"""
cleanast(ex) = rmlines(ex) |>
    x -> binarize!(x, :(+)) |>
    x -> binarize!(x, :(*))

interp_dol(ex::Expr, mod::Module) =
    Meta.isexpr(ex, :$) ? mod.eval(ex.args[1]) : ex
interp_dol(any, mod::Module) = any

function interpolate_dollar(ex, mod::Module)
    df_walk(interp_dol, ex, mod)
end


text:
export score_driven_recursion, fitted_mean, simulate_recursion

"""
    score_driven_recursion(sd_model::SDM, observations::Vector{T}) where T

start with the stationary params for a
"""
function score_driven_recursion(gas::Model{D, T}, observations::Vector{T};
                                 initial_params::Matrix{T} = stationary_initial_params(gas)) where {D, T}
    @assert gas.scaling in SCALINGS
    # Allocations
    n = length(observations)
    n_params = num_params(D)
    param = Matrix{T}(undef, n + 1, n_params)
    param_tilde = Matrix{T}(undef, n + 1, n_params)
    scores_tilde = Matrix{T}(undef, n, n_params)
    aux = AuxiliaryLinAlg{T}(n_params)

    # Query the biggest lag
    biggest_lag = number_of_lags(gas)

    # initial_values
    @inbounds for i in 1:biggest_lag
        for p in 1:n_params
            param[i, p] = initial_params[i, p]
        end
        link!(param_tilde, D, param, i)
        score_tilde!(scores_tilde, observations[i], D, param, aux, gas.scaling, i)
    end

    update_param_tilde!(param_tilde, gas.ω, gas.A, gas.B, scores_tilde, biggest_lag)

    @inbounds for i in biggest_lag + 1:n
        univariate_score_driven_update!(param, param_tilde, scores_tilde, observations[i], aux, gas, i)
    end
    update_param!(param, param_tilde, D, n + 1)

    return param
end

function simulate_recursion(gas::Model{D, T}, n::Int;
                            initial_params::Matrix{T} = stationary_initial_params(gas),
                            update::Function = sample_observation) where {D, T}
    # Allocations
    serie = zeros(n)
    n_params = num_params(D)
    param = Matrix{T}(undef, n, n_params)
    param_tilde = Matrix{T}(undef, n, n_params)
    scores_tilde = Matrix{T}(undef, n, n_params)

    aux = AuxiliaryLinAlg{T}(n_params)

    # Auxiliary Allocation
    param_dist = zeros(T, 1, n_params)

    biggest_lag = number_of_lags(gas)

    # initial_values
    for t in 1:biggest_lag
        for p in 1:n_params
            @inbounds param[t, p] = initial_params[t, p]
        end
        link!(param_tilde, D, param, t)
        # Sample
        updated_dist = update_dist(D, param, t)
        serie[t] = update(updated_dist)
        score_tilde!(scores_tilde, serie[t], D, param, aux, gas.scaling, t)
    end

    update_param_tilde!(param_tilde, gas.ω, gas.A, gas.B, scores_tilde, biggest_lag)
    unlink!(param, D, param_tilde, biggest_lag + 1)
    updated_dist = update_dist(D, param, biggest_lag + 1)
    serie[biggest_lag + 1] = update(updated_dist)

    for i in biggest_lag + 1:n-1
        # update step
        univariate_score_driven_update!(param, param_tilde, scores_tilde, serie[i], aux, gas, i)
        # Sample from the updated distribution
        unlink!(param, D, param_tilde, i + 1)
        updated_dist = update_dist(D, param, i + 1)
        serie[i + 1] = update(updated_dist)
    end
    update_param!(param, param_tilde, D, n)

    return serie, param
end

function univariate_score_driven_update!(param::Matrix{T}, param_tilde::Matrix{T},
                                         scores_tilde::Matrix{T},
                                         observation::T, aux::AuxiliaryLinAlg{T},
                                         gas::Model{D, T}, i::Int) where {D <: Distribution, T <: AbstractFloat}
    # update param
    update_param!(param, param_tilde, D, i)
    # evaluate score
    score_tilde!(scores_tilde, observation, D, param, aux, gas.scaling, i)
    # update param_tilde
    update_param_tilde!(param_tilde, gas.ω, gas.A, gas.B, scores_tilde, i)
    return
end

function update_param!(param::Matrix{T}, param_tilde::Matrix{T}, D::Type{<:Distribution}, i::Int) where T
    unlink!(param, D, param_tilde, i)
    # Some treatments
    NaN2zero!(param, i)
    big_threshold!(param, BIG_NUM, i)
    small_threshold!(param, SMALL_NUM, i)
    return
end

function update_param_tilde!(param_tilde::Matrix{T}, ω::Vector{T}, A::Dict{Int, Matrix{T}},
                             B::Dict{Int, Matrix{T}}, scores_tilde::Matrix{T}, i::Int) where T
    for p in eachindex(ω)
        @inbounds param_tilde[i + 1, p] = ω[p]
    end
    for (lag, mat) in A
        for p in axes(mat, 1)
            @inbounds param_tilde[i + 1, p] += mat[p, p] * scores_tilde[i - lag + 1, p]
        end
    end
    for (lag, mat) in B
        for p in axes(mat, 1)
            @inbounds param_tilde[i + 1, p] += mat[p, p] * param_tilde[i - lag + 1, p]
        end
    end
    return
end

"""
    fitted_mean(gas::Model{D, T}, observations::Vector{T}, initial_params::Vector{T}) where {D, T}

return the fitted mean.. #TODO
"""
function fitted_mean(gas::Model{D, T}, observations::Vector{T};
                     initial_params::Matrix{T} = stationary_initial_params(gas)) where {D, T}
    params_fitted = score_driven_recursion(gas, observations; initial_params = initial_params)

    # Discard last step of the update
    params_fitted = params_fitted[1:end-1, :]

    n = size(params_fitted, 1)
    fitted_mean = Vector{T}(undef, n)

    for t in axes(params_fitted, 1)
        sdm_dist = update_dist(D, params_fitted, t)
        fitted_mean[t] = mean(sdm_dist)
    end

    return fitted_mean
end


text:
using Documenter
using EzXML: AttributeNode, ElementNode, HTMLDocumentNode, TextNode, link!, prettyprint
using Kelpie
using Suppressor
using Test

# Set up doctests
DocMeta.setdocmeta!(Kelpie, :DocTestSetup, :(using Kelpie); recursive=true)

prettystring(xml) = @capture_out prettyprint(xml)

# Set up the Bootstrap 5.1.3 starter template in EzXML nodes
function bootstrap_starter()
    # Create the nodes
    doc_node = HTMLDocumentNode("about:legacy-compat", nothing)
    html_node = ElementNode("html")
    head_node = ElementNode("head")
    meta_charset_node = ElementNode("meta")
    meta_viewport_node = ElementNode("meta")
    css_link_node = ElementNode("link")
    title_node = ElementNode("title")
    body_node = ElementNode("body")
    h1_node = ElementNode("h1")
    script_node = ElementNode("script")

    # Add attributes
    link!(html_node, AttributeNode("lang", "en"))
    link!(meta_charset_node, AttributeNode("charset", "utf-8"))
    link!(meta_viewport_node, AttributeNode("name", "viewport"))
    link!(
        meta_viewport_node, AttributeNode("content", "width=device-width, initial-scale=1")
    )
    link!(
        css_link_node,
        AttributeNode(
            "href",
            "https://cdn.jsdelivr.net/npm/bootstrap@5.1.3/dist/css/bootstrap.min.css",
        ),
    )
    link!(css_link_node, AttributeNode("rel", "stylesheet"))
    link!(
        css_link_node,
        AttributeNode(
            "integrity",
            "sha384-1BmE4kWBq78iYhFldvKuhfTAU6auU8tT94WrHftjDbrCEXSU1oBoqyl2QvZ6jIW3",
        ),
    )
    link!(css_link_node, AttributeNode("crossorigin", "anonymous"))
    link!(
        script_node,
        AttributeNode(
            "src",
            "https://cdn.jsdelivr.net/npm/bootstrap@5.1.3/dist/js/bootstrap.bundle.min.js",
        ),
    )
    link!(
        script_node,
        AttributeNode(
            "integrity",
            "sha384-ka7Sk0Gln4gmtz2MlQnikT1wXgYsOg+OMhuP+IlRH9sENBO0LRn5q+8nbTov4+1p",
        ),
    )
    link!(script_node, AttributeNode("crossorigin", "anonymous"))

    # Add content
    link!(title_node, TextNode("Hello, world!"))
    link!(h1_node, TextNode("Hello, world!"))
    link!(script_node, TextNode(""))

    # Link the tree together
    link!(head_node, meta_charset_node)
    link!(head_node, meta_viewport_node)
    link!(head_node, css_link_node)
    link!(head_node, title_node)
    link!(body_node, h1_node)
    link!(body_node, script_node)
    link!(html_node, head_node)
    link!(html_node, body_node)
    link!(doc_node, html_node)

    return prettystring(doc_node)
end #function

# Set up the Zurb Foundation 6.7.4 starter template in EzXML nodes
function foundation_template()
    # Create the nodes
    doc_node = HTMLDocumentNode("about:legacy-compat", nothing)
    html_node = ElementNode("html")
    head_node = ElementNode("head")
    meta_charset_node = ElementNode("meta")
    meta_http_equiv_node = ElementNode("meta")
    meta_viewport_node = ElementNode("meta")
    title_node = ElementNode("title")
    css_link_node = ElementNode("link")
    body_node = ElementNode("body")
    h1_node = ElementNode("h1")
    script_jquery_node = ElementNode("script")
    script_whatinput_node = ElementNode("script")
    script_foundation_node = ElementNode("script")
    script_launch_node = ElementNode("script")

    # Add attributes
    link!(html_node, AttributeNode("class", "no-js"))
    link!(html_node, AttributeNode("lang", "en"))
    link!(meta_charset_node, AttributeNode("charset", "utf-8"))
    link!(meta_http_equiv_node, AttributeNode("http-equiv", "x-ua-compatible"))
    link!(meta_http_equiv_node, AttributeNode("content", "ie=edge"))
    link!(meta_viewport_node, AttributeNode("name", "viewport"))
    link!(
        meta_viewport_node,
        AttributeNode("content", "width=device-width, initial-scale=1.0"),
    )
    link!(css_link_node, AttributeNode("rel", "stylesheet"))
    link!(
        css_link_node,
        AttributeNode(
            "href",
            "https://cdn.jsdelivr.net/npm/foundation-sites@6.7.4/dist/css/foundation.min.css",
        ),
    )
    link!(css_link_node, AttributeNode("crossorigin", "anonymous"))
    link!(
        script_jquery_node,
        AttributeNode(
            "src", "https://cdn.jsdelivr.net/npm/jquery@3.6.0/dist/jquery.min.js"
        ),
    )
    link!(script_jquery_node, AttributeNode("crossorigin", "anonymous"))
    link!(
        script_whatinput_node,
        AttributeNode(
            "src", "https://cdn.jsdelivr.net/npm/what-input@5.2.10/dist/what-input.min.js"
        ),
    )
    link!(script_whatinput_node, AttributeNode("crossorigin", "anonymous"))
    link!(
        script_foundation_node,
        AttributeNode(
            "src",
            "https://cdn.jsdelivr.net/npm/foundation-sites@6.7.4/dist/js/foundation.min.js",
        ),
    )
    link!(script_foundation_node, AttributeNode("crossorigin", "anonymous"))

    # Add content
    link!(title_node, TextNode("Foundation Starter Template"))
    link!(h1_node, TextNode("Hello, world!"))
    link!(script_jquery_node, TextNode(""))
    link!(script_whatinput_node, TextNode(""))
    link!(script_foundation_node, TextNode(""))
    link!(script_launch_node, TextNode("\$(document).foundation();"))

    # Link the tree together
    link!(head_node, meta_charset_node)
    link!(head_node, meta_http_equiv_node)
    link!(head_node, meta_viewport_node)
    link!(head_node, title_node)
    link!(head_node, css_link_node)
    link!(body_node, h1_node)
    link!(body_node, script_jquery_node)
    link!(body_node, script_whatinput_node)
    link!(body_node, script_foundation_node)
    link!(body_node, script_launch_node)
    link!(html_node, head_node)
    link!(html_node, body_node)
    link!(doc_node, html_node)

    return prettystring(doc_node)
end #function

@testset "Kelpie.jl" begin
    @testset "Doctests" begin
        doctest(Kelpie)
    end #testset

    @testset "Frameworks" begin
        bootstrap_kelpie = html(
            head(
                meta(; charset="utf-8"),
                meta(; name="viewport", content="width=device-width, initial-scale=1"),
                link(;
                    href="https://cdn.jsdelivr.net/npm/bootstrap@5.1.3/dist/css/bootstrap.min.css",
                    rel="stylesheet",
                    integrity="sha384-1BmE4kWBq78iYhFldvKuhfTAU6auU8tT94WrHftjDbrCEXSU1oBoqyl2QvZ6jIW3",
                    crossorigin="anonymous",
                ),
                title("Hello, world!"),
            ),
            body(
                h1("Hello, world!"),
                script(
                    "";
                    src="https://cdn.jsdelivr.net/npm/bootstrap@5.1.3/dist/js/bootstrap.bundle.min.js",
                    integrity="sha384-ka7Sk0Gln4gmtz2MlQnikT1wXgYsOg+OMhuP+IlRH9sENBO0LRn5q+8nbTov4+1p",
                    crossorigin="anonymous",
                ),
            );
            lang="en",
        )
        foundation_kelpie = html(
            head(
                meta(; charset="utf-8"),
                meta(; http_equiv="x-ua-compatible", content="ie=edge"),
                meta(; name="viewport", content="width=device-width, initial-scale=1.0"),
                title("Foundation Starter Template"),
                link(;
                    rel="stylesheet",
                    href="https://cdn.jsdelivr.net/npm/foundation-sites@6.7.4/dist/css/foundation.min.css",
                    crossorigin="anonymous",
                ),
            ),
            body(
                h1("Hello, world!"),
                script(
                    "";
                    src="https://cdn.jsdelivr.net/npm/jquery@3.6.0/dist/jquery.min.js",
                    crossorigin="anonymous",
                ),
                script(
                    "";
                    src="https://cdn.jsdelivr.net/npm/what-input@5.2.10/dist/what-input.min.js",
                    crossorigin="anonymous",
                ),
                script(
                    "";
                    src="https://cdn.jsdelivr.net/npm/foundation-sites@6.7.4/dist/js/foundation.min.js",
                    crossorigin="anonymous",
                ),
                script("\$(document).foundation();"),
            );
            class="no-js",
            lang="en",
        )

        @test prettystring(bootstrap_kelpie) == bootstrap_starter()
        @test prettystring(foundation_kelpie) == foundation_template()
    end #testset

    @testset "Behaviors" begin
        # Test that attributes with underscores become attributes with dashes
        @test occursin("", prettystring(i(; data_tooltip="foo")))

        # Test that edge case functions can slurp their arguments
        @test occursin(
            """

foo

bar

""",
            prettystring(html_div(p("foo"), p("bar"))),
        )
    end #testset
end #testset


text:
using SparseInverseProblems
import SparseInverseProblems: lmo, phi, solveFiniteDimProblem, localDescent,
getStartingPoint, parameterBounds, computeGradient
using SparseInverseProblems.Util
const min_theta = 0.001
const max_theta = pi
const min_r = 1E-3
const max_r = 0.999
const max_x = 1.0
const max_b = 1.0
struct LinearSysID <: BoxConstrainedDifferentiableModel
  u :: Vector{Float64}  #u is the input time series
  grid
  normalization :: Function # (r,theta) -> n(r,theta)
  dndr :: Function # (r,theta) -> dn/dr
  dndt :: Function # (r,theta) -> dn/dt
  LinearSysID(u, gr,gt,n,dndr,dndt) = new(u, grid(gr,gt,u),n,dndr,dndt)
end

function getStartingPoint(model :: LinearSysID, v :: Vector{Float64})
  (gridR, gridT, l_x, l_b)  = model.grid

  scores = zeros(length(gridR), length(gridT))
  for i = 1:length(gridR)
    for j = 1:length(gridT)
      r,theta = gridR[i], gridT[j]
      nf = model.normalization(r,theta)
      scores[i,j] = nf*(max_x*sumabs(v'*l_x[i,j]) + max_b*sumabs(v'*l_b[i,j]))
    end
  end
  ind = indmax(scores)
  (rind,tind) = ind2sub((length(gridR), length(gridT)),ind)
  best_x = -max_x * sign(v'*l_x[ind])
  best_b = -max_b * sign(v'*l_b[ind])
  gridP = [gridR[rind]; gridT[tind]; vec(best_x); vec(best_b)]

  return gridP
end

parameterBounds(model :: LinearSysID) =
  ([min_r, min_theta,-max_x,-max_x, -max_b, -max_b],
  [max_r, max_theta,max_x,max_x, max_b, max_b])

function computeGradient(model :: LinearSysID, weights :: Vector{Float64},
  thetas :: Matrix{Float64}, v :: Vector{Float64})
  gradient = zeros(thetas)
  for i = 1:size(thetas,2)
    gradient[:,i] = weights[i]*jacobian(model, vec(thetas[:,i]))* v
  end
  return gradient
end

function phi(sim :: LinearSysID,parameters :: Matrix{Float64}, weights :: Vector{Float64})
  return runSim(parameters,sim.u,sim.normalization)*weights
end

function runSim(atoms,u,n)
  nAtoms = size(atoms,2)
  results = zeros(length(u),nAtoms)
  for i = 1:nAtoms
    w = vec(atoms[:,i])
    r,t,x,b = unpackParameters(w)
    l_x,l_b = getLinearOps(r,t,u)
    results[:,i] = (l_x*x + l_b*b)*n(r,t)
  end
  return results
end

function getLinearOps(r,theta,u; l_x :: Matrix{Float64} = zeros(length(u),2),l_b :: Matrix{Float64} = zeros(length(u),2) )
  n = length(u)
  rt = r.^(1:n)
  ttheta = (1:n)*theta
  ct = cos(ttheta)
  st = sin(ttheta)
  @inbounds for t = 1:n
    l_x[t,1] = rt[t]*ct[t]
    l_x[t,2] = -rt[t]*st[t]
     @inbounds for tau = 1:t-2 #u[0] = 0
      flipped_t = t-1-tau
      c = u[tau]*rt[flipped_t]
      l_b[t,1] += c*ct[flipped_t]
      l_b[t,2] -= c*st[flipped_t]
    end
  end
  @inbounds for t = 2:n
    l_b[t,1] += u[t-1] # flipped_t = 0
  end
  return (l_x, l_b)
end

function jacobian(s :: LinearSysID, w :: Vector{Float64}; j_r :: Vector{Float64} = zeros(length(s.u)), j_theta :: Vector{Float64} = zeros(length(s.u)), l_x :: Matrix{Float64} = zeros(length(s.u),2), l_b :: Matrix{Float64} = zeros(length(s.u),2))
  r,theta,x,b = unpackParameters(w)
  n = length(s.u)
  u :: Vector{Float64} = s.u
  rt = r.^(1:n)
  ttheta = (1:n)*theta
  ct = cos(ttheta)
  st = sin(ttheta)

  @inbounds for t = 1:n
    l_x[t,1] = rt[t]*ct[t]
    l_x[t,2] = -rt[t]*st[t]
    j_theta[t] -= t*rt[t]*st[t]*x[1]
    j_theta[t] -= t*rt[t]*ct[t]*x[2]
     for tau = 1:t-2 #u[0] = 0
      flipped_t = t-1-tau
      c = u[tau]*rt[flipped_t]
      l_b[t,1] += c*ct[flipped_t]
      l_b[t,2] -= c*st[flipped_t]
      j_theta[t] -= u[tau]*flipped_t*rt[flipped_t]*st[flipped_t]*b[1]
      j_theta[t] -= u[tau]*flipped_t*rt[flipped_t]*ct[flipped_t]*b[2]
    end
     for tau = 1:t-3 #u[0] = 0
      flipped_t = t-1-tau
      j_r[t] += u[tau]*flipped_t*rt[flipped_t-1]*ct[flipped_t]*b[1]
      j_r[t] -= u[tau]*flipped_t*rt[flipped_t-1]*st[flipped_t]*b[2]
    end
  end
  @inbounds for t = 2:n
    j_r[t] += t*rt[t-1]*ct[t]*x[1]
    j_r[t] -= t*rt[t-1]*st[t]*x[2]
    l_b[t,1] += u[t-1] # flipped_t = 0
  end
  #normalization!
  nf = s.normalization(r,theta)
  o = l_x*x + l_b*b
  j_r = j_r*nf + o*s.dndr(r,theta)
  j_theta = j_theta*nf + o*s.dndt(r,theta)
  return [j_r j_theta l_x*nf l_b*nf]'
end

function grid(nr,nt,u)
  gridR = max_r * sqrt(linspace(min_r,1.0,nr))
  gridT = linspace(min_theta,max_theta,nt)
  l_x = Array(Matrix{Float64},nr,nt)
  l_b = Array(Matrix{Float64},nr,nt)
  for i = 1:nr, j = 1:nt
    r = gridR[i]
    theta = gridT[j]
    (l_x[i,j],l_b[i,j]) = getLinearOps(r,theta,u)
  end
  return (gridR, gridT, l_x,l_b)
end

unpackParameters(p) = (p[1],p[2],p[3:4],p[5:6]) #r, t, x_0, B

function solveFiniteDimProblem(sim :: LinearSysID, lossFn :: LSLoss, thetas :: Matrix{Float64}, y :: Vector{Float64}, tau :: Float64)
  A = runSim(thetas,sim.u,sim.normalization)
  return lasso_kernel(A'*A,A'*y,tau)
end


text:
"""
Compute a scattered grid of M points randomly distributed in `Ω`, using
the uniform probability measure on `Ω`.
"""
function randomgrid(Ω::Domain, M::Int)
    points = [randompoint(Ω, boundingbox(Ω)) for m in 1:M]
    ScatteredGrid(sort(points))
end

"Generate a single random point inside the given box, with `eltype` `T`."
function randompoint(dom::ProductDomain)
    convert(eltype(dom),map(randompoint,elements(dom)))
end

# Don't return an SVector in 1d, just a value
function randompoint(dom::AbstractInterval)
    T = float(eltype(dom))
    val = rand(T)
    convert(T,val * infimum(dom) + (1-val) * supremum(dom))
end

"""
Generate a single random point inside the given domain, with `eltype` `T`.
Random points are generated inside the given box, until one is inside the domain.
"""
function randompoint(Ω::Domain, box = boundingbox(Ω))
    local point
    indomain = false
    while ~indomain
        point = randompoint(box)
        indomain = point ∈ Ω
    end
    point
end