text:
"""
    store(prop::Properties, filename, comment::Union{String,Nothing})

Store the property-dictionary in a text file called `filename`.
Additional comment can be added, that will be written to the beginning of the file.
"""
function store(prop::Properties, filename, comment::Union{String,Nothing})
    open(abspath(filename), "w") do io
        # add comment
        if comment != nothing
            # if multiline comment, place # before every line
            comment = replace(comment, "\n" => "\n#")

            print(io, "#")
            println(io, comment)
        end
        # add date every time
        println(io, "#", now())
        # write every key-value pair
        for v in keys(prop.properties)
            println(io, v, "=" , prop.properties[v])
        end
    end
    nothing
end

"""
    store(prop::Properties, filename)

Store the property-dictionary in a text file called `filename`.
"""
function store(prop::Properties, filename)
    store(prop, filename, nothing)
    nothing
end

"""
    load(filename)

Load `filename` propertyfile into a property-dictionary.
"""
function load(filename)
    p = Dict{String,String}()
    lines = readlines(abspath(filename))
    for line in lines
        # skip comments
        startswith(line, "#") && continue
        #split by "="
        splitted = split(line, "=")
        p[splitted[1]] = splitted[2]
    end
    Properties(p)
end


text:
doc = """PCL visualizer

Usage:
  pcl_visualizer.jl <pcd_file> [options]
  pcl_visualizer.jl -h | --help
  pcl_visualizer.jl --version

Options:
  -h --help       Show this screen.
  --version       Show version.
  --point_type=T  Point type [default: PointXYZ]
"""

using DocOpt
using PCL
using Cxx

let
    args = docopt(doc, version=v"0.0.1")
    pcd_file = args["<pcd_file>"]

    # e.g. pcl.PointXYZ
    T = eval(Expr(:., :pcl, QuoteNode(Symbol(args["--point_type"]))))

    cloud = pcl.PointCloud{T}(pcd_file)

    viewer = pcl.PCLVisualizer("pcl visualizer")
    pcl.addPointCloud(viewer, cloud)

    pcl.run(viewer)
end

ext: jl


text:
export DMD_SPG_Options, solveDMD_withSPG

using StatsBase, Printf

# SPG solver for large scale DMD problem
mutable struct DMD_SPG_Options
    tau::Integer	# sample size
    eta::AbstractFloat # step size
    dms::Integer # dimishing speed
    itm::Integer
    tol::AbstractFloat
    ptf::Integer
    true_obj::Bool
    prox::Function
end

# constructor of the SPG options
function DMD_SPG_Options(tau, eta; dms=500, itm=1000, tol=1e-5, ptf=100,
                         true_obj=false,
                         prox=prox_null::Function)

    return DMD_SPG_Options(tau, eta, dms, itm, tol, ptf, true_obj, prox)
end

# SPG solver
function solveDMD_withSPG(params, opts)
    # load all the variables
    tau    = opts.tau;
    eta    = opts.eta;
    dms    = opts.dms;
    prox = opts.prox;
    n    = params.n;
    k    = params.k;
    T    = typeof(real(params.X[1]));

    ar   = params.ar;
    arold = copy(ar);
    prox(ar);
    dar  = zeros(T, 2 * k);
    ind  = collect(1:n);

    fars = zeros(T, n)
    Gars = zeros(T, 2 * k, n)
    gars = Array{typeof(view(Gars, :, 1))}(undef, n)
    for i = 1:n
        gars[i] = view(Gars, :, i)
    end

    tfar = T(0.0); tgar = zeros(T, 2 * k);
    dfar = T(0.0); dgar = zeros(T, 2 * k);
    rfar = T(0.0); rgar = zeros(T, 2 * k);

    itm  = opts.itm;
    tol  = opts.tol;
    ptf  = opts.ptf;

    obj_his = zeros(T, itm);
    err_his = zeros(T, itm);

    # first full step
    # initialize the gradient (might take a while)
    for id = 1:n
        fars[id] = abFunc(params, id);
        abGrad(gars[id], params, id);
    end
    tfar = sum(fars);
    sum!(tgar, Gars);
    # update alpha
    copyto!(dar, tgar); dar .*= eta / n;
    copyto!(arold, ar);
    ar .-= dar;

    prox(ar);
    BLAS.axpy!(-1.0, ar, arold)
    err  = sqrt(sum(abs2, arold));
    noi  = 0;

    ind2 = collect(1:tau);

    while err >= tol
        # random sample columns
        sample!(ind, ind2, replace=false); fill!(tgar, T(0.0)); dfar = T(0.0);
        for i = 1:tau
            id = ind2[i];
            # calculate the objecitve
            rfar = abFunc(params, id);
            dfar = dfar + rfar - fars[id];
            fars[id] = rfar;
            # calculate the gradient
            abGrad(gars[id], params, id)
            tgar .+= gars[id];
        end
        # update alpha
        copyto!(dar, tgar); dar .*= eta / (tau * (floor(noi / dms) + 1));
        copyto!(arold, ar); ar .-= dar; prox(ar);
        # update tfar
        if opts.true_obj
            tfar = BFunc(params);
        else
            tfar  += dfar;
        end

        # update information
        BLAS.axpy!(-1.0, ar, arold)
        err = sqrt(sum(abs2, arold));
        noi = noi + 1;
        obj_his[noi] = tfar;
        err_his[noi] = err;

        # print information
        noi % ptf == 0 && @printf("iter %5d, obj %1.2e, err %1.2e\n",
                                  noi, tfar, err);
        noi >= itm && break;
    end
    return obj_his[1:noi], err_his[1:noi]
end


text:
#! /usr/bin/julia

# Rosetta Code, URL parser

using URIParser
const FIELDS = names(URI)

function detailview(uri::URI, indentlen::Int=4)
    indent = " "^indentlen
    s = String[]
    for f in FIELDS
        d = string(getfield(uri, f))
        !isempty(d) || continue
        f != :port || d != "0" || continue
        push!(s, @sprintf("%s%s:  %s", indent, string(f), d))
    end
    join(s, "\n")
end

test = ["foo://example.com:8042/over/there?name=ferret#nose",
        "urn:example:animal:ferret:nose",
        "jdbc:mysql://test_user:ouupppssss@localhost:3306/sakila?profileSQL=true",
        "ftp://ftp.is.co.za/rfc/rfc1808.txt",
        "http://www.ietf.org/rfc/rfc2396.txt#header1",
        "ldap://[2001:db8::7]/c=GB?objectClass=one&objectClass=two",
        "mailto:John.Doe@example.com",
        "news:comp.infosystems.www.servers.unix",
        "tel:+1-816-555-1212",
        "telnet://192.0.2.16:80/",
        "urn:oasis:names:specification:docbook:dtd:xml:4.1.2",
        "This is not a URI!",
        "ssh://alice@example.com",
        "https://bob:pass@example.com/place",
        "http://example.com/?a=1&b=2+2&c=3&c=4&d=%65%6e%63%6F%64%65%64"]

isfirst = true
for st in test
    if isfirst
        isfirst = false
    else
        println()
    end
    println("Attempting to parse\n  \"", st, "\" as a URI:")
    uri = try
        URI(st)
    catch
        println("URIParser failed to parse this URI, is it OK?")
        continue
    end
    print("This URI is parsable ")
    if isvalid(uri)
        println("and appears to be valid.")
    else
        println("but may be invalid.")
    end
    println(detailview(uri))
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
module human
    include("./fixed_parameters.jl")
    using .fixed_parameters

    mutable struct Human
        radius::Float64
        position::Tuple
        intern::Bool
        random_infection_rate::Float64
        infected::Bool
        infected_time::Int64
        show_infection_sign::Bool
        in_quarantine::Bool
        quarantine_time::Int64
        immunized::Bool
        alive::Bool
        destination::Any
        step::Float64
        break_rules::Float64

        contact::Any
        get_override_permission::Any
        move::Any

        function Human(spawnpoint::Tuple{Float64, Float64}, break_rules::Float64, intern::Bool, random_infection_rate::Float64)
            new(radius_without_masks, spawnpoint, intern, random_infection_rate, false, 0, false, false, 0, false, true, nothing, move_speed, break_rules, contact, get_override_permission, move)
        end
    end

    function contact(self::Human, human::Any)
        distance = sqrt((self.position[1]-human.position[1])^2+(self.position[2]-human.position[2])^2)
        if distance>=radius_without_masks
            return false
        end
        infection_distance = get_infection_distance(self.radius, human.radius)
        if distance < infection_distance
            infection_rate = get_infection_rate(distance, infection_distance)
            infection_active_rate = get_infection_active_rate(self.infected_time/day_tick)
            if rand()<=infection_rate*infection_active_rate
                human.infected = true
                return true
            end
        end
        return false
    end

    function get_override_permission(self::Human, fakewalls::Array)
        for wall in lundwalls
            if !(wall in fakewalls) && wall[1][1]>min(self.position[1], self.destination[1]) && wall[1][1]<max(self.position[1], self.destination[1])
                height = ((self.destination[2]-self.position[2])/(self.destination[1]-self.position[1]))*(wall[1][1]-self.position[1])+self.position[2]
                if height>=wall[2][1] && height<=wall[2][2]
                    return true
                end
            end
        end
        return false
    end

    function move(self::Human, fakewalls::Array)
        if !self.alive || self.in_quarantine
            return nothing
        end
        if self.destination==nothing || rand()<=move_random_rate
            angle = rand(1:360)
            x = cos(angle)*self.step
            y = sin(angle)*self.step
        else
            destination_override = self.get_override_permission(self, fakewalls)
            if destination_override
                temp = self.destination
                if self.position[2]<=cwh || self.position[2]>=lund_size[2]-hwh
                    self.destination = (self.position[1], rand(cwh:0.01:lund_size[2]-hwh))
                else
                    self.destination = (self.destination[1], rand(cwh:0.01:lund_size[2]-hwh))
                end
            end
            length_x = self.destination[1]-self.position[1]
            length_y = self.destination[2]-self.position[2]
            x = (self.step*length_x)/sqrt(length_x^2+length_y^2)
            y = (self.step*length_y)/sqrt(length_x^2+length_y^2)
            if destination_override
                self.destination = temp
            end
        end
        new_position = [self.position[1]+x, self.position[2]+y]

        for wall in lundwalls
            if !(wall in fakewalls) && wall[1][1]<max(new_position[1], self.position[1]) && wall[1][1]>min(new_position[1], self.position[1])
                height = (y*(wall[1][1]-self.position[1]))/x+self.position[2]
                if height>=wall[2][1] && height<=wall[2][2]
                    new_position[1] = wall[1][1]-(x/abs(x))*d
                end
                break
            end
        end
        for i in 1:2
            if new_position[i] < d
                new_position[i] = d
            elseif new_position[i] > lund_size[i]-d
                new_position[i] = lund_size[i]-d
            end
        end
        for wall in lundwalls
            if !(wall in fakewalls) && abs(new_position[1]-wall[1][1]) < d
                new_position[1] = wall[1][1]+((new_position[1]-wall[1][1])*d)/abs(new_position[1]-wall[1][1])
            end
        end
        self.position = Tuple(new_position)
    end

    export Human

end


text:
# What if we have multiple solution and we know the solution we want is sparse?
# We are going to ues matrix.

using LassoPath


text:
"""
XCB window manager wrapper
This should go into an own package, best with a highlevel interface abstracting away os
specifics.
This will be pretty much GLFW in Julia!
"""
module XCB

@enum(xcb_cw_t,
 	CW_BACK_PIXMAP = 1,
 	CW_BACK_PIXEL = 2,
 	CW_BORDER_PIXMAP = 4,
 	CW_BORDER_PIXEL = 8,
 	CW_BIT_GRAVITY = 16,
 	CW_WIN_GRAVITY = 32,
 	CW_BACKING_STORE = 64,
 	CW_BACKING_PLANES = 128,
 	CW_BACKING_PIXEL = 256,
 	CW_OVERRIDE_REDIRECT = 512,
 	CW_SAVE_UNDER = 1024,
	CW_EVENT_MASK = 2048,
	CW_DONT_PROPAGATE = 4096,
	CW_COLORMAP = 8192,
	CW_CURSOR = 16384
)

@enum(xcb_event_mask_t,
	EVENT_MASK_NO_EVENT = 0,
	EVENT_MASK_KEY_PRESS = 1,
	EVENT_MASK_KEY_RELEASE = 2,
	EVENT_MASK_BUTTON_PRESS = 4,
	EVENT_MASK_BUTTON_RELEASE = 8,
	EVENT_MASK_ENTER_WINDOW = 16,
	EVENT_MASK_LEAVE_WINDOW = 32,
	EVENT_MASK_POINTER_MOTION = 64,
	EVENT_MASK_POINTER_MOTION_HINT = 128,
	EVENT_MASK_BUTTON_1_MOTION = 256,
	EVENT_MASK_BUTTON_2_MOTION = 512,
	EVENT_MASK_BUTTON_3_MOTION = 1024,
	EVENT_MASK_BUTTON_4_MOTION = 2048,
	EVENT_MASK_BUTTON_5_MOTION = 4096,
	EVENT_MASK_BUTTON_MOTION = 8192,
	EVENT_MASK_KEYMAP_STATE = 16384,
	EVENT_MASK_EXPOSURE = 32768,
	EVENT_MASK_VISIBILITY_CHANGE = 65536,
	EVENT_MASK_STRUCTURE_NOTIFY = 131072,
	EVENT_MASK_RESIZE_REDIRECT = 262144,
	EVENT_MASK_SUBSTRUCTURE_NOTIFY = 524288,
	EVENT_MASK_SUBSTRUCTURE_REDIRECT = 1048576,
	EVENT_MASK_FOCUS_CHANGE = 2097152,
	EVENT_MASK_PROPERTY_CHANGE = 4194304,
	EVENT_MASK_COLOR_MAP_CHANGE = 8388608,
	EVENT_MASK_OWNER_GRAB_BUTTON = 16777216
)

@enum(xcb_window_class_t,
	WINDOW_CLASS_COPY_FROM_PARENT = 0,
	WINDOW_CLASS_INPUT_OUTPUT = 1,
	WINDOW_CLASS_INPUT_ONLY = 2
)

const xcb_connection_t = Ptr{Void} # TODO create correct composite type
const xcb_window_t = UInt32
const xcb_colormap_t = UInt32
const xcb_visualid_t = UInt32
const xcb_keycode_t = UInt8

struct xcb_setup_t
	status::UInt8
	pad0::UInt8
	uprotocol_major_version::UInt16
	uprotocol_minor_version::UInt16
	length::UInt16
	release_number::UInt32
	resource_id_base::UInt32
	resource_id_mask::UInt32
	motion_buffer_size::UInt32
	vendor_len::UInt16
	umaximum_request_length::UInt16
	roots_len::UInt8
	pixmap_formats_len::UInt8
	image_byte_order::UInt8
	bitmap_format_bit_order::UInt8
	bitmap_format_scanline_unit::UInt8
	bitmap_format_scanline_pad::UInt8
	min_keycode::xcb_keycode_t
	max_keycode::xcb_keycode_t
	pad1::NTuple{4, UInt8}
end

struct xcb_screen_t
	root::xcb_window_t
	default_colormap::xcb_colormap_t
	white_pixel::UInt32
	black_pixel::UInt32
	current_input_masks::UInt32
	width_in_pixels::UInt16
	height_in_pixels::UInt16
	width_in_millimeters::UInt16
	height_in_millimeters::UInt16
	min_installed_maps::UInt16
	max_installed_maps::UInt16
	root_visual::xcb_visualid_t
	backing_stores::UInt8
	save_unders::UInt8
	root_depth::UInt8
	allowed_depths_len::UInt8
end
const COPY_FROM_PARENT = Clong(0)
struct xcb_screen_iterator_t
	data::Ptr{xcb_screen_t}
	rem::Cint
	index::Cint
end
struct xcb_void_cookie_t
	sequence::Cuint
end
function xcb_connect(displayname, screenp)
	ccall((:xcb_connect, "libxcb"), xcb_connection_t, (Ptr{Cchar}, Ptr{Cint}), displayname, screenp)
end
function connect(displayname)
	err = Ref{Cint}(0)
	if isa(displayname, AbstractString)
		displayname = ascii(displayname)
	end
	connection = xcb_connect(displayname, err)
	if connection == C_NULL
		error("Cannot find a compatible Vulkan installable client driver (ICD).\nExiting ...\n")
	end
	connection, err[]
end

function xcb_get_setup(connection)
	ccall((:xcb_get_setup, "libxcb"), Ptr{xcb_setup_t}, (xcb_connection_t,), connection)
end
function get_setup(connection)
	setup = xcb_get_setup(connection)
	if setup == C_NULL
		error("Couldn't get setup")
	end
	setup
end

function xcb_setup_roots_iterator(setup)
	ccall((:xcb_setup_roots_iterator, "libxcb"), xcb_screen_iterator_t, (Ptr{xcb_setup_t},), setup)
end
setup_roots_iterator(setup) = xcb_setup_roots_iterator(setup)
function xcb_screen_next(iter)
	ccall((:xcb_screen_next, "libxcb"), Void, (Ptr{xcb_screen_iterator_t},), iter)
end
screen_next(iter) = xcb_screen_next(iter)

function xcb_generate_id(connection)
	ccall((:xcb_generate_id, "libxcb"), Cuint, (xcb_connection_t,), connection)
end
generate_id(connection) = xcb_generate_id(connection)

function create_window(
		connection,
		depth,
		wid,
		parent,
		x,
		y,
		width,
		height,
		border_width,
		_class,
		visual,
		value_mask,
		value_list
	)
	ccall(
		(:xcb_create_window, "libxcb"), xcb_void_cookie_t,
		(
			xcb_connection_t,
			UInt8           ,
			xcb_window_t    ,
			xcb_window_t    ,
			Int16           ,
			Int16           ,
			UInt16          ,
			UInt16          ,
			UInt16          ,
			UInt16          ,
			xcb_visualid_t  ,
			UInt32          ,
			Ptr{UInt32}
		),
		connection,
		depth,
		wid,
		parent,
		x,
		y,
		width,
		height,
		border_width,
		_class,
		visual,
		value_mask,
		value_list
	)
end

end


text:
using JuAFEM
using JuAFEM:Vec
using Tensors
using SparseArrays
using LinearAlgebra
include("../src/FEMRefinements.jl")
using .FEMRefinements

const ∇ = Tensors.gradient
const Δ = Tensors.hessian;

v1 = Vec((2.0, 0.0))
v2 = Vec((0.0, 2.0))
grid = generate_grid2D(40, 40, v1, v2)

dim = 2
ip = Lagrange{dim,RefCube,1}()
qr = QuadratureRule{dim,RefCube}(2)
qr_face = QuadratureRule{dim - 1,RefCube}(2)
cellvalues = CellScalarValues(qr, ip);
facevalues = FaceScalarValues(qr_face, ip);

dh = DofHandler(grid)
push!(dh, :u, 1)
close!(dh)

function u_ana(x::Vec{2,T}) where {T}
    xs = (Vec{2}((0.5,  1.5)),
          Vec{2}((1.5, 1.5)),
          Vec{2}((0.5,  0.5)))
    σ = 1 / 8
    s = zero(eltype(x))
    for i in 1:3
        s += exp(- norm(x - xs[i])^2 / σ^2)
    end
    return max(1e-15 * one(T), s) # Denormals, be gone
end;

dbcs = ConstraintHandler(dh)
dbc = Dirichlet(:u, union(getfaceset(grid, "top"), getfaceset(grid, "right")), (x, t) -> u_ana(x))
add!(dbcs, dbc)
close!(dbcs)
update!(dbcs, 0.0)

K = create_sparsity_pattern(dh);

function doassemble(cellvalues::CellScalarValues{dim}, facevalues::FaceScalarValues{dim},
                         K::SparseMatrixCSC, dh::DofHandler) where {dim}
    b = 1.0
    f = zeros(ndofs(dh))
    assembler = start_assemble(K, f)

    n_basefuncs = getnbasefunctions(cellvalues)
    global_dofs = zeros(Int, ndofs_per_cell(dh))

    fe = zeros(n_basefuncs) # Local force vector
    Ke = zeros(n_basefuncs, n_basefuncs) # Local stiffness mastrix

    @inbounds for (cellcount, cell) in enumerate(CellIterator(dh))
        fill!(Ke, 0)
        fill!(fe, 0)
        coords = getcoordinates(cell)

        reinit!(cellvalues, cell)
        for q_point in 1:getnquadpoints(cellvalues)
            dΩ = getdetJdV(cellvalues, q_point)
            coords_qp = spatial_coordinate(cellvalues, q_point, coords)
            f_true = -LinearAlgebra.tr(hessian(u_ana, coords_qp)) + u_ana(coords_qp)
            for i in 1:n_basefuncs
                δu = shape_value(cellvalues, q_point, i)
                ∇δu = shape_gradient(cellvalues, q_point, i)
                fe[i] += (δu * f_true) * dΩ
                for j in 1:n_basefuncs
                    u = shape_value(cellvalues, q_point, j)
                    ∇u = shape_gradient(cellvalues, q_point, j)
                    Ke[i, j] += (∇δu ⋅ ∇u + δu * u) * dΩ
                end
            end
        end
        for face in 1:nfaces(cell)
            if (cellcount, face) ∈ getfaceset(grid, "left") ||
                    (cellcount, face) ∈ getfaceset(grid, "bottom")
                reinit!(facevalues, cell, face)
                for q_point in 1:getnquadpoints(facevalues)
                    coords_qp = spatial_coordinate(facevalues, q_point, coords)
                    n = getnormal(facevalues, q_point)
                    g = gradient(u_ana, coords_qp) ⋅ n
                    dΓ = getdetJdV(facevalues, q_point)
                    for i in 1:n_basefuncs
                        δu = shape_value(facevalues, q_point, i)
                        fe[i] += -(δu * g) * dΓ
                        for j in 1:n_basefuncs
                            ∇u = shape_gradient(cellvalues, q_point, j)
                            Ke[i, j] += (δu * ∇u ⋅ n) * dΓ
                        end
                    end
                end
            end
        end

        celldofs!(global_dofs, cell)
        assemble!(assembler, global_dofs, fe, Ke)
    end
    return K, f
end;

K, f = doassemble(cellvalues, facevalues, K, dh);
apply!(K, f, dbcs)
u = Symmetric(K) \ f;

AdaptiveMesh.plot_solution3D(dh, u; plotgrid=true)


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests of filters.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Filters" begin
    expected = """

      Row
      Col. 1
      Col. 3

      2
      2
      2.0

      4
      4
      4.0

      6
      6
      6.0

"""

    result = pretty_table(
        String, data;
        backend = Val(:html),
        filters_row = ((data,i) -> i%2 == 0,),
        filters_col = ((data,i) -> i%2 == 1,),
        formatters = ft_printf("%.1f",3),
        standalone = false,
        show_row_number = true
    )

    @test result == expected

    expected = """

      Row
      Col. 1
      Col. 3

      2
      2
      2.0

      4
      4
      4.0

      6
      6
      6.0

"""

    result = pretty_table(
        String, data;
        alignment = [:c, :l, :l, :c],
        backend = Val(:html),
        filters_row = ((data, i) -> i % 2 == 0,),
        filters_col = ((data ,i) -> i % 2 == 1,),
        formatters = ft_printf("%.1f", 3),
        standalone = false,
        show_row_number = true
    )

    @test result == expected
end


text:
using ADMPS
using ADMPS:leftenv,rightenv,norm_FL,norm_FR
using CUDA
using KrylovKit
using LinearAlgebra
using Random
using Test
using OMEinsum
CUDA.allowscalar(false)

@testset "eigsolve with $atype{$dtype}" for atype in [Array, CuArray], dtype in [Float64]
    Random.seed!(100)
    D,d = 10,2
    A = atype(rand(D,d,D))
    工 = ein"asc,bsd -> abcd"(A,conj(A))
    λLs, Ls, info = eigsolve(L -> ein"ab,abcd -> cd"(L,工), atype(rand(D,D)), 1, :LM)
    λL, L = λLs[1], Ls[1]
    @test imag(λL) ≈ 0
    @test Array(ein"ab,ab -> "(L,L))[] ≈ 1
    @test λL * L ≈ ein"ab,abcd -> cd"(L,工)

    λRs, Rs, info = eigsolve(R -> ein"abcd,cd -> ab"(工,R), atype(rand(D,D)), 1, :LM)
    λR, R = λRs[1], Rs[1]
    @test imag(λR) ≈ 0
    @test Array(ein"ab,ab -> "(R,R))[] ≈ 1
    @test λR * R ≈ ein"abcd,cd -> ab"(工,R)
    @test λL ≈ λR
end

@testset "leftenv and rightenv with $atype{$dtype}" for atype in [Array, CuArray], dtype in [Float64, ComplexF64]
    Random.seed!(100)
    d = 2
    D = 10

    β = rand()
    Au = atype(rand(dtype,D,d,D))
    Ad = atype(rand(dtype,D,d,D))
    M = atype(model_tensor(Ising(β)))

    λL,FL = leftenv(Au, Ad, M)
    @test λL * FL ≈ ein"((adf,abc),dgeb),fgh -> ceh"(FL,Au,M,Ad)

    λR,FR = rightenv(Au, Ad, M)
    @test λR * FR ≈ ein"((abc,ceh),dgeb),fgh -> adf"(Au,FR,M,Ad)
end

@testset "normalization leftenv and rightenv with $atype{$dtype}" for atype in [Array, CuArray], dtype in [Float64, ComplexF64]
    Random.seed!(100)
    d = 2
    D = 10

    Au = atype(rand(dtype,D,d,D))
    Ad = atype(rand(dtype,D,d,D))

    λL,FL = norm_FL(Au, Ad)
    @test λL * FL ≈ ein"(ad,acb), dce -> be"(FL,Au,Ad)
    λR,FR = norm_FR(Au, Ad)
    @test λR * FR ≈ ein"(be,acb), dce -> ad"(FR,Au,Ad)
end


text:
nbdir = ENV["HOME"] * "/.julia/dev/DeltaRCWA/notebooks/"
include("$(nbdir)NystromMethodQP.jl")
import FFTW: fftfreq
import LinearAlgebra: I, Diagonal

"""
    compute_BIE_method(η::Function, μ::Function, L::Float64, N::Int, k::Float64, i::Int)

Compute the fields using Carlos' Boundary Integral Equation (BIE) method.
The output fields are on a grid compatible with the RCWA method.

Arguments:
η:: A periodic function returning the metasurface magnetic conductance
μ:: A periodic function returning the metasurface electric resistance
L:: the length of the unit cell (some multiple of the period of η, μ)
N:: the number of discretization points along the unit cell
k:: the wavenumber/frequency of the incident planewave (assumed TM polarization)
i:: the index of the incident planewave in the Fourier basis

Keyword arguments:
method:: the Nystrom method to be used (choose default "MK" or "Alpert)
NC:: the number of cells in a windowed sum for Green's function (default 5)
NP:: the number of points in BIE discretization per wavelength (default 20)
s:: a scale factor to stretch the y coordinate of the output grid

Returns:
NamedTuple with fields:
I₁:: an array of the incident fields in the position basis for y<0
O₁:: an array of the outgoing fields in the position basis for y<0
I₂:: an array of the incident fields in the position basis for y>0
O₂:: an array of the outgoing fields in the position basis for y>0
"""
function compute_BIE_method(
    η::Function, μ::Function, L::Float64, N::Int, k::Float64, i::Int;
    method="MK", NC=5, NP=20, s=1
)
    dx = L/N
    k⃗ = 2π*fftfreq(N, 1/dx)
    β⃗ = @. sqrt(Complex(k^2 - k⃗^2))
    α = k⃗[i]
    β = β⃗[i]
    uInc(x,y)= exp(im*α*x+im*β*y)  # incident planewave
    # Defining geometry
    # Vertices that define the curve Γ
    ver = [
        L/2   0;   #1
        -L/2   0   #2
    ]
    # In this case, Γ is just a line segment from (L/2,0) to (-L/2,0).
    # This matrix defines the orientation of the curve
            #starting node  #ending node  #domain on the right #domain on the left
    info = [1               2             1                    0]
    # The parameter M defines the number of boundary points to used in the discretization of the BIE
    λ = 2π/k
    M  = Int(round(max(NP,NP/λ))) # ~20 pts per unit length
    # Γ contains all the information of the discretized curve that is needed for the numerical solution of the problem
    Γ  = SkeletonCctor(ver,info,[M])

    # Solving the boundary integral equations
    # method::Selection of the Nystrom method to be used
    # NC::Number of unit cells considered in the windowed sum (A = NC*L)
    mat = matricesSkeleton(1,Γ,k,L,α,NC,method) # Construction of the Nystrom matrices

    x1 = Γ.parts[1].x[:,1] # x-coordinate of the discretization points on Γ

    f1(x) = -2*im*k*η(x)*exp(im*α*x)
    A⁺ = -0.5*I+((im*k)*Diagonal(η.(x1)))*mat.S  # System matrix for Ω⁺
    b⁺ = f1.(x1) # boundary data for Ω⁺
    ψ⁺ = A⁺\b⁺

    f2(x) = -2*im*β*exp(im*α*x)
    A⁻ = -0.5*I+((im*k)*Diagonal(μ.(x1)))*mat.S # System matrix for Ω⁻
    b⁻ =  f2.(x1) # boundary data for Ω⁻
    ψ⁻ = A⁻\b⁻

    φ⁺ = 0.5*(ψ⁺+ψ⁻)
    φ⁻ = 0.5*(ψ⁺-ψ⁻)

    x = range(0, step=dx, length=N)
    y = range(-s*L, s*L, length=2N)
    y₁ = y[y .< 0] # y points in domain below sheet
    y₂ = y[y .> 0] # y points in domain above sheet

    ver⁺ = [L 0;-L 0;-L 2s*L;L 2s*L]
    Ω⁺ = isin(x,y,ver⁺)

    ver⁻ = [L 0;-L 0;-L -2s*L;L -2s*L]
    Ω⁻ = isin(x,y,ver⁻)

    UInc = reshape(uInc.(Ω⁺.pts[:,1],Ω⁺.pts[:,2]),Ω⁺.Ny,Ω⁺.Nx)

    pot⁺ = potentialsSkeleton(Ω⁺.pts,Γ,k,L,α,NC);
    U⁺scat = reshape(pot⁺.SL*φ⁺,Ω⁺.Ny,Ω⁺.Nx)
    U⁺   =  Ω⁺.In .* (U⁺scat .+ UInc)

    pot⁻= potentialsSkeleton(Ω⁻.pts,Γ,k,L,α,NC);
    U⁻scat = reshape(pot⁻.SL*φ⁻,Ω⁻.Ny,Ω⁻.Nx)
    U⁻   =  Ω⁻.In .* (U⁻scat .+ UInc)

    Ω⁺mask = isfinite.(Ω⁺.In)
    I₂ = reshape(UInc[Ω⁺mask] .- UInc[Ω⁺mask], length.((y₂, x)))
    O₂ = reshape(U⁺[Ω⁺mask], length.((y₂, x)))

    Ω⁻mask = isfinite.(Ω⁻.In)
    I₁ = reshape(UInc[Ω⁻mask], length.((y₁, x)))
    O₁ = reshape(U⁻[Ω⁻mask] .- UInc[Ω⁻mask], length.((y₁, x)))
    (x=x, y₁=y₁, y₂=y₂, I₁=I₁, O₁=O₁, I₂=I₂, O₂=O₂)
end


text:
module TraitWrappers

export AbstractTraitWrapper, object

"""
	TraitWrapper{T}

Abtract type wrapper
"""
abstract type AbstractTraitWrapper end

"""
	object(t::TraitWrapper)

Return the object wrapped in the trait
"""
object(t::AbstractTraitWrapper) = t.object

include("IterableTraitWrapper.jl")

end # module


text:
function cable_transform(y, z)
    v1 = [0.0, 0.0, 1.0]
    # if norm(z) > norm(y)
    #     v2 = y[1:3,1] - z[1:3,1]
    # else
    #     v2 = z[1:3,1] - y[1:3,1]
    # end
    v2 = y[1:3,1] - z[1:3,1]
    normalize!(v2)
    ax = cross(v1, v2)
    ang = acos(v1'*v2)
    R = AngleAxis(ang, ax...)

    if any(isnan.(R))
        R = I
    else
        nothing
    end

    compose(Translation(z), LinearMap(R))
end

function ModifiedMeshFileObject(obj_path::String, material_path::String;
        scale::T = 0.1) where {T}
    obj = MeshFileObject(obj_path)
    rescaled_contents = rescale_contents(obj_path, scale = scale)
    material = select_material(material_path)
    mod_obj = MeshFileObject(
        rescaled_contents,
        obj.format,
        material,
        obj.resources,
        )
    return mod_obj
end

function rescale_contents(obj_path::String; scale::T = 0.1) where T
    lines = readlines(obj_path)
    rescaled_lines = copy(lines)
    for (k,line) in enumerate(lines)
        if length(line) >= 2
            if line[1] == 'v'
                stringvec = split(line, " ")
                vals = map(x -> parse(Float64, x), stringvec[2:end])
                rescaled_vals = vals .* scale
                rescaled_lines[k] = join([stringvec[1]; string.(rescaled_vals)], " ")
            end
        end
    end
    rescaled_contents = join(rescaled_lines, "\r\n")
    return rescaled_contents
end

function select_material(material_path::String)
    mtl_file = open(material_path)
    mtl = read(mtl_file, String)
    return mtl
end

function default_background!(vis)
    setvisible!(vis["/Background"], true)
    setprop!(vis["/Background"], "top_color", RGBA(1.0, 1.0, 1.0, 1.0))
    setprop!(vis["/Background"], "bottom_color", RGBA(1.0, 1.0, 1.0, 1.0))
    setvisible!(vis["/Axes"], false)
end

function pad_trajectory(x, shift, T_shift)
	[[x[1] + shift for i = 1:T_shift]..., [_x + shift for _x in x]..., [x[end] + shift for i = 1:T_shift]...]
end


text:
# Classic Lax-Friedrichs Method
# Reference:
# R. Leveque. Finite Volume Methods for Hyperbolic Problems.Cambridge University
# Press. New York 2002

struct LaxFriedrichsScheme <: AbstractFVAlgorithm end

struct LocalLaxFriedrichsScheme <: AbstractFVAlgorithm end

mutable struct GlobalLaxFriedrichsScheme{T} <: AbstractFVAlgorithm
  αf :: Function #viscosity coefficient
  α :: T
end

function update_dt(alg::GlobalLaxFriedrichsScheme{T1},u,Flux,
    CFL,mesh) where {T1}
  alg.α = alg.αf(u, Flux, mesh)
  @assert (abs(alg.α) > eps(T1))
  dx = cell_volume(mesh, 1)
  CFL*dx/alg.α
end

function GlobalLaxFriedrichsScheme(;αf = nothing)
    if αf == nothing
        αf = maxfluxρ
    end
    GlobalLaxFriedrichsScheme(αf,0.0)
end

function update_flux_value(uold,node_idx,dt,dx,mesh,Flux,alg::LaxFriedrichsScheme)
    # Local speeds of propagation
    @inbounds ul=cellval_at_left(node_idx,uold,mesh)
    @inbounds ur=cellval_at_right(node_idx,uold,mesh)
    # Numerical Fluxes
    return 0.5*(Flux(ul)+Flux(ur))-dx/(2*dt)*(ur-ul)
end

"""
compute_fluxes!(hh, Flux, u, mesh, dt, M, alg::LocalLaxFriedrichsScheme, ::Type{Val{true}})
Numerical flux of local lax friedrichs algorithm in 1D
"""
function compute_fluxes!(fluxes, Flux, u, mesh, dt, alg::LocalLaxFriedrichsScheme, nonscalar::Bool,::Type{Val{false}})
    dx = cell_volume(mesh, 1)
    ul=cellval_at_left(1,u,mesh)
    αl = fluxρ(ul, Flux)
    #update vector
    for j in node_indices(mesh)
        # Local speeds of propagation
        @inbounds ul=cellval_at_left(j,u,mesh)
        @inbounds ur=cellval_at_right(j,u,mesh)
        αr = fluxρ(ur, Flux)
        αk = max(αl, αr)
        # Numerical Fluxes
        if nonscalar
            fluxes[:,j] .= 0.5*(Flux(ul)+Flux(ur))-αk*(ur-ul)
        else
            fluxes[j] = 0.5*(Flux(ul)+Flux(ur))-αk*(ur-ul)
        end
        αl = αr
    end
end

function update_flux_value(uold,node_idx,dt,dx,mesh,Flux, alg::GlobalLaxFriedrichsScheme)
    # Local speeds of propagation
    @inbounds ul=cellval_at_left(node_idx,uold,mesh)
    @inbounds ur=cellval_at_right(node_idx,uold,mesh)
    # update numerical flux
    return 0.5*(Flux(ul)+Flux(ur))-alg.α*(ur-ul)
end


text:
__precompile__(false)

module POMDPModelChecking

# package code goes here
using POMDPs
using POMDPModelTools
using POMDPPolicies
using LightGraphs
using Parameters
using Random
using DiscreteValueIteration
using LinearAlgebra
using Distributions
using Spot

export
    ReachabilitySolver,
    ReachabilityPolicy,
    ReachabilityMDP,
    ReachabilityPOMDP

include("reachability.jl")

export
    ProductState,
    ProductMDP,
    ProductPOMDP,
    labels

include("product.jl")

export
    # graph analysis
    maximal_end_components,
    mdp_to_graph,
    sub_mdp,
    accepting_states!

include("end_component.jl")

export
    ModelCheckingSolver,
    ModelCheckingPolicy

include("model_checking_solver.jl")

end # POMDPModelChecking module


text:
module GaussianRV
using FFTW

include("fft_sample.jl")
export sample_field

end # module


text:
module LoadingSpinner

using ..Ahorn
using Gtk, Gtk.ShortNames
using Cairo

const baseWidth = ceil(Int, width(Ahorn.Assets.crowIdle))
const baseHeight = ceil(Int, height(Ahorn.Assets.crowIdle))

const startFrame = Ahorn.scalePixbufSimple(Ahorn.pixbufFromSurface(Ahorn.Assets.crowIdle), baseWidth * 4, baseHeight * 4)
const idleFrame = Ahorn.scalePixbufSimple(Ahorn.pixbufFromSurface(Ahorn.Assets.crowIdle), baseWidth * 4, baseHeight * 4)
const peck1Frame = Ahorn.scalePixbufSimple(Ahorn.pixbufFromSurface(Ahorn.Assets.crowPeck1), baseWidth * 4, baseHeight * 4)
const peck2Frame = Ahorn.scalePixbufSimple(Ahorn.pixbufFromSurface(Ahorn.Assets.crowPeck2), baseWidth * 4, baseHeight * 4)

# Peck twice, return to idle position
const frames = Tuple{Gtk.Pixbuf, Number}[
    (peck1Frame, 0.08),
    (peck2Frame, 0.08),
    (idleFrame, 0.08),
    (peck1Frame, 0.08),
    (peck2Frame, 0.08),
    (idleFrame, 0.08),
]

function animate(dlg::Gtk.GObject)
    for (frame, delay) in frames
        Ahorn.setProgressDialogPixbuf!(dlg, frame)
        sleep(delay)
    end
end

end


text:
#---------------------------------------------------------
# # [Simple Examples](@id 1_Examples)
#---------------------------------------------------------

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/1_Examples.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/1_Examples.ipynb)

#md # !!! note
#md #     These examples are available as Jupyter notebooks.
#md #     You can execute them online with [binder](https://mybinder.org/) or just view them with [nbviewer](https://nbviewer.jupyter.org/) by clicking on the badges above!

# These examples show what UnitfulRecipes is all about.

# First we need to tell Julia we are using Plots, Unitful, and UnitfulRecipes

using Plots, Unitful, UnitfulRecipes

# ## Simplest plot

# This is the most basic example

y = randn(10)*u"kg"
plot(y)

# Add some more plots, and it will be aware of the units you used previously (note `y2` is about 10 times smaller than `y1`)

y2 = 100randn(10)*u"g"
plot!(y2)

# UnitfulRecipes will not allow you to plot with different unit-dimensions, so
# ```julia
# plot!(rand(10)*u"m")
# ```
# won't work here.
#
# But you can add inset subplots with different axes that have different dimensions

plot!(rand(10)*u"m", inset=bbox(0.5, 0.5, 0.3, 0.3), subplot=2)

# ## Axis label

# If you specify an axis label, the unit will be appended to it.

plot(y, ylabel="mass")

# Unless you want it untouched, in which case you can use a "protected" string using the `@P_str` macro.

plot(y, ylabel=P"mass in kilograms")

# Just like with the `label` keyword for legends, no axis label is added if you specify the axis label to be an empty string.

plot(y, ylabel="")

# ## Axis unit

# You can use the axis-specific keyword arguments to convert units on the fly

plot(y, yunit=u"g")

# ## Axis limits

# Setting the axis limits can be done with units

plot(y, ylims=(-1000u"g",2000u"g"))

# or without

plot(y, ylims=(-1,2))

# ## Multiple series

# You can plot multiple series as 2D arrays

x, y = rand(10,3)*u"m", rand(10,3)*u"g"
plot(x, y)

# Or vectors of vectors (of potnetially different lengths)

x, y = [rand(10), rand(15), rand(20)]*u"m", [rand(10), rand(15), rand(20)]*u"g"
plot(x, y)

# ## 3D

# It works in 3D

x, y = rand(10)*u"km", rand(10)*u"hr"
z = x ./ y
plot(x, y, z)

# ## Scatter plots

# You can do scatter plots

scatter(x, y, zcolor=z, clims=(5,20).*unit(eltype(z)))

# and 3D scatter plots too

scatter(x, y, z, zcolor=z)

# ## Contour plots

# for contours plots

x, y = (1:0.01:2)*u"m", (1:0.02:2)*u"s"
z = x' ./ y
contour(x, y, z)

# and filled contours

contourf(x, y, z)


text:
type CompositeMap{T}<:AbstractLinearMap{T}
#    t::Type{T} # just stores T, trick to make CompositeMap{T} actually depend on T
    maps::Vector{AbstractLinearMap} # stored in order of application to vector
    function CompositeMap(maps::Vector{AbstractLinearMap})
        N=length(maps)
        for n=2:N
            size(maps[n],2)==size(maps[n-1],1) || throw(DimensionMismatch("CompositeMap"))
        end
        for n=1:N
            promote_type(T,eltype(maps[n]))==T || throw(InexactError())
        end
        new(maps)
    end
end

# basic methods
Base.size(A::CompositeMap,n)=(n==1 ? size(A.maps[end],1) : (n==2 ? size(A.maps[1],2) : error("AbstractLinearMap objects have only 2 dimensions")))
Base.size(A::CompositeMap)=(size(A,1),size(A,2))
Base.isreal(A::CompositeMap)=all(isreal,A.maps) # sufficient but not necessary
# the following rules are sufficient but not necessary
function Base.issym(A::CompositeMap)
    N=length(A.maps)
    if isodd(N)
        issym(A.maps[div(N+1,2)]) || return false
    end
    for n=1:div(N,2)
        A.maps[n]==transpose(A.maps[N-n+1]) || return false
    end
    return true
end
function Base.ishermitian(A::CompositeMap)
    N=length(A.maps)
    if isodd(N)
        ishermitian(A.maps[div(N+1,2)]) || return false
    end
    for n=1:div(N,2)
        A.maps[n]==ctranspose(A.maps[N-n+1]) || return false
    end
    return true
end
function Base.isposdef(A::CompositeMap)
    N=length(A.maps)
    if isodd(N)
        isposdef(A.maps[div(N+1,2)]) || return false
    end
    for n=1:div(N,2)
        A.maps[n]==ctranspose(A.maps[N-n+1]) || return false
    end
    return true
end

# composition of linear maps
function *(A1::CompositeMap,A2::CompositeMap)
    size(A1,2)==size(A2,1) || throw(DimensionMismatch("*"))
    T=promote_type(eltype(A1),eltype(A2))
    return CompositeMap{T}(AbstractLinearMap[A2.maps...,A1.maps...])
end
function *(A1::AbstractLinearMap,A2::CompositeMap)
    size(A1,2)==size(A2,1) || throw(DimensionMismatch("*"))
    T=promote_type(eltype(A1),eltype(A2))
    return CompositeMap{T}(AbstractLinearMap[A2.maps...,A1])
end
function *(A1::CompositeMap,A2::AbstractLinearMap)
    size(A1,2)==size(A2,1) || throw(DimensionMismatch("*"))
    T=promote_type(eltype(A1),eltype(A2))
    return CompositeMap{T}(AbstractLinearMap[A2,A1.maps...])
end
function *(A1::AbstractLinearMap,A2::AbstractLinearMap)
    size(A1,2)==size(A2,1) || throw(DimensionMismatch("*"))
    T=promote_type(eltype(A1),eltype(A2))
    return CompositeMap{T}(AbstractLinearMap[A2,A1])
end

# comparison of CompositeMap objects
==(A::CompositeMap,B::CompositeMap)=(eltype(A)==eltype(B) && A.maps==B.maps)

# special transposition behavior
Base.transpose(A::CompositeMap)=CompositeMap{eltype(A)}(AbstractLinearMap[transpose(M) for M in reverse(A.maps)])
Base.ctranspose(A::CompositeMap)=CompositeMap{eltype(A)}(AbstractLinearMap[ctranspose(M) for M in reverse(A.maps)])

# multiplication with vectors
function Base.A_mul_B!(y::AbstractVector,A::CompositeMap,x::AbstractVector)
    # no size checking, will be done by individual maps
    N=length(A.maps)
    if N==1
        Base.A_mul_B!(y,A.maps[1],x)
    else
        T=promote_type(eltype(A),eltype(x))
        dest=Array(T,size(A.maps[1],1))
        Base.A_mul_B!(dest,A.maps[1],x)
        source=dest
        if N>2
            dest=Array(T,size(A.maps[2],1))
        end
        for n=2:N-1
            resize!(dest,size(A.maps[n],1))
            Base.A_mul_B!(dest,A.maps[n],source)
            dest,source=source,dest # alternate dest and source
        end
        Base.A_mul_B!(y,A.maps[N],source)
    end
    return y
end
function Base.At_mul_B!(y::AbstractVector,A::CompositeMap,x::AbstractVector)
    # no size checking, will be done by individual maps
    N=length(A.maps)
    if N==1
        Base.At_mul_B!(y,A.maps[1],x)
    else
        T=promote_type(eltype(A),eltype(x))
        dest=Array(T,size(A.maps[N],2))
        Base.At_mul_B!(dest,A.maps[N],x)
        source=dest
        if N>2
            dest=Array(T,size(A.maps[N-1],2))
        end
        for n=N-1:-1:2
            resize!(dest,size(A.maps[n],2))
            Base.At_mul_B!(dest,A.maps[n],source)
            dest,source=source,dest # alternate dest and source
        end
        Base.At_mul_B!(y,A.maps[1],source)
    end
    return y
end
function Base.Ac_mul_B!(y::AbstractVector,A::CompositeMap,x::AbstractVector)
    # no size checking, will be done by individual maps
    N=length(A.maps)
    if N==1
        Base.Ac_mul_B!(y,A.maps[1],x)
    else
        T=promote_type(eltype(A),eltype(x))
        dest=Array(T,size(A.maps[N],2))
        Base.Ac_mul_B!(dest,A.maps[N],x)
        source=dest
        if N>2
            dest=Array(T,size(A.maps[N-1],2))
        end
        for n=N-1:-1:2
            resize!(dest,size(A.maps[n],2))
            Base.Ac_mul_B!(dest,A.maps[n],source)
            dest,source=source,dest # alternate dest and source
        end
        Base.Ac_mul_B!(y,A.maps[1],source)
    end
    return y
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Trilinos_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Trilinos")
JLLWrappers.@generate_main_file("Trilinos", UUID("b6fd3212-6f87-5999-b9ea-021e9cd21b17"))
end  # module Trilinos_jll


text:
"""
    ThermalLoad

Thermal load.
"""
struct ThermalLoad <: Resource
    index::Int
    num_timesteps::Int

    temp_min::Vector{Float64}  # Minimum temperature at each time-step
    temp_max::Vector{Float64}  # Maximum temperature at each time-step
    temp_ext::Vector{Float64}  # Outside temperature at each time-step
    temp_init::Float64         # Initial temperature
    pwr_min::Vector{Float64}  # Minimum power when device is on
    pwr_max::Vector{Float64}  # Maximum power when device is on
    C::Float64  # Heat capacity
    η::Float64  # Power efficiency of the heating unit
    μ::Float64  # Thermal conduction coefficient

    binflag::Bool  # Hard/soft on-off flag

    function ThermalLoad(;
        index=0, T=0,
        temp_min=Float64[], temp_max=Float64[], temp_ext=Float64[], temp_init=0.0,
        pwr_min=Float64[], pwr_max=Float64[],
        C=1.0, η=1.0, μ=1.0,
        binflag=true
    )

        # Sanity checks
        C  > 0.0 || throw(DomainError("C  must be positive"))
        T == length(temp_min) || throw(DimensionMismatch(
            "T=$T but temp_min has size $(length(temp_min))"
        ))
        T == length(temp_max) || throw(DimensionMismatch(
            "T=$T but temp_max has size $(length(temp_max))"
        ))
        T == length(temp_ext) || throw(DimensionMismatch(
            "T=$T but temp_ext has size $(length(temp_ext))"
        ))
        T == length(pwr_min) || throw(DimensionMismatch(
            "T=$T but pwr_min has size $(length(pwr_min))"
        ))
        T == length(pwr_max) || throw(DimensionMismatch(
            "T=$T but pwr_max has size $(length(pwr_max))"
        ))

        return new(
            index, T,
            temp_min, temp_max, temp_ext, temp_init,
            pwr_min, pwr_max,
            C, η, μ,
            binflag
        )
    end
end

function add_resource_to_model!(
    m::MOI.ModelLike,
    h::House,
    l::ThermalLoad,
    var2idx, con2idx
)
    T = l.num_timesteps

    # ==========================================
    #    I. Add local variables
    # ==========================================
    pnet = MOI.add_variables(m, T)  # net load
    temp = MOI.add_variables(m, T)  # net load
    uind = MOI.add_variables(m, T)  # net load

    # Variable bounds
    for t in 1:T
        # Net load
        MOI.add_constraint(m, MOI.SingleVariable(pnet[t]),
            MOI.Interval(l.pwr_min[t], l.pwr_max[t])
        )
        # Inside temperature
        MOI.add_constraint(m, MOI.SingleVariable(temp[t]),
            MOI.Interval(l.temp_min[t], l.temp_max[t])
        )
        # On-off indicator
        if l.binflag
            MOI.add_constraint(m, MOI.SingleVariable(uind[t]),
                MOI.ZeroOne()
            )
        else
            MOI.add_constraint(m, MOI.SingleVariable(uind[t]),
                MOI.Interval(0.0, 1.0)
            )
        end
    end
    # Update variable indices
    for t in 1:T
        var2idx[(:thermal, l.index, :pnet, t)] = pnet[t]
        var2idx[(:thermal, l.index, :temp, t)] = temp[t]
        var2idx[(:thermal, l.index, :uind, t)] = uind[t]
    end

    # ==========================================
    #    II. Update linking constraints
    # ==========================================
    for t in 1:T
        cidx = con2idx[(:house, h.index, :link, t)]
        MOI.modify(m, cidx, MOI.ScalarCoefficientChange(pnet[t], -1.0))
    end

    # ==========================================
    #    III. Add local constraints
    # ==========================================
    for t in 1:T
        # On-off
        # Pmin[t] * u[t] <= p[t], i.e., p[t] - Pmin[t] * u[t] >= 0
        MOI.add_constraint(m,
            MOI.ScalarAffineFunction{Float64}([
                MOI.ScalarAffineTerm(1.0, pnet[t]),
                MOI.ScalarAffineTerm(-l.pwr_min[t], uind[t])
            ], 0.0),
            MOI.GreaterThan(0.0)
        )
        # p[t] <= Pmax[t] * u[t], i.e., p[t] - Pmax[t] * u[t] <= 0
        MOI.add_constraint(m,
            MOI.ScalarAffineFunction{Float64}([
                MOI.ScalarAffineTerm(1.0, pnet[t]),
                MOI.ScalarAffineTerm(-l.pwr_max[t], uind[t])
            ], 0.0),
            MOI.LessThan(0.0)
        )
    end

    # Heat transfer, t = 1
    MOI.add_constraint(m,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm((l.η / l.C), pnet[1]),
            MOI.ScalarAffineTerm(-1.0, temp[1])
        ], 0.0),
        MOI.EqualTo(-(l.μ / l.C) * l.temp_ext[1] + (-1.0 + l.μ / l.C) * l.temp_init)
    )

    # Heat transfer, t > 1
    for t in 2:T
        MOI.add_constraint(m,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm((l.η / l.C), pnet[t]),
            MOI.ScalarAffineTerm(-1.0, temp[t]),
            MOI.ScalarAffineTerm(1.0 - (l.μ / l.C), temp[t-1])
        ], 0.0),
        MOI.EqualTo(-(l.μ / l.C) * l.temp_ext[t])
    )
    end

    # ==========================================
    #    IV. Add sub-resources to current model
    # ==========================================
    # (None here)

    return nothing
end


text:
using GLPK
using JuMP
function scan_maker(A)
    m= Model(GLPK.Optimizer)

    level = size(A, 2)
    v = zeros(Int, level)
    ub = zeros(Int, level)
    lb = zeros(Int, level)

    @variable(m, x[1:level])
    @constraint(m, con, A*x.>=0)

    function setc(a)
        set_normalized_rhs.(con, float.(a))
    end

    function scan(c::Channel)
        i = 1
        init = 1
        while i > 0
            if i >= init
                @objective(m, Max, x[i])
                res = JuMP.optimize!(m)
                if termination_status(m) == MOI.OPTIMAL || termination_status(m) == MOI.DUAL_INFEASIBLE
                    ub[i] = round(Int, getvalue(x[i]))
                    set_objective_sense(m, MOI.MIN_SENSE)
                    res = JuMP.optimize!(m)
                    @assert termination_status(m) == MOI.OPTIMAL || termination_status(m) == MOI.DUAL_INFEASIBLE
                    lb[i] = round(Int, getvalue(x[i]))

                    v[i] = lb[i]
                    init += 1
                else
                    @assert termination_status(m) == MOI.INFEASIBLE
                    i -= 1
                    continue
                end
            elseif v[i] < ub[i]
                v[i] += 1
            else
                set_upper_bound(x[i], Inf)
                set_lower_bound(x[i], -Inf)
                init -= 1
                i -= 1
                continue
            end

            if i >= level
                put!(c, v)
                continue
            else
                set_upper_bound(x[i], v[i])
                set_lower_bound(x[i], v[i])
                i += 1
            end
        end
        close(c)
    end

    return setc, scan
end


text:
################################################################################
#
#   (10,3)d
#
################################################################################

# Implementation
# - implementation 1
# - labels <: Any
function getUnitcell_10_3_d(
			unitcell_type  :: Type{U},
			implementation :: Val{1}
		) :: U where {LS,LB,S<:AbstractSite{LS,3},B<:AbstractBond{LB,3},U<:AbstractUnitcell{S,B}}

	# define some constants
	a = 0.25*(2 - sqrt(2))
	c = 0.5
	# return a new Unitcell
	return newUnitcell(
		# Type of the unitcell
		U,
		# lattice vectors
		Vector{Float64}[
			Float64[0.5, -0.5, 0],
			Float64[0.5, 0.5, 0.0],
			Float64[0.0, 0.0, c],
		],
		# sites
		S[
			newSite(S, Float64[0.0, -a, 0.75*c], getDefaultLabelN(LS,1)),
			newSite(S, Float64[-a, 0.0, 0.5*c], getDefaultLabelN(LS,2)),
			newSite(S, Float64[0, a, 0.25*c], getDefaultLabelN(LS,3)),
			newSite(S, Float64[a, 0.0, 0.0], getDefaultLabelN(LS,4)),
			newSite(S, Float64[-a, -0.5, 0.25*c], getDefaultLabelN(LS,5)),
			newSite(S, Float64[0, a - 0.5, 0.5*c], getDefaultLabelN(LS,6)),
			newSite(S, Float64[a, -0.5, 0.75*c], getDefaultLabelN(LS,7)),
			newSite(S, Float64[0, -a - 0.5, 0.0], getDefaultLabelN(LS,8)),
		],
		# bonds
		B[
			newBond(B, 1, 2, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 3, 4, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 5, 6, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 2, 1, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 4, 3, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 6, 5, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 2, 3, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 3, 2, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 6, 7, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 7, 6, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 1, 6, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 6, 1, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 4, 5, getDefaultLabel(LB), (0, 1, 0)),
			newBond(B, 5, 4, getDefaultLabel(LB), (0, -1, 0)),
			newBond(B, 8, 5, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 5, 8, getDefaultLabel(LB), (0, 0, 0)),
			newBond(B, 8, 7, getDefaultLabel(LB), (0, 0, -1)),
			newBond(B, 7, 8, getDefaultLabel(LB), (0, 0, 1)),
			newBond(B, 2, 7, getDefaultLabel(LB), (-1, 0, 0)),
			newBond(B, 7, 2, getDefaultLabel(LB), (1, 0, 0)),
			newBond(B, 3, 8, getDefaultLabel(LB), (-1, 1, 0)),
			newBond(B, 8, 3, getDefaultLabel(LB), (1, -1, 0)),
			newBond(B, 1, 4, getDefaultLabel(LB), (0, 0, 1)),
			newBond(B, 4, 1, getDefaultLabel(LB), (0, 0, -1)),
		]
	)
end

# Implementation
# - implementation 4 (Kitaev)
# - labels <: Any
function getUnitcell_10_3_d(
			unitcell_type  :: Type{U},
			implementation :: Val{4}
		) :: U where {LS,LB,S<:AbstractSite{LS,3},B<:AbstractBond{LB,3},U<:AbstractUnitcell{S,B}}

	# define some constants
	a = 0.25*(2 - sqrt(2))
	c = 0.5
	# return a new Unitcell
	return newUnitcell(
		# Type of the unitcell
		U,
		# lattice vectors
		Vector{Float64}[
			Float64[0.5, -0.5, 0],
			Float64[0.5, 0.5, 0.0],
			Float64[0.0, 0.0, c],
		],
		# sites
		S[
			newSite(S, Float64[0.0, -a, 0.75*c], getDefaultLabelN(LS,1)),
			newSite(S, Float64[-a, 0.0, 0.5*c], getDefaultLabelN(LS,2)),
			newSite(S, Float64[0, a, 0.25*c], getDefaultLabelN(LS,3)),
			newSite(S, Float64[a, 0.0, 0.0], getDefaultLabelN(LS,4)),
			newSite(S, Float64[-a, -0.5, 0.25*c], getDefaultLabelN(LS,5)),
			newSite(S, Float64[0, a - 0.5, 0.5*c], getDefaultLabelN(LS,6)),
			newSite(S, Float64[a, -0.5, 0.75*c], getDefaultLabelN(LS,7)),
			newSite(S, Float64[0, -a - 0.5, 0.0], getDefaultLabelN(LS,8)),
		],
		# bonds
		B[
			newBond(B, 1, 2, getDefaultLabelX(LB), (0, 0, 0)),
			newBond(B, 3, 4, getDefaultLabelX(LB), (0, 0, 0)),
			newBond(B, 5, 6, getDefaultLabelY(LB), (0, 0, 0)),
			newBond(B, 2, 1, getDefaultLabelX(LB), (0, 0, 0)),
			newBond(B, 4, 3, getDefaultLabelX(LB), (0, 0, 0)),
			newBond(B, 6, 5, getDefaultLabelY(LB), (0, 0, 0)),
			newBond(B, 2, 3, getDefaultLabelY(LB), (0, 0, 0)),
			newBond(B, 3, 2, getDefaultLabelY(LB), (0, 0, 0)),
			newBond(B, 6, 7, getDefaultLabelX(LB), (0, 0, 0)),
			newBond(B, 7, 6, getDefaultLabelX(LB), (0, 0, 0)),
			newBond(B, 1, 6, getDefaultLabelZ(LB), (0, 0, 0)),
			newBond(B, 6, 1, getDefaultLabelZ(LB), (0, 0, 0)),
			newBond(B, 4, 5, getDefaultLabelZ(LB), (0, 1, 0)),
			newBond(B, 5, 4, getDefaultLabelZ(LB), (0, -1, 0)),
			newBond(B, 8, 5, getDefaultLabelX(LB), (0, 0, 0)),
			newBond(B, 5, 8, getDefaultLabelX(LB), (0, 0, 0)),
			newBond(B, 8, 7, getDefaultLabelY(LB), (0, 0, -1)),
			newBond(B, 7, 8, getDefaultLabelY(LB), (0, 0, 1)),
			newBond(B, 2, 7, getDefaultLabelZ(LB), (-1, 0, 0)),
			newBond(B, 7, 2, getDefaultLabelZ(LB), (1, 0, 0)),
			newBond(B, 3, 8, getDefaultLabelZ(LB), (-1, 1, 0)),
			newBond(B, 8, 3, getDefaultLabelZ(LB), (1, -1, 0)),
			newBond(B, 1, 4, getDefaultLabelY(LB), (0, 0, 1)),
			newBond(B, 4, 1, getDefaultLabelY(LB), (0, 0, -1)),
		]
	)
end


text:
function prepare_sequence(sequence)
    N_list = [1]
    append!(N_list, [get_length(i) for i in sequence.pulses])
    totN = Int(sum(N_list) - 1)
    N_list = cumsum(N_list)
    ΩArr = zeros(totN)
    ϕArr = zeros(totN)
    dtArr = zeros(totN)

    for (i, p) in enumerate(sequence.pulses)
        pulse_t = range(0, p.T, step=p.dt)
        N = get_length(p)
        D = p.(pulse_t)

        @views ΩArr[N_list[i]:N_list[i + 1] - 1] = [D[i][1] for i = 1:N]
        @views ϕArr[N_list[i]:N_list[i + 1] - 1] = [D[i][2] for i = 1:N]
        @views dtArr[N_list[i]:N_list[i + 1] - 1] = repeat([p.dt], N)
    end
    (ΩArr, ϕArr, dtArr)
end


text:
using PairDistribution
using Test
using Random
using DelimitedFiles

function readlfdem(fname)
    pnum, vf, lx, ly, lz = open(fname, "r") do io
        readline(io)
        parse(Int, split(readline(io), " ")[end]),
        parse(Float32, split(readline(io), " ")[end]),
        parse(Float32, split(readline(io), " ")[end]),
        parse(Float32, split(readline(io), " ")[end]),
        parse(Float32, split(readline(io), " ")[end])
    end
    box = Float32[lx,ly,lz]
    X = readdlm(fname, ' ', Float32, skipstart=23, comments=true)
    X, box, pnum, vf
end

function frames(fname)
    X, box, pnum, _  = readlfdem(fname)
    # make frames
    tot = size(X,1)
    fnum = Int(tot/pnum) # number of frames
    F = Array{Float32,3}(undef, pnum, 4, fnum)
    for i in 1:fnum
        F[:,:,i] .= X[(i-1)*pnum+1:i*pnum,2:5]
    end
    F, box
end

@testset "In-Box Distnace" begin
    vbox = [2;2]
    @test PairDistribution.distpbc([0;0],[0;1],vbox) == [0;1]
    @test PairDistribution.distpbc([0;-1],[0;1],vbox) == [0;0]
end

max_stars_repo_name: wildart/PairDistribution

max_issues_repo_name: wildart/PairDistribution

max_forks_repo_name: wildart/PairDistribution


text:
#using Distributed
#addprocs(Sys.CPU_THREADS)

import PyPlot; const plt = PyPlot

import KirchMig

dz, dx = 15, 15
x = 0:dx:1000; nx = length(x)
z = 0:dz:1000; nz = length(z)

rho = 1000ones(nz, nx)
rho[div(330,dx):end,:] .+= 1000
rho[div(670,dz):end,:] .-= 500

vel = 2000 .+ zero(rho)

imp = 1e-3rho.*vel;

refl = [zeros(nx)'
    (imp[2:end,:] .- imp[1:end-1,:])./(imp[2:end,:] .+ imp[1:end-1,:])]

refl[div(nz, 2)-1, div(nx, 2)] += 0.2
refl[div(nz, 2)+1, div(nx, 2)] += 0.2;

fig, ax = plt.subplots(1, 3, figsize=(10,2.5))

# Density
cax = ax[1].imshow(rho, extent=[x[1], x[end], z[end], z[1]], aspect="equal", cmap="Blues", interpolation=nothing)
ax[1].set(xlabel="Position [m]", ylabel="Depth [m]")
cbar = fig.colorbar(cax, ax=ax[1])
cbar.ax.set(ylabel="Density [kg/m³]")

# Velocity
cax = ax[2].imshow(vel, extent=[x[1], x[end], z[end], z[1]], aspect="equal", cmap="Blues", interpolation=nothing)
ax[2].set(xlabel="Position [m]", yticks=[])
cbar = fig.colorbar(cax, ax=ax[2])
cbar.ax.set(ylabel="Velocity [m/s]");

cax = ax[3].imshow(imp, extent=[x[1], x[end], z[end], z[1]], aspect="equal", cmap="Blues", interpolation=nothing)
ax[3].set(xlabel="Position [m]", yticks=[])
cbar = fig.colorbar(cax, ax=ax[3])
cbar.ax.set(ylabel="Impedance [kPa·s/m]");
fig.tight_layout()

xran = x[end]-x[1]
mod = reshape(refl, nz, nx)
vmin, vmax = maximum(abs.(mod))*[-1,1]
fig, ax = plt.subplots(figsize=(6,4))
cax = ax.imshow(mod, extent=[x[1], x[end], z[end], z[1]],
    vmin=vmin, vmax=vmax, aspect="equal", cmap="gray", interpolation=nothing)
cbar = fig.colorbar(cax, ax=ax)
cbar.ax.set(ylabel="Amplitude");
ax.plot(xran*mod[:,div(end,2)]./(4vmax) .+ xran./2, z, color="#d62728");
ax.set(xlim=(x[1], x[end]), ylim=(z[end-1], z[1]), xlabel="Position [m]", ylabel="Depth [m]",
    title="Reflectivity")
fig.tight_layout()

ns = 50

src_z = zeros(ns)
src_x = 0:20:(ns-1)*20

trav = KirchMig.eikonal_const_vel([src_z src_x], z, x, vel[1]);

fig, ax = plt.subplots(figsize=(6,4))
cax = ax.imshow(trav[:,:,div(ns,3)], extent=[x[1], x[end], z[end], z[1]], aspect="equal", cmap="Paired",
    interpolation="bilinear", vmin=0, vmax=maximum(trav[:,:,div(ns,3)]))
cbar = fig.colorbar(cax, ax=ax)
cbar.ax.set(ylabel="Traveltime [s]")
cbar.ax.invert_yaxis()
ax.scatter(src_x[1:4:end], src_z[1:4:end], color="#1f77b4", marker="v", s=100, clip_on=false, zorder=100)
ax.scatter(src_x[3:4:end], src_z[3:4:end], color="#d62728", alpha=0.4, marker="*", s=200, clip_on=false, zorder=100)
ax.scatter(src_x[div(ns,3)], src_z[div(ns,3)], color="#d62728", marker="*", s=200, clip_on=false, zorder=100)
con = ax.contour(trav[:,:,div(ns,3)], extent=[x[1], x[end], z[end], z[1]], origin="upper", colors="k")
plt.clabel(con, fontsize=9, inline=1, fmt="%1.2fs");
ax.set(xlabel="Position [m]", ylabel="Depth [m]", xlim=(x[1], x[end]), ylim=(z[end-1], z[1]))
fig.tight_layout()

t = 0:0.008:1; nt = length(t)
G = KirchMig.KirchMap(t, trav);

@time Gm = G*refl[:]; # refl is 67×67 a array, refl[:] is a 4489-element vector

ricker(t0, f) = @. (1 - 2pi^2 * f^2 * t0^2) * exp(-pi^2 * f^2 * t0^2)
rick_dtt = ricker(t .- t[div(nt,5)], 15);
@views rick_dtt[2:end-1] = (rick_dtt[1:end-2] - 2.0*rick_dtt[2:end-1] + rick_dtt[3:end])/(t[2] - t[1])^2;
rick_dtt /= maximum(abs.(rick_dtt));

W = KirchMig.ConvMap(rick_dtt, ns, ns, nt);

d = W*Gm;

data = reshape(d, ns, ns, nt)

rows = 10
fig, ax = plt.subplots(rows, div(ns, rows), figsize=(12,4rows))
idx = reshape(1:rows*div(ns, rows), div(ns, rows), rows)'
for (i, axi) in enumerate(ax)
    csg = Array(data[:,idx[i],:]')
    cax = axi.imshow(csg, extent=[src_x[1], src_x[end], t[end], t[1]],
        vmin=-0.5maximum(abs.(csg)), vmax=0.5maximum(abs.(csg)), aspect="auto", cmap="gray", interpolation="none")
    axi.scatter(src_x[idx[i]], t[1], color="#d62728", marker="*", s=200, clip_on=false, zorder=100)
    axi.set(xlim=(src_x[1],src_x[end]), ylim=(t[end], t[1]), xlabel="Position [m]",
        xticks = range(minimum(src_x), stop=maximum(src_x), length=5))
    if idx[i] != 1 && idx[i] != div(ns,rows)+1
        axi.set(yticks=[])
    else
        axi.set(ylabel="Time [s]")
    end
end
fig.tight_layout()

L = W*G
@time m_mig = L'd;

xran = x[end]-x[1]
mod_mig = reshape(m_mig, nz, nx)
vmin, vmax = maximum(abs.(mod_mig))*[-1,1]
fig, ax = plt.subplots(figsize=(6,4))
cax = ax.imshow(mod_mig, extent=[x[1], x[end], z[end], z[1]],
    vmin=vmin, vmax=vmax, aspect="equal", cmap="gray", interpolation=nothing)
cbar = fig.colorbar(cax, ax=ax)
cbar.ax.set(ylabel="Amplitude");
ax.plot(xran*mod_mig[:,div(end,2)]./(4vmax) .+ xran./2, z, color="#d62728");
ax.set(xlim=(x[1], x[end]), ylim=(z[end], z[1]), xlabel="Position [m]", ylabel="Depth [m]",
    title="Migrated image")
fig.tight_layout()

@time m_lsm, hist_m, hist_r = KirchMig.cg(L'L, m_mig, maxiter=15, log=true);

xran = x[end]-x[1]
mod_lsm = reshape(m_lsm, nz, nx)
vmin, vmax = maximum(abs.(mod_lsm))*[-1,1]
fig, ax = plt.subplots(figsize=(6,4))
cax = ax.imshow(mod_lsm, extent=[x[1], x[end], z[end], z[1]],
    vmin=vmin, vmax=vmax, aspect="equal", cmap="gray", interpolation=nothing)
cbar = fig.colorbar(cax, ax=ax)
cbar.ax.set(ylabel="Amplitude");
ax.plot(xran*mod_lsm[:,div(end,2)]./(4vmax) .+ xran./2, z, color="#d62728");
ax.set(xlim=(x[1], x[end]), ylim=(z[end], z[1]), xlabel="Position [m]", ylabel="Depth [m]",
    title="Migrated image")
fig.tight_layout()

import LinearAlgebra: norm
model_residuals = [ norm(m - refl[:])/norm(refl) for m in hist_m];

its = 0:(length(hist_r)-1)

fig, ax = plt.subplots(2, 2, figsize=(8,6), sharex=true)
ax[1].plot(its, hist_r, "k")
ax[1].set(ylabel="Data residual")
ax[3].semilogy(its, hist_r, "k")
ax[3].set(ylabel="Data residual")
ax[2].plot(its, model_residuals, "k")
ax[2].set(xlabel="Iteration", ylabel="Model residual")
ax[4].semilogy(its, model_residuals, "k")
ax[4].set(xlim=(its[1], its[end]), xlabel="Iteration", ylabel="Model residual")
ax[4].xaxis.set_major_locator(plt.matplotlib.ticker.MaxNLocator(integer=true))
fig.suptitle("Residuals", y=1.)
fig.tight_layout();

err_mig = refl/maximum(refl) - mod_mig/maximum(mod_mig)
err_lsm = refl/maximum(refl) - mod_lsm/maximum(mod_lsm)

fig, ax = plt.subplots(1, 2, figsize=(12,4))
pclip=.5
cax = ax[1].imshow(err_mig, extent=[x[1], x[end], z[end], z[1]],
    vmin=-pclip, vmax=pclip, aspect="equal", cmap="gray", interpolation="none")
ax[1].set(xlabel="Position [m]", ylabel="Depth [m]", title="Normalized migration error")
cbar = fig.colorbar(cax, ax=ax[1])
cbar.ax.set(ylabel="Amplitude");

vmin, vmax = pclip*maximum(abs.(mod_lsm))*[-1,1]
cax = ax[2].imshow(err_lsm, extent=[x[1], x[end], z[end], z[1]],
    vmin=-pclip, vmax=pclip, aspect="equal", cmap="gray", interpolation="none")
ax[2].set(xlabel="Position [m]",title="Normalized LSM error", yticks=[]);
cbar = fig.colorbar(cax, ax=ax[2])
cbar.ax.set(ylabel="Amplitude");

fig, ax = plt.subplots(2, 1, figsize=(8,4))
ax[1].plot(mod_mig[:,div(end,2)], "k", linewidth=2)
ax[1].set(ylabel="Mig. Amplitude", ylim=[-1.1,1.1]*maximum(abs.(mod_mig[:,div(end,2)])))
ax2 = ax[1].twinx()
ax2.plot(mod_lsm[:,div(end,2)], "--", color="#1f77b4", linewidth=2)
ax2.set(ylim=[-1.1,1.1]*maximum(abs.(mod_lsm[:,div(end,2)])))
ax2.set_ylabel("LSM Amplitude", color="#1f77b4")

ax[2].plot(mod_mig[:,1], "k", linewidth=2)
ax[2].set(ylabel="Mig. Amplitude", ylim=[-1.1,1.1]*maximum(abs.(mod_mig[:,1])))
ax2 = ax[2].twinx()
ax2.plot(mod_lsm[:,1], "--", color="#1f77b4", linewidth=2)
ax2.set(ylim=[-1.1,1.1]*maximum(abs.(mod_lsm[:,1])))
ax2.set_ylabel("LSM Amplitude", color="#1f77b4");


text:
# Weave readme
using Pkg
cd("c:/git/CuCountMap")
Pkg.activate("c:/git/CuCountMap/readme-env")

using Weave

weave("README.jmd", out_path=:pwd, doctype="github")

if false
    tangle("README.jmd")
end


text:
#  Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at https://mozilla.org/MPL/2.0/.

# With `Iterators.ProductIterator`, everything works except when the
# `IteratorSize` of an iterator is not `Base.HasShape{1}` or `Base.HasLength`
# which is notably the case for scalars for it is `Base.HasShape{0}` and
# multidimensional arrays for which it is `Base.HasShape{N}` where `N` is the
# dimension. For instance:
# ```julia
# julia> collect(Iterators.product(2, 3))
# 0-dimensional Array{Tuple{Int64,Int64},0}:
# (2, 3)
# ```
# while we could like it to be a `3`-dimensional array of size `(1, 1)`.
# When the user does `@container([2, 3], 1)`, a `DenseAxisArray` of size
# `(1, 1)`. Another example:
# ```julia
# julia> collect(Iterators.product([1 2; 3 4]))
# 2×2 Array{Tuple{Int64},2}:
# (1,)  (2,)
# (3,)  (4,)
# ```
# while we need the size to be `(4,)`, not `(2, 2)` when the user does
# `@container([i = [1, 2; 3 4]], i^2)`.
# Long story short, we want to tried everything as a interator without shape
# while `Iterators.ProductIterator` does care about preserving the shape
# when doing the cartesian product.
"""
    struct VectorizedProductIterator{T}
        prod::Iterators.ProductIterator{T}
    end

Cartesian product of the iterators `prod.iterators`. It is the same iterator as
`Base.Iterators.ProductIterator` except that it is independent of the
`IteratorSize` of the elements of `prod.iterators`.
For instance:
* `size(Iterators.product(1, 2))` is `tuple()` while
  `size(VectorizedProductIterator(1, 2))` is `(1, 1)`.
* `size(Iterators.product(ones(2, 3)))` is `(2, 3)` while
  `size(VectorizedProductIterator(ones(2, 3)))` is `(1, 1)`.
"""
struct VectorizedProductIterator{T}
    prod::Iterators.ProductIterator{T}
end
function vectorized_product(iterators...)
    return VectorizedProductIterator(Iterators.product(iterators...))
end
function Base.IteratorSize(
    ::Type{<:VectorizedProductIterator{<:Tuple{Vararg{Any,N}}}},
) where {N}
    return Base.HasShape{N}()
end
Base.IteratorEltype(::Type{<:VectorizedProductIterator}) = Base.EltypeUnknown()
Base.size(it::VectorizedProductIterator) = _prod_size(it.prod.iterators)
_prod_size(::Tuple{}) = ()
_prod_size(t::Tuple) = (length(t[1]), _prod_size(Base.tail(t))...)
Base.axes(it::VectorizedProductIterator) = _prod_indices(it.prod.iterators)
_prod_indices(::Tuple{}) = ()
function _prod_indices(t::Tuple)
    return (Base.OneTo(length(t[1])), _prod_indices(Base.tail(t))...)
end
Base.ndims(it::VectorizedProductIterator) = length(axes(it))
Base.length(it::VectorizedProductIterator) = prod(size(it))
Base.iterate(it::VectorizedProductIterator, args...) = iterate(it.prod, args...)
_eltype_or_any(it::VectorizedProductIterator) = eltype(it.prod)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using YaoBase, YaoArrayRegister, Random
using BitBasis
export Measure,
    MeasureAndReset, AllLocs, ComputationalBasis, chmeasureoperator, num_measured

"""
    Measure{D,K, OT, LT, PT, RNG} <: PrimitiveBlock{D}
    Measure(n::Int; rng=Random.GLOBAL_RNG, operator=ComputationalBasis(), locs=1:n, resetto=nothing, remove=false, nlevel=2)

Measure operator, currently only qudits are supported.
"""
mutable struct Measure{D,K,OT,LT<:Union{NTuple{K,Int},AllLocs},PT<:PostProcess,RNG} <:
               PrimitiveBlock{D}
    n::Int
    rng::RNG
    operator::OT
    locations::LT
    postprocess::PT
    results::Any
    function Measure{D,K,OT,LT,PT,RNG}(n::Int,
        rng::RNG,
        operator::OT,
        locations::LT,
        postprocess::PT,
    ) where {RNG,D,K,OT,LT<:Union{NTuple{K,Int},AllLocs},PT<:PostProcess}
        locations isa AllLocs || @assert_locs_safe n locations
        new{D,K,OT,LT,PT,RNG}(n, rng, operator, locations, postprocess)
    end
end
nqudits(m::Measure) = m.n

function Measure{D}(n::Int,
    rng::RNG,
    operator::OT,
    locations::LT,
    postprocess::PT,
) where {D,RNG,OT,LT<:Union{NTuple{K,Int} where K,AllLocs},PT<:PostProcess}
    k = locations isa AllLocs ? n : length(locations)
    Measure{D,k,OT,LT,PT,RNG}(n, rng, operator, locations, postprocess)
end

const MeasureAndReset{D,K,OT,LT,RNG} = Measure{D,K,OT,LT,ResetTo{BitStr64{K}},RNG}
function MeasureAndReset(
    N::Int,
    resetto = 0;
    operator = ComputationalBasis,
    rng = Random.GLOBAL_RNG,
    nlevel=2,
)
    Measure{nlevel}(N; nlevel=nlevel, postprocess = ResetTo(resetto), operator = operator, rng = rng, locs = locs)
end

"""
    chmeasureoperator(m::Measure, op::AbstractBlock)

change the measuring `operator`. It will also discard existing measuring results.
"""
function chmeasureoperator(m::Measure{D}, op::AbstractBlock) where D
    Measure{D}(m.n, m.rng, op, m.locations, m.postprocess)
end

function Base.:(==)(m1::Measure, m2::Measure)
    res =
        m1.rng == m2.rng &&
        m1.operator == m2.operator &&
        m1.locations == m2.locations &&
        m1.postprocess == m2.postprocess
    res = res && isdefined(m1, :results) == isdefined(m2, :results)
    res && (!isdefined(m1, :results) || m1.results == m2.results)
end

num_measured(::Measure{D,K}) where {D,K} = K

"""
    Measure(n::Int; rng=Random.GLOBAL_RNG, operator=ComputationalBasis(), locs=AllLocs(), resetto=nothing, remove=false)

Create a `Measure` block with number of qudits `n`.

# Example

You can create a `Measure` block on given basis (default is the computational basis).

```jldoctest; setup=:(using YaoBlocks)
julia> Measure(4)
Measure(4)
```

Or you could specify which qudits you are going to measure

```jldoctest; setup=:(using YaoBlocks)
julia> Measure(4; locs=1:3)
Measure(4;locs=(1, 2, 3))
```

by default this will collapse the current register to measure results.

```jldoctest; setup=:(using YaoBlocks, YaoArrayRegister, Random; Random.seed!(2))
julia> r = normalize!(arrayreg(bit"000") + arrayreg(bit"111"))
ArrayReg{2, ComplexF64, Array...}
    active qudits: 3/3
    nlevel: 2

julia> state(r)
8×1 Matrix{ComplexF64}:
 0.7071067811865475 + 0.0im
                0.0 + 0.0im
                0.0 + 0.0im
                0.0 + 0.0im
                0.0 + 0.0im
                0.0 + 0.0im
                0.0 + 0.0im
 0.7071067811865475 + 0.0im

julia> r |> Measure(3)
ArrayReg{2, ComplexF64, Array...}
    active qudits: 3/3
    nlevel: 2

julia> state(r)
8×1 Matrix{ComplexF64}:
 1.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
```

But you can also specify the target bit configuration you want to collapse to with keyword `resetto`.

```jldoctest; setup=:(using YaoBlocks; using BitBasis)
julia> m = Measure(4; resetto=bit"0101")
Measure(4;postprocess=ResetTo{BitStr{4,Int64}}(0101 ₍₂₎))

julia> m.postprocess
ResetTo{BitStr{4,Int64}}(0101 ₍₂₎)```
"""
function Measure(
    n::Int;
    rng::RNG = Random.GLOBAL_RNG,
    operator::OT = ComputationalBasis(),
    locs = AllLocs(),
    resetto = nothing,
    remove = false,
    nlevel = 2,
) where {OT,RNG}
    if resetto !== nothing
        if remove
            error(
                "invalid keyword combination, expect resetto or remove, got (resetto=$resetto, remove=true)",
            )
        else
            postprocess = ResetTo(BitStr64{locs isa AllLocs ? n : length(locs)}(resetto))
        end
    else
        if remove
            postprocess = RemoveMeasured()
        else
            postprocess = NoPostProcess()
        end
    end
    if locs isa AllLocs
        K = n
    else
        locs = (locs...,)
        K = length(locs)
    end
    Measure{nlevel,K,OT,typeof(locs),typeof(postprocess),RNG}(n,rng, operator, locs, postprocess)
end

Measure(;
    rng = Random.GLOBAL_RNG,
    locs = AllLocs(),
    operator = ComputationalBasis(),
    resetto = nothing,
    remove = false,
    nlevel = 2,
) = @λ(
    n -> Measure(
        n;
        rng = rng,
        locs = locs,
        operator = operator,
        resetto = resetto,
        remove = remove,
        nlevel = nlevel,
    )
)
mat(x::Measure) = error("use BlockMap to get its matrix.")

function _apply!(r::AbstractRegister{D}, m::Measure{D}) where {D}
    m.results = measure!(m.postprocess, m.operator, r, m.locations; rng = m.rng)
    return r
end

occupied_locs(m::Measure) = m.locations isa AllLocs ? (1:nqudits(m)...,) : m.locations


text:
# BBoxSet

export BBoxSet
abstract type BBoxSet{D} end


text:
import LibGit2
import PkgBenchmark
import Statistics

# include("./utils/github/_httpjson_github_api_unauthenticated.jl")
# include("./utils/github/_httpjson_github_api_authenticated.jl")

struct AllowedToIgnoreThisError end

_single_line_travis_ignore_errors(x::AbstractString) = _single_line_travis_ignore_errors(convert(String, x))

function _single_line_travis_ignore_errors(line::String)::Bool
    _line::String = strip(line)
    _regex_ignore_errors = r"^\d*: \[IGNORE_ERRORS_BENCHMARK_STAGE\]"
    ignore_errors::Bool = occursin(_regex_ignore_errors, _line)
    return ignore_errors
end

travis_ignore_errors(x::AbstractString) = travis_ignore_errors(convert(String, x))

function travis_ignore_errors(commit_message::String)::Bool
    lines::Vector{String} = split(strip(commit_message), "\n")
    vector_ignore_errors::Vector{Bool} = Vector{Bool}(undef, 0)
    for line in lines
        _line = strip(line)
        if isempty(_line)
        elseif startswith(_line, "Merge #")
        elseif startswith(_line, "Try #")
        elseif startswith(_line, "Co-authored-by:")
        else
            line_ignore_errors = _single_line_travis_ignore_errors(_line)
            push!(vector_ignore_errors, line_ignore_errors)
        end
    end
    ignore_errors = _all_and_notempty(vector_ignore_errors)
    return ignore_errors
end

function get_travis_git_commit_message(a::AbstractDict = ENV)::String
    result::String = strip(get(a, "TRAVIS_COMMIT_MESSAGE", ""))
    return result
end

_single_line_travis_allow_regressions(x::AbstractString) = _single_line_travis_allow_regressions(convert(String, x))

function _single_line_travis_allow_regressions(line::String)::Tuple{Bool, Bool}
    _line::String = strip(line)
    _regex_allow_onlytime_regressions = r"^\d*: \[ALLOW_TIME_REGRESSIONS\]"
    _regex_allow_onlymemory_regressions = r"^\d*: \[ALLOW_MEMORY_REGRESSIONS\]"
    _regex_allow_bothtimeandmemory_regressions = r"^\d*: \[ALLOW_TIME\+MEMORY_REGRESSIONS\]"
    _allow_onlytime_regressions::Bool = occursin(_regex_allow_onlytime_regressions, _line)
    _allow_onlymemory_regressions::Bool = occursin(_regex_allow_onlymemory_regressions, _line)
    _allow_bothtimeandmemory_regressions::Bool = occursin(_regex_allow_bothtimeandmemory_regressions, _line)
    allow_time_regressions::Bool = _allow_onlytime_regressions || _allow_bothtimeandmemory_regressions
    allow_memory_regressions::Bool = _allow_onlymemory_regressions || _allow_bothtimeandmemory_regressions
    return allow_time_regressions, allow_memory_regressions
end

travis_allow_regressions(x::AbstractString) = travis_allow_regressions(convert(String, x))

_all_and_notempty(x) = !isempty(x) && all(x)

function travis_allow_regressions(commit_message::String)::Tuple{Bool, Bool}
    lines::Vector{String} = split(strip(commit_message), "\n")
    vector_allow_time_regressions::Vector{Bool} = Vector{Bool}(undef, 0)
    vector_allow_memory_regressions::Vector{Bool} = Vector{Bool}(undef, 0)
    for line in lines
        _line = strip(line)
        if isempty(_line)
        elseif startswith(_line, "Merge #")
        elseif startswith(_line, "Try #")
        elseif startswith(_line, "Co-authored-by:")
        else
            line_allow_time_regressions, line_allow_memory_regressions = _single_line_travis_allow_regressions(_line)
            push!(vector_allow_time_regressions, line_allow_time_regressions)
            push!(vector_allow_memory_regressions, line_allow_memory_regressions)
        end
    end
    allow_time_regressions = _all_and_notempty(vector_allow_time_regressions)
    allow_memory_regressions = _all_and_notempty(vector_allow_memory_regressions)
    return allow_time_regressions, allow_memory_regressions
end

function benchmarkpkg_do_not_ignore_errors(pkg_name, git_identifier)
    result = PkgBenchmark.benchmarkpkg(pkg_name, git_identifier; retune = true)
    return result
end

function benchmarkpkg_ignore_errors(pkg_name, git_identifier)
    try
        result = PkgBenchmark.benchmarkpkg(pkg_name, git_identifier; retune = true)
        return result
    catch ex
        showerror(stderr, ex)
        Base.show_backtrace(stderr, catch_backtrace())
        println(stderr)
        return AllowedToIgnoreThisError()
    end
end

function judge_do_not_ignore_errors(target_run, baseline_run, f)
    result = PkgBenchmark.judge(target_run, baseline_run, f)
    return result
end

function judge_ignore_errors(target_run, baseline_run, f)
    try
        result = PkgBenchmark.judge(target_run, baseline_run, f)
        return result
    catch ex
        showerror(stderr, ex)
        Base.show_backtrace(stderr, catch_backtrace())
        println(stderr)
        return AllowedToIgnoreThisError()
    end
end

function all_versions(repo_path::AbstractString)
    _repo_path::String = strip(repo_path)
    _repo = LibGit2.GitRepoExt(_repo_path)
    _ref_list = LibGit2.ref_list(_repo)
    _all_versions = Vector{VersionNumber}(undef, 0)
    for _ref in _ref_list
        if occursin(r"^refs\/tags\/v(\d*).(\d*).(\d*)$", _ref)
            _m = match(r"^refs\/tags\/v(\d*).(\d*).(\d*)$", _ref)
            _v = VersionNumber("$(_m[1]).$(_m[2]).$(_m[3])")
            push!(_all_versions, _v)
        end
    end
    unique!(_all_versions)
    sort!(_all_versions)
    return _all_versions
end

latest_semver_version(repo_path) = maximum(all_versions(repo_path))

function _run_benchmarks(
        ;
        target::Union{String, PkgBenchmark.BenchmarkConfig},
        baseline::Union{String, PkgBenchmark.BenchmarkConfig},
        )::Nothing
    # git_commit_message::String = _httpjson_get_github_pull_request_title_unauthenticated()
    # git_commit_message::String = _httpjson_get_github_pull_request_title_authenticated()
    git_commit_message::String = get_travis_git_commit_message()

    ignore_errors = travis_ignore_errors(git_commit_message)
    allow_time_regressions, allow_memory_regressions = travis_allow_regressions(git_commit_message)

    @info("Ignore errors during the benchmark suite: $(ignore_errors)")
    @info("Allow time regressions: $(allow_time_regressions)")
    @info("Allow memory regressions: $(allow_memory_regressions)")
    @info("Target: $(target)")
    @info("Baseline: $(baseline)")

    project_root = dirname(dirname(@__FILE__))

    proof_of_concept_dataframes = joinpath(project_root, "test", "integration-tests", "test-proof-of-concept-dataframes.jl")
    proof_of_concept_linearmodel = joinpath(project_root, "test", "integration-tests", "test-proof-of-concept-linearmodel.jl")
    proof_of_concept_mlj = joinpath(project_root, "test", "integration-tests", "test-proof-of-concept-mlj.jl")

    # run each of the scripts once to force compilation of the functions
    include(proof_of_concept_dataframes)
    include(proof_of_concept_linearmodel)
    include(proof_of_concept_mlj)

    if ignore_errors
        results_of_run_on_target = benchmarkpkg_ignore_errors("ModelSanitizer", target)
        results_of_run_on_baseline = benchmarkpkg_ignore_errors("ModelSanitizer", baseline)
    else
        results_of_run_on_target = benchmarkpkg_do_not_ignore_errors("ModelSanitizer", target)
        results_of_run_on_baseline = benchmarkpkg_do_not_ignore_errors("ModelSanitizer", baseline)
    end
    if ( ignore_errors ) && ( (results_of_run_on_target isa AllowedToIgnoreThisError) || (results_of_run_on_baseline isa AllowedToIgnoreThisError) )
        @error("One or more errors occurred that prevented us from continuing.")
    else
        if ignore_errors
            judgement_minimum = judge_ignore_errors(results_of_run_on_target, results_of_run_on_baseline, minimum)
            judgement_median = judge_ignore_errors(results_of_run_on_target, results_of_run_on_baseline, Statistics.median)
        else
            judgement_minimum = judge_do_not_ignore_errors(results_of_run_on_target, results_of_run_on_baseline, minimum)
            judgement_median = judge_do_not_ignore_errors(results_of_run_on_target, results_of_run_on_baseline, Statistics.median)
        end
        if ( ignore_errors ) && ( (judgement_minimum isa AllowedToIgnoreThisError) || (judgement_median isa AllowedToIgnoreThisError) )
            @error("One or more errors occurred that prevented us from continuing.")
        else
            println(stdout)
            println(stdout, "<!--- Package: `ModelSanitizer`. Target: `$(target)`. Baseline: `$(baseline)`. Estimator function: `f = minimum`. Invariant rows are not shown. -->")
            println(stdout, "<!--- BEGIN MARKDOWN -->")
            PkgBenchmark.export_markdown(stdout, judgement_minimum; export_invariants = true)
            println(stdout)
            println(stdout, "<!--- END MARKDOWN -->")
            println(stdout, "<!--- Package: `ModelSanitizer`. Target: `$(target)`. Baseline: `$(baseline)`. Estimator function: `f = median`. Invariant rows are not shown. -->")
            println(stdout, "<!--- BEGIN MARKDOWN -->")
            PkgBenchmark.export_markdown(stdout, judgement_median; export_invariants = true)
            println(stdout)
            println(stdout, "<!--- END MARKDOWN -->")
            println(stdout)

            this_judgement_was_failed_for_time = false
            this_judgement_was_failed_for_memory = false

            for i in ["integration-tests"]
                for j in ["proof-of-concept-dataframes", "proof-of-concept-linearmodel", "proof-of-concept-mlj"]
                    trial_judgement = PkgBenchmark.benchmarkgroup(judgement_minimum).data[i].data[j]
                    if PkgBenchmark.time(trial_judgement) == :regression
                        if allow_time_regressions
                            @error("Time regression (allowed) detected in $(i)/$(j)", trial_judgement)
                        else
                            this_judgement_was_failed_for_time = true
                            @error("Time regression detected in $(i)/$(j)", trial_judgement)
                        end
                    end
                    if PkgBenchmark.memory(trial_judgement) == :regression
                        if allow_memory_regressions
                            @error("Memory regression (allowed) detected in $(i)/$(j)", trial_judgement)
                        else
                            this_judgement_was_failed_for_memory = true
                            @error("Memory regression regression detected in $(i)/$(j)", trial_judgement)
                        end
                    end
                end
            end

            if this_judgement_was_failed_for_time || this_judgement_was_failed_for_memory
                error_message = string("FAILURE: ",
                                       "One or more fatal performance ",
                                       "performance regressions were detected.\n",
                                       "To ignore only time regressions, ",
                                       "begin your pull request title with ",
                                       "\"",
                                       "[ALLOW_TIME_REGRESSIONS]",
                                       "\" (without the quotation marks).\n",
                                       "To ignore only memory regressions, ",
                                       "begin your pull request title with ",
                                       "\"",
                                       "[ALLOW_MEMORY_REGRESSIONS]",
                                       "\".\n",
                                       "To ignore both time and memory regressions, ",
                                       "begin your pull request title with ",
                                       "\"",
                                       "[ALLOW_TIME+MEMORY_REGRESSIONS]",
                                       "\".\n")
                travis_branch = lowercase(strip(get(ENV, "TRAVIS_BRANCH", "")))
                travis_pull_request = lowercase(strip(get(ENV, "TRAVIS_PULL_REQUEST", "")))
                if travis_branch == "trying" && travis_pull_request == "false"
                    @error(error_message)
                else
                    error(error_message)
                end
            else
                @info("SUCCESS: No fatal performance regressions were detected.")
            end
        end
    end
    return nothing
end

function run_benchmarks(repo_path = pwd())::Nothing
    travis_branch::String = lowercase(strip(get(ENV, "TRAVIS_BRANCH", "")))
    travis_pull_request::String = lowercase(strip(get(ENV, "TRAVIS_PULL_REQUEST", "")))
    latest_version::String = string("v", latest_semver_version(repo_path))
    target::String = "HEAD"
    if travis_branch == "master" && travis_pull_request == "false"
        baseline = latest_version
    else
        baseline = "master"
    end
    _run_benchmarks(; target = target, baseline = baseline)
    return nothing
end

run_benchmarks()


text:
s = """
s = "a" + "b"
"""
msgs = lintstr(s)
@test msgs[1].code == :E422
@test contains(msgs[1].message, "string uses * to concatenate")

s = """
function f(x)
    Dict("a" + "b" => x)
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E422
@test contains(msgs[1].message, "string uses * to concatenate")

@test messageset(lintstr("""
s = String(1)
""")) == Set([:E539])

s = """
b = string(12)
s = "a" + b
"""
msgs = lintstr(s)
@test msgs[1].code == :E422
@test contains(msgs[1].message, "string uses * to concatenate")

s = """
function f()
    b = repeat(" ", 10)
    @lintpragma("Info type b")
    b
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271
@test contains(msgs[1].message, "typeof(b) == $(Compat.String)")

u = """
안녕하세요 = "Hello World"

Hello
World
"""
msgs = lintstr(u)
@test msgs[1].code == :E321
@test msgs[1].variable == "Hello"
@test msgs[2].code == :E321
@test msgs[2].variable == "World"


text:
struct Ellipse{Dim,T} <: Primitive{Dim,T}
    center::Point{Dim,T}
    weights::Point{Dim,T}
    radius::T
end

Ellipse(center::NTuple{Dim,T}, weights::NTuple{Dim,T}, radius::T) where {Dim,T} = Ellipse(Point(center), Point(weights), radius)

struct HalfSpace{Dim,T}
    p::Vec{Dim,T}
    offset::T
end

norm2(x) = sum(map(*, x, x))

sdf(h::HalfSpace) = TrueSDF(x -> sum(coordinates(x) .* h.p) - h.offset)
sdf(e::Ellipse) = TrueSDF(x -> hypot(((e.center - x) ./ coordinates(e.weights))...) - e.radius)
sdf(s::Sphere) = TrueSDF(x -> hypot((s.center - x)...) - s.radius)

function sdf(b::Box)
    function _sdf(x)
        c = center(b)
        radii = b.max - c
        p = x - c
        z = zero(typeof(p))
        q = abs.(p) - radii
        hypot(maximum.(zip(q, z))...) + min(maximum(q), zero(eltype(p)))
    end
    TrueSDF(_sdf)
end


text:
#######################################################################
# BoxIterator
#######################################################################

"""
    BoxIterator(seed::SVector{N,Int}; maxiterations = missing)

Cartesian iterator `iter` over `SVector{N,Int}`s (`cell`s) that starts at `seed` and
grows outwards in the form of a box of increasing sides (not necesarily equal) until
it encompasses a certain N-dimensional region. To signal that a cell is in the desired
region the user calls `acceptcell!(iter, cell)`.
"""
struct BoxIterator{N,MI<:Union{Int,Missing}}
    seed::SVector{N,Int}
    maxiter::MI
    dimdir::MVector{2,Int}
    nmoves::MVector{N,Bool}
    pmoves::MVector{N,Bool}
    npos::MVector{N,Int}
    ppos::MVector{N,Int}
end

Base.IteratorSize(::BoxIterator) = Base.SizeUnknown()

Base.IteratorEltype(::BoxIterator) = Base.HasEltype()

Base.eltype(::BoxIterator{N}) where {N} = SVector{N,Int}

Base.CartesianIndices(b::BoxIterator) =
    CartesianIndices(UnitRange.(Tuple(b.npos), Tuple(b.ppos)))

function BoxIterator(seed::SVector{N}; maxiterations = missing) where {N}
    BoxIterator(seed, maxiterations, MVector(1, 2),
        ones(MVector{N,Bool}), ones(MVector{N,Bool}), MVector{N,Int}(seed), MVector{N,Int}(seed))
end

function iteratorreset!(b::BoxIterator{N}) where {N}
    b.dimdir[1] = 1
    b.dimdir[2] = 2
    b.nmoves .= ones(MVector{N,Bool})
    b.pmoves .= ones(MVector{N,Bool})
    b.npos   .= MVector(b.seed)
    b.ppos   .= MVector(b.seed)
    return nothing
end

struct BoxIteratorState{N}
    range::CartesianIndices{N, NTuple{N,UnitRange{Int}}}
    rangestate::CartesianIndex{N}
    iteration::Int
end

Base.iterate(b::BoxIterator{0}) = (SVector{0,Int}(), nothing)
Base.iterate(b::BoxIterator{0}, state) = nothing

function Base.iterate(b::BoxIterator{N}) where {N}
    range = CartesianIndices(ntuple(i -> b.seed[i]:b.seed[i], Val(N)))
    itrange = iterate(range)
    if itrange === nothing
        return nothing
    else
        (cell, rangestate) = itrange
        return (SVector(Tuple(cell)), BoxIteratorState(range, rangestate, 1))
    end
end

function Base.iterate(b::BoxIterator{N}, s::BoxIteratorState{N}) where {N}
    itrange = iterate(s.range, s.rangestate)
    facedone = itrange === nothing
    if facedone
        alldone = !any(b.pmoves) && !any(b.nmoves) || isless(b.maxiter, s.iteration)
        if alldone  # Last shells in all directions were empty, trim from boundingboxcorners
            b.npos .+= 1
            b.ppos .-= 1
            return nothing
        else
            newrange = nextface!(b)
            # newrange === nothing && return nothing
            itrange = iterate(newrange)
            # itrange === nothing && return nothing
            (cell, rangestate) = itrange
            return (SVector(Tuple(cell)), BoxIteratorState(newrange, rangestate, s.iteration + 1))
        end
    else
        (cell, rangestate) = itrange
        return (SVector(Tuple(cell)), BoxIteratorState(s.range, rangestate, s.iteration + 1))
    end
end

function nextface!(b::BoxIterator{N}) where {N}
    @inbounds for i in 1:2N
        nextdimdir!(b)
        newdim, newdir = Tuple(b.dimdir)
        if newdir == 1
            if b.nmoves[newdim]
                b.npos[newdim] -= 1
                b.nmoves[newdim] = false
                return newrangeneg(b, newdim)
            end
        else
            if b.pmoves[newdim]
                b.ppos[newdim] += 1
                b.pmoves[newdim] = false
                return newrangepos(b, newdim)
            end
        end
    end
    return nothing
end

function nextdimdir!(b::BoxIterator{N}) where {N}
    dim, dir = Tuple(b.dimdir)
    if dim < N
        dim += 1
    else
        dim = 1
        dir = ifelse(dir == 1, 2, 1)
    end
    b.dimdir[1] = dim
    b.dimdir[2] = dir
    return nothing
end

@inline function newrangeneg(b::BoxIterator{N}, dim) where {N}
    return CartesianIndices(ntuple(
        i -> b.npos[i]:(i == dim ? b.npos[i] : b.ppos[i]),
        Val(N)))
end

@inline function newrangepos(b::BoxIterator{N}, dim) where {N}
    return CartesianIndices(ntuple(
        i -> (i == dim ? b.ppos[i] : b.npos[i]):b.ppos[i],
        Val(N)))
end

function acceptcell!(b::BoxIterator{N}, cell) where {N}
    dim, dir = Tuple(b.dimdir)
    if dir == 1
        @inbounds for i in 1:N
            (cell[i] == b.ppos[i]) && (b.pmoves[i] = true)
            (i == dim || cell[i] == b.npos[i]) && (b.nmoves[i] = true)
        end
    else
        @inbounds for i in 1:N
            (i == dim || cell[i] == b.ppos[i]) && (b.pmoves[i] = true)
            (cell[i] == b.npos[i]) && (b.nmoves[i] = true)
        end
    end
    return nothing
end

# Fallback for non-BoxIterators
acceptcell!(b, cell) = nothing

#######################################################################
# CoSort
#######################################################################

struct CoSortTup{T,Tv}
    x::T
    y::Tv
end

mutable struct CoSort{T,Tv,S<:AbstractVector{T},C<:AbstractVector{Tv}} <: AbstractVector{CoSortTup{T,Tv}}
    sortvector::S
    covector::C
    offset::Int
    function CoSort{T,Tv,S,C}(sortvector, covector, offset) where {T,Tv,S,C}
        length(covector) >= length(sortvector) ? new(sortvector, covector, offset) :
            throw(DimensionMismatch("Coarray length should exceed sorting array"))
    end
end

CoSort(sortvector::S, covector::C) where {T,Tv,S<:AbstractVector{T},C<:AbstractVector{Tv}} =
    CoSort{T,Tv,S,C}(sortvector, covector, 0)

Base.size(c::CoSort) = (size(c.sortvector, 1) - c.offset,)

Base.getindex(c::CoSort, i) =
    CoSortTup(getindex(c.sortvector, i + c.offset), getindex(c.covector, i + c.offset))

Base.setindex!(c::CoSort, t::CoSortTup, i) =
    (setindex!(c.sortvector, t.x, i + c.offset); setindex!(c.covector, t.y, i + c.offset); c)

Base.isless(a::CoSortTup, b::CoSortTup) = isless(a.x, b.x)

Base.Sort.defalg(v::C) where {T<:Union{Number, Missing}, C<:CoSort{T}} = Base.DEFAULT_UNSTABLE

isgrowing(c::CoSort) = isgrowing(c.sortvector, c.offset + 1)

#######################################################################
# SparseMatrixBuilder
#######################################################################

mutable struct SparseMatrixBuilder{T,S<:AbstractSparseMatrix{T}}
    matrix::S
    colcounter::Int
    rowvalcounter::Int
    cosorter::CoSort{Int,T,Vector{Int},Vector{T}}
end

struct UnfinalizedSparseMatrixCSC{Tv,Ti} <: AbstractSparseMatrix{Tv,Ti}
    m::Int
    n::Int
    colptr::Vector{Ti}
    rowval::Vector{Ti}
    nzval::Vector{Tv}
end

function SparseMatrixBuilder{Tv}(m, n, nnzguess = missing) where {Tv}
    colptr = [1]
    rowval = Int[]
    nzval = Tv[]
    matrix = UnfinalizedSparseMatrixCSC(m, n, colptr, rowval, nzval)
    builder = SparseMatrixBuilder(matrix, 1, 0, CoSort(rowval, nzval))
    nnzguess === missing || sizehint!(builder, nnzguess)
    return builder
end

# Unspecified size constructor
SparseMatrixBuilder{Tv}() where Tv = SparseMatrixBuilder{Tv}(0, 0)

function SparseMatrixBuilder(s::SparseMatrixCSC{Tv,Int}) where {Tv}
    colptr = getcolptr(s)
    rowval = rowvals(s)
    nzval = nonzeros(s)
    nnzguess = length(nzval)
    resize!(rowval, 0)
    resize!(nzval, 0)
    resize!(colptr, 1)
    colptr[1] = 1
    builder = SparseMatrixBuilder(s, 1, 0, CoSort(rowval, nzval))
    sizehint!(builder, nnzguess)
    return builder
end

SparseArrays.getcolptr(s::UnfinalizedSparseMatrixCSC) = s.colptr
SparseArrays.nonzeros(s::UnfinalizedSparseMatrixCSC) = s.nzval
SparseArrays.rowvals(s::UnfinalizedSparseMatrixCSC) = s.rowval
Base.size(s::UnfinalizedSparseMatrixCSC) = (s.m, s.n)
Base.size(s::UnfinalizedSparseMatrixCSC, k) = size(s)[k]

function Base.sizehint!(s::SparseMatrixBuilder, n)
    sizehint!(getcolptr(s.matrix), n + 1)
    sizehint!(nonzeros(s.matrix), n)
    sizehint!(rowvals(s.matrix), n)
    return s
end

function pushtocolumn!(s::SparseMatrixBuilder, row::Int, x, skipdupcheck::Bool = true)
    nrows = size(s.matrix, 1)
    nrows == 0 || 1 <= row <= size(s.matrix, 1) || throw(ArgumentError("tried adding a row $row out of bounds ($(size(s.matrix, 1)))"))
    if skipdupcheck || !isintail(row, rowvals(s.matrix), getcolptr(s.matrix)[s.colcounter])
        push!(rowvals(s.matrix), row)
        push!(nonzeros(s.matrix), x)
        s.rowvalcounter += 1
    end
    return s
end

function isintail(element, container, start::Int)
    for i in start:length(container)
        container[i] == element && return true
    end
    return false
end

function finalizecolumn!(s::SparseMatrixBuilder, sortcol::Bool = true)
    size(s.matrix, 2) > 0 && s.colcounter > size(s.matrix, 2) && throw(DimensionMismatch("Pushed too many columns to matrix"))
    if sortcol
        s.cosorter.offset = getcolptr(s.matrix)[s.colcounter] - 1
        sort!(s.cosorter)
        isgrowing(s.cosorter) || throw(error("Internal error: repeated rows"))
    end
    s.colcounter += 1
    push!(getcolptr(s.matrix), s.rowvalcounter + 1)
    return nothing
end

function finalizecolumn!(s::SparseMatrixBuilder, ncols::Int)
    for _ in 1:ncols
        finalizecolumn!(s)
    end
    return nothing
end

function SparseArrays.sparse(s::SparseMatrixBuilder{<:Any,<:SparseMatrixCSC})
    completecolptr!(getcolptr(s.matrix), size(s.matrix, 2), s.rowvalcounter)
    return s.matrix
end

function completecolptr!(colptr, cols, lastrowptr)
    colcounter = length(colptr)
    if colcounter < cols + 1
        resize!(colptr, cols + 1)
        for col in (colcounter + 1):(cols + 1)
            colptr[col] = lastrowptr + 1
        end
    end
    return colptr
end

function SparseArrays.sparse(s::SparseMatrixBuilder{<:Any,<:UnfinalizedSparseMatrixCSC})
    m, n = size(s.matrix)
    rowval = rowvals(s.matrix)
    colptr = getcolptr(s.matrix)
    nzval = nonzeros(s.matrix)
    if m != 0 && n != 0
        completecolptr!(colptr, n, s.rowvalcounter)
    else # determine size of matrix after the fact
        m, n = isempty(rowval) ? 0 : maximum(rowval), s.colcounter - 1
    end
    return SparseMatrixCSC(m, n, colptr, rowval, nzval)
end

#######################################################################
# SparseMatrixReader
#######################################################################

struct SparseMatrixReader{Tv,Ti}
    matrix::SparseMatrixCSC{Tv,Ti}
end

Base.IteratorSize(::SparseMatrixReader) = Base.HasLength()

Base.IteratorEltype(::SparseMatrixReader) = Base.HasEltype()

Base.eltype(::SparseMatrixReader{Tv,Ti}) where {Tv,Ti} = Tuple{Ti,Ti,Tv}

Base.length(s::SparseMatrixReader) = nnz(s.matrix)

function Base.iterate(s::SparseMatrixReader, state = (1, 1))
    (ptr, col) = state
    ptr > length(s) && return nothing
    @inbounds while ptr > s.matrix.colptr[col + 1] - 1
         col += 1
    end
    return (s.matrix.rowval[ptr], col, s.matrix.nzval[ptr], ptr), (ptr + 1, col)
end

enumerate_sparse(s::SparseMatrixCSC) = SparseMatrixReader(s)