text:
JSONRPC.parse_params(::Type{Val{Symbol("\$/cancelRequest")}}, params) = CancelParams(params)
function process(r::JSONRPC.Request{Val{Symbol("\$/cancelRequest")},CancelParams}, server) end

JSONRPC.parse_params(::Type{Val{Symbol("\$/setTraceNotification")}}, params) = params
function process(r::JSONRPC.Request{Val{Symbol("\$/setTraceNotification")},Dict{String,Any}}, server) end

function JSONRPC.parse_params(::Type{Val{Symbol("julia/lint-package")}}, params) end
function process(r::JSONRPC.Request{Val{Symbol("julia/lint-package")},Nothing}, server) end

function JSONRPC.parse_params(::Type{Val{Symbol("julia/reload-modules")}}, params) end
function process(r::JSONRPC.Request{Val{Symbol("julia/reload-modules")},Nothing}, server) end

JSONRPC.parse_params(::Type{Val{Symbol("julia/toggleFileLint")}}, params) = params
function process(r::JSONRPC.Request{Val{Symbol("julia/toggleFileLint")}}, server) end

function JSONRPC.parse_params(::Type{Val{Symbol("julia/toggle-log")}}, params) end
function process(r::JSONRPC.Request{Val{Symbol("julia/toggle-log")},Nothing}, server)
    server.debug_mode = !server.debug_mode
end

JSONRPC.parse_params(::Type{Val{Symbol("julia/getCurrentBlockRange")}}, params) = TextDocumentPositionParams(params)
function process(r::JSONRPC.Request{Val{Symbol("julia/getCurrentBlockRange")},TextDocumentPositionParams}, server)
    tdpp = r.params
    doc = getdocument(server, URI2(tdpp.textDocument.uri))
    offset = get_offset(doc, tdpp.position)
    x = getcst(doc)
    loc = 0
    p1, p2, p3 = 0, x.span, x.fullspan
    if typof(x) === CSTParser.FileH
        (offset > x.fullspan || x.args === nothing) && return Position(get_position_at(doc, p1)...), Position(get_position_at(doc, p2)...), Position(get_position_at(doc, p3)...)
        for a in x.args
            if loc <= offset < loc + a.fullspan
                if typof(a) === CSTParser.ModuleH
                    if loc + a.args[1].fullspan + a.args[2].fullspan < offset < loc + a.args[1].fullspan + a.args[2].fullspan + a.args[3].fullspan
                        loc0 = loc +  a.args[1].fullspan + a.args[2].fullspan
                        loc += a.args[1].fullspan + a.args[2].fullspan
                        for b in a.args[3].args
                            if loc <= offset < loc + b.fullspan
                                p1, p2, p3 = loc, loc + b.span, loc + b.fullspan
                                break
                            end
                            loc += b.fullspan
                        end
                    else
                        p1, p2, p3 = loc, loc + a.span, loc + a.fullspan
                    end
                elseif typof(a) === CSTParser.TopLevel
                    p1, p2, p3 = loc, loc + a.span, loc + a.fullspan
                    for b in a.args
                        if loc <= offset < loc + b.fullspan
                            p1, p2, p3 = loc, loc + b.span, loc + b.fullspan
                        end
                        loc += b.fullspan
                    end
                else
                    p1, p2, p3 = loc, loc + a.span, loc + a.fullspan
                end
            end
            loc += a.fullspan
        end
    end
    return Position(get_position_at(doc, p1)...), Position(get_position_at(doc, p2)...), Position(get_position_at(doc, p3)...)
end

JSONRPC.parse_params(::Type{Val{Symbol("julia/activateenvironment")}}, params) = params
function process(r::JSONRPC.Request{Val{Symbol("julia/activateenvironment")}}, server)
    server.env_path = r.params

    trigger_symbolstore_reload(server)
end


text:
module IfSimpleLoadBearing

include("./../../constraints/geometric/constants.jl")
using JuMP
using .Constants

export if_simple_lb

function if_simple_lb(m)

    #Basically just check whether box j is located right above box i

    # Zj >= Zi + hobi
    @variable(m, if1_zj_geq_zi_add_hobi_ij[i=1:length(boxes), j=1:length(boxes); i!=j], Bin)
    @constraint(m, c_if1_zj_geq_zi_add_hobi_ij1[i=1:length(boxes), j=1:length(boxes); i!=j],
        H * (1 - m[:if1_zj_geq_zi_add_hobi_ij][i, j]) >=
        (m[:flbz_front_left_bot_z][i] + m[:hob_height_of_box][i]) - m[:flbz_front_left_bot_z][j]
    )

    @constraint(m, c_if1_zj_geq_zi_add_hobi_ij2[i=1:length(boxes), j=1:length(boxes); i!=j],
        H * m[:if1_zj_geq_zi_add_hobi_ij][i, j] >=
        m[:flbz_front_left_bot_z][j] - (m[:flbz_front_left_bot_z][i] + m[:hob_height_of_box][i]) + 0.001
    )

    # Xj < Xi + Xcoveri
    @variable(m, if2_xj_l_xi_xcoveri_ij[i=1:length(boxes), j=1:length(boxes); i!=j], Bin)
    @constraint(m, c_if2_xj_l_xi_xcoveri_ij1[i=1:length(boxes), j=1:length(boxes); i!=j],
        L * (1 - m[:if2_xj_l_xi_xcoveri_ij][i, j]) >=
        m[:flbx_front_left_bot_x][j] + 0.001 - (m[:flbx_front_left_bot_x][i] + m[:box_x_axis_cover][i])
    )

    @constraint(m, c_if2_xj_l_xi_xcoveri_ij2[i=1:length(boxes), j=1:length(boxes); i!=j],
        L * m[:if2_xj_l_xi_xcoveri_ij][i, j] >=
        (m[:flbx_front_left_bot_x][i] + m[:box_x_axis_cover][i]) - m[:flbx_front_left_bot_x][j] - 0.001
    )

    # Xj + Xcoverj < Xi
    @variable(m, if3_xj_xcoverj_g_xi_ij[i=1:length(boxes), j=1:length(boxes); i!=j], Bin)
    @constraint(m, c_if3_xj_xcoverj_g_xi_ij1[i=1:length(boxes), j=1:length(boxes); i!=j],
        L * (1 - m[:if3_xj_xcoverj_g_xi_ij][i, j]) >=
        m[:flbx_front_left_bot_x][i] + 0.001 - (m[:flbx_front_left_bot_x][j] + m[:box_x_axis_cover][j])
    )

    @constraint(m, c_if3_xj_xcoverj_g_xi_ij2[i=1:length(boxes), j=1:length(boxes); i!=j],
        L * m[:if3_xj_xcoverj_g_xi_ij][i, j] >=
        (m[:flbx_front_left_bot_x][j] + m[:box_x_axis_cover][j]) - m[:flbx_front_left_bot_x][i] - 0.001
    )

    # Yj < Yi + Ycoveri
    @variable(m, if4_yj_l_yi_ycoveri_ij[i=1:length(boxes), j=1:length(boxes); i!=j], Bin)
    @constraint(m, c_if4_yj_l_yi_ycoveri_ij1[i=1:length(boxes), j=1:length(boxes); i!=j],
        W * (1 - m[:if4_yj_l_yi_ycoveri_ij][i, j]) >=
        m[:flby_front_left_bot_y][j] + 0.001 - (m[:flby_front_left_bot_y][i] + m[:box_y_axis_cover][i])
    )

    @constraint(m, c_if4_yj_l_yi_ycoveri_ij2[i=1:length(boxes), j=1:length(boxes); i!=j],
        W * m[:if4_yj_l_yi_ycoveri_ij][i, j] >=
        (m[:flby_front_left_bot_y][i] + m[:box_y_axis_cover][i]) - m[:flby_front_left_bot_y][j] - 0.001
    )

    # Yj < Yi + Ycoveri
    @variable(m, if5_yj_ycoverj_g_yi_ij[i=1:length(boxes), j=1:length(boxes); i!=j], Bin)
    @constraint(m, c_if5_yj_ycoverj_g_yi_ij1[i=1:length(boxes), j=1:length(boxes); i!=j],
        W * (1 - m[:if5_yj_ycoverj_g_yi_ij][i, j]) >=
        m[:flby_front_left_bot_y][i] + 0.001 - (m[:flby_front_left_bot_y][j] + m[:box_y_axis_cover][j])
    )

    @constraint(m, c_if5_yj_ycoverj_g_yi_ij2[i=1:length(boxes), j=1:length(boxes); i!=j],
        W * m[:if5_yj_ycoverj_g_yi_ij][i, j] >=
        (m[:flby_front_left_bot_y][j] + m[:box_y_axis_cover][j]) - m[:flby_front_left_bot_y][i] - 0.001
    )

    #attempt at advanced AND
    @variable(m, if_simple_lb_var[i=1:length(boxes), j=1:length(boxes); i!=j], Bin)
    @constraint(m, c_if_simple_lb1[i=1:length(boxes), j=1:length(boxes); i!=j],
        m[:if_simple_lb_var][i, j] >=
        m[:if1_zj_geq_zi_add_hobi_ij][i, j] +
        m[:if2_xj_l_xi_xcoveri_ij][i, j] +
        m[:if3_xj_xcoverj_g_xi_ij][i, j] +
        m[:if4_yj_l_yi_ycoveri_ij][i, j] +
        m[:if5_yj_ycoverj_g_yi_ij][i, j] - 4
    )

    #@expression(m, if_simple_lb_objective[i=1:length(boxes), j=1:length(boxes); i!=j],
    #    (1 - m[:if_simple_lb_var][i, j])
    #)
    return m
end

end

ext: jl


text:
module HiddenMarkovModel

using Distributions
using StatsBase:sample,WeightVec

# HMM constructors
export HMM, generate
include("HMM.jl")

# Forward/Backward Algorithm and Viterbi estimation of state sequence
export forward_backward, viterbi, smoothed_forward_backward
include("forward_backward.jl")

# Fitting algorithms (under development)
export fit!
include("fit.jl")

# Discrete/toy HMM -- include("dHMM.jl")

end


text:
using Pkg
using Base: UUID

# Pkg.PackageSpec return different types in different Julia versions so...
const PkgSpec = typeof(Pkg.PackageSpec(name="dummy"))

# If we don't have `stdlib_version` from Pkg, recreate it ourselves
if !isdefined(Pkg.Types, :stdlib_version)
    function stdlib_version(uuid::Base.UUID, julia_version::Union{VersionNumber,Nothing})::Union{VersionNumber,Nothing}
        last_stdlibs = Pkg.Types.get_last_stdlibs(julia_version)
        if !(uuid in keys(last_stdlibs))
            return nothing
        end
        return last_stdlibs[uuid][2]
    end
else
    const stdlib_version = Pkg.Types.stdlib_version
end

if isdefined(Pkg, :Registry) && isdefined(Pkg.Registry, :registry_info)
    const registry_info = Pkg.Registry.registry_info
elseif isdefined(Pkg, :RegistryHandling) && isdefined(Pkg.RegistryHandling, :registry_info)
    const registry_info = Pkg.RegistryHandling.registry_info
end

"""
    collect_jll_uuids(manifest::Pkg.Types.Manifest, dependencies::Set{UUID})

Return a `Set` of all JLL packages in the `manifest` with `dependencies` being
the list of direct dependencies of the environment.
"""
function collect_jll_uuids(manifest::Pkg.Types.Manifest, dependencies::Set{UUID})
    jlls = copy(dependencies)
    for (uuid, pkg_entry) in manifest
        if uuid in jlls
            for (dep_name, dep_uuid) in pkg_entry.deps
                if endswith(dep_name, "_jll")
                    push!(jlls, dep_uuid)
                end
            end
        end
    end
    if jlls == dependencies
        return jlls
    else
        return collect_jll_uuids(manifest, jlls)
    end
end

"""
    get_addable_spec(name::AbstractString, version::VersionNumber)

Given a JLL name and registered version, return a `PackageSpec` that, when passed as a
`Dependency`, ensures that exactly that version will be installed.  Example usage:

    dependencies = [
        BuildDependency(get_addable_spec("LLVM_jll", v"9.0.1+0")),
    ]
"""
function get_addable_spec(name::AbstractString, version::VersionNumber;
                          ctx = Pkg.Types.Context(), verbose::Bool = false)
    # First, resolve the UUID
    uuid = first(Pkg.Types.registry_resolve!(ctx.registries, Pkg.Types.PackageSpec(;name))).uuid

    # Next, determine the tree hash from the registry
    repo_urls = Set{String}()
    tree_hashes = Set{Base.SHA1}()
    for reg in ctx.registries
        if !haskey(reg, uuid)
            continue
        end

        pkg_info = registry_info(reg[uuid])
        if pkg_info.repo !== nothing
            push!(repo_urls, pkg_info.repo)
        end
        if pkg_info.version_info !== nothing
            if haskey(pkg_info.version_info, version)
                version_info = pkg_info.version_info[version]
                push!(tree_hashes, version_info.git_tree_sha1)
            end
        end
    end

    if isempty(tree_hashes)
        @error("Unable to find dependency!",
            name,
            version,
            registries=ctx.registries,
        )
        error("Unable to find dependency!")
    end
    if length(tree_hashes) != 1
        @error("Multiple treehashes found!",
            name,
            version,
            tree_hashes,
            registries=ctx.registries,
        )
        error("Multiple treehashes found!")
    end

    tree_hash_sha1 = first(tree_hashes)

    # Once we have a tree hash, turn that into a git commit sha
    git_commit_sha = nothing
    valid_url = nothing
    for url in repo_urls
        git_commit_sha = get_commit_sha(url, tree_hash_sha1; verbose)
        # Stop searching urls as soon as we find one
        if git_commit_sha !== nothing
            valid_url = url
            break
        end
    end

    if git_commit_sha === nothing
        @error("Unable to find revision for specified dependency!",
            name,
            version,
            tree_hash = bytes2hex(tree_hash_sha1.bytes),
            repo_urls,
        )
        error("Unable to find revision for specified dependency!")
    end

    return Pkg.Types.PackageSpec(
        name=name,
        uuid=uuid,
        #version=version,
        tree_hash=tree_hash_sha1,
        repo=Pkg.Types.GitRepo(rev=git_commit_sha, source=valid_url),
    )
end

# We only want to update the registry once per run
registry_updated = false
function update_registry(outs = stdout, )
    global registry_updated
    if !registry_updated
        Pkg.Registry.update(
            [Pkg.RegistrySpec(uuid = "23338594-aafe-5451-b93e-139f81909106")];
            io=outs,
        )
        registry_updated = true
    end
end


text:
#=
Created on Tuesday 11 August 2020
Last update: Tuesday 18 August 2020

@author: Michiel Stock
michielfmstock@gmail.com

Linear models with zero-"norm" regularization.
Useful for feature selection.
=#

module ZeroNorm

export momentum!, flip, findzeronorm, SimulatedAnnealing, L0Solver, HillClimbing

using Zygote, LinearAlgebra

"""
momentum!(w, s, dl!, dw;
		α=0.01, β=0.5, niter=100, ϵ=1e-4)

Gradient descent with momentum to fit sparse linear models. Fits a parameter
vector `w` with sparsity structure `s`. Requires a function `dl!` (in-place
gradient of the loss function) and a vector `dw` to allocate the gradients.
"""
function momentum!(w::AbstractVecOrMat,
					s::AbstractVecOrMat{Bool},
					dl!,
					dw::AbstractVecOrMat;
					α=0.01, β=0.5, niter=100, ϵ=1e-4)
	@assert size(w) == size(s) == size(dw)
	dw .= 0.0
	for i in 1:niter
		# update gradient
		dl!(w, dw, β)
		dw .*= s
		norm(dw) < ϵ && break
		w .-= α .* dw
	end
	return w
end

abstract type L0Solver end

# SIMULATED ANNEALING

struct SimulatedAnnealing <: L0Solver
	p₀
	p₁
	r
	kT
	Tmin
	Tmax
end

"""
	SimulatedAnnealing(;p₀=0.05, p₁=0.05, r=0.8, kT=1000, Tmin=1e-3, Tmax=1e3)

Solver for L0-norm minimization. Uses parameters:

	- `p₀` : probability of true -> false in `s`
	- `p₁`: probability of false -> true in `s`
	- `r` : cooling rate
	- `kT` : number of repetitions per temperature
	- `Tmax` : starting temperature
	- `Tmin` : stopping temperature
"""
SimulatedAnnealing(;p₀=0.05, p₁=0.05, r=0.8, kT=1000,
			Tmin=1e-3, Tmax=1e3) = SimulatedAnnealing(p₀, p₁, r, kT, Tmin, Tmax)

struct HillClimbing <: L0Solver end

"""
Turns a `true` to a `false` with probability `p₀` and a `false` to a `true`
with probability `p₁`.
"""
flip(bit::Bool, p₀, p₁) = bit ? rand() > p₀ : rand() < p₁

function findzeronormSA(w, loss, λ, s, p₀, p₁, r, kT,
			Tmin, Tmax; verbose::Bool=false,
			momentumpars...)
	n = length(w)
	dw = similar(w)
	# compute gradient for momentum
	dl!(w, dw, β) = dw .= (1.0 - β) .* loss'(w) .+ β .* dw
	snew = copy(s)
	sbest = copy(s)
	j(w, s) = loss(w .* s) + λ * sum(s)
	T = Tmax
	best_obj = Inf
	current_obj = Inf
	# create structure to store the objective through the iterations
	obj_vals = Vector{Float64}(undef, ceil(Int, (log(Tmin) - log(Tmax)) / log(r)))
	iter = 0
	while T > Tmin
		iter += 1
		for rep in 1:kT
			# perturbate s
			snew .= flip.(s, p₀, p₁)
			# update w
			momentum!(w, snew, dl!, dw; momentumpars...)
			obj = j(w, snew)
			if obj < current_obj || rand() < exp((current_obj - obj)/T)
				current_obj = obj
				s .= snew
			end
			if obj < best_obj
				best_obj = obj
				sbest .= s
			end
		end
		verbose && println("T=$T : obj=$best_obj, nfeatures=$(sum(sbest))")
		obj_vals[iter] = best_obj
		T *= r
		end
	return sbest
end

# HILL CLIMBING

function findzeronormHC(w, loss, λ, s; verbose::Bool=false,
		momentumpars...)
	n = length(w)
	dw = similar(w)
	# compute gradient for momentum
	dl!(w, dw, β) = dw .= (1.0 - β) .* loss'(w) .+ β .* dw
	j(w, s) = loss(w .* s) + λ * sum(s)
	best_obj = Inf
	improved = true
	iter = 0
	while improved
		improved = false
		iter += 1
		best_flip = 0
		for i in 1:n
			# perturbate s
			s[i] = !s[i]
			# update w
			momentum!(w, s, dl!, dw; momentumpars...)
			obj = j(w, s)
			if obj < best_obj
				best_obj = obj
				improved = true
				best_flip = i
			end
			# change s back
			s[i] = !s[i]
		end
	# update s
	improved && (s[best_flip] = !s[best_flip])
	verbose && improved && println("iteration $iter, $(s[best_flip] ? "added $best_flip" : "removed $best_flip"), obj=$best_obj, nfeatures=$(sum(s))")
	end
return s, w
end

function findzeronorm(w0, loss, λ, s, sa::SimulatedAnnealing; kwargs...)
	return findzeronormSA(w0, loss, λ, s, sa.p₀, sa.p₁, sa.r, sa.kT, sa.Tmin,
				sa.Tmax; kwargs...), w0
end

findzeronorm(w0, loss, λ, s, ::HillClimbing; kwargs...) = findzeronormHC(w0, loss, λ, s; kwargs...)

findzeronorm(w0, loss, λ, sol::L0Solver; kwargs...) = findzeronorm(w0, loss, λ, ones(Bool, length(w0)), sol; kwargs...)

"""
	findzeronorm(w0, loss, λ[, s0], sol; verbose=true, momentumpars...)

Finds the zero norm for a model with loss function regularized with a zero-"norm":

	min_w loss(w) + ||w||_0

Uses Simulated Annealing or Hill Climbing and returns `s` a binary vector containing the sparsity
and a vector with the objective value throughout the iterations.

For SA, use `sol=SimulatedAnnealing()` (with optional parameters), for HC, use `sol=HillClimbing()`.

Parameters:
	- `w0` : starting value for the parameter vector
	- `loss` : a differentiable loss function (we use autodiff)
	- `λ` : regularization parameter
	- `s0` : initial guess for `s`, optional, uses all features initially by default
	- `verbose` : flag to print the convergence during fitting
	- `momentumpars` : parameters to determine behaviour for momentum
"""
findzeronorm(w0, loss, λ, s, sol::L0Solver; kwargs...) = nothing

end  # module ZeroNorm


text:
# getfunction.jl
#
# Makro to ease access to C functions; mainly adopted from the GetC-Package
# for Julia. What is desired is e.g. the following:
#
# INPUT:
# void mgl_surf(HMGL gr, HCDT z, char *sch, char *opt);
# mgl_surf(gr::HMGL, z::HCDT, sch::Ptr{Cchar}, opt::Ptr{Cchar})::Void
#
# OUTPUT:
# function mgl_surf(gr, z, sch, opt)
#     ccall(("libmgl", :mgl_surf), Void, (HMGL, HCDT, Ptr{Cchar}, Ptr{Cchar}), gr, z, sch, opt)
# end

macro getf(lib, fun)
    # Name of the function
    cname     = string(fun.args[1].args[1])                # :mgl_surf
    jname     = Symbol(replace(string(cname), "mgl_", "")) # :surf
    arguments = fun.args[1].args[2:end]                    # [:(gr::HMGL), :(z::HCDT), ...]

    # Get info out of arguments of 'fun'
    inputC   = [ arg.args[2] for arg in arguments ]   # [ :HMGL, :HCDT, ...]
    varnames = [ arg.args[1] for arg in arguments ]   # [ :gr, :z, :sch, :opt ]
    rettype  = fun.args[2]                            # :Void

    # Construct the result.
    csym = Expr(:quote, cname)
    body = Expr(:ccall, :($csym, $lib), rettype, Expr(:tuple, inputC...),
                varnames...)
    ret = Expr(:function, Expr(:call, jname, varnames...), body)
    return esc(ret)
end


text:
#=
# 2D Poisson, Dirichlet bc
# Uses Matlab imported as a custom gen target
=#
if !@isdefined(Femshop)
    include("../Femshop.jl");
    using .Femshop
end
init_femshop("poisson2dcustomnew");

# Try making an optional log
useLog("poisson2dcustomnewlog")

# Generate for the target in customtarget.jl
generateFor("target_matlab_cg.jl")

# Set up the configuration (order doesn't matter)
domain(2)
functionSpace(order=2)

# Specify the problem
mesh(QUADMESH, elsperdim=30)

u = variable("u")

testSymbol("v")

boundary(u, 1, DIRICHLET, "sin(3*pi*x)")

# Write the weak form
coefficient("f", "-2*pi*pi*sin(pi*x)*sin(pi*y)")
weakForm(u, "-dot(grad(u),grad(v)) - f*v")

solve(u);

finalize_femshop()


text:
module ReferenceTests

# Write your package code here.
using Test
using MeshIO
using FileIO
using MacroTools
using Makie
using Makie: Record, Stepper, Axis
using Makie.FFMPEG
using Printf
using ghr_jll
using Tar
using Downloads
using Pkg.TOML
using Statistics
using ImageShow

basedir(files...) = normpath(joinpath(@__DIR__, "..", files...))
loadasset(files...) = FileIO.load(assetpath(files...))

# The version in Images.jl throws an error... whyyyyy!?
# TODO look into error!
using Images, FixedPointNumbers, Colors, ColorTypes

include("database.jl")
include("stable_rng.jl")
include("runtests.jl")
include("image_download.jl")
include("html_rendering.jl")

end


text:
using ProgressBars, Plots, LaTeXStrings, JLD

choices = ((1,0),(0,1),(-1,0),(0,-1))

function RecurseSAW(R, Network, N)
    if N == 0
        return 1
    end
    count_ = 0
    Network[R...] = true
    for step in choices
        Rₙ = R .+ step
        if !Network[Rₙ...]
            count_ += RecurseSAW(Rₙ, Network, N-1)
        end
    end
    Network[R...] = false
    return count_
end

function maincal(N)
    dim = 2N+1
    Network = falses(dim,dim)
    R = [N+1,N+1]
    return RecurseSAW(R, Network, N)
end

TotList = [maincal(n) for n ∈ ProgressBar(0:17)]
save("../../Data/Q3/Q3-SAW.jld", "Data", TotList)

P1 = begin
    plot(0:17,log.(TotList), c = :black, label = nothing)
    scatter!(0:17,log.(TotList), c = :steelblue,framestyle = :box,
        ylabel = L"ln\ (Number\ of\ Paths)", xlabel = L"N",
        label = L"Self-Avoiding\ Walker", legend = 140)
        plot!(0:17, log.(4 .^(0:17)), c = :black, label = nothing)
        scatter!(0:17, log.(4 .^(0:17)), c = :purple, label = L"Random\ Walker")
end

P2 = begin
    plot(0:17,TotList, c = :black, label = nothing)
    scatter!(0:17,TotList, c = :steelblue,framestyle = :box,
        ylabel = L"Number\ of\ Paths", xlabel = L"N",
        label = L"Self-Avoiding\ Walker", legend = 140)
        plot!(0:13,4 .^(0:13), c = :black, label = nothing)
        scatter!(0:13, 4 .^(0:13), c = :purple, label = L"Random\ Walker")
end

savefig(P1,"../../Figs/Q3/Q3-1.pdf")
savefig(P2,"../../Figs/Q3/Q3-2.pdf")

P3 = begin
    plot(0:17, TotList./ 4 .^(0:17), c = :black, label = nothing)
    scatter!(0:17, TotList./ 4 .^(0:17), c = :steelblue,framestyle = :box,
        ylabel = L"Ratio\ of\ Number\ of\ Paths", xlabel = L"N", legend = nothing,
            title = L"Ratio\ of\ Number\ of\ Paths\ of\ SAW\ and\ RW")
end

savefig(P3,"../../Figs/Q3/Q3-3.pdf")


text:
# # SNAP example using GaN
#
# In this example we extract the bispectrum from SNAP.
# This example demonstrates how to install and use LAMMPS.jl
#
# ## Install dependencies
#
# First let's make sure we have all required packages installed

# ```julia
# using Pkg
# pkg"add LAMMPS"
# ```

# LAMMPS.jl will automatically download and install a pre-built version of LAMMPS.
#
# We start off by importing `LAMMPS`.

using LAMMPS

const DATA = joinpath(dirname(pathof(LAMMPS)), "..", "examples", "example_GaN_data")

function run_snap(lmp, path, rcut, twojmax)
    read_data_str = "read_data " * path

    command(lmp, "log none")
    command(lmp, "units metal")
    command(lmp, "boundary p p p")
    command(lmp, "atom_style atomic")
    command(lmp, "atom_modify map array")
    command(lmp, read_data_str)
    command(lmp, "pair_style zero $rcut")
    command(lmp, "pair_coeff * *")
    command(lmp, "compute PE all pe")
    command(lmp, "compute S all pressure thermo_temp")
    command(lmp, "compute SNA all sna/atom $rcut 0.99363 $twojmax 0.5 0.5 1.0 0.5 rmin0 0.0 bzeroflag 0 quadraticflag 0 switchflag 1")
    command(lmp, "compute SNAD all snad/atom $rcut 0.99363 $twojmax 0.5 0.5 1.0 0.5 rmin0 0.0 bzeroflag 0 quadraticflag 0 switchflag 1")
    command(lmp, "compute SNAV all snav/atom $rcut 0.99363 $twojmax 0.5 0.5 1.0 0.5 rmin0 0.0 bzeroflag 0 quadraticflag 0 switchflag 1")
    command(lmp, "thermo_style custom pe")
    command(lmp, "run 0")

    ## Extract bispectrum
    bs = extract_compute(lmp, "SNA", LAMMPS.API.LMP_STYLE_ATOM,
                                     LAMMPS.API.LMP_TYPE_ARRAY)
    return bs
end

function calculate_snap_bispectrum(path, rcut, twojmax, M, N1, N2)
    J = twojmax / 2
    ncoeff = round(Int, (J+1)*(J+2)*((J+(1.5))/3) + 1)

    A = LMP(["-screen","none"]) do lmp
        A = Array{Float64}(undef, M, 2*ncoeff) # bispectrum is 2*(ncoeff - 1) + 2
        for m in 1:M
            data = joinpath(path, string(m), "DATA")
            bs = run_snap(lmp, data, rcut, twojmax)

            ## Make bispectrum sum vector
            row = Float64[]

            push!(row, N1)

            for k in 1:(ncoeff-1)
                acc = 0.0
                for n in 1:N1
                    acc += bs[k, n]
                end
                push!(row, acc)
            end

            push!(row, N2)
            for k in 1:(ncoeff-1)
                acc = 0.0
                for n in N1 .+ (1:N2)
                    acc += bs[k, n]
                end
                push!(row, acc)
            end

            A[m, :] = row

            command(lmp, "clear")
        end
        return A
    end
    return A, ncoeff
end

const rcut = 3.5
const twojmax = 6
const M  = 48 # number of input files
const N1 = 96 # number of atoms of the first type
const N2 = 96 # number of atoms of the second type
A, ncoeff = calculate_snap_bispectrum(DATA, rcut, twojmax, M, N1, N2)


text:
function rms{T<:Real}(s::Array{T,1})
    # TODO: add argument valdiation
    m = zero(T)
    for si in s
        m += si * si
    end
    sqrt(m/length(s))
end

function rms{T<:Real}(s::Array{Complex{T},1})
    # TODO: add argument valdiation
    m = zero(T)
    for si in s
        m += abs2(si)
    end
    sqrt(m/length(s))
end


text:
# parser trial

PUNCTS = [',']
ispunc(x) = x in PUNCTS

DELIMS = ['.','(',')']
isdelim(x) = x in DELIMS
delimDict=Dict('.' => :dot, '('=>:lpar, ')'=>:rpar)

islpar(x) = x == '('
isrpar(x) = x == ')'
isdot(x)  = x == '.'

LSYMS = ['∀','∃','∨','∧','⇒','¬','≡']
islsym(x) = x in LSYMS

lsymDict=Dict(
'¬'=>:not,
'∀'=>:all,
'∃'=>:some,
'∨'=>:or,
'∧'=>:and,
'⇒'=>:imply,
'≡'=>:eqv
)

isforall(x) = x=="∀"
isexist(x)  = x=="∃"
isor(x)     = x=="∨"
isand(x)    = x=="∧"
isthen(x)   = x=="⇒"
isiff(x)    = x=="≡"
isnot(x)    = x=="¬"

function breaksym(str)
 ix = 1
 while ix <= lastindex(str)
  println(str[ix])
  ix = nextind(str, ix)
 end
end

function findsym(str, six)
  ix = six
  while ix<=lastindex(str)    &&
    	!isspace(str[ix]) &&
	!isdelim(str[ix]) &&
	!ispunc(str[ix])  &&
	!islsym(str[ix])
   ix = nextind(str, ix)
  end
  if ix == six; error("usage") end
  return(str[six:ix-1], ix-1)
end

# sample code for tokenizer
function tokenizer(str)
 ix = 1
 toklist = []
 nx = 0
 while ix <= lastindex(str)
   c = str[ix]
   if isspace(c);
   elseif ispunc(c);
   elseif isdelim(c)
     push!(toklist,delimDict[c])
   elseif islsym(c)
     push!(toklist,lsymDict[c])
   else
     sym,ix = findsym(str, ix)
     push!(toklist,Symbol(sym))
   end

   ix = nextind(str, ix)

 end
 return toklist
end

function parser(str)
 strn = str
 ix = 1

 toklist = []
 while ix < length(str)
   nx = findnext(SEPRS, str, ix)
   if nx == 0
     tok=strip(str[ix:end])
     println("$(tok)")
     push!(toklist,Symbol(tok))
     return toklist,length(str)
   end
   tok=strip(str[ix:(nx-1)])
   if !all(isspace, tok)
     println("$(tok)")
     push!(toklist,Symbol(tok))
   end
   sep=str[nx]
   if islpar(sep);
     println("(")
     block,nnx = parser(str[nx+1:end])
     nx = nx + nnx + 1
     push!(toklist, block)
   end

   if isrpar(sep)
     println(")")
     return toklist, nx
   end

#   println("sep=$(str[nx])")
   ix = nx+1
 end
 return toklist, ix
end

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
"""
    plot(lattice)

Creates a 2D plot of lattice using PyPlot.
"""
function plot(lattice::T, args...; kwargs...) where {T <: AbstractGraph}
    D = dims(lattice.unitcell)
    if D == 2
        plot2D(lattice, args...; kwargs...)
    elseif D == 3
        plot3D(lattice, args...; kwargs...)
    else
        error("No method available to plot $D-dimensional lattice.")
    end
end

function plot2D(lattice::T) where {T <: AbstractGraph}
    @assert dims(lattice.unitcell) == 2 "plot2D only works for 2D lattices."

    # Recover node positions
    positions = [
        get_pos(
            lattice.unitcell,
            node.uvw[1],
            node.uvw[2:end]
        ) for node in lattice.nodes
    ]

    # Plot nodes/sites
    PyPlot.plot(
        [p[1] for p in positions[:]],
        [p[2] for p in positions[:]],
        "ko"
    )

    # plot edges/bonds
    for level in lattice.edges
        for edge in level
            PyPlot.plot(
                [positions[edge.node1.uvw...][1], positions[edge.node2.uvw...][1]],
                [positions[edge.node1.uvw...][2], positions[edge.node2.uvw...][2]],
                "k-"
            )
        end
    end

    PyPlot.axis("equal")
    lattice
end


text:
import MPI

include("types.jl")
include("init.jl")
include("get_sim_params.jl")
include("update.jl")
include("io.jl")
include("get_sim_params.jl")

function main()

#temp1_init: temperature init on borders
temp1_init = 10.0
#temp2_init: temperature init inside domain
temp2_init = -10.0
#diffusion coefficient
k0 = Float64(1)

#define rank of root process
root = 0

#MPI Initialization
MPI.Init()
comm = MPI.COMM_WORLD
my_id = MPI.Comm_rank(comm)
nproc = MPI.Comm_size(comm)

#simulation parameters
params_int = Array{Int64}(undef, 5)
params_double = Array{Float64}(undef, 2)
if (my_id == root)
    params_int, params_double, output_path = get_sim_params(ARGS)
end

MPI.Barrier(comm)
MPI.Bcast!(params_int, 5, 0, comm)
MPI.Bcast!(params_double, 2, 0, comm)

size_x    = params_int[1]
size_y    = params_int[2]
nx_domains = params_int[3]
ny_domains = params_int[4]
maxStep   = params_int[5]
dt1       = params_double[1]
epsilon   = params_double[2]

#Warning message if dimensions and number of processes don't match
if ((my_id == root) && (nproc != (nx_domains * ny_domains)))
    println("ERROR - Number of processes not equal to number of subdomains")
end

#Various other variables
size_global_x = size_x + 2
size_global_y = size_y + 2
hx = Float64(1.0 / size_global_x)
hy = Float64(1.0 / size_global_y)
dt2 = 0.25 * min(hx, hy)^2 / k0 #fraction of CFL condition
size_total_x = size_x + 2 * nx_domains + 2 #including ghost cells
size_total_y = size_y + 2 * ny_domains + 2 #including ghost cells

#Take a right time step for convergence
if (dt1 >= dt2)
    if (my_id == 0)
        println()
        println("Time step too large ==> Taking convergence criterion")
    end
    dt = dt2
else
    dt = dt1
end

#2D solution including ghost cells
u0 = twoDArray(Float64, size_total_x, size_total_y)
u = twoDArray(Float64, size_total_x, size_total_y)

#Allocate coordinates of processes (start cell, end cell)
xs = oneDArray(Int, nproc)
xe = oneDArray(Int, nproc)
ys = oneDArray(Int, nproc)
ye = oneDArray(Int, nproc)

#Size of each physical domain
xcell = Int(size_x / nx_domains)
ycell = Int(size_y / ny_domains)

#allocate flattened (1D) local physical solution (i.e., relative to sub-domain)
u_local = oneDArray(Float64, xcell * ycell)
#allocate flattened (1D) global physical solution
u_global = oneDArray(Float64, size_x * size_y)

#find processes surrounding mine (if they exist)
my_neighbors = neighbors(my_id, nproc, nx_domains, ny_domains)

#compute coordinates of processes for each sub-domain
process_coordinates!(xs, ys, xe, ye, xcell, ycell, nx_domains, ny_domains, nproc)

#initialize domain
init_values(u0, size_total_x, size_total_y, temp1_init, temp2_init)

#update ghost cells
updateBound!(u0, size_total_x, size_total_y, my_neighbors, comm,
            my_id, xs, ys, xe, ye, xcell, ycell, nproc)

#Initialize step, time and convergence boolean
step = 0
t = 0.0
converged = false

#Starting time
if (my_id==0)
    time_init = time()
end

#Main loop : until convergence
while (!converged)
    #increment step and time
    step += 1
    t += dt
    #perform one step of the explicit scheme
    local_diff = computeNext!(u0, u, size_total_x, size_total_y, dt, hx, hy,
        my_id, xs, ys, xe, ye, nproc, k0)

    #update ghost cells
    updateBound!(u0, size_total_x, size_total_y, my_neighbors, comm,
        my_id, xs, ys, xe, ye, xcell, ycell, nproc)

    MPI.Barrier(comm)

    #sum local_diff to get global difference
    global_diff = MPI.Allreduce(local_diff, MPI.SUM, comm)
    global_diff = sqrt(global_diff)
    #break if convergence reached or step greater than maxStep
    if ((global_diff <= epsilon) || (step >= maxStep))
        converged = true
    end
end

#get ending time
if (my_id == 0)
    time_final = time()
    #elapsed time
    elapsed_time = time_final - time_init
    println("Elapsed time = ", elapsed_time)
    println("Steps = ", step)
end

#find solution on my sub-domain (as a 1-dimensional array)
i = 1
for j = ys[my_id+1]:ye[my_id+1]
    u_local[(i-1)*xcell+1:i*xcell] = u0[xs[my_id+1]:xe[my_id+1],j]
    i = i+1
end

#gather local solution to global solution (as a 1-dimensional array)
u_global = MPI.Gather(u_local, root, comm)

#write to disk
if (my_id == 0)
    write_to_disk(u_global, nx_domains, ny_domains, xcell, ycell, temp1_init, output_path)
end

MPI.Finalize()

end

main()


text:
import Observables
import AbstractPlotting
import FileIO

scene = AbstractPlotting.Scene()
circle = U1(rand() * 2pi - pi)
s3rotation = Quaternion(rand() * 2pi - pi, ℝ³(rand(3)))
config = Biquaternion(Quaternion(rand() * 2pi - pi, ℝ³(rand(3))), ℝ³(rand(3)))
segments = rand(5:10)
color = FileIO.load("data/basemap90grid.png")
transparency = rand(1:2) == 1 ? true : false
frame = Frame(scene,
              circle,
              color,
              s3rotation = s3rotation,
              config = config,
              segments = segments,
              transparency = transparency)

## update circle

circle2 = U1(rand() * 2pi - pi)
update(frame, circle2)

@test isapprox(frame.circle, circle2)

## update S³ rotation

s3rotation2 = Quaternion(rand() * 2pi - pi, ℝ³(rand(3)))
update(frame, s3rotation2)

@test isapprox(frame.s3rotation, s3rotation2)

## update configuration

config2 = config * Biquaternion(Quaternion(rand() * 2pi - pi, ℝ³(rand(3))), ℝ³(rand(3)))
value1 = getsurface(frame.observable, frame.segments)
update(frame, config2)
value2 = getsurface(frame.observable, frame.segments)

@test isapprox(frame.config, config2)
@test isapprox(value1, value2) == false


text:
include("../pwgrid_02/PWGrid.jl")

function latvec_hexagonal(a; ca=1.0)
  LL = zeros(3,3)
  LL[1,:] = [1.0, 0.0, 0.0]
  LL[2,:] = [cos(pi/3.0), sin(pi/3.0), 0.0]
  LL[3,:] = [0.0, 0.0, ca]
  return a*LL
end

function ewald()

  Ns = [60, 60, 60]
  const alat = 16.0
  LatVecs = alat*diamg(ones(3))
  pw = PWGrid( Ns, LatVecs )

  atpos = [0.0, 0.0, 0.0]

  nat = 1
  ntyp = 1
  ityp = [1]  # ityp(nat)
  zv = [1.0]     # zv(ntyp)  # need to conform with `ntypx` in PWSCF ?
  at = LatVecs'/alat
  bg = inv(at')
  #bg = pw.RecVecs'/(2*pi) # ???
  tau = atpos/alat
  omega = pw.Ω
  g  = pw.G
  gg = pw.G2
  ngm = prod(Ns) # for the moment
  gcutm =

  println( at )
  println( bg )
  println( at * bg' )

  type_args = (
               Ref{Float64},    # alat
               Ref{Int64},      # nat
               Ref{Int64},      # ntyp
               Ptr{Int64},      # ityp
               Ptr{Float64},    # zv
               Ptr{Float64},    # at
               Ptr{Float64},    # bg
               Ptr{Float64},    # tau
               Ref{Float64},    # omega
               Ptr{Float64},    # g
               Ptr{Float64},    # gg
               Ref{Int64},      # ngm
               Ref{Float64},    # gcutm
               Ref{Int64},      # gstart
               Ref{Int32},      # gamma_only
               Ptr{Complex128}, # strf
               )

  #ccall( (:ewald_, "libqe.so"), Float64,
  #       type_args,
  #       alat, nat, ntyp, ityp, zv, at, bg, tau, omega, g, gg,
  #       ngm, gcutm, gstart, gamma_only, strf )
end

ewald()


text:
@testset "emission.jl" begin
    println("             > emission.jl")
    m_out  = 16
    m_in   = 16

    state = BinomialState(128, m_out, m_in, :cpu)
    model = BinomialModel(128, m_out, :cpu)
    @test_throws ErrorException propagate_emit!(state, model)

    state = ScalarBinomialState(10, 2)
    model = ScalarBinomialModel(10, 0.85, 1.0, 0.2, 0.2)
    @test propagate_emit!(state, model).dt > 0
    @test propagate_emit!(state, model, FixedTimestep(0.3)).dt == 0.3
    @test propagate_emit!(state, model, RandomTimestep(Exponential(0.5))).dt > 0
end


text: asset_path(files...) = joinpath(@__DIR__, "..", "assets", files...)


text:
using Test
import Pluto

@testset "Existing notebooks " begin
    server = Pluto.ServerSession()

    project_dir = joinpath(@__DIR__, "..", "notebook")

    for path in readdir(project_dir)
        @testset "$path" begin
            path = normpath(joinpath(project_dir, path))

            nb = @test_nowarn Pluto.load_notebook(path)

            client = Pluto.ClientSession(Symbol("client", rand(UInt16)), nothing)
            client.connected_notebook = nb

            server.connected_clients[client.id] = client

            Pluto.update_run!(server, nb, nb.cells)
            for c in nb.cells
                @test c.errored == false
            end
        end
    end
end

max_stars_repo_name: MattAmos/Example-Pluto-CI

max_issues_repo_name: MattAmos/Example-Pluto-CI

max_forks_repo_name: MattAmos/Example-Pluto-CI


text:
using RecipesBase

@recipe function f(img::AstroImage)
    seriestype   := :heatmap
    aspect_ratio := :equal
    # Right now we only support single frame images,
    # gray scale is a good choice.
    color        := :grays
    img.data
end


text:
let
	startcalibration = 2008		# originally 1960
	endcalibration = 2010

	calibrationtime = 1:iyear(endcalibration)
	rcplist = ["RCP3PD", "RCP45", "RCP6", "RCP85"]
	cache_constantforcing =
		Dict(r => sum(forcing_conc_RCP[r][g][calibrationtime] for g in GAS) + otherforcing[r][calibrationtime] for r in rcplist)

	global function forcingerror(aerosolforcingfactor, p::ClimateParams, rcp)
		TotalRadiativeForcing_annual = cache_constantforcing[rcp] + aerosolforcingfactor * aerosolforcing[rcp][calibrationtime]

		s = ClimateState()
		init_temperatures!(s)
		Error_Temp = 0.0

		for t in YEARS[1]:p.timestep:endcalibration
			s.TotalRadiativeForcing = interpolate(t, TotalRadiativeForcing_annual)
			temperatures!(s, p)
			if t >= startcalibration
				Error_Temp += (interpolate(t, histTempGISS) - s.Temp_global)^2
			end
		end

		return Error_Temp
	end

	global function calibrateforcing!(p::ClimateParams, rcp)
		result = optimize(aerosolforcingfactor -> forcingerror(aerosolforcingfactor,p,rcp), -2.0, 5.0)
		aerosolforcingfactor = Optim.minimizer(result) 	# 1.3925
		@pack! p = aerosolforcingfactor
	end

	function concentrationerror(fertilization, annualemissions, p::ClimateParams, rcp)
		@pack! p = fertilization
		results::Vector{ClimateState} = solveclimate(annualemissions, p, 1765, endcalibration, false, rcp)

		error = 0.0
		for i=iyear(startcalibration):iyear(endcalibration)
			error += (results[i].Concentration[:CO2] - conc_RCP[rcp][:CO2][i])^2
		end

		return error
	end

	global function calibratefertilization!(annualemissions, p::ClimateParams, rcp)
		result = optimize(fert -> concentrationerror(fert, annualemissions, p, rcp), 0.2, 0.8, rel_tol=1e-4)

		fertilization::Float64 = Optim.minimizer(result) 	# 0.5886
		#println("mean concentration error: ",sqrt(Optim.minimum(result))/(endcalibration-startcalibration+1))
		@pack! p = fertilization
	end
end

max_stars_repo_name: niclasmattsson/ClimateCalculator

max_issues_repo_name: niclasmattsson/ClimateCalculator

max_forks_repo_name: niclasmattsson/ClimateCalculator


text:
module Tracker

export TrackedArray, param, back!

data(x) = x
istracked(x) = false

struct Call{F,As<:Tuple}
  func::F
  args::As
end

Call(f, args...) = Call{typeof(f),typeof(args)}(f, args)

(c::Call)() = c.func(data.(c.args)...)

mutable struct TrackedArray{T,N,A} <: AbstractArray{T,N}
  ref::UInt32
  f::Call
  data::A
  grad::A
  TrackedArray{T,N,A}(f::Call, data::A) where {T,N,A} = new(0, f, data)
  TrackedArray{T,N,A}(f::Call, data::A, grad::A) where {T,N,A} = new(0, f, data, grad)
end

TrackedScalar{T,A} = TrackedArray{T,0,A}
TrackedVector{T,A} = TrackedArray{T,1,A}
TrackedMatrix{T,A} = TrackedArray{T,2,A}
TrackedVecOrMat{T,A} = Union{TrackedVector{T,A},TrackedMatrix{T,A}}

TrackedArray(c::Call, x::A) where A <: AbstractArray =
  TrackedArray{eltype(A),ndims(A),A}(c, x)

TrackedArray(c::Call, x::A, Δ::A) where A <: AbstractArray =
  TrackedArray{eltype(A),ndims(A),A}(c, x, Δ)

TrackedArray(c::Call) = TrackedArray(c, c())

TrackedArray(x::AbstractArray) = TrackedArray(Call(nothing), x, zeros(x))

param(xs) = TrackedArray(AbstractFloat.(xs))
param(xs::Real) = param(fill(xs))

istracked(x::TrackedArray) = true
data(x::TrackedArray) = x.data
grad(x::TrackedArray) = x.grad

# Fallthrough methods

for f in :[Base.size, Base.ndims].args
  @eval @inline $f(x::TrackedArray, a...) = $f(data(x), a...)
end

Base.similar(x::TrackedArray, dims::Union{AbstractUnitRange,Integer}...) =
  similar(data(x), dims...)

Base.similar(x::TrackedArray, T::Type) = similar(data(x), T)

value(x) = x
value(x::TrackedArray) = data(x)
value(x::TrackedScalar) = data(x)[]

Base.:(==)(x::TrackedArray, y) = value(x) == y
Base.:(==)(y, x::TrackedArray) = y == value(x)
Base.:(==)(x::TrackedArray, y::TrackedArray) = value(x) == value(x)

Base.isless(x::TrackedScalar, y) = isless(value(x), y)
Base.isless(x, y::TrackedScalar) = isless(x, value(y))
Base.isless(x::TrackedScalar, y::TrackedScalar) = isless(value(x), value(y))

Base.show(io::IO, ::Type{TrackedArray{T,N,A}}) where {T,N,A<:AbstractArray{T,N}} =
  print(io, "TrackedArray{…,$A}")

function Base.showarray(io::IO, X::TrackedArray, repr::Bool = true; header = true)
  if repr
    print(io, "param(")
    Base.showarray(io, data(X), true)
    print(io, ")")
  else
    header && print(io, "Tracked ")
    Base.showarray(io, data(X), false, header = header)
  end
end

Base.setindex!(xs::TrackedArray, v, i...) =
  error("Can't differentiate `setindex!`")

include("back.jl")
include("lib.jl")
include("numeric.jl")

import NNlib.adapt

adapt(T, xs::TrackedArray) = TrackedArray(xs.f, adapt(T, xs.data), adapt(T, xs.grad))

end


text:
__precompile__()

module CatViews

using Base.Cartesian
import Iterators: chain, repeated

export CatView, splitview

include("catview.jl")
include("splitview.jl")

# end module
end


text:
mutable struct URLTest
    name::String
    url::String
    isconnect::Bool
    offsets::NTuple{7, HTTP.URIs.Offset}
    shouldthrow::Bool
end

URLTest(nm::String, url::String, isconnect::Bool, shouldthrow::Bool) = URLTest(nm, url, isconnect, ntuple(x->HTTP.URIs.Offset(), 7), shouldthrow)

@testset "HTTP.URI" begin
    # constructor
    @test string(HTTP.URI("")) == ""
    @test HTTP.URI(hostname="google.com") == HTTP.URI("http://google.com")
    @test HTTP.URI(hostname="google.com", path="/") == HTTP.URI("http://google.com/")
    @test HTTP.URI(hostname="google.com", userinfo="user") == HTTP.URI("http://user@google.com")
    @test HTTP.URI(hostname="google.com", path="user") == HTTP.URI("http://google.com/user")
    @test HTTP.URI(hostname="google.com", query=Dict("key"=>"value")) == HTTP.URI("http://google.com?key=value")
    @test HTTP.URI(hostname="google.com", fragment="user") == HTTP.URI("http://google.com/#user")

    urls = [("hdfs://user:password@hdfshost:9000/root/folder/file.csv#frag", ["root", "folder", "file.csv"]),
            ("https://user:password@httphost:9000/path1/path2;paramstring?q=a&p=r#frag", ["path1", "path2;paramstring"]),
            ("https://user:password@httphost:9000/path1/path2?q=a&p=r#frag", ["path1","path2"]),
            ("https://user:password@httphost:9000/path1/path2;paramstring#frag", ["path1","path2;paramstring"]),
            ("https://user:password@httphost:9000/path1/path2#frag", ["path1","path2"]),
            ("file:///path/to/file/with%3fshould%3dwork%23fine", ["path","to","file","with%3fshould%3dwork%23fine"]),
            ("ftp://ftp.is.co.za/rfc/rfc1808.txt", ["rfc","rfc1808.txt"]),
            ("http://www.ietf.org/rfc/rfc2396.txt", ["rfc","rfc2396.txt"]),
            ("ldap://[2001:db8::7]/c=GB?objectClass?one", ["c=GB"]),
            ("mailto:John.Doe@example.com", ["John.Doe@example.com"]),
            ("news:comp.infosystems.www.servers.unix", ["comp.infosystems.www.servers.unix"]),
            ("tel:+1-816-555-1212", ["+1-816-555-1212"]),
            ("telnet://192.0.2.16:80/", String[]),
            ("urn:oasis:names:specification:docbook:dtd:xml:4.1.2", ["oasis:names:specification:docbook:dtd:xml:4.1.2"])
            ]

    for (url, splpath) in urls
        u = parse(HTTP.URI, url)
        @test string(u) == url
        @test isvalid(u)
        @test HTTP.splitpath(u) == splpath
    end

    @test parse(HTTP.URI, "hdfs://user:password@hdfshost:9000/root/folder/file.csv") == HTTP.URI(hostname="hdfshost", path="/root/folder/file.csv", scheme="hdfs", port=9000, userinfo="user:password")
    @test parse(HTTP.URI, "http://google.com:80/some/path") == HTTP.URI(hostname="google.com", path="/some/path")

    @test isempty(HTTP.URIs.Offset())
    @test HTTP.lower(UInt8('A')) == UInt8('a')
    @test HTTP.escape(Char(1)) == "%01"

    @test HTTP.escape(Dict("key1"=>"value1", "key2"=>["value2", "value3"])) == "key2=value2&key2=value3&key1=value1"

    @test HTTP.escape("abcdef αβ 1234-=~!@#\$()_+{}|[]a;") == "abcdef%20%CE%B1%CE%B2%201234-%3D%7E%21%40%23%24%28%29_%2B%7B%7D%7C%5B%5Da%3B"
    @test HTTP.unescape(HTTP.escape("abcdef 1234-=~!@#\$()_+{}|[]a;")) == "abcdef 1234-=~!@#\$()_+{}|[]a;"
    @test HTTP.unescape(HTTP.escape("👽")) == "👽"

    @test HTTP.escape([("foo", "bar"), (1, 2)]) == "foo=bar&1=2"
    @test HTTP.escape(Dict(["foo" => "bar", 1 => 2])) in ("1=2&foo=bar", "foo=bar&1=2")
    @test HTTP.escape(["foo" => "bar", 1 => 2]) == "foo=bar&1=2"

    @test "user:password" == HTTP.userinfo(parse(HTTP.URI, "https://user:password@httphost:9000/path1/path2;paramstring?q=a&p=r#frag"))

    @test false == isvalid(parse(HTTP.URI, "file:///path/to/file/with?should=work#fine"))
    @test true == isvalid( parse(HTTP.URI, "file:///path/to/file/with%3fshould%3dwork%23fine"))

    @test parse(HTTP.URI, "s3://bucket/key") == HTTP.URI(hostname="bucket", path="/key", scheme="s3")

    @test sprint(show, parse(HTTP.URI, "http://google.com")) == "HTTP.URI(\"http://google.com\")"

    # Error paths
    # Non-ASCII characters
    @test_throws HTTP.URIs.URLParsingError parse(HTTP.URI, "http://🍕.com")
    # Unexpected start of URL
    @test_throws HTTP.URIs.URLParsingError parse(HTTP.URI, ".google.com")
    # Unexpected character after scheme
    @test_throws HTTP.URIs.URLParsingError parse(HTTP.URI, "ht!tp://google.com")

    #  Issue #27
    @test HTTP.escape("t est\n") == "t%20est%0A"

    @testset "HTTP.parse(HTTP.URI, str)" begin

        urltests = URLTest[
        URLTest("proxy request"
         ,"http://hostname/"
         ,false
             ,(HTTP.URIs.Offset(1, 4) # UF_SCHEMA
             ,HTTP.URIs.Offset(8, 8) # UF_HOST
             ,HTTP.URIs.Offset(0, 0) # UF_PORT
             ,HTTP.URIs.Offset(16, 1) # UF_PATH
             ,HTTP.URIs.Offset(0, 0) # UF_QUERY
             ,HTTP.URIs.Offset(0, 0) # UF_FRAGMENT
             ,HTTP.URIs.Offset(0, 0) # UF_USERINFO
             )
         ,false
         ), URLTest("proxy request with port"
         ,"http://hostname:444/"
         ,false
             ,(HTTP.URIs.Offset(1, 4) # UF_SCHEMA
             ,HTTP.URIs.Offset(8, 8) # UF_HOST
             ,HTTP.URIs.Offset(17, 3) # UF_PORT
             ,HTTP.URIs.Offset(20, 1) # UF_PATH
             ,HTTP.URIs.Offset(0, 0) # UF_QUERY
             ,HTTP.URIs.Offset(0, 0) # UF_FRAGMENT
             ,HTTP.URIs.Offset(0, 0) # UF_USERINFO
             )
         ,false
         ), URLTest("CONNECT request"
         ,"hostname:443"
         ,true
             ,(HTTP.URIs.Offset(0, 0) # UF_SCHEMA
             ,HTTP.URIs.Offset(1, 8) # UF_HOST
             ,HTTP.URIs.Offset(10, 3) # UF_PORT
             ,HTTP.URIs.Offset(0, 0) # UF_PATH
             ,HTTP.URIs.Offset(0, 0) # UF_QUERY
             ,HTTP.URIs.Offset(0, 0) # UF_FRAGMENT
             ,HTTP.URIs.Offset(0, 0) # UF_USERINFO
             )
         ,false
         ), URLTest("proxy ipv6 request"
         ,"http://[1:2::3:4]/"
         ,false
             ,(HTTP.URIs.Offset(1, 4) # UF_SCHEMA
             ,HTTP.URIs.Offset(9, 8) # UF_HOST
             ,HTTP.URIs.Offset(0, 0) # UF_PORT
             ,HTTP.URIs.Offset(18, 1) # UF_PATH
             ,HTTP.URIs.Offset(0, 0) # UF_QUERY
             ,HTTP.URIs.Offset(0, 0) # UF_FRAGMENT
             ,HTTP.URIs.Offset(0, 0) # UF_USERINFO
             )
         ,false
         ), URLTest("proxy ipv6 request with port"
         ,"http://[1:2::3:4]:67/"
         ,false
             ,(HTTP.URIs.Offset(1, 4) # UF_SCHEMA
             ,HTTP.URIs.Offset(9, 8) # UF_HOST
             ,HTTP.URIs.Offset(19, 2) # UF_PORT
             ,HTTP.URIs.Offset(21, 1) # UF_PATH
             ,HTTP.URIs.Offset(0, 0) # UF_QUERY
             ,HTTP.URIs.Offset(0, 0) # UF_FRAGMENT
             ,HTTP.URIs.Offset(0, 0) # UF_USERINFO
             )
         ,false
         ), URLTest("CONNECT ipv6 address"
         ,"[1:2::3:4]:443"
         ,true
             ,(HTTP.URIs.Offset(0, 0) # UF_SCHEMA
             ,HTTP.URIs.Offset(2, 8) # UF_HOST
             ,HTTP.URIs.Offset(12, 3) # UF_PORT
             ,HTTP.URIs.Offset(0, 0) # UF_PATH
             ,HTTP.URIs.Offset(0, 0) # UF_QUERY
             ,HTTP.URIs.Offset(0, 0) # UF_FRAGMENT
             ,HTTP.URIs.Offset(0, 0) # UF_USERINFO
             )
         ,false
         ), URLTest("ipv4 in ipv6 address"
         ,"http://[2001:0000:0000:0000:0000:0000:1.9.1.1]/"
         ,false
             ,(HTTP.URIs.Offset(1, 4) # UF_SCHEMA
             ,HTTP.URIs.Offset(9,37) # UF_HOST
             ,HTTP.URIs.Offset(0, 0) # UF_PORT
             ,HTTP.URIs.Offset(47, 1) # UF_PATH
             ,HTTP.URIs.Offset(0, 0) # UF_QUERY
             ,HTTP.URIs.Offset(0, 0) # UF_FRAGMENT
             ,HTTP.URIs.Offset(0, 0) # UF_USERINFO
             )
         ,false
         ), URLTest("extra ? in query string"
         ,"http://a.tbcdn.cn/p/fp/2010c/??fp-header-min.css,fp-base-min.css,fp-channel-min.css,fp-product-min.css,fp-mall-min.css,fp-category-min.css,fp-sub-min.css,fp-gdp4p-min.css,fp-css3-min.css,fp-misc-min.css?t=20101022.css"
         ,false
             ,(HTTP.URIs.Offset(1, 4) # UF_SCHEMA
             ,HTTP.URIs.Offset(8,10) # UF_HOST
             ,HTTP.URIs.Offset(0, 0) # UF_PORT
             ,HTTP.URIs.Offset(18,12) # UF_PATH
             ,HTTP.URIs.Offset(31,187) # UF_QUERY
             ,HTTP.URIs.Offset(0, 0) # UF_FRAGMENT
             ,HTTP.URIs.Offset(0, 0) # UF_USERINFO
             )
         ,false
         ), URLTest("space URL encoded"
         ,"/toto.html?toto=a%20b"
         ,false
             ,(HTTP.URIs.Offset(0, 0) # UF_SCHEMA
             ,HTTP.URIs.Offset(0, 0) # UF_HOST
             ,HTTP.URIs.Offset(0, 0) # UF_PORT
             ,HTTP.URIs.Offset(1,10) # UF_PATH
             ,HTTP.URIs.Offset(12,10) # UF_QUERY
             ,HTTP.URIs.Offset(0, 0) # UF_FRAGMENT
             ,HTTP.URIs.Offset(0, 0) # UF_USERINFO
             )
         ,false
         ), URLTest("URL fragment"
         ,"/toto.html#titi"
         ,false
             ,(HTTP.URIs.Offset(0, 0) # UF_SCHEMA
             ,HTTP.URIs.Offset(0, 0) # UF_HOST
             ,HTTP.URIs.Offset(0, 0) # UF_PORT
             ,HTTP.URIs.Offset(1,10) # UF_PATH
             ,HTTP.URIs.Offset(0, 0) # UF_QUERY
             ,HTTP.URIs.Offset(12, 4) # UF_FRAGMENT
             ,HTTP.URIs.Offset(0, 0) # UF_USERINFO
             )
         ,false
         ), URLTest("complex URL fragment"
         ,"http://www.webmasterworld.com/r.cgi?f=21&d=8405&url=http://www.example.com/index.html?foo=bar&hello=world#midpage"
         ,false
         ,(HTTP.URIs.Offset(  1,  4) # UF_SCHEMA
          ,HTTP.URIs.Offset(  8, 22) # UF_HOST
          ,HTTP.URIs.Offset(  0,  0) # UF_PORT
          ,HTTP.URIs.Offset( 30,  6) # UF_PATH
          ,HTTP.URIs.Offset( 37, 69) # UF_QUERY
          ,HTTP.URIs.Offset(107,  7) # UF_FRAGMENT
          ,HTTP.URIs.Offset(  0,  0) # UF_USERINFO
          )
         ,false
         ), URLTest("complex URL from node js url parser doc"
         ,"http://host.com:8080/p/a/t/h?query=string#hash"
         ,false
         ,(   HTTP.URIs.Offset(1, 4) # UF_SCHEMA
             ,HTTP.URIs.Offset(8, 8) # UF_HOST
             ,HTTP.URIs.Offset(17, 4) # UF_PORT
             ,HTTP.URIs.Offset(21, 8) # UF_PATH
             ,HTTP.URIs.Offset(30,12) # UF_QUERY
             ,HTTP.URIs.Offset(43, 4) # UF_FRAGMENT
             ,HTTP.URIs.Offset(0, 0) # UF_USERINFO
             )
             ,false
         ), URLTest("complex URL with basic auth from node js url parser doc"
         ,"http://a:b@host.com:8080/p/a/t/h?query=string#hash"
         ,false
         ,(   HTTP.URIs.Offset(1, 4) # UF_SCHEMA
             ,HTTP.URIs.Offset(12, 8) # UF_HOST
             ,HTTP.URIs.Offset(21, 4) # UF_PORT
             ,HTTP.URIs.Offset(25, 8) # UF_PATH
             ,HTTP.URIs.Offset(34,12) # UF_QUERY
             ,HTTP.URIs.Offset(47, 4) # UF_FRAGMENT
             ,HTTP.URIs.Offset(8, 3) # UF_USERINFO
             )
            ,false
         ), URLTest("double @"
         ,"http://a:b@@hostname:443/"
         ,false
         ,true
         ), URLTest("proxy empty host"
         ,"http://:443/"
         ,false
         ,true
         ), URLTest("proxy empty port"
         ,"http://hostname:/"
         ,false
         ,true
         ), URLTest("CONNECT with basic auth"
         ,"a:b@hostname:443"
         ,true
         ,true
         ), URLTest("CONNECT empty host"
         ,":443"
         ,true
         ,true
         ), URLTest("CONNECT empty port"
         ,"hostname:"
         ,true
         ,true
         ), URLTest("CONNECT with extra bits"
         ,"hostname:443/"
         ,true
         ,true
         ), URLTest("space in URL"
         ,"/foo bar/"
         ,false
         ,true # s_dead
         ), URLTest("proxy basic auth with space url encoded"
         ,"http://a%20:b@host.com/"
         ,false
             ,(HTTP.URIs.Offset(1, 4) # UF_SCHEMA
              ,HTTP.URIs.Offset(15, 8) # UF_HOST
              ,HTTP.URIs.Offset(0, 0) # UF_PORT
              ,HTTP.URIs.Offset(23, 1) # UF_PATH
              ,HTTP.URIs.Offset(0, 0) # UF_QUERY
              ,HTTP.URIs.Offset(0, 0) # UF_FRAGMENT
              ,HTTP.URIs.Offset(8, 6) # UF_USERINFO
             )
         ,false
         ), URLTest("carriage return in URL"
         ,"/foo\rbar/"
         ,false
         ,true # s_dead
         ), URLTest("proxy double : in URL"
         ,"http://hostname::443/"
         ,false
         ,true # s_dead
         ), URLTest("proxy basic auth with double :"
         ,"http://a::b@host.com/"
         ,false
             ,(HTTP.URIs.Offset(1, 4) # UF_SCHEMA
             ,HTTP.URIs.Offset(13, 8) # UF_HOST
             ,HTTP.URIs.Offset(0, 0) # UF_PORT
             ,HTTP.URIs.Offset(21, 1) # UF_PATH
             ,HTTP.URIs.Offset(0, 0) # UF_QUERY
             ,HTTP.URIs.Offset(0, 0) # UF_FRAGMENT
             ,HTTP.URIs.Offset(8, 4) # UF_USERINFO
             )
         ,false
         ), URLTest("line feed in URL"
         ,"/foo\nbar/"
         ,false
         ,true # s_dead
         ), URLTest("proxy empty basic auth"
         ,"http://@hostname/fo"
         ,false
             ,(HTTP.URIs.Offset(1, 4) # UF_SCHEMA
             ,HTTP.URIs.Offset(9, 8) # UF_HOST
             ,HTTP.URIs.Offset(0, 0) # UF_PORT
             ,HTTP.URIs.Offset(17, 3) # UF_PATH
             ,HTTP.URIs.Offset(0, 0) # UF_QUERY
             ,HTTP.URIs.Offset(0, 0) # UF_FRAGMENT
             ,HTTP.URIs.Offset(0, 0) # UF_USERINFO
             )
         ,false
         ), URLTest("proxy line feed in hostname"
         ,"http://host\name/fo"
         ,false
         ,true # s_dead
         ), URLTest("proxy % in hostname"
         ,"http://host%name/fo"
         ,false
         ,true # s_dead
         ), URLTest("proxy ; in hostname"
         ,"http://host;ame/fo"
         ,false
         ,true # s_dead
         ), URLTest("proxy basic auth with unreservedchars"
         ,"http://a!;-_!=+\$@host.com/"
         ,false
             ,(HTTP.URIs.Offset(1, 4) # UF_SCHEMA
             ,HTTP.URIs.Offset(18, 8) # UF_HOST
             ,HTTP.URIs.Offset(0, 0) # UF_PORT
             ,HTTP.URIs.Offset(26, 1) # UF_PATH
             ,HTTP.URIs.Offset(0, 0) # UF_QUERY
             ,HTTP.URIs.Offset(0, 0) # UF_FRAGMENT
             ,HTTP.URIs.Offset(8, 9) # UF_USERINFO
             )
         ,false
         ), URLTest("proxy only empty basic auth"
         ,"http://@/fo"
         ,false
         ,true # s_dead
         ), URLTest("proxy only basic auth"
         ,"http://toto@/fo"
         ,false
         ,true # s_dead
         ), URLTest("proxy = in URL"
         ,"http://host=ame/fo"
         ,false
         ,true # s_dead
         ), URLTest("ipv6 address with Zone ID"
         ,"http://[fe80::a%25eth0]/"
         ,false
             ,(HTTP.URIs.Offset(1, 4) # UF_SCHEMA
             ,HTTP.URIs.Offset(9,14) # UF_HOST
             ,HTTP.URIs.Offset(0, 0) # UF_PORT
             ,HTTP.URIs.Offset(24, 1) # UF_PATH
             ,HTTP.URIs.Offset(0, 0) # UF_QUERY
             ,HTTP.URIs.Offset(0, 0) # UF_FRAGMENT
             ,HTTP.URIs.Offset(0, 0) # UF_USERINFO
             )
         ,false
         ), URLTest("ipv6 address with Zone ID, but '%' is not percent-encoded"
         ,"http://[fe80::a%eth0]/"
         ,false
             ,(HTTP.URIs.Offset(1, 4) # UF_SCHEMA
             ,HTTP.URIs.Offset(9,12) # UF_HOST
             ,HTTP.URIs.Offset(0, 0) # UF_PORT
             ,HTTP.URIs.Offset(22, 1) # UF_PATH
             ,HTTP.URIs.Offset(0, 0) # UF_QUERY
             ,HTTP.URIs.Offset(0, 0) # UF_FRAGMENT
             ,HTTP.URIs.Offset(0, 0) # UF_USERINFO
             )
         ,false
         ), URLTest("ipv6 address ending with '%'"
         ,"http://[fe80::a%]/"
         ,false
         ,true # s_dead
         ), URLTest("ipv6 address with Zone ID including bad character"
         ,"http://[fe80::a%\$HOME]/"
         ,false
         ,true # s_dead
         ), URLTest("just ipv6 Zone ID"
         ,"http://[%eth0]/"
         ,false
         ,true # s_dead
         ), URLTest("tab in URL"
         ,"/foo\tbar/"
         ,false
             ,(HTTP.URIs.Offset(0, 0) # UF_SCHEMA
             ,HTTP.URIs.Offset(0, 0) # UF_HOST
             ,HTTP.URIs.Offset(0, 0) # UF_PORT
             ,HTTP.URIs.Offset(1, 9) # UF_PATH
             ,HTTP.URIs.Offset(0, 0) # UF_QUERY
             ,HTTP.URIs.Offset(0, 0) # UF_FRAGMENT
             ,HTTP.URIs.Offset(0, 0) # UF_USERINFO
             )
         ,false
         ), URLTest("form feed in URL"
         ,"/foo\fbar/"
         ,false
             ,(HTTP.URIs.Offset(0, 0) # UF_SCHEMA
             ,HTTP.URIs.Offset(0, 0) # UF_HOST
             ,HTTP.URIs.Offset(0, 0) # UF_PORT
             ,HTTP.URIs.Offset(1, 9) # UF_PATH
             ,HTTP.URIs.Offset(0, 0) # UF_QUERY
             ,HTTP.URIs.Offset(0, 0) # UF_FRAGMENT
             ,HTTP.URIs.Offset(0, 0) # UF_USERINFO
             )
         ,false
         )
        ]

        for u in urltests
            println("TEST - uri.jl: $(u.name)")
            if u.shouldthrow
                @test_throws HTTP.URIs.URLParsingError parse(HTTP.URI, u.url; isconnect=u.isconnect)
            else
                url = parse(HTTP.URI, u.url; isconnect=u.isconnect)
                @test u.offsets == url.offsets
            end
        end
    end
end; # @testset


text:
abstract TOFilter

type ThresHoldDeviance <: TOFilter
    thres::Float64

    function ThresHoldDeviance(thres::Float64)
        if thres < 0.000001
            error("Specified threshold must be larger than zero.")
        end
        return new(thres)
    end
end

####################
## old code below ##
####################

################################
## filters to reduce turnover ##
################################

## all filters should be functions that receive an Investments object
## of intended portfolio weights as well as an associated Timematr
## object containing asset returns (in order to be able to derive
## automatic portfolio weight evolutions)
##
## function filterName(invs::Investments, tm::Timematr; opts="frequency")
##
## 	return filtWeights # invs object
## end

## various filters:
## - based on ranges (if deviation from desired wgts is too large,
## 	automatic trading will occur to bring weights back to
## 	permitted ranges)
## - based on rebalancing frequencies: rebalance each n days
## - based on overall distance to desired weights
## - combination of all

## function regularRB(invs::Investments, discRet::Timematr; freq=30)
##     ## rebalance each "freq" days

##     ## check inputs
##     AssetMgmt.chkMatchInvData(invs,discRet)

##     ## extract core values
##     wgts = AssetMgmt.core(invs)
##     rets = core(discRet)
##     (nObs, nAss) = size(wgts)

##     ## invRets = invRetCore(wgts, rets)

##     filteredWgts = zeros(nObs, nAss)
##     filteredWgts[1, :] = wgts[1, :]
##     for ii=2:nObs

##         ## is it time to rebalance?
##         timeToRebalance = false
##         if mod(ii, freq) == 0
##             timeToRebalance = true
##         end

##         if timeToRebalance
##             filteredWgts[ii, :] = wgts[ii, :]
##         else
##             ## get portfolio return
##             pfRet = sum(filteredWgts[ii-1, :] .* rets[ii-1, :])

##             filteredWgts[ii, :] = filteredWgts[ii-1, :] .*
##             (1 .+ rets[ii-1, :]) ./ (1 .+ pfRet)
##         end
##     end
##     wgtsDf = composeDataFrame(filteredWgts, names(discRet))

##     return Investments(wgtsDf, idx(invs))
## end

## function rangeRB(invs::Investments, discRet::Timematr; freq=30)

##     return Investments(wgtsDf, idx(invs))
## end

## function distThreshold

## end


text:
"""
    contourf(C::ClimGrid; region::String="auto", level, caxis, titlestr::String, cm::String="", ncolors::Int, center_cs::Bool, filename::String, cs_label::String)

Maps the time-mean average of ClimGrid C. If a filename is provided, the figure is saved in a png format.

Optional keyworkd includes precribed regions (keyword *region*, see list below). For 4D data, keyword *level* is used to map a given level (defaults to 1). *caxis* is used to limit the colorscale. *cm* is used to manually set the colorscale (see Python documentation for native colorscale keyword), *ncolors* is used to set the number of color classes (defaults to 12). Set *center_cs* to true to center the colorscale (useful for divergent results, such as anomalies, positive/negative temprature). *cs_label* is used for custom colorscale label.

## Arguments for keyword *region* (and shortcuts)
- Europe ("EU")
- NorthAmerica ("NA")
- Canada ("CA")
- Quebec, QuebecNSP ("QC", "QCNSP")
- Americas ("Ams")
- World, WorldAz, WorldEck4 ("W", "Waz", "Weck4")
- Greenwich ("Gr")

"""
function contourf(C::ClimGrid; region::String="auto", states::Bool=false,  rivers::Bool=false, borders::Bool=false, lakes::Bool=false, level=1, caxis=[], titlestr::String="", cm::String="", ncolors::Int=12, center_cs::Bool=false, filename::String="", cs_label::String="", dpi=200)

    # Build basic plot information
    PI = ClimatePlots.build_PlotInfo(C, region, :contourf, caxis, cm, states, rivers, borders, lakes, center_cs, ncolors, level)

    # Empty-map generator
    status, ax, projection = mapclimgrid(PI)

    cs = ax.contourf(C.longrid, C.latgrid, PI.data, transform=cartopy.crs.PlateCarree(), cmap=PI.cm, vmin=PI.vmin, vmax=PI.vmax)

    # Colorbar
    if isempty(cs_label)
        cs_label = getunitslabel(C)
    end
    cbar = colorbar(cs, orientation = "vertical", shrink = 0.7, label=cs_label)

    # Title
    if isempty(titlestr)
        titlestr = titledef(C)
    end
    title(titlestr)

    # Save to "filename" if not empty
    if !isempty(filename)
        PyPlot.savefig(filename, dpi=dpi)
    end

    return true, ax, projection, cbar
end

"""
    contour(C::ClimGrid; region::String="auto", level, caxis, titlestr::String, cm::String="", ncolors::Int, center_cs::Bool, filename::String, cs_label::String)

Maps the time-mean average of ClimGrid C. If a filename is provided, the figure is saved in a png format.

Optional keyworkd includes precribed regions (keyword *region*, see list below). For 4D data, keyword *level* is used to map a given level (defaults to 1). *caxis* is used to limit the colorscale. *cm* is used to manually set the colorscale (see Python documentation for native colorscale keyword), *ncolors* is used to set the number of color classes (defaults to 12). Set *center_cs* to true to center the colorscale (useful for divergent results, such as anomalies, positive/negative temprature). *cs_label* is used for custom colorscale label.

## Arguments for keyword *region* (and shortcuts)
- Europe ("EU")
- NorthAmerica ("NA")
- Canada ("CA")
- Quebec, QuebecNSP ("QC", "QCNSP")
- Americas ("Ams")
- World, WorldAz, WorldEck4 ("W", "Waz", "Weck4")
- Greenwich ("Gr")

"""
function contour(C::ClimGrid; region::String="auto", states::Bool=false,  rivers::Bool=false, borders::Bool=false, lakes::Bool=false, level=1, caxis=[], titlestr::String="", cm::String="", ncolors::Int=12, center_cs::Bool=false, filename::String="", cs_label::String="", dpi=200)

    # Build basic plot information
    PI = ClimatePlots.build_PlotInfo(C, region, :contour, caxis, cm, states, rivers, borders, lakes, center_cs, ncolors, level)

    # Empty-map generator
    status, ax, projection = mapclimgrid(PI)

    cs = ax.contour(C.longrid, C.latgrid, PI.data, transform=cartopy.crs.PlateCarree(), cmap=PI.cm, vmin=PI.vmin, vmax=PI.vmax)

    # Colorbar
    if isempty(cs_label)
        cs_label = getunitslabel(C)
    end
    cbar = colorbar(cs, orientation = "vertical", shrink = 0.7, label=cs_label)

    # Title
    if isempty(titlestr)
        titlestr = titledef(C)
    end
    title(titlestr)

    # Save to "filename" if not empty
    if !isempty(filename)
        PyPlot.savefig(filename, dpi=dpi)
    end

    return true, ax, projection, cbar

end

"""
    pcolormesh(C::ClimGrid; region::String="auto", level, caxis, titlestr::String, cm::String="", ncolors::Int, center_cs::Bool, filename::String, cs_label::String)

Maps the time-mean average of ClimGrid C. If a filename is provided, the figure is saved in a png format.

Optional keyworkd includes precribed regions (keyword *region*, see list below). For 4D data, keyword *level* is used to map a given level (defaults to 1). *caxis* is used to limit the colorscale. *cm* is used to manually set the colorscale (see Python documentation for native colorscale keyword), *ncolors* is used to set the number of color classes (defaults to 12). Set *center_cs* to true to center the colorscale (useful for divergent results, such as anomalies, positive/negative temprature). *cs_label* is used for custom colorscale label.

## Arguments for keyword *region* (and shortcuts)
- Europe ("EU")
- NorthAmerica ("NA")
- Canada ("CA")
- Quebec, QuebecNSP ("QC", "QCNSP")
- Americas ("Ams")
- World, WorldAz, WorldEck4 ("W", "Waz", "Weck4")
- Greenwich ("Gr")

"""
function pcolormesh(C::ClimGrid; region::String="auto", states::Bool=false,  rivers::Bool=false, borders::Bool=false, lakes::Bool=false, level=1, caxis=[], titlestr::String="", cm::String="", ncolors::Int=12, center_cs::Bool=false, filename::String="", cs_label::String="", dpi=200)

    # Build basic plot information
    PI = ClimatePlots.build_PlotInfo(C, region, :pcolormesh, caxis, cm, states, rivers, borders, lakes, center_cs, ncolors, level)

    # Empty-map generator
    status, ax, projection = mapclimgrid(PI)

    cs = ax.pcolormesh(C.longrid, C.latgrid, PI.data, transform=cartopy.crs.PlateCarree(), cmap=PI.cm, vmin=PI.vmin, vmax=PI.vmax)

    # Colorbar
    if isempty(cs_label)
        cs_label = getunitslabel(C)
    end
    cbar = colorbar(cs, orientation = "vertical", shrink = 0.9, ticklocation="auto", label=cs_label)

    # Title
    if isempty(titlestr)
        titlestr = titledef(C)
    end
    title(titlestr)

    # Save to "filename" if not empty
    if !isempty(filename)
        PyPlot.savefig(filename, dpi=dpi)
    end

    return true, ax, projection, cbar
end


text:
using Test
x = GRASP.QAProgram(GRASP.binary_rnucleus, ["y","n"])

GRASP.run!(x)


text:
export se, SE

# Tr0, R0 = SA[0.0; 0], SA[1.0 0; 0 1]
# p0 = ProductRepr(Tr0, R0)

# t0, w0 = [0.0; 0],  get_vector(M_R, R0, 0, B)
# x0 = ProductRepr(t0, w0)
# p0_ = exp(G, p0, x0)  # or,  exp(M, p0, x0)

# t1a, w1a = [1.0; 0],  get_vector(M_R, R0, 0  , B)
# x1a = ProductRepr(t1a, w1a)   # Lie algebra element
# p1a = exp(G, p0, x1a)         # Lie group element

"""
    $TYPEDEF

Lie Group representation for `Manifolds.SpecialEuclidean(N)`

Related

[`se`](@ref)
"""
struct SE{N, T}
  p::T # Group
end

function SE(R::Matrix{T}, t::Vector{T}) where T<:Real
    N = length(t)
    @assert size(R,1) == size(R,2) == N

    _t = SVector{N}(t)
    _R = SMatrix{N,N}(R)

    p = ProductRepr(_t, _R)

    return SE{N, typeof(P)}(p)
end

manifold(::SE{N,T}) where {N,T} = SpecialEuclidean(N)

function _identity_SE(N::Int, ::Type{T}=Float64) where {T<:Real}
  p1 = zeros(SVector{N, T})
  p2 = SMatrix{N,N,T}(one(T)I)
  o = ProductRepr(p1, p2)
  SE{N, typeof(o)}(o)
end

"""
    $TYPEDEF

Lie algebra representation for `Manifolds.SpecialEuclidean(N)`

Related

[`SE`](@ref)
"""
struct se{N, T}
  X::T # Tangent Vector (Lie algebra, skew simetric matrix)
end

function se(X::Vector)
    N = length(X)

    for D in 2:3
        if manifold_dimension(SpecialEuclidean(D)) == N
            M = SpecialEuclidean(D)
            e = _identity_SE(D, eltype(X)).p
            _X = hat(M, e, X)

            return se{D, typeof(_X)}(_X)
        end
    end
    error("Only se(2) and se(3) suppored")
end

manifold(::se{N,T}) where {N,T} = SpecialEuclidean(N)

# p0 = SE([1.0 0.0; 0 1.0], [0.0, 0.0])

function ManifoldsBase.vee(X::se)
    M = manifold(se)
    return vee(M, Identity(M, X.X), X.X)
end

# we must use
# seVec = Manifolds.hat(M, coords)

function ⊕(p::SE, X::ProductRepr)
    M = manifold(p)
    e = identity(M, p.p)
    typeof(p)(compose(M, p.p, exp(M, e, X)))
end

⊕(X::ProductRepr, p::SE) =  typeof(p)(compose(manifold(p), exp(manifold(p), I_SO2, X), p.p))

# I_SO2 = ProductRepr(SA[0.0, 0.0], SA[1.0 0.0; 0 1.0]);

# M = SpecialEuclidean(2)
# X1 = hat(M, I_SO2, SA[1., 0, pi/2]);
# X2 = hat(M, I_SO2, SA[1., 0, -pi/2]);

# p1 = p0 ⊕ X1

# p2 = p1 ⊕ X

##

#


text:
module TestDiscreteRv

using QuantEcon
using Base.Test
using FactCheck
using DataStructures

# set up
srand(42)
n = 10
x = rand(n)
x ./= sum(x)
drv = DiscreteRV(x)

facts("Testing discrete_rv.jl") do
    # test Q sums to 1
    @fact drv.Q[end] --> roughly(1.0)

    # test lln
    draws = draw(drv, 100000)
    c = counter(draws)
    counts = Array(Float64, n)
    for i=1:n
        counts[i] = c[i]
    end
    counts ./= sum(counts)

    @fact Base.maxabs(counts - drv.q) --> roughly(0.0; atol=1e-2)

end  # facts
end  # module


text:
module Test_pscf

using DescriptorSystems
using LinearAlgebra
using Polynomials
using Test

println("Test_pscf")
@testset "grcf/gnrcf and glcf/gnlcf" begin

@testset "grcf/gnrcf" begin

sys = rdss(0,0,0);
@time sysn, sysm = grcf(sys);
@test gnrank(sys*sysm-sysn) == 0   &&   #  G(s)*M(s)-N(s) = 0
      isproper(sysm) & isproper(sysn) && # checking properness of factors
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7)))

a = -(1. +2im)*ones(1,1); e = (1. -im)*ones(1,1); b = (2. +im)*ones(1,1);
c = (-1. +im)*ones(1,1); d = complex(0)*ones(1,1);
sys = dss(a,e,b,c,d);
@time sysn, sysm = grcf(sys,atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysm-sysn,atol=1.e-7) == 0 &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7)))

a = -(1. +2im)*ones(1,1); e = (1. -im)*ones(1,1); b = (2. +im)*ones(1,1);
c = (-1. +im)*ones(1,1); d = complex(0)*ones(1,1);
sys = dss(a,e,b,c,d,Ts=1);
@time sysn, sysm = grcf(sys,atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysm-sysn,atol=1.e-7) == 0 &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7)))

s = Polynomial([0, 1],'s'); z = Polynomial([0, 1],'z');     # define the complex variables s and z

# Gc = [s^2 s/(s+1); 0 1/s]     # define the 2-by-2 improper Gc(s)
Nc = [s^2 s; 0 1]; Dc = [1 (s+1); 1 s];
sysc = dss(Nc,Dc);             # build continuous-time descriptor system realization
@time sysn, sysm = grcf(sysc, evals = [-2,-3,-4], sdeg = -0.99, mindeg = true, mininf = true);
@test gnrank(sysc*sysm-sysn,atol=1.e-7) == 0 &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7))) &&
      sort(gpole(sysm)) ≈ sort([-2,-3,-4])

#test pole assignment feature
sysc = dss(A = diagm(rand(3)),B = rand(3,2), C = rand(2,3));
@time sysn, sysm = grcf(sysc, evals = [-2+im,-4,-2-im], mindeg = true, mininf = true);
@test gnrank(sysc*sysm-sysn,atol=1.e-7) == 0 &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7))) &&
      sort(real(gpole(sysm))) ≈ sort([-2,-2,-4]) &&
      sort(imag(gpole(sysm))) ≈ sort([1,-1,0])

sysc = sysn'; gpole(sysc)
@time sysn, sysm = grcf(sysc, evals = [-2+im,-4,-2-im], mindeg = true, mininf = true);
@test gnrank(sysc*sysm-sysn,atol=1.e-7) == 0 &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7))) &&
      sort(real(gpole(sysm))) ≈ sort([-2,-2,-4]) &&
      sort(imag(gpole(sysm))) ≈ sort([1,-1,0])

e = rand(3,3)
sysc = dss(A = e*diagm(rand(3)), E = e, B = rand(3,2), C = rand(2,3));

@time sysn, sysm = grcf(sysc, evals = [-2+im,-4,-2-im], mindeg = true, mininf = true);
@test gnrank(sysc*sysm-sysn,atol=1.e-7) == 0 &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7))) &&
      sort(real(gpole(sysm))) ≈ sort([-2,-2,-4]) &&
      sort(imag(gpole(sysm))) ≈ sort([1,-1,0])

sysc = sysn';
@time sysn, sysm = grcf(sysc, evals = [-2+im,-4,-2-im], mindeg = true, mininf = true);
@test gnrank(sysc*sysm-sysn,atol=1.e-7) == 0 &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7))) &&
      sort(real(gpole(sysm))) ≈ sort([-2,-2,-4]) &&
      sort(imag(gpole(sysm))) ≈ sort([1,-1,0])

# enforce interchange of 1x1 and 2x2 blocks
sysc = dss(A = [1 100 0 1;0 2 3 1; 0 -3 2 0; 0 0 0 3], B = rand(4,2), C = rand(2,4) );
@time sysn, sysm = grcf(sysc, evals = [-2+im,-4+2im,-4-2im,-2-im], mindeg = true, mininf = true);
@test gnrank(sysc*sysm-sysn,atol=1.e-7) == 0 &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7))) &&
      sort(real(gpole(sysm))) ≈ sort([-2,-2,-4,-4]) &&
      sort(imag(gpole(sysm))) ≈ sort([1,-1,2,-2])

sysc = dss(A = [1 100 0 1;0 2 3 1; 0 -3 2 0; 0 0 0 3], E = 2*eye(4),  B = rand(4,2), C = rand(2,4) );
@time sysn, sysm = grcf(sysc, evals = [-2+im,-4+2im,-4-2im,-2-im], mindeg = true, mininf = true);
@test gnrank(sysc*sysm-sysn,atol=1.e-7) == 0 &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7))) &&
      sort(real(gpole(sysm))) ≈ sort([-2,-2,-4,-4]) &&
      sort(imag(gpole(sysm))) ≈ sort([1,-1,2,-2])

# Gd = [z^2 z/(z-2); 0 1/z]     # define the 2-by-2 improper Gd(z)
Nd = [z^2 z; 0 1]; Dd = [1 z-2; 1 z];
# build LTI (nonminimal) descriptor realizations of Gc(s) and Gd(z)
sysd = dss(Nd,Dd,Ts = 1);      # build discrete-time descriptor system realization

@time sysn, sysm = grcf(sysd);
@test gnrank(sysd*sysm-sysn,atol=1.e-7) == 0  &&    #  Gd(z)*M(z)-N(z) = 0
      isproper(sysm) & isproper(sysn)  &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7)))

s = Polynomial([0, 1],'s');
gn = [(s-1) s 1;
    0 (s-2) (s-2);
    (s-1) (s^2+2*s-2) (2*s-1)];
gd = [(s+2) (s+2) (s+2);
    1 (s+1)^2 (s+1)^2;
    (s+2) (s+1)*(s+2) (s+1)*(s+2)];
sys = dss(gn,gd,minimal = true, atol = 1.e-7);

@time sysn, sysm = grcf(sys,mindeg=true);
@test gnrank(sys*sysm-sysn,atol=1.e-7) == 0  &&
      isproper(sysm) & isproper(sysn)  &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7))) &&
      order(sysm) == 0

@time sysn, sysm = grcf(sys',sdeg=-0.5);
@test gnrank(sys'*sysm-sysn,atol=1.e-7) == 0  &&
      isproper(sysm) & isproper(sysn)  &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7))) &&
      order(sysm) == 4

s = Polynomial([0, 1],'s');
gn = [(s-1) s 1;
    0 (s-2) (s-2);
    (s-1) (s^2+2*s-2) (2*s-1)];
gd = [(s-2) (s+2) (s+2);
    1 (s+1)^2 (s+1)^2;
    (s+2) (s-1)*(s+2) (s+1)*(s-2)];
sys = dss(gn,gd,minimal = true, atol = 1.e-7);

@time sysn, sysm = grcf(sys,mindeg = true, atol = 1.e-7);
@test gnrank(sys*sysm-sysn,atol=1.e-7) == 0  &&
      isproper(sysm) & isproper(sysn)  &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7))) &&
      order(sysm) == 3

@time sysn, sysm = gnrcf(sys,atol = 1.e-7);
@test gnrank(sys*sysm-sysn,atol=1.e-7) == 0  &&
      isproper(sysm) & isproper(sysn)  &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7))) &&
      iszero(sysn'*sysn+sysm'*sysm-I,atol=1.e-7)

z = Polynomial([0, 1],'z')
g = [z^2+z+1 4*z^2+3*z+2 2*z^2-2;
    z 4*z-1 2*z-2;
    z^2 4*z^2-z 2*z^2-2*z];
sys = dss(g,minimal = true, atol = 1.e-7,Ts = 1);

sysn, sysm = grcf(sys,mindeg = true,atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysm-sysn,atol=1.e-7) == 0  &&    #  Gd(z)*M(z)-N(z) = 0
      isproper(sysm) & isproper(sysn)  &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7)))

@time sysn, sysm = gnrcf(sys, ss = true, atol = 1.e-7);
@test gnrank(sys*sysm-sysn,atol=1.e-7) == 0  &&
      isproper(sysm) & isproper(sysn)  &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7))) &&
      iszero(sysn'*sysn+sysm'*sysm-I,atol=1.e-7)

z = Polynomial([0, 1],'z')
# Eigenvalue(s) on the unit circle
gn = 1
gd = z-1
sys = dss(gn,gd,Ts = 1)
sysn, sysm = grcf(sys)
@test gnrank(sys*sysm-sysn,atol=1.e-7) == 0  &&    #  Gd(z)*M(z)-N(z) = 0
      isproper(sysm) & isproper(sysn)  &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7)))

fast = true; Ty = Complex{Float64}; Ty = Float64
n = 5; p = 3; m = 2;

for fast in (true, false)
# random examples
for Ty in (Float64, Complex{Float64})

# continuous, standard
sys = rss(n,p,m,T = Ty,disc=false);
@time sysn, sysm = grcf(sys, smarg = -1, fast = fast,atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysm-sysn,atol=1.e-7) == 0   &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,-0.99, atol=1.e-7) &&
      isstable(sysn,-0.99,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7),atol=1.e-7))

#fail
sys = rss(n,p,m,T = Ty,disc=false);
@time sysn, sysm = grcf(sys, smarg = -1, evals = -Vector(1:n),fast = fast,atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysm-sysn,atol=1.e-7) == 0   &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,-0.99, atol=1.e-7) &&
      isstable(sysn,-0.99,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7),atol=1.e-7))

# discrete, standard
sys = rss(n,p,m,T = Ty,disc=true);
@time sysn, sysm = grcf(sys, smarg = 0.8, sdeg = 0.5, fast = fast,atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysm-sysn,atol=1.e-7) == 0   &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,0.81,atol=1.e-7) &&
      isstable(sysn,0.81,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7),atol=1.e-7))

# discrete, standard
sys = rss(n,p,m,T = Ty,disc=true);
@time sysn, sysm = grcf(sys, evals = -Vector(1:n)/10, fast = fast,atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysm-sysn,atol=1.e-7) == 0   &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7),atol=1.e-7))

# continuous, standard, uncontrollable
sys = rss(n,p,m,T = Ty,disc=false, nuc = 5);
@time sysn, sysm = grcf(sys, fast = fast,atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysm-sysn,atol=1.e-7) == 0   &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7),atol=1.e-7))

# discrete, standard, uncontrollable
sys = rss(n,p,m,T = Ty,disc=true, nuc = 5);
@time sysn, sysm = grcf(sys, fast = fast,atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysm-sysn,atol=1.e-7) == 0   &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7),atol=1.e-7))

# continuous, descriptor, no infinite eigenvalues
sys = rdss(n,p,m,T = Ty,disc=false);
@time sysn, sysm = grcf(sys, fast = fast,atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysm-sysn,atol=1.e-7) == 0   &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7)))

# discrete, descriptor, no infinite eigenvalues
sys = rdss(n,p,m,T = Ty, disc=true);
@time sysn, sysm = grcf(sys, fast = fast,atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysm-sysn,atol=1.e-7) == 0   &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7)))

# continuous, descriptor, proper
sys = rdss(n,p,m,T = Ty,disc=false,id=ones(Int,3));
@time sysn, sysm = grcf(sys, fast = fast,atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysm-sysn,atol=1.e-7) == 0   &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7)))

# discrete, descriptor, proper
sys = rdss(n,p,m,T = Ty, disc=true);
@time sysn, sysm = grcf(sys, fast = fast,atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysm-sysn,atol=1.e-7) == 0   &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7)))

# continuous, descriptor, proper
sys = rdss(n,p,m,T = Ty,disc=false,id=ones(Int,3));
@time sysn, sysm = grcf(sys, fast = fast, mindeg = true, mininf = true, atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysm-sysn,atol=1.e-7) == 0   &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7)))

# discrete, descriptor, proper
sys = rdss(n,p,m,T = Ty, disc=true,id=ones(Int,3));
@time sysn, sysm = grcf(sys, fast = fast,mindeg = true, mininf = true, atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysm-sysn,atol=1.e-7) == 0   &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7)))

# continuous, descriptor, infinite poles
sys = rdss(n,p,m,T = Ty,disc=false,id=[3*ones(Int,1);2*ones(Int,1)]);
@time sysn, sysm = grcf(sys, fast = fast, mindeg = true, mininf = true, atol=1.e-7,atol3=1.e-7);
#@time sysn, sysm = grcf(sys, fast = fast, mindeg = false, mininf = false, atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysm-sysn,atol=1.e-7) == 0   &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7)))

# discrete, descriptor, infinite poles
sys = rdss(n,p,m,T = Ty, disc=true,id=[3*ones(Int,1);2*ones(Int,1)]);
@time sysn, sysm = grcf(sys, fast = fast,mindeg = true, mininf = true, atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysm-sysn,atol=1.e-7) == 0   &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7)))

# continuous, descriptor, proper, uncontrollable infinite eigenvalues
sys = rdss(n,p,m,T = Ty,disc=false,iduc=[3*ones(Int,1);2*ones(Int,1)]);
@time sysn, sysm = grcf(sys, fast = fast, mindeg = true, mininf = true, atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysm-sysn,atol=1.e-7) == 0   &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7)))

# discrete, descriptor, proper, uncontrollable infinite eigenvalues
sys = rdss(n,p,m,T = Ty,disc=true,iduc=[3*ones(Int,1);2*ones(Int,1)]);
@time sysn, sysm = grcf(sys, fast = fast, mindeg = true, mininf = true, atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysm-sysn,atol=1.e-7) == 0   &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7)))

# continuous, descriptor, proper, uncontrollable finite eigenvalues
sys = rdss(n,p,m,T = Ty,disc=false, nfuc = 5);
@time sysn, sysm = grcf(sys, fast = fast, mindeg = true, mininf = true, atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysm-sysn,atol=1.e-7) == 0   &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7)))

# discrete, descriptor, proper, uncontrollable finite eigenvalues
sys = rdss(n,p,m,T = Ty,disc=true, nfuc = 5);
@time sysn, sysm = grcf(sys, fast = fast, mindeg = true, mininf = true, atol=1.e-7,atol3=1.e-7);
@test gnrank(sys*sysm-sysn,atol=1.e-7) == 0   &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7)))

end
end

end # grcf/gnrcf

@testset "glcf/gnlcf" begin

sys = rdss(0,0,0);
@time sysn, sysm = glcf(sys);
@test gnrank(sysm*sys-sysn) == 0   &&   #  M(s)*G(s)-N(s) = 0
      isproper(sysm) & isproper(sysn) && # checking properness of factors
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn sysm])))

s = Polynomial([0, 1],'s');
gn = [(s-1) s 1;
    0 (s-2) (s-2);
    (s-1) (s^2+2*s-2) (2*s-1)];
gd = [(s-2) (s+2) (s+2);
    1 (s+1)^2 (s+1)^2;
    (s+2) (s-1)*(s+2) (s+1)*(s-2)];
sys = dss(gn,gd,minimal = true, atol = 1.e-7);

@time sysn, sysm = glcf(sys,mindeg = true, atol = 1.e-7);
@test gnrank(sysm*sys-sysn,atol=1.e-7) == 0  &&
      isproper(sysm) & isproper(sysn)  &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn sysm],atol=1.e-7))) &&
      order(sysm) == 3

@time sysn, sysm = gnlcf(sys,atol = 1.e-7);
@test gnrank(sysm*sys-sysn,atol=1.e-7) == 0  &&
      isproper(sysm) & isproper(sysn)  &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn sysm],atol=1.e-7))) &&
      iszero(sysn*sysn'+sysm*sysm'-I,atol=1.e-7)

z = Polynomial([0, 1],'z')
g = [z^2+z+1 4*z^2+3*z+2 2*z^2-2;
    z 4*z-1 2*z-2;
    z^2 4*z^2-z 2*z^2-2*z];
sys = dss(g,minimal = true, atol = 1.e-7,Ts = 1);

sysn, sysm = glcf(sys,mindeg = true,atol=1.e-7,atol3=1.e-7);
@test gnrank(sysm*sys-sysn,atol=1.e-7) == 0  &&    #  Gd(z)*M(z)-N(z) = 0
      isproper(sysm) & isproper(sysn)  &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn sysm],atol=1.e-7)))

@time sysn, sysm = gnlcf(sys, ss = true, atol = 1.e-7);
@test gnrank(sysm*sys-sysn,atol=1.e-7) == 0  &&
      isproper(sysm) & isproper(sysn)  &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn sysm],atol=1.e-7))) &&
      iszero(sysn*sysn'+sysm*sysm'-I,atol=1.e-7)

z = Polynomial([0, 1],'z')
# Eigenvalue(s) on the unit circle
gn = 1
gd = z-1
sys = dss(gn,gd,Ts = 1)
sysn, sysm = grcf(sys)
@test gnrank(sys*sysm-sysn,atol=1.e-7) == 0  &&    #  Gd(z)*M(z)-N(z) = 0
      isproper(sysm) & isproper(sysn)  &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn;sysm],atol=1.e-7)))

fast = true; Ty = Complex{Float64}; #Ty = Float64
n = 5; p = 3; m = 2;
for fast in (true, false)
# random examples
for Ty in (Float64, Complex{Float64})

# continuous, standard
sys = rss(n,p,m,T = Ty,disc=false);
@time sysn, sysm = glcf(sys, fast = fast,atol=1.e-7,atol3=1.e-7);
@test gnrank(sysm*sys-sysn,atol=1.e-7) == 0   &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn sysm],atol=1.e-7),atol=1.e-7))

# discrete, standard
sys = rss(n,p,m,T = Ty,disc=true);
@time sysn, sysm = glcf(sys, fast = fast,atol=1.e-7,atol3=1.e-7);
@test gnrank(sysm*sys-sysn,atol=1.e-7) == 0   &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn sysm],atol=1.e-7),atol=1.e-7))

# continuous, standard, unobservable
sys = rss(n,p,m,T = Ty,disc=false, nuo = 5);
@time sysn, sysm = glcf(sys, fast = fast,atol=1.e-7,atol3=1.e-7);
@test gnrank(sysm*sys-sysn,atol=1.e-7) == 0   &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn sysm],atol=1.e-7),atol=1.e-7))

# discrete, standard, unobservable
sys = rss(n,p,m,T = Ty,disc=true, nuo = 5);
@time sysn, sysm = glcf(sys, fast = fast,atol=1.e-7,atol3=1.e-7);
@test gnrank(sysm*sys-sysn,atol=1.e-7) == 0   &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn sysm],atol=1.e-7),atol=1.e-7))

# continuous, descriptor, proper, unobservable infinite eigenvalues
sys = rdss(n,p,m,T = Ty,disc=false,iduo=[3*ones(Int,1);2*ones(Int,1)]);
@time sysn, sysm = glcf(sys, fast = fast, mindeg = true, mininf = true, atol=1.e-7,atol3=1.e-7);
@test gnrank(sysm*sys-sysn,atol=1.e-7) == 0   &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn sysm])))

# discrete, descriptor, proper, unobservable infinite eigenvalues
sys = rdss(n,p,m,T = Ty,disc=true,iduo=[3*ones(Int,1);2*ones(Int,1)]);
@time sysn, sysm = glcf(sys, fast = fast, mindeg = true, mininf = true, atol=1.e-7,atol3=1.e-7);
@test gnrank(sysm*sys-sysn,atol=1.e-7) == 0   &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn sysm])))

# continuous, descriptor, proper, unobservable finite eigenvalues
sys = rdss(n,p,m,T = Ty,disc=false, nfuo = 5);
@time sysn, sysm = glcf(sys, fast = fast, mindeg = true, mininf = true, atol=1.e-7,atol3=1.e-7);
@test gnrank(sysm*sys-sysn,atol=1.e-7) == 0   &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn sysm])))

# discrete, descriptor, proper, unobservable finite eigenvalues
sys = rdss(n,p,m,T = Ty,disc=true, nfuo = 5);
@time sysn, sysm = glcf(sys, fast = fast, mindeg = true, mininf = true, atol=1.e-7,atol3=1.e-7);
@test gnrank(sysm*sys-sysn,atol=1.e-7) == 0   &&
      isproper(sysm) && isproper(sysn) &&
      isstable(sysm,atol=1.e-7) &&
      isstable(sysn,atol=1.e-7)  &&
      isempty(gzero(gminreal([sysn sysm])))

end
end

end # glcf/gnlcf

end #test

end #module


text:
using Documenter, BackedUpImmutable

makedocs(sitename = "BackedUpImmutable Module Documentation", format = Documenter.HTML(prettyurls = false))


text:
using Plots

function Plots.plot(m::marginal{Ω, <: SmoothCDF} where Ω)
  plot(x -> pdf(m, x), quantile(m, .01), quantile(m, 0.99))
end


text:
export simple

function simple(a::Any, depth=0)
    a
end

function simple(arr::AbstractArray, depth=0)
    out = string("[", join([simple(x, depth+1) for x in arr], ", "), "]")
    0 == depth ?  println(out) : out
end

function simple(dict::Dict, depth=0)
    output = []
    for (k, v) in dict
        push!(output, string(simple(k, simple), ": ", simple(v, depth+1)))
    end
    out = string("{", join(output, ", "), "}")
    0 == depth ?  println(out) : out
end


text:
using Distributions
using Test, Random

# Test the consistency between the recursive and nonrecursive computation of the pdf
# of the Binomial distribution
Random.seed!(1234)
for (p, n) in [(0.6, 10), (0.8, 6), (0.5, 40), (0.04, 20), (1., 100), (0., 10), (0.999999, 1000), (1e-7, 1000)]
    local p

    d = Binomial(n, p)

    a = pdf.(d, 0:n)
    for t=0:n
        @test pdf(d, t) ≈ a[1+t]
    end

    li = rand(0:n, 2)
    rng = minimum(li):maximum(li)
    b = pdf.(d, rng)
    for t in rng
        @test pdf(d, t) ≈ b[t - first(rng) + 1]
    end

end

# Test calculation of expectation value for Binomial distribution
@test Distributions.expectation(Binomial(6), identity) ≈ 3.0
@test Distributions.expectation(Binomial(10, 0.2), x->-x) ≈ -2.0