text:
export total_degree, bezout_number, total_degree_start_solutions

"""
    total_degree(
        F::System;
        parameters = nothing,
        gamma = cis(2π * rand()),
        tracker_options = TrackerOptions(),
        endgame_options = EndgameOptions(),
    )

Solve the system `F` using a total degree homotopy.
This returns a path tracker ([`EndgameTracker`](@ref) or [`OverdeterminedTracker`](@ref)) and an iterator to compute the start solutions.
If the system `F` has declared `variable_groups` then a multi-homogeneous
a start system following [^Wam93] will be constructed.

[^Wam93]: An efficient start system for multi-homogeneous polynomial continuation, Wampler, C.W. Numer. Math. (1993) 66: 517. https://doi.org/10.1007/BF01385710
"""
function total_degree(
    F::Union{System,AbstractSystem};
    compile::Union{Bool,Symbol} = COMPILE_DEFAULT[],
    kwargs...,
)
    if isa(F, AbstractSystem) || isnothing(variable_groups(F))
        return total_degree_variables(F; compile = compile, kwargs...)
    else
        return total_degree_variable_groups(F; compile = compile, kwargs...)
    end
end

function total_degree_variables(
    F::Union{System,AbstractSystem};
    target_parameters = nothing,
    γ = cis(2π * rand()),
    gamma = γ,
    tracker_options = TrackerOptions(),
    endgame_options = EndgameOptions(),
    compile::Union{Bool,Symbol} = COMPILE_DEFAULT[],
    kwargs...,
)
    unsupported_kwargs(kwargs)
    m, n = size(F)

    @unique_var x[1:n] t s[1:n]
    support, coeffs = support_coefficients(System(F(x, target_parameters), x))

    homogeneous = true
    D = zeros(Int, length(support))
    for (k, A) in enumerate(support)
        d = 0
        for j = 1:size(A, 2)
            dⱼ = 0
            for i = 1:size(A, 1)
                dⱼ += A[i, j]
            end
            homogeneous = homogeneous && (j == 1 || dⱼ == d)
            d = max(d, dⱼ)
        end
        D[k] = d
    end
    scaling = maximum.(abs ∘ float, coeffs)

    m ≥ (n - homogeneous) || throw(FiniteException(n - homogeneous - m))

    overdetermined = m > n - homogeneous
    # sort by descending degree for overdetermined systems
    if overdetermined && F isa System
        F = deepcopy(F)
        perm = sortperm(D; rev = true)
        permute!(F.expressions, perm)
        permute!(D, perm)
        permute!(scaling, perm)
    end

    if F isa System
        F = fixed(F; compile = compile)
    end
    if target_parameters !== nothing
        F = fix_parameters(F, target_parameters; compile = compile)
    end
    # if homogeneous put on on an affine chart
    if homogeneous && overdetermined
        F = on_affine_chart(F)
        push!(D, 1)
        push!(scaling, 1.0)
        homogeneous = false
        F = square_up(F)
        scaling = [LA.I F.A] * scaling
        D = D[1:n]
    elseif overdetermined
        F = square_up(F)
        scaling = [LA.I F.A] * scaling
        D = D[1:n]
    end
    if homogeneous
        G = fixed(
            System(s[1:n-1] .* (x[1:n-1] .^ D[1:n-1] .- x[end] .^ D[1:n-1]), x, s[1:n-1]);
            compile = compile,
        )
    else
        G = fixed(
            System(s .* (x .^ D .- 1), x, s);
            compile = compile,
            optimizations = false,
        )
    end
    G = fix_parameters(G, scaling)

    H = StraightLineHomotopy(G, F; γ = gamma)
    if homogeneous
        H = on_affine_chart(H)
    end
    T = EndgameTracker(H, tracker_options = tracker_options, options = endgame_options)
    if overdetermined
        T = OverdeterminedTracker(T, F)
    end
    starts = total_degree_start_solutions(D; homogeneous = homogeneous)

    T, starts
end

function total_degree_variable_groups(
    F::System;
    γ = cis(2π * rand()),
    gamma = γ,
    target_parameters = nothing,
    tracker_options = TrackerOptions(),
    endgame_options = EndgameOptions(),
    compile::Union{Bool,Symbol} = COMPILE_DEFAULT[],
    kwargs...,
)
    unsupported_kwargs(kwargs)
    m, n = size(F)
    homogeneous = is_homogeneous(F)
    D = multi_degrees(F)
    vargroups = variable_groups(F)
    projective_dims = length.(vargroups) .- homogeneous
    M = length(projective_dims)

    m ≥ (n - M * homogeneous) || throw(FiniteException(n - M * homogeneous - m))

    overdetermined = m > n - M * homogeneous

    F = fixed(F; compile = compile)
    if target_parameters !== nothing
        F = fix_parameters(F, target_parameters; compile = compile)
    end

    if homogeneous && overdetermined
        F = on_affine_chart(F, projective_dims)
        D = [D LA.I]
        homogeneous = false
        F = square_up(F)
        D = max.(D[:, 1:n], maximum(D[:, n+1:end], dims = 2))
        projective_dims .+= 1
    elseif overdetermined
        F = square_up(F)
        D = max.(D[:, 1:n], maximum(D[:, n+1:end], dims = 2))
    end

    g, C = multi_homogeneous_system(D, vargroups; homogeneous = homogeneous)
    P = parameters(g)
    p = randn(length(P))
    G = fix_parameters(fixed(g; compile = compile), p)
    starts = MultiBezoutSolutionsIterator(
        D,
        C(P => p),
        projective_dims;
        homogeneous = homogeneous,
    )
    H = StraightLineHomotopy(G, F; gamma = gamma)
    if homogeneous
        H = on_affine_chart(H, projective_dims)
    end
    T = EndgameTracker(H, tracker_options = tracker_options, options = endgame_options)
    if overdetermined
        T = OverdeterminedTracker(T, F)
    end
    T, starts
end

function multi_homogeneous_system(D, vargroups; homogeneous::Bool)
    M, n = size(D)
    projective_dims = length.(vargroups) .- homogeneous
    @unique_var c[1:n+M, 1:n]
    C = zeros(Expression, n + M, n)
    for i = 1:n
        j = 1
        for kⱼ in projective_dims
            for l = 1:kⱼ
                C[j, i] = begin
                    if (i == l && i ≤ kⱼ)
                        1
                    elseif i ≠ l && i ≤ kⱼ
                        0
                    else
                        c[j, i]
                    end
                end
                j += 1
            end
            C[j, i] = 1
            j += 1
        end
    end

    P = variables(C)

    if homogeneous
        Z = vargroups
    else
        Z = vcat.(vargroups, 1)
    end

    m, n = size(D)
    g = map(1:n) do i
        s = Ref(1)
        prod(1:m) do j
            dᵢⱼ = D[j, i]
            kⱼ = length(Z[j]) - 1
            if dᵢⱼ == 0
                s[] += kⱼ + 1
                return 1
            end
            bᵢⱼ = sum(1:kⱼ) do l
                cz = C[s[], i] * Z[j][l]
                s[] += 1
                cz
            end
            s[] += 1
            bᵢⱼ^dᵢⱼ - Z[j][end]^dᵢⱼ
        end
    end
    System(g; variable_groups = vargroups, parameters = P), C
end

struct TotalDegreeStartSolutionsIterator{Iter}
    degrees::Vector{Int}
    homogeneous::Bool
    iterator::Iter
end
function TotalDegreeStartSolutionsIterator(degrees, homogeneous::Bool)
    iterator = Iterators.product(map(d -> 0:(d-1), degrees)...)
    TotalDegreeStartSolutionsIterator(degrees, homogeneous, iterator)
end
function Base.show(io::IO, iter::TotalDegreeStartSolutionsIterator)
    print(io, "$(length(iter)) total degree start solutions for degrees $(iter.degrees)")
end

function Base.iterate(iter::TotalDegreeStartSolutionsIterator)
    indices, state = iterate(iter.iterator)
    _value(iter, indices), state
end
function Base.iterate(iter::TotalDegreeStartSolutionsIterator, state)
    it = iterate(iter.iterator, state)
    it === nothing && return nothing
    _value(iter, first(it)), last(it)
end

function _value(iter::TotalDegreeStartSolutionsIterator, indices)
    x = map((k, d) -> cis(2π * k / d), indices, iter.degrees)
    iter.homogeneous && push!(x, 1)
    x
end

Base.length(iter::TotalDegreeStartSolutionsIterator) = length(iter.iterator)
Base.eltype(iter::Type{<:TotalDegreeStartSolutionsIterator}) = Vector{Complex{Float64}}

"""
    total_degree_start_solutions(degrees)

Returns an iterator of the start solutions of the system
```math
\\begin{array}{rl}
    z_1^{d_1} - 1 &= 0 \\\\
    z_2^{d_2} - 1 &= 0 \\\\
    &\\vdots \\\\
    z_n^{d_n} - 1 &= 0 \\\\
\\end{array}
```
where ``d_i`` is `degrees[i]`.

## Example
```julia-repl
julia> iter = total_degree_start_solutions([2, 2])
4 total degree start solutions for degrees [2, 2]

julia> collect(iter)
4-element Array{Array{Complex{Float64},1},1}:
 [1.0 + 0.0im, 1.0 + 0.0im]
 [-1.0 + 1.2246467991473532e-16im, 1.0 + 0.0im]
 [1.0 + 0.0im, -1.0 + 1.2246467991473532e-16im]
 [-1.0 + 1.2246467991473532e-16im, -1.0 + 1.2246467991473532e-16im]
```
"""
total_degree_start_solutions(
    degrees::AbstractVector{<:Integer};
    homogeneous::Bool = false,
) = TotalDegreeStartSolutionsIterator(degrees, homogeneous)

###
### Multihomogeneous
###

struct MultiBezoutIndicesIterator
    D::Matrix{Int} # Degree matrix
    k::Vector{Int} # (projective) dimensions

    function MultiBezoutIndicesIterator(D::Matrix{Int}, k::Vector{Int})
        @assert size(D, 1) == length(k)
        new(D, k)
    end
end
Base.IteratorSize(::Type{MultiBezoutIndicesIterator}) = Base.SizeUnknown()
Base.eltype(::Type{MultiBezoutIndicesIterator}) = Tuple{Int,Vector{Int}}
function Base.iterate(
    iter::MultiBezoutIndicesIterator,
    state = [i for i = 1:length(iter.k) for j = 1:iter.k[i]],
)
    d = 0
    m, n = size(iter.D)
    p = state
    while d == 0
        state[1] > n && return nothing
        p, state = nextpermutation(1:m, state)
        let p = p, D = iter.D
            d = prod(i -> D[p[i], i], 1:n)::Int
        end
    end

    (d, p), state
end

# Adopted from Combinatorics.jl:
# https://github.com/JuliaMath/Combinatorics.jl/blob/8d9571402319799b29da2005a65b627e8771c1e4/src/permutations.jl#L47
function nextpermutation(m, state)
    n = length(state)
    perm = [m[state[i]] for i = 1:n]
    s = copy(state)
    i = n - 1
    while i >= 1 && s[i] >= s[i+1]
        i -= 1
    end
    if i > 0
        j = n
        while j > i && s[i] >= s[j]
            j -= 1
        end
        s[i], s[j] = s[j], s[i]
        reverse!(s, i + 1)
    else
        s[1] = n + 1
    end
    return (perm, s)
end

struct MultiBezoutSolutionsIterator
    indices::MultiBezoutIndicesIterator
    coeffs::Matrix{Float64}
    homogeneous::Bool
    # precomputed stuff and cache
    roots_of_unity::Matrix{ComplexF64} # Lower Triangular matrix storing the precomputed roots
    A::Vector{Matrix{Float64}}
    b::Vector{Vector{ComplexF64}}
    ranges::Vector{UnitRange{Int}}
end

function MultiBezoutSolutionsIterator(
    indices::MultiBezoutIndicesIterator,
    coeffs::Matrix;
    homogeneous::Bool,
)
    d_max = maximum(indices.D)
    roots_of_unity = zeros(ComplexF64, d_max, d_max)
    for i = 1:d_max, j = 0:i-1
        roots_of_unity[i, j+1] = cis(2π * j / i)
    end
    A = map(kⱼ -> zeros(kⱼ, kⱼ), indices.k)
    b = map(kⱼ -> zeros(ComplexF64, kⱼ), indices.k)
    r_stop = Ref(1)
    ranges = map(indices.k) do kⱼ
        range = r_stop[]:(kⱼ+r_stop[]-1)
        r_stop[] += kⱼ + 1
        range
    end
    MultiBezoutSolutionsIterator(indices, coeffs, homogeneous, roots_of_unity, A, b, ranges)
end
function MultiBezoutSolutionsIterator(
    D::Matrix,
    C::Matrix,
    proj_dims::Vector{Int};
    kwargs...,
)
    MultiBezoutSolutionsIterator(MultiBezoutIndicesIterator(D, proj_dims), C; kwargs...)
end

function Base.show(io::IO, iter::MultiBezoutSolutionsIterator)
    print(io, "Solutions iterator for a multi-homogeneous start system")
end
Base.show(io::IO, ::MIME"application/prs.juno.inline", x::MultiBezoutSolutionsIterator) = x
Base.IteratorSize(::Type{<:MultiBezoutSolutionsIterator}) = Base.SizeUnknown()
Base.eltype(::Type{MultiBezoutSolutionsIterator}) = Vector{ComplexF64}

function Base.iterate(iter::MultiBezoutSolutionsIterator)
    (_, perm), indices_state = iterate(iter.indices)
    dᵢ = ntuple(i -> iter.indices.D[perm[i], i], length(perm))
    d_iter = Iterators.product(map(dᵢ -> 1:dᵢ, dᵢ)...)
    q, d_state = iterate(d_iter)

    x = compute_solution(iter, perm, q, dᵢ)
    x, (perm, indices_state, dᵢ, d_iter, d_state)
end

function Base.iterate(iter::MultiBezoutSolutionsIterator, state)
    perm, indices_state, dᵢ, d_iter, d_state = state

    next_d_state = iterate(d_iter, d_state)
    if next_d_state !== nothing
        q, d_state = next_d_state
        x = compute_solution(iter, perm, q, dᵢ)
        return x, (perm, indices_state, dᵢ, d_iter, d_state)
    end

    # We exhausted our current batch. move indices forward
    indices_newstate = iterate(iter.indices, indices_state)
    # check whether we are completely done
    indices_newstate === nothing && return nothing

    (_, perm), indices_state = indices_newstate
    # Create new d_iter
    new_dᵢ::typeof(dᵢ) = let D = iter.indices.D, perm = perm
        ntuple(i -> D[perm[i], i], length(perm))
    end
    new_d_iter::typeof(d_iter) = Iterators.product(map(dᵢⱼ -> 1:dᵢⱼ, new_dᵢ)...)
    q, d_state = iterate(new_d_iter)
    x = compute_solution(iter, perm, q, new_dᵢ)
    return x, (perm, indices_state, new_dᵢ, new_d_iter, d_state)
end

function compute_solution(iter::MultiBezoutSolutionsIterator, perm, q, dᵢ)
    m, n = size(iter.indices.D)
    # for each variable group we have to solve a linear system
    t = 1
    if iter.homogeneous
        data = zeros(ComplexF64, n + m)
    else
        data = zeros(ComplexF64, n)
    end
    for j = 1:m
        kⱼ = iter.indices.k[j]
        Aⱼ = iter.A[j]
        bⱼ = iter.b[j]
        s = 1
        for i = 1:n
            perm[i] == j || continue
            range = iter.ranges[j]
            for (t, l) in enumerate(range)
                Aⱼ[s, t] = iter.coeffs[l, i]
            end
            bⱼ[s] = iter.roots_of_unity[dᵢ[i], q[i]]
            s += 1
        end

        LinearAlgebra.ldiv!(LinearAlgebra.lu!(Aⱼ), bⱼ)

        data[t:(t+kⱼ-1)] .= bⱼ
        if iter.homogeneous
            data[t+kⱼ] = 1
            t += kⱼ + 1
        else
            t += kⱼ
        end
    end

    data
end

function paths_to_track(f, ::Val{:total_degree})
    target_parameters = nparameters(f) > 0 ? zeros(nparameters(f)) : nothing
    _, starts = total_degree(f; target_parameters = target_parameters)
    if Base.IteratorSize(typeof(starts)) == Base.SizeUnknown()
        k = 0
        for s in starts
            k += 1
        end
        k
    else
        length(starts)
    end
end
@deprecate bezout_number(f::Union{System,AbstractSystem}) paths_to_track(
    f;
    start_system = :total_degree,
)


text:
# using Base.Test
# using MIToS.MSA

print("""

Tests for Raw Multiple Sequence Alignment a.k.a Matrix{Residue}
================================================================
""")

print("""

Parse Raw
---------
""")

let raw = read(joinpath(pwd(), "data", "gaps.txt"), Raw),
  raw_string = """THAYQAIHQV
       THAYQAIHQ-
       THAYQAIH--
       THAYQAI---
       THAYQA----
       THAYQ-----
       THAY------
       THA-------
       TH--------
       T---------
       """

  @test asciisequence(raw, 1)  == "THAYQAIHQV"
  @test asciisequence(raw, 10) == "T---------"

  @test parse(raw_string, Raw) == raw

  print("""

Print Raw
---------
""")

  let io = IOBuffer()
    print(io, raw, Raw)
    @test takebuf_string(io) == raw_string
  end

  print("""

Test %
------
""")

  @test gapfraction(raw) == 0.45
  @test gapfraction(raw, 1) == [0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9]
  @test gapfraction(raw, 2) == [0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9]

  @test residuefraction(raw) == 0.55
  @test residuefraction(raw, 1) == [1.0,0.9,0.8,0.7,0.6,0.5,0.4,0.3,0.2,0.1]
  @test residuefraction(raw, 2) == [1.0,0.9,0.8,0.7,0.6,0.5,0.4,0.3,0.2,0.1]

  @test coverage(raw) == [1.0,0.9,0.8,0.7,0.6,0.5,0.4,0.3,0.2,0.1]

  @test columngapfraction(raw) == [0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9]

  print("""

Reference and Gapstrip
----------------------
""")

  let gs = gapstrip(raw,  coveragelimit=0.5, gaplimit=0.5)
    @test getsequence(gs, 1) == res"THAYQAIH"
    @test ncolumns(gs) == 8
    @test nsequences(gs) == 6
  end

  let ref = setreference!(copy(raw), 2)
    @test getsequence(ref, 1) == res"THAYQAIHQ-"
    ref = adjustreference(ref)
    @test getsequence(ref, 1) == res"THAYQAIHQ"
  end

end

ext: jl


text:
# result reference counter
# -------------------------
# Used to keep track if a result should be kept around because it is needed by another node in the DAG.
# If not, then results can be cleaned to free up resources.

"""
Simple reference counting done by the scheduler master
"""
module SimpleRC

import ..DagScheduler
import ..DagScheduler: TaskIdType

const REFCOUNT = Dict{TaskIdType,Int}()

refcount(id::TaskIdType, cnt::Int) = remotecall_fetch(_refcount, 1, id, cnt)
clean_refcount() = remotecall_fetch(_clean_refcount, 1)
nz_refcounts() = remotecall_fetch(_nz_refcounts, 1)

function _refcount(id::TaskIdType, cnt::Int)
    ref = REFCOUNT[id] = (get!(REFCOUNT, id, 0) + cnt)
    if ref == 0
        delete!(REFCOUNT, id)
    end
    ref
end

function _clean_refcount()
    empty!(REFCOUNT)
    nothing
end

function _nz_refcounts()
    filter((n,v)->(v > 0), REFCOUNT)
end

end # SimpleRC

"""
No reference counting done by the scheduler
"""
module NoRC

import ..DagScheduler
import ..DagScheduler: TaskIdType

refcount(id::TaskIdType, cnt::Int) = 2
clean_refcount() = nothing
nz_refcounts() = Dict{TaskIdType,Int}()

end # NoRC

"""
Maintain reference counts on foundation db
"""
module FdbRC

import ..DagScheduler
import ..DagScheduler: TaskIdType
using FoundationDB

const RefCountRoot = UInt8[0]
const fdb_db = Ref{Union{Void,FDBDatabase}}(nothing)

function set_fdb_db()
    fdb_db[] = DagScheduler.FdbMeta.fdb_db[]
    nothing
end

function refcount(id::TaskIdType, cnt::Int)
    key = vcat(RefCountRoot, reinterpret(UInt8, [id]))
    open(FDBTransaction(fdb_db[])) do tran
        atomic_add(tran, key, cnt)
    end
    ref = open(FDBTransaction(fdb_db[])) do tran
        ref = atomic_integer(Int, getval(tran, key))
        if ref == 0
            clearkey(tran, key)
        end
        ref
    end
    ref
end

function clean_refcount()
    start_key = vcat(RefCountRoot, UInt8[0])
    end_key = vcat(RefCountRoot, UInt8[0xff])
    open(FDBTransaction(fdb_db[])) do tran
        clearkeyrange(tran, start_key, end_key)
    end
    nothing
end

function nz_refcounts()
    start_key = vcat(RefCountRoot, UInt8[0])
    end_key = vcat(RefCountRoot, UInt8[0xff])
    kvs = open(FDBTransaction(fdb_db[])) do tran
        kvs, more = getrange(tran, keysel(FDBKeySel.first_greater_than, start_key), keysel(FDBKeySel.first_greater_than, end_key))
        kvs
    end
    res = Dict{TaskIdType,Int}()
    if kvs !== nothing
        for kv in kvs
            key,val = kv
            id = reinterpret(TaskIdType, key[(end-sizeof(TaskIdType)+1):end])[1]
            cnt = atomic_integer(Int, val)
            res[id] = cnt
        end
    end
    res
end

end # FdbRC

const RefCounter = Ref{Module}(NoRC)


text:
module SymbolicMath

import Base: +, -, *, /, ^, cos, sin, log, exp, abs, real, imag, convert
import Base.==

abstract type AbstractTerm end
struct ExprTerm <: AbstractTerm
    ex
end
convert(::Type{ExprTerm}, ex) = ExprTerm(ex)

macro ExprTerm(ex)
    ExprTerm(ex)
end

function createsymbol(s::Symbol, m::Module)
    @eval m $s = $(Meta.quot(s)) |> ExprTerm
end
function createsymbol(s::Symbol)
    createsymbol(s, Main)
end

macro symbols(args::Symbol...)
    createsymbol.(args, Ref(__module__))
end
macro symbol(s::Symbol)
    createsymbol(s, __module__)
end

# euqality comparison
==(t1::ExprTerm, t2::ExprTerm) = (t1.ex == t2.ex)
==(t::ExprTerm, other) = (t.ex == other)
==(other, t::ExprTerm) = (t == other)

# operator definitions
binary_opeators = [:+, :-, :*, :/, :^]
unary_operators = [:cos, :sin, :log, :exp, :abs, :real, :imag]

for op in binary_opeators
    @eval $op(t1::ExprTerm, t2::ExprTerm) = Expr(:call, $(Meta.quot(op)), t1.ex, t2.ex) |> ExprTerm
    @eval $op(t::ExprTerm, other) = Expr(:call, $(Meta.quot(op)), t.ex, other) |> ExprTerm
    @eval $op(other, t::ExprTerm) = Expr(:call, $(Meta.quot(op)), other, t.ex) |> ExprTerm
end
for op in unary_operators
    @eval $op(t::ExprTerm) = Expr(:call, $(Meta.quot(op)), t.ex) |> ExprTerm
end

export @symbol, @symbols, @ExprTerm, ExprTerm

end # module SymbolicMath

max_stars_repo_name: timkittel/SymbolicMath


text:
using Gadfly
using Cairo
using Fontconfig

type Coordinate
    x :: Float64
    y :: Float64
end

function input_coordinates(filename::String)
    coords = Array(Coordinate, 0)
    open(filename, "r") do file
        for line in eachline(file)
            coord_ary = split(strip(line, ['\n', ',']), ' ')
            coord = Coordinate(float(coord_ary[1]), float(coord_ary[2]))
            push!(coords, coord)
        end
    end

    return coords
end

function calc_distance{N <: Integer}(path::Array{N}, coords::Array{Coordinate})
    distance = 0
    for i in 1:(length(coords) - 1)
        present_idx, next_idx = path[i], path[i + 1]
        distance += sqrt((coords[next_idx].x - coords[present_idx].x)^2
                         + (coords[next_idx].y - coords[present_idx].y)^2)
    end

    return distance
end

function random_path{N <: Integer}(n::N)
    path = collect(1:(n - 1))
    path = shuffle(path)
    push!(path, path[1])
    return path
end

function random_path(coords::Array{Coordinate})
    n = length(coords)
    random_path(n)
end

function metropolis{N <: Integer, R <: Real}(current_path::Array{N},
                                             coords::Array{Coordinate},
                                             temp::R)
    current_distance = calc_distance(current_path, coords)
    candidate_path = copy(current_path)
    n = length(current_path)

    swap_cities_indices = Array(Any, 2)
    while true
        swap_cities_indices = rand(2:(length(current_path) - 1), 2)
        if swap_cities_indices[1] != swap_cities_indices[2]
            break
        end
    end
    candidate_path[swap_cities_indices[1]], candidate_path[swap_cities_indices[2]] =
        candidate_path[swap_cities_indices[2]], candidate_path[swap_cities_indices[1]]
    candidate_distance = calc_distance(candidate_path, coords)

    Δd = candidate_distance - current_distance
    if Δd < 0
        return candidate_path, candidate_distance
    end

    prob = exp(-Δd / temp)
    if rand() <= prob
        return candidate_path, candidate_distance
    else
        return current_path, current_distance
    end
end

function anneal{N <: Integer}(path::Array{N}, coords::Array{Coordinate}; n_iter::N=100000,
                              pmelt=0.7, tgt=0.1, stagfactor=0.05, procplt::Bool=false)
    n_cities = length(path)
    initial_distance = calc_distance(path, coords)
    min_distance, max_distance = initial_distance, initial_distance

    optimized_distances = Array(Float64, 0)
    distances = Array(Float64, 0)
    push!(optimized_distances, initial_distance)
    push!(distances, initial_distance)
    for i in 1:(max(0.01 * n_cities, 2))
        path_, disntance = metropolis(path, coords, typemax(Int64))
        if disntance < min_distance
            min_distance = disntance
        end
        if max_distance < disntance
            max_distance = disntance
        end
    end

    range = (max_distance - min_distance) * pmelt
    temp = tgt ^ (1 / n_iter)
    optimized_distance = initial_distance
    optimized_step = 1
    optimized_path = copy(path)
    path_ = copy(path)

    for i in 1:n_iter
        print("$(i) / $(n_iter) processing...\r")
        flush(STDOUT)

        dtemp = range * (temp ^ i)
        path_, distance = metropolis(path_, coords, dtemp)

        if distance < optimized_distance
            optimized_distance = distance
            optimized_path = copy(path_)
            optimized_step = i
        end
        push!(optimized_distances, optimized_distance)
        push!(distances, distance)

        # Reheat
        if i - optimized_step == stagfactor * n_iter
            temp = temp ^ (0.05 * i / n_iter)
        end
    end
    print("\n")

    if procplt
        white_panel = Theme(background_color="white")
        p = plot(layer(x=1:length(distances), y=distances, Geom.line,
                       Theme(default_color="gray")),
                 layer(x=1:length(distances), y=optimized_distances, Geom.line,
                       Theme(default_color="orange", line_width=2px)),
                 white_panel,
                 Coord.Cartesian(xmin=0, xmax=length(distances)))
        img = PNG("result.png", 1280px, 720px)
        draw(img, p)
    end

    return optimized_path
end

coords = input_coordinates("prefs.in")
init_path = random_path(coords)
path = anneal(init_path, coords, n_iter=100000, procplt=true)


text:
include("code_needed.jl")
myedge_list,names_data = run_script()
# reformat data to smat
writeToSMAT(myedge_list,"batter_pitcher_edges.smat")
writedlm("batter_pitcher_names.labels",names_data)

max_stars_repo_name: empet/more-graphs

max_issues_repo_name: empet/more-graphs

max_forks_repo_name: empet/more-graphs


text:
using NURBS
using Base.Test

@testset "Order Increasing" begin

    @testset "From order 2 to 3" begin
        x = knot(4,2)
        b = [1. 2 3 4;1 2 2 1;1 1 1 1]
        @test typeof(raise23(b,x,4)) == Tuple{Array{Float64,2},Int64}
        @test raise23(b,x,4)[1] == [1 1.5 2 2.5 3 3.5 4; 1 1.5 2 2 2 1.5 1; 1 1 1 1 1 1 1]
        @test raise23(b,x,4)[2] == 7
        x = [0, 0., 1, 1, 2, 2]
        @test typeof(raise23(b,x,4)) == Tuple{Array{Float64,2},Int64}
        @test raise23(b,x,4)[1] == [1 1.5 2 3 3.5 4; 1 1.5 2 2 1.5 1; 1 1 1 1 1 1]
        @test raise23(b,x,4)[2] == 6
        x = knotu(4,2)
        @test typeof(raise23(b,x,4)) == Tuple{Array{Float64,2},Int64}
        @test raise23(b,x,4)[1] == [0.5 1 1.5 2 2.5 3 3.5 4 2;0.5 1 1.5 2 2 2 1.5 1 0.5;0.5 1 1 1 1 1 1 1 0.5]
        @test raise23(b,x,4)[2] == 9
    end

    @testset "From order 4 to 5" begin
        x = knot(4,4)
        b = [1. 2 3 4 ; 1 2 2 1 ; 1 1 1 1]
        @test typeof(raise45(b,x,4)) == Tuple{Array{Float64,2},Int64}
        @test raise45(b,x,4)[1] == [1 1.75 2.5 3.25 4; 1 1.75 2 1.75 1; 1 1 1 1 1]
        @test raise45(b,x,4)[2] == 5
        x = knot(5,4)
        b = [1. 2 3 4 5; 1 2 2 1 0; 1 1 1 1 1]
        @test typeof(raise45(b,x,5)) == Tuple{Array{Float64,2},Int64}
        @test raise45(b,x,5)[1] == [1 1.75 2.25 3 3.75 4.25 5; 1 1.75 2 1.875 1.25 0.75 0.0; 1 1 1 1 1 1 1]
        @test raise45(b,x,5)[2] == 7
        x = knot(6,4)
        b = [1. 2 3 4 5 6; 1 2 2 1 0 0; 1 1 1 1 1 1]
        @test typeof(raise45(b,x,6)) == Tuple{Array{Float64,2},Int64}
        @test isapprox(raise45(b,x,6)[1], [1 1.75 2.25 2.95833 3.5 4.04167 4.75 5.25 6; 1 1.75 2 1.91667 1.5 0.958333 0.25 0.0 0.0; 1 1 1 1 1 1 1 1 1], atol=5e-4)
        @test raise45(b,x,6)[2] == 9
        x = knot(7,4)
        b = [1. 2 3 4 5 6 7; 1 2 2 1 0 0 1; 1 1 1 1 1 1 1]
        @test typeof(raise45(b,x,7)) == Tuple{Array{Float64,2},Int64}
        @test isapprox(raise45(b,x,7)[1], [1 1.75 2.25 2.95833 3.5 4 4.5 5.04167 5.75 6.25 7; 1 1.75 2 1.91667 1.5 1 0.5 0.08333 0 0.25 1; 1 1 1 1 1 1 1 1 1 1 1]; atol=1e-5)
        @test raise45(b,x,7)[2] == 11
    end
end


text:
####################################
## Demo script for Bayesian neural network.

using PyPlot, PyCall

push!(LOAD_PATH, ".")
import BayesNeuralNet

"""
Sample neural nets from prior.
"""
function sample_test()
    # model parameters
    D = 1 # output
    K = 3 # hidden
    M = 2 # input
    sigma2_w = 10.0
    sigma2_y = 0.1

    xmin = -5
    xmax = 5
    N_lin = 1000
    X_lin = ones(M, N_lin)
    X_lin[1,:] = linspace(xmin, xmax, N_lin)
    X_lin[2,:] = 1 # bias

    # visualize
    num_samples = 5
    figure("Function samples")
    clf()
    for i in 1 : num_samples
        _, Y_true, _, _ = BayesNeuralNet.sample_data_from_prior(X_lin, sigma2_w, sigma2_y, D, K)
        plot(X_lin[1,:], Y_true)
        xlim([xmin, xmax])
    end
    ratey = (ylim()[2] - ylim()[1]) * 0.1
    ratex = (xlim()[2] - xlim()[1]) * 0.1
    text(xlim()[1] + ratex, ylim()[2] - ratey, @sprintf("K=%d", K), fontsize=18)
    show()
end

"""
Run a test script of variational inference for Bayesian neural net.
"""
function test()
    #################
    # prepara data

    # data size
    D = 1 # output
    M = 2 # input

    # function setting
    xmin = -2
    xmax = 4
    N_lin = 1000
    X_lin = ones(M, N_lin)
    X_lin[1,:] = linspace(xmin, xmax, N_lin)
    X_lin[2,:] = 1 # bias

    # training data
    N = 50 # data size
    X = 2*rand(M, N) - 0.0 # input
    X[2,:] = 1.0 # bias
    Y = 0.5*sin.(2*pi * X[1,:]/3) + 0.05 * randn(N)

    # model parameters
    K = 5
    sigma2_w = 10.0
    sigma2_y = 0.01

    ################
    # inference
    alpha = 1.0e-5
    max_iter = 100000
    mu1, rho1, mu2, rho2 = BayesNeuralNet.VI(Y, X, sigma2_w, sigma2_y, K, alpha, max_iter)
    Y_mean = [mu2'* tanh.(mu1'X_lin[:,n]) for n in 1 : N_lin]

    ################
    # visualize
    figure("result")
    clf()
    Y_list = []
    num_samples = 100
    for i in 1 : num_samples
        Y_est, _ = BayesNeuralNet.sample_data_from_posterior(X_lin, mu1, rho1, mu2, rho2, sigma2_y, D)
        push!(Y_list, Y_est)
        plot(X_lin[1,:], Y_est, "-c", alpha=0.25)
    end
    plot(X[1,:], Y, "ok")
    plot(X_lin[1,:], Y_mean, "b-")
    xlim([xmin, xmax])
    xlabel("x")
    ylabel("y")
    show()
end

#sample_test()
test()


text:
X = randn(100,10)
truebeta = rand(10)
y = X*truebeta+randn(100)

function LinAlg.linreg(x::Array{Float64,2}, y::Array{Float64,1})
    n, p = size(x)
    beta = [ones(n) x]\y
    return Dict(:intercept => beta[1], :beta => beta[2:(p+1)])
end

out = linreg(X, y)

out[:intercept]
out[:beta]

using Plots
plot(zeros(10), label = "", ylims = (-1, 1), size=(1000,500))
plot!((truebeta-out[:beta]), label = "")
scatter!((truebeta-out[:beta]), label = "")
title!("Bias of the beta parameters")
savefig("linreg")

# Also does ridge regression

X2 = randn(100, 200)
truebeta = rand(200) - 0.5
y2 = X2*truebeta+randn(100)
out2 = linreg(X2, y2)

plot(zeros(200), label = "", ylims = (-1, 1), size=(1000,500))
plot!((truebeta-out2[:beta]), label = "")
scatter!((truebeta-out2[:beta]), label = "")
title!("Bias of the beta parameters")
savefig("linregridge")

max_stars_repo_name: vankesteren/julia-tryout

max_issues_repo_name: vankesteren/julia-tryout

max_forks_repo_name: vankesteren/julia-tryout


text:
# ---
# title: 542. 01 Matrix
# id: problem542
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Depth-first Search, Breadth-first Search
# link: <https://leetcode.com/problems/01-matrix/description/>
# hidden: true
# ---
#
# Given a matrix consists of 0 and 1, find the distance of the nearest 0 for
# each cell.
#
# The distance between two adjacent cells is 1.
#
#
#
# **Example 1:**
#
#
#
#     Input:
#     [[0,0,0],
#      [0,1,0],
#      [0,0,0]]
#
#     Output:
#     [[0,0,0],
#      [0,1,0],
#      [0,0,0]]
#
#
# **Example 2:**
#
#
#
#     Input:
#     [[0,0,0],
#      [0,1,0],
#      [1,1,1]]
#
#     Output:
#     [[0,0,0],
#      [0,1,0],
#      [1,2,1]]
#
#
#
#
# **Note:**
#
#   1. The number of elements of the given matrix will not exceed 10,000.
#   2. There are at least one 0 in the given matrix.
#   3. The cells are adjacent in only four directions: up, down, left and right.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
"""
$(TYPEDEF)

Tokenstream allows to read tokenized data from file without keeping
the file ocntent in memory.

$(TYPEDFIELDS)
"""
mutable struct TokenStream
    """
    Input stream
    """
    input::IOStream

    """
    Array of current tokens kept in memory.
    """
    tokens::Array{SubString{String},1}

    """
    Position of actual token in tokens array
    """
    itoken::Int

    """
    Line number in IOStream
    """
    lineno::Int

    """
    Comment character
    """
    comment::Char

    """
    Function telling if given character is a delimiter.
    """
    dlm::Function
    TokenStream()=new()
end

"""
$(TYPEDSIGNATURES)

Tokenstream destructor should close input
"""
destruct!(tks::TokenStream)=close(tks.input)

"""
$(TYPEDEF)

Error thrown when the token expected  in expect!  is not there.

$(TYPEDFIELDS)
"""
struct UnexpectedTokenError <: Exception
    found::String
    expected::String
    lineno::Int
end

Base.showerror(io::IO, e::UnexpectedTokenError) = print(io,"Unexpected token in line $(e.lineno): $(e.found) (expected $(e.expected))")

"""
$(TYPEDSIGNATURES)

Create Tokenstream with file name argument.
"""
TokenStream(filename::String;comment='#', dlm=isspace)=TokenStream(open(filename),comment=comment,dlm=dlm)

"""
$(TYPEDSIGNATURES)

Create Tokenstream with IOStream argument.
"""
function TokenStream(input::IOStream; comment='#', dlm=isspace)
    tks=TokenStream()
    tks.input=input
    tks.lineno=0
    tks.comment=comment
    tks.dlm=dlm
    tks.tokens=split("")
    tks.itoken=1
    _fetch!(tks)
    return tks
end

"""
$(TYPEDSIGNATURES)

Check if all tokens have been consumed.
"""
Base.eof(tks::TokenStream)=!_fetch!(tks)

#=
Read next line with tokens and split it.
Return if tokens are left
=#
function _fetch!(tks::TokenStream)
    if (tks.itoken<=length(tks.tokens))
        return true
    end
    hasline=false
    buffer=""
    while !hasline && !eof(tks.input)
        buffer=readline(tks.input)
        tks.lineno+=1
        if (length(buffer)>0) && (buffer[1]!=tks.comment)
            hasline=true
        end
    end
    # TODO: replace this by something non-allocating
    tks.tokens=split(buffer,tks.dlm,keepempty=false)
    tks.itoken=1
    return (tks.itoken<=length(tks.tokens))
end

"""
$(TYPEDSIGNATURES)

Get next token from tokenstream.
"""
function  gettoken(tks::TokenStream)
    if _fetch!(tks)
        token=tks.tokens[tks.itoken]
        tks.itoken+=1
        return token
    end
    nothing
end

"""
$(TYPEDSIGNATURES)

Expect keyword token.

If token is missing, an UnexpectedTokenError is thrown
If the token  has been found, reading will continue  at the position
after the token found.
"""
function expecttoken(tks::TokenStream,expected::String)
    token=gettoken(tks)
    if token!=expected
        throw(UnexpectedTokenError(token,expected,tks.lineno))
    end
    true
end

"""
$(TYPEDSIGNATURES)

Try for keyword token.

It token is missing, the token read is put back into stream,
a value of false is returned and the next try/gettoken command continues at the same position,

Otherwise, true is returned, and reading continues after the token found.
"""
function trytoken(tks::TokenStream,expected::String)
    token=gettoken(tks)
    if token!=expected
        tks.itoken=tks.itoken-1
        return false
    end
    return true
end


text:
@block SimonDanisch ["2d"] begin

    @cell "Test heatmap + image overlap" [image, heatmap, transparency] begin
        heatmap(rand(32, 32))
        image!(map(x->RGBAf0(x,0.5, 0.5, 0.8), rand(32,32)))
    end

    @cell "Interaction" [scatter, linesegment, record] begin
        scene = Scene()
        f(t, v, s) = (sin(v + t) * s, cos(v + t) * s)
        time_node = Node(0.0)
        p1 = scatter!(scene, lift(t-> f.(t, range(0, stop = 2pi, length = 50), 1), time_node))[end]
        p2 = scatter!(scene, lift(t-> f.(t * 2.0, range(0, stop = 2pi, length = 50), 1.5), time_node))[end]
        points = lift(p1[1], p2[1]) do pos1, pos2
            map((a, b)-> (a, b), pos1, pos2)
        end
        linesegments!(scene, points)
        N = 150
        record(scene, @outputfile(mp4), range(0, stop = 10, length = N)) do i
            push!(time_node, i)
        end
    end
    @cell "barplot" [barplot] begin
        # barplot(1:10, rand(10))
        # barplot(rand(10))
        barplot(rand(10), color = rand(10))
        # barplot(rand(3), color = [:red, :blue, :green])
    end
    @cell "quiver" [quiver, arrows, vectorfield, gradient] begin
        using ImageFiltering
        x = range(-2, stop = 2, length = 21)
        y = x
        z = x .* exp.(-x .^ 2 .- (y') .^ 2)
        scene = contour(x, y, z, levels = 10, linewidth = 3)
        u, v = ImageFiltering.imgradients(z, KernelFactors.ando3)
        arrows!(x, y, u, v, arrowsize = 0.05)
    end
    @cell "image" [image] begin
        AbstractPlotting.hbox(
            image(Makie.logo(), scale_plot = false),
            image(rand(100, 500), scale_plot = false),
        )
    end
    @cell "scatter colormap" [scatter, colormap] begin
        scatter(rand(10), rand(10), color = rand(10))
    end
    @cell "FEM polygon 2D" [fem, poly] begin
        coordinates = [
            0.0 0.0;
            0.5 0.0;
            1.0 0.0;
            0.0 0.5;
            0.5 0.5;
            1.0 0.5;
            0.0 1.0;
            0.5 1.0;
            1.0 1.0;
        ]
        connectivity = [
            1 2 5;
            1 4 5;
            2 3 6;
            2 5 6;
            4 5 8;
            4 7 8;
            5 6 9;
            5 8 9;
        ]
        color = [0.0, 0.0, 0.0, 0.0, -0.375, 0.0, 0.0, 0.0, 0.0]
        poly(coordinates, connectivity, color = color, strokecolor = (:black, 0.6), strokewidth = 4)
    end
    @cell "FEM mesh 2D" [fem, mesh] begin
        coordinates = [
            0.0 0.0;
            0.5 0.0;
            1.0 0.0;
            0.0 0.5;
            0.5 0.5;
            1.0 0.5;
            0.0 1.0;
            0.5 1.0;
            1.0 1.0;
        ]
        connectivity = [
            1 2 5;
            1 4 5;
            2 3 6;
            2 5 6;
            4 5 8;
            4 7 8;
            5 6 9;
            5 8 9;
        ]
        color = [0.0, 0.0, 0.0, 0.0, -0.375, 0.0, 0.0, 0.0, 0.0]
        scene = mesh(coordinates, connectivity, color = color, shading = false)
        wireframe!(scene[end][1], color = (:black, 0.6), linewidth = 3)
    end
    @cell "colored triangle" [mesh, polygon] begin
        mesh(
            [(0.0, 0.0), (0.5, 1.0), (1.0, 0.0)], color = [:red, :green, :blue],
            shading = false
        )
    end
    @cell "heatmap interpolation" [heatmap, interpolate, subscene, theme] begin
        using AbstractPlotting: hbox, vbox
        data = rand(50, 100)
        p1 = heatmap(data, interpolate = true)
        p2 = heatmap(data, interpolate = false)
        t = Theme(align = (:left, :bottom), raw = true, camera = campixel!)
        title1 = text(t, "Interpolate = true")
        title2 = text(t, "Interpolate = false")
        s = vbox(
            hbox(p1, title1),
            hbox(p2, title2),
        )
    end
    @cell "colored triangle" [polygon] begin
        poly(
            [(0.0, 0.0), (0.5, 1.0), (1.0, 0.0)],
            color = [:red, :green, :blue],
            strokecolor = :black, strokewidth = 2
        )
    end
    @cell "Subscenes" [image, scatter, subscene] begin
        img = rand(RGBAf0, 100, 100)
        scene = image(img, show_axis = false)
        subscene = Scene(scene, IRect(100, 100, 300, 300))
        scatter!(subscene, rand(100) * 200, rand(100) * 200, markersize = 4)
        scene
    end

    @cell "Polygons" [poly, polygon, linesegments] begin
        using GeometryTypes
        scene = Scene(resolution = (500, 500))
        points = decompose(Point2f0, Circle(Point2f0(50), 50f0))
        pol = poly!(scene, points, color = :gray, strokewidth = 10, strokecolor = :red)
        # Optimized forms
        poly!(scene, [Circle(Point2f0(50+300), 50f0)], color = :gray, strokewidth = 10, strokecolor = :red)
        poly!(scene, [Circle(Point2f0(50+i, 50+i), 10f0) for i = 1:100:400], color = :red)
        poly!(scene, [Rectangle{Float32}(50+i, 50+i, 20, 20) for i = 1:100:400], strokewidth = 2, strokecolor = :green)
        linesegments!(scene,
            [Point2f0(50 + i, 50 + i) => Point2f0(i + 70, i + 70) for i = 1:100:400], linewidth = 8, color = :purple
        )
    end

    @cell "Contour Function" [contour] begin
        r = range(-10, stop = 10, length = 512)
        z = ((x, y)-> sin(x) + cos(y)).(r, r')
        contour(r, r, z, levels = 5, color = :viridis, linewidth = 3)
    end
    @cell "Hbox" [lines, scatter, hbox] begin
        t = range(-122277.93103448274, stop=-14798.035304081845, length=29542)
        x = -42 .- randn(length(t))
        sc1 = scatter(t, x, color=:black, markersize=sqrt(length(t)/20))
        sc2 = lines(t[1:end-1], diff(x), color = :blue)
        hbox(sc2, sc1)
    end
    @cell "Customize Axes" [lines, axis] begin
        x = LinRange(0,3pi,200); y = sin.(x)
        lin = lines(x, y, padding = (0.0, 0.0), axis = (
            names = (axisnames = ("", ""),),
            grid = (linewidth = (0, 0),),
        ))
    end
    @cell "contour" [contour] begin
        y = range(-0.997669, stop = 0.997669, length = 23)
        contour(range(-0.99, stop = 0.99, length = 23), y, rand(23, 23), levels = 10)
    end

    @cell "Heatmap" [heatmap] begin
        heatmap(rand(32, 32))
    end

    @cell "Animated Scatter" [animated, scatter, updating, record] begin
        N = 10
        r = [(rand(7, 2) .- 0.5) .* 25 for i = 1:N]
        scene = scatter(r[1][:, 1], r[1][:, 2], markersize = 1, limits = FRect(-25/2, -25/2, 25, 25))
        s = scene[end] # last plot in scene
        record(scene, @outputfile(mp4), r) do m
            s[1] = m[:, 1]
            s[2] = m[:, 2]
        end
    end

    @cell "Text Annotation" [text, align, annotation] begin
        text(
            ". This is an annotation!",
            position = (300, 200),
            align = (:center,  :center),
            textsize = 60,
            font = "Blackchancery"
        )
    end

    @cell "Text rotation" [text, rotation] begin
        scene = Scene()
        pos = (500, 500)
        posis = Point2f0[]
        for r in range(0, stop = 2pi, length = 20)
            global pos, posis
            p = pos .+ (sin(r)*100.0, cos(r) * 100)
            push!(posis, p)
            t = text!(
                scene, "test",
                position = p,
                textsize = 50,
                rotation = 1.5pi - r,
                align = (:center, :center)
            )
        end
        scatter!(scene, posis, markersize = 10)
    end

    @cell "The famous iris example" [RDatasets, DataFrames, scatter, axis] begin
        using DataFrames, RDatasets # do Pkg.add.(["DataFrames", "RDatasets"]) if you don't have these packages installed
        iris = dataset("datasets", "iris")

        x = iris[:SepalWidth]
        y = iris[:SepalLength]

        scene = Scene()
        colors = [:red, :green, :blue]
        i = 1 #color incrementer
        for sp in unique(iris[:Species])
            idx = iris[:Species] .== sp
            sel = iris[idx, [:SepalWidth, :SepalLength]]
            scatter!(scene, sel[:,1], sel[:,2], color = colors[i], limits = FRect(1.5, 4.0, 3.0, 4.0))
            global i = i+1
        end
        scene
        axis = scene[Axis] # get axis
        axis[:names][:axisnames] = ("Sepal width", "Sepal length")
        scene
    end

    # @cell "Chess Game" [heatmap, scatter, interactive] begin
    #     using Base.Iterators: repeated
    #     r = 1:8
    #     board = isodd.(r .+ r')
    #     scene = Scene(resolution = (1000, 1000))
    #     heatmap!(scene, board, scale_plot = false, show_axis = false)
    #     white = ['♕', '♔', '♖', '♖', '♗', '♗', '♘', '♘', repeated('♙', 8)...]
    #     wx_positions = [4, 5, 1, 8, 3, 6, 2, 7, (1:8)...]
    #     wy_positions = [repeated(1, 8)..., repeated(2, 8)...]
    #     w_positions = Point2.(wx_positions, wy_positions)
    #     white_game = scatter!(
    #         scene, w_positions, marker = white,
    #         scale_plot = false, show_axis = false,
    #         markersize = 0.5, marker_offset = Vec2f0(-0.7)
    #     )[end]
    #     black = Char.(Int.(white) .+ 6)
    #     b_positions = Point2f0.(wx_positions, [repeated(8, 8)..., repeated(7, 8)...])
    #     black_game = scatter!(
    #         scene, b_positions, marker = black,
    #         scale_plot = false, show_axis = false,
    #         markersize = 0.5, marker_offset = Vec2f0(-0.7)
    #     )[end]
    #
    #     function move_fig!(color, figure, target)
    #         game = color == :white ? white_game : black_game
    #         game[1][][figure] = target
    #         game[1][] = game[1][]
    #     end
    #     display(scene)
    #     move_fig!(:white, 9, (1, 4))
    # end
end


text:
# --------------------------------------------
# -   Plotting script for the dubins example -
# --------------------------------------------

# Python plotting with matplotlib
using PyCall, LaTeXStrings
import PyPlot; const plt = PyPlot

include("../Models/dubins_free.jl")
include("../SCP/scp_problem.jl")

# --------------------
# Plotting with PyPlot
# --------------------
function plt_circle(ax, pos, radius; color="k", alpha=1., label="None")
    # Filled circle
    circle = plt.matplotlib.patches.Circle(pos, radius=radius,
                    color=color, fill=true, alpha=alpha)
    ax.add_patch(circle)
    # Edge of circle, with alpha 1.
    if label=="None"
        c_edge = plt.matplotlib.patches.Circle(pos, radius=radius,
                        color=color, alpha=1, fill=false)
    else
        c_edge = plt.matplotlib.patches.Circle(pos, radius=radius,
                        color=color, alpha=1, fill=false, label=label)
    end
    ax.add_patch(c_edge)
    return ax
end

function plt_solutions(scp_problem::SCPProblem, model, X_all, U_all;
                        xlims=[-0.5,3.], ylims=[0.0,6.], figsize=(8,6), B_plot_labels=true)
    N = length(X_all)
    idx = [1,2]

    fig = plt.figure(figsize=figsize)
    ax  = plt.gca()

    # Plot SCP solutions
    plt.plot(X_all[1][idx[1],:], X_all[1][idx[2],:], "--",
                    label="Initializer", linewidth=3, color="b")
    for iter = 2:length(X_all)#-3
        X = X_all[iter]
        ax.plot(X[idx[1],:], X[idx[2],:], "o--",
                    label="Iterate $(iter - 1)", linewidth=2, markersize=7)
    end

    # initial / final conditions

    plt.scatter(model.x_init[idx[1]], model.x_init[idx[2]], color="black")
    plt.scatter(model.x_final[idx[1]], model.x_final[idx[2]], color="black")

    # Plot obstacles
    for obs_i = 1:length(model.obstacles)
        p_obs, obs_radius = model.obstacles[obs_i][1], model.obstacles[obs_i][2]
        plt_circle(ax, p_obs[idx], obs_radius; color="r", alpha=0.3)
    end

    # Settings / Style / Parameters
    PyPlot.rc("text", usetex=true)
    rcParams = PyDict(plt.matplotlib["rcParams"])
    rcParams["font.size"] = 30
    rcParams["font.family"] = "Helvetica"
    plt.xlim(xlims)
    plt.ylim(ylims)
    if B_plot_labels
        plt.legend(loc="upper left", fontsize=27.5,
                                      labelspacing=0.1)
    end
    plt.grid(alpha=0.3)
    plt.xticks([0,1,2,3])
    plt.yticks([0,1,2,3])

    plt.draw()

    return fig
end

function plt_final_solution(scp_problem::SCPProblem, model, X, U)
    N = length(X_all)
    idx = [1,2]

    fig = plt.figure(figsize=(4.5, 7.5))
    ax  = plt.gca()

    # Plot final solution
    plt.plot(X[idx[1],:], X[idx[2],:], "bo-",
                linewidth=2, markersize=6)
    plt.plot(Inf*[1,1],Inf*[1,1], "b-", label="Trajectory") # for legend

    # Plot Thrust
    for k = 1:(length(X[1,:])-1)
        xk, uk =  X[:,k], U[:,k]

        uMax, mag = 23.2, 1.5

        plt.arrow(xk[idx[1]], xk[idx[2]],
                    mag*(uk[idx[1]]/uMax), mag*(uk[idx[2]]/uMax),
                    color="g")
    end
    plt.plot(Inf*[1,1],Inf*[1,1], "g-", label="Thrust") # for legend

    # Plot obstacles
    for obs_i = 1:length(model.obstacles)
        p_obs, obs_radius = model.obstacles[obs_i][1], model.obstacles[obs_i][2]
        plt_circle(ax, p_obs[idx], obs_radius; color="r", alpha=0.4)
    end
    plt.plot(Inf*[1,1],Inf*[1,1], "r-", label="Obstacle") # for legend

    # Settings / Style / Parameters
    PyPlot.rc("text", usetex=true)
    rcParams = PyDict(plt.matplotlib["rcParams"])
    rcParams["font.size"] = 20
    rcParams["font.family"] = "Helvetica"
    plt.title("Final Dubins Trajectory", pad=15)
    plt.xlim([-0.5,3.])
    plt.ylim([ 0.0,6.])
    plt.xlabel("E")
    plt.ylabel("N")
    plt.legend(loc="lower right", fontsize=18)
    plt.grid(alpha=0.3)

    plt.draw()

    return fig
end

function plt_final_angle_accel(scp_problem::SCPProblem, model, X, U)
    N = length(X_all)
    t_max = 2.

    times = collect(range(0,stop=(SCPproblem.N-1)*SCPproblem.dt,length=SCPproblem.N))[1:SCPproblem.N-1]
    norms_U = sqrt.(U[1,:].^2+U[2,:].^2+U[3,:].^2)

    fig = plt.figure(figsize=(5.5, 7.5))

    # -------------
    # Tilt Angle
    plt.subplot(2,1,1)
    tilt_angles = U[3,:]./norms_U
    tilt_angles = (180. / pi) * tilt_angles
    plt.plot(times, tilt_angles, "bo-",
                linewidth=1, markersize=4)
    # max tilt angle
    theta_max = (180/pi) * (pi/3.0)
    plt.plot([0,t_max], theta_max*ones(2),
                    color="r", linestyle="dashed", linewidth=2)
    plt.fill_between([0,t_max], theta_max*ones(2), 90,
                        color="r", alpha=0.2)
    # Params
    plt.title("Tilt Angle", pad=10)
    plt.xlim([0, t_max])
    plt.ylim([0, 65   ])
    # plt.xlabel("Time [s]")
    plt.ylabel(L"$\theta$ [deg]")
    plt.grid(alpha=0.3)
    plt.draw()

    # -------------
    # Acceleration
    plt.subplot(2,1,2)
    fig.tight_layout(pad=2.0)

    plt.plot(times, norms_U, "bo-",
                linewidth=1, markersize=4)

    # max/min acceleration
    a_min, a_max = 0.6, 23.2
    plt.plot([0,t_max], a_max*ones(2),
                    color="r", linestyle="dashed", linewidth=2)
    plt.fill_between([0,t_max], a_max*ones(2), 90,
                        color="r", alpha=0.2)
    plt.plot([0,t_max], a_min*ones(2),
                    color="r", linestyle="dashed", linewidth=2)
    plt.fill_between([0,t_max], a_min*ones(2), -5,
                        color="r", alpha=0.2)

    # Parameters / Settings / Labels
    PyPlot.rc("text", usetex=true)
    rcParams = PyDict(plt.matplotlib["rcParams"])
    rcParams["font.size"] = 20
    rcParams["font.family"] = "Helvetica"
    plt.title("Cmd. Acceleration", pad=10)
    plt.xlim([0, t_max])
    plt.ylim([-1, 25  ])
    plt.xlabel("Time [s]")
    plt.ylabel(L"$u [m/s^2]$")
    plt.grid(alpha=0.3)
    plt.draw()

    return fig
end

max_stars_repo_name: StanfordASL/jlSCP

max_issues_repo_name: StanfordASL/jlSCP

max_forks_repo_name: StanfordASL/jlSCP


text:
## Two electron integrals

function coulomb(aexpn::Float64,ax::Float64,ay::Float64,az::Float64,
                 aI::Int64,aJ::Int64,aK::Int64,
                 bexpn::Float64,bx::Float64,by::Float64,bz::Float64,
                 bI::Int64,bJ::Int64,bK::Int64,
                 cexpn::Float64,cx::Float64,cy::Float64,cz::Float64,
                 cI::Int64,cJ::Int64,cK::Int64,
                 dexpn::Float64,dx::Float64,dy::Float64,dz::Float64,
                 dI::Int64,dJ::Int64,dK::Int64)
    # This is the slow method of computing integrals from Huzinaga et al.
    # Use the HRR/VRR scheme from Head-Gordon & Pople instead

    rab2 = dist2(ax-bx,ay-by,az-bz)
    rcd2 = dist2(cx-dx,cy-dy,cz-dz)

    px,py,pz = gaussian_product_center(aexpn,ax,ay,az,bexpn,bx,by,bz)
    qx,qy,qz = gaussian_product_center(cexpn,cx,cy,cz,dexpn,dx,dy,dz)
    rpq2 = dist2(px-qx,py-qy,pz-qz)

    g1 = aexpn+bexpn
    g2 = cexpn+dexpn
    delta = 0.25*(1/g1+1/g2)

    Bx = Barray(aI,bI,cI,dI,px,ax,bx,qx,cx,dx,g1,g2,delta)
    By = Barray(aJ,bJ,cJ,dJ,py,ay,by,qy,cy,dy,g1,g2,delta)
    Bz = Barray(aK,bK,cK,dK,pz,az,bz,qz,cz,dz,g1,g2,delta)

    s = 0
    #println("$(aI+bI+cI+dI),$(aJ+bJ+cJ+dJ),$(aK+bK+cK+dK)")
    for I in 0:(aI+bI+cI+dI)
        for J in 0:(aJ+bJ+cJ+dJ)
            for K in 0:(aK+bK+cK+dK)
                #println("coul: $I,$J,$K,$(Bx[I+1]),$(By[J+1]),$(Bz[K+1])")
                s += Bx[I+1]*By[J+1]*Bz[K+1]*Fgamma(I+J+K,0.25*rpq2/delta)
            end
        end
    end
    return 2pi^(2.5)/(g1*g2*sqrt(g1+g2))*exp(-aexpn*bexpn*rab2/g1)*exp(-cexpn*dexpn*rcd2/g2)*s
end

function coulomb(a::PGBF,b::PGBF,c::PGBF,d::PGBF)
    return a.norm*b.norm*c.norm*d.norm*coulomb(a.expn,a.x,a.y,a.z,a.I,a.J,a.K,
        b.expn,b.x,b.y,b.z,b.I,b.J,b.K,
        c.expn,c.x,c.y,c.z,c.I,c.J,c.K,
        d.expn,d.x,d.y,d.z,d.I,d.J,d.K)
end

fB(i::Int64,l1::Int64,l2::Int64,p::Float64,a::Float64,b::Float64,r::Int64,g::Float64) = binomial_prefactor(i,l1,l2,p-a,p-b)*B0(i,r,g)
B0(i::Int64,r::Int64,g::Float64) = fact_ratio2(i,r)*(4g)^(r-i)
fact_ratio2(a::Int64,b::Int64) = factorial(a)factorial(b)/factorial(a-2b)

function Bterm(i1::Int64,i2::Int64,r1::Int64,r2::Int64,u::Int64,
               l1::Int64,l2::Int64,l3::Int64,l4::Int64,
               Px::Float64,Ax::Float64,Bx::Float64,Qx::Float64,Cx::Float64,Dx::Float64,
               gamma1::Float64,gamma2::Float64,delta::Float64)
    # THO eq. 2.22
    #print("Bterm($i1,$i2,$r1,$r2,$u,$l1,$l2,$l3,$l4,$Px,$Ax,$Bx,$Qx,$Cx,$Dx,$gamma1,$gamma2,$delta)=")
    val = (-1)^(i2+u)*fB(i1,l1,l2,Px,Ax,Bx,r1,gamma1)*fB(i2,l3,l4,Qx,Cx,Dx,r2,gamma2)*(
            fact_ratio2(i1+i2-2*(r1+r2),u)*(Qx-Px)^(i1+i2-2*(r1+r2)-2*u)/delta^(i1+i2-2*(r1+r2)-u))
    #println("$val")
    return val
end

function Barray(l1::Int64,l2::Int64,l3::Int64,l4::Int64,p::Float64,a::Float64,b::Float64,
                q::Float64,c::Float64,d::Float64,g1::Float64,g2::Float64,delta::Float64)
    Imax = l1+l2+l3+l4+1
    B = zeros(Float64,Imax)
    for i1 in 0:(l1+l2)
        for i2 in 0:(l3+l4)
            for r1 in 0:div(i1,2)
                for r2 in 0:div(i2,2)
                    for u in 0:(div(i1+i2,2)-r1-r2)
                        I = i1+i2-2*(r1+r2)-u
                        B[I+1] += Bterm(i1,i2,r1,r2,u,l1,l2,l3,l4,p,a,b,q,c,d,g1,g2,delta)
                    end
                end
            end
        end
    end
    return B
end

coulomb(a::CGBF,b::CGBF,c::CGBF,d::CGBF) = contract(coulomb,a,b,c,d)

function test_two_terms()
    @testset "test_two_terms" begin
    @test fB(0,0,0,0.0,0.0,0.0,0,2.0) == 1
    @test fB(0,0,0,1.0,1.0,1.0,0,2.0) == 1
    @test fB(0,0,0,0.0,0.0,0.0,0,2.0 ) == 1
    @test fB(1,0,1,0.0,0.0,0.0,0,2.0 ) == 0.125
    @test B0(0,0,2.0) == 1
    @test fact_ratio2(0,0) == 1
    @test Bterm(0,0,0,0,0,0,0,0,0,0.0,0.0,0.0,0.0,0.0,0.0,2.0,2.0,0.25)==1
    @test Bterm(0,1,0,0,0,0,0,0,1,0.0,0.0,0.0,0.0,0.0,0.0,2.0,2.0,0.25)==0
    end #testset
end

function test_coul1()
    @testset "test_coul1" begin
    s = pgbf(1.0)
    px = pgbf(1.0,0,0,0,1,0,0)
    @test coulomb(s,s,s,px)==0 # 0
    @test coulomb(s,s,px,px) ≈ 0.9403159725793305
    end #testset
end

function coulomb_hgp(a::PGBF,b::PGBF,c::PGBF,d::PGBF)
    return a.norm*b.norm*c.norm*d.norm*hrr(a.expn,a.x,a.y,a.z,a.I,a.J,a.K,
        b.expn,b.x,b.y,b.z,b.I,b.J,b.K,
        c.expn,c.x,c.y,c.z,c.I,c.J,c.K,
        d.expn,d.x,d.y,d.z,d.I,d.J,d.K)
end
coulomb_hgp(a::CGBF,b::CGBF,c::CGBF,d::CGBF) = contract(coulomb_hgp,a,b,c,d)


text:
#####
# Synthetic Bench
#####
#
# For a Bench-like glacier as used in SHMIP
# https://shmip.bitbucket.io/instructions.html#sec-2-2

function LoadPara(gid::SyntheticGlacier; pl_kwargs...)
    pl_kwargs = deepcopy(pl_kwargs)
    pl = LoadPara(; gid=gid, pl_kwargs...)
    merge!((a,b)->a, pl.dataset_opts[ParaData][:pm],
           Dict(
               :bandsize => 30,
               :window_dem_smooth=>0,
               :window_width_smooth=>0, # can make a noisy tau
               :window=>100.0)
           )
    for DK in datakinds()
        pl.dataset_LOADERS[DK] = SyntheticBenchLoader
    end
    return pl
end

# surface functions
_bench_g(y) = 0.5e-6 * abs(y)^3
_bench_ginv(x) = (x/0.5e-6).^(1/3)
_bench_surface(x,y) = x<0 ? zero(x) : 100(x+200)^(1/4) + 1/60*x - 2e10^(1/4) + 1

# bed function from SHMIP:
const para_bench = 0.05  # gives depth of overdeepening
_bench_h(x, para) = (-4.5*x/6e3 + 5) * (_bench_surface(x,0)-_bench_f(x, para)) /
    (_bench_surface(x,0)-_bench_f(x, para_bench)+eps())
_bench_f(x,para) = (_bench_surface(6e3,0) - para*6e3)/6e3^2 * x^2 + para*x
_bench_bed(x,y, para) = _bench_f(x,para) + _bench_g(y) * _bench_h(x,para)

# outline, upper branch
_bench_outline(x) = _bench_ginv( (_bench_surface(x,0)-_bench_f(x,0.05))/(_bench_h(x,0.05)+eps()) )

# bed function from forward model

function load!(::DataSet{OutlineData, SyntheticBenchLoader}, gid::ANY)
    x = 0:10:6e3
    _outline = _bench_outline.(x)
    x = vcat(0, x[1:end-1], reverse(x), 0)
    y = vcat(0, -_outline[1:end-1], reverse(_outline), 0)
    outline = hcat(x, y)'
    iscalving = false
    isicecap = false
    proj = "n/a"
    outline = Outline(gid, outline)
    return outline, iscalving::Bool, isicecap::Bool, proj::String
end
const _dx = 25.0
const _buf = 10*_dx
function load!(ds::DataSet{DEMData, SyntheticBenchLoader}, outline,
               pp::Phys)::Tuple{Gridded,Matrix{Bool},Matrix{Bool},Matrix{F},F,Tuple,Any}
    opts = ds.opts
    grid = get(opts, :grid, (-_buf:_dx:6000+_buf, -500-_buf:_dx:500+_buf))
    cropbox = get(opts, :cropbox, (minimum(grid[1]), maximum(grid[1]), minimum(grid[2]), maximum(grid[1])))
    dist_land_maxdist = get(opts, :dist_land_maxdist, -10)

    dem = Gridded(grid[1], grid[2],_bench_surface.(grid[1], grid[2]'))
    glaciermask, landmask = maskit(dem, outline)
    dist_land = dist_to_mask(dem, landmask, glaciermask, dist_land_maxdist)
    proj = "n/a"
    return dem, glaciermask, landmask, dist_land, dist_land_maxdist, cropbox, proj
end

"""
opts:
- use_forward_model -- if true use fwd to generate IV field
  - if this is set then also pass in the model parameters:
    opts[:fwd_para] = (gl, pp, pm, pn)
  - as a hack it returns the function of the bed elevation inside
    the opts-struct: opts[:fwd_iv_fn] and opts[:fwd_iv_err_fn]
- iv_grid -- grid to return iv on
- err_sigma -- add a random error with this std
- err_cor_length -- and this correlation length
"""
function load!(ds::DataSet{IVData, SyntheticBenchLoader}, dem, glaciermask, target_proj, outline, cropbox,
               pp::Phys)::Tuple{SomeData,Matrix{Bool}}
    opts = ds.opts
    @get! use_forward_model = opts | false
    if use_forward_model
        gl, pp, pm, pn = opts[:fwd_para]
        @get! iv_grid = opts | (gl.dem.x[1:3:end], gl.dem.y[1:3:end])
        @get! err_cor_length = opts | 200.0
        @get! err_sigma = opts | 2.0

        # run forward model
        iv1d, iv2d = fwdm(gl, pp, pm, pn)[[3,6]]
        itp = Interp.interpolate((gl.dem.x, gl.dem.y), iv2d,
                                 Interp.Gridded(Interp.Linear()) )
        ivfn = (x,y) -> itp[x,y]
        opts[:fwd_iv_fn] = ivfn # return the "true" thickness function also
        opts[:fwd_true_iv1d] = iv1d

        # add some error (poor-man's spatially correlated errors)

        err_grid = (linspace(gl.dem.x[1], gl.dem.x[end], max(2, (gl.dem.x[end]-gl.dem.x[1])÷err_cor_length+1)),
                    linspace(gl.dem.y[1], gl.dem.y[end], max(2, (gl.dem.y[end]-gl.dem.y[1])÷err_cor_length+1)))

        itp_err = Interp.interpolate(err_grid, randn(length.(err_grid)...)*err_sigma,
                                 Interp.Gridded(Interp.Linear()) )
        err_fn = (x,y) -> itp_err[x,y]
        opts[:fwd_iv_err_fn] = err_fn

        # make it
        ivs = zeros(F, length.(iv_grid)...)
        ivmask = falses(length.(iv_grid)...)
        for (j,y) in enumerate(iv_grid[2])
            for (i,x) in enumerate(iv_grid[1])
                if 0<=x<=6e3 && _bench_outline(x)>=abs(y)
                    ivs[i,j] = ivfn(x,y) + err_fn(x,y)
                    # avoid negative values:
                    ivs[i,j] = max(ivs[i,j], 0)
                    ivmask[i,j] = true
                else
                    ivs[i,j] = 0
                    ivmask[i,j] = false
                end
            end
        end
        iv = SomeData(vals=Gridded(iv_grid..., ivs), sigma=err_sigma)
    else
        # only one point with velocity
        iv = SomeData(vals=Gridded(2000.0:2000, 0.0:0, 30.0+zeros(1,1)), sigma=5)
        ivmask = trues(1,1)
    end
    return iv, ivmask
end

"""
opts:
- use_forward_model=false -- if true, use the forward-model to generate
  the ice thickness
  - if this is set then also pass in the model parameters:
    opts[:fwd_para] = (gl, pp, pm, pn)
  - as a hack it returns the function of the bed elevation inside
    the opts-struct:  opts[:fwd_thick_fn] and opts[:fwd_thick_err_fn]
- para_bench=0.05 -- how much of an overdeepeing there is
- profile_xlocs=[2000.0,4000.0] -- locations of cross-profiles
- profile_ylocs=[50.0] -- locations of along-profiles
- profile_margin_dist=100.0 -- distance from margin at profile start
- profile_step=10.0 -- step between profile points
- err_sigma -- add a random error with this std
- err_cor_length -- and this correlation length
"""
function load!(ds::DataSet{ThicknessData, SyntheticBenchLoader}, dem,
               glaciermask, target_proj, outline, cropbox,
               pp::Phys)
    opts = ds.opts
    @get! use_forward_model = opts | false
    if use_forward_model
        gl, pp, pm, pn = opts[:fwd_para]
        h1d, hs2d = fwdm(gl, pp, pm, pn)[[1,4]]
        itp = Interp.interpolate((gl.dem.x, gl.dem.y), hs2d,
                                 Interp.Gridded(Interp.Linear()) )
        hfn = (x,y) -> itp[x,y]
        opts[:fwd_thick_fn] = hfn # return the "true" thickness function also
        opts[:fwd_true_h1d] = h1d
    else
        @get! para_bench = opts | 0.05
        # TODO add length-profiles, add arbitrary trajectories
        hfn = (x,y) -> _bench_surface(x,y) - _bench_bed(x,y,para_bench)
    end

    # produce data along profiles:
    @get! profile_xlocs = opts | [850, 2000.0, 4000.0]
    @get! profile_ylocs = opts | [-100.0, 200.0]
    @get! profile_margin_dist = opts | 100.0
    @get! profile_step = opts | 10.0
    @get! err_sigma = opts | 10.0
    @get! err_cor_length = opts | 100.0
    @assert err_cor_length<1080 "err_cor_length too big"

    x = F[]
    y = F[]
    thick = F[]
    surf = F[]
    splits = UnitRange{Int}[]
    hfns_err = []
    for xl in profile_xlocs
        w = _bench_outline(xl)-profile_margin_dist
        ys = -w:profile_step:w
        xs = ys*0 + xl

        # make an error
        ys_err = linspace( ys[1], ys[end], max(2, (ys[end]-ys[1])÷err_cor_length+1))
        hfn_err = let
            tmp = Interp.interpolate((ys_err,), randn(length(ys_err))*err_sigma,
                                     Interp.Gridded(Interp.Linear()) )
            (x,y) -> tmp[y]
        end

        push!(hfns_err, hfn_err)

        s = _bench_surface.(xs,ys)
        push!(splits, length(thick)+1:length(thick)+length(xs))
        append!(surf, s)
        append!(thick, hfn.(xs,ys) + hfn_err.(xs,ys))
        append!(x,xs)
        append!(y,ys)
    end

    for yl in profile_ylocs
        # find starting point
        if abs(yl)<75
            # the outline has \inf slope at 0,6e3:
            x0, x1 = 1.0, 6e3-5
        else
            rootfn = x -> abs(yl) - _bench_outline(x)
            x0 = Roots.fzero(rootfn, 1e-3, 2e3)
            x1 = Roots.fzero(rootfn, 4e3, 6e3-1e-3)
        end
        xs = x0+profile_margin_dist:profile_step:x1-profile_margin_dist
        ys = 0*xs + yl

        # make an error
        xs_err = linspace( xs[1], xs[end], max(2, (xs[end]-xs[1])÷err_cor_length+1))
        hfn_err = let
            tmp = Interp.interpolate((xs_err,), randn(length(xs_err))*err_sigma,
                                     Interp.Gridded(Interp.Linear()) )
            (x,y) -> tmp[y]
        end
        push!(hfns_err, hfn_err)

        s = _bench_surface.(xs,ys)
        push!(splits, length(thick)+1:length(thick)+length(xs))
        append!(surf, s)
        append!(thick, hfn.(xs,ys) + hfn_err.(xs,ys))
        append!(x,xs)
        append!(y,ys)
    end
    opts[:fwd_thick_err_fn] = hfns_err

    thickness = SomeData(vals=Traj(x,y, thick, surf, splits),
                         sigma=err_sigma)
    return thickness::SomeData
end

"""
Construct a plausible SMB field:

- ELA set to median elevation (2D) (can be set)
- gradient 0.008 (range 0.002 (continental) to 0.02 (maritime))
- gradient reduction in accumulation area: 0.55 (range 0.1 to 1)
"""
function load!(ds::DataSet{BdotData, SyntheticBenchLoader}, dem, glaciermask, target_proj, outline, cropbox,
              pp::Phys)
    opts = ds.opts
    ela_elevation = get(opts, :ela_elevation, median(dem.v[glaciermask]))
    gradient = get(opts, :smb_gradient, 0.008)
    gradient_corr_accumulation = get(opts, :gradient_corr_accumulation, 0.55)
    bd = dem.v-ela_elevation
    for i in eachindex(bd)
        b = bd[i]
        bd[i] = b>0 ? b*gradient*gradient_corr_accumulation : b*gradient
    end
    # scale to m ice / a
    scale_m_ice!(bd, pp)

    bdot = SomeData(vals=Gridded(dem.x, dem.y, bd),
                    sigma=get(opts, :sigma, 2.0))
    return bdot::SomeData
end

"""
Makes a 2D dh/dt field such that: first dh/dt= 0 from top to ELA
elevation, then such that flux at terminus is given.

NOTE: this can give a glacier of zero thickness in the middle as the
flux can go negative, as long as it comes back to the given terminus
flux.

- ELA is taken where bdot==0.  Or can be specified in ds.opts[:ela_elevation].
  Note that if the ELA is inconsistent with bdot, then this produces a inconsistent
  dhdt in the sense that flux will go negative for some of the glacier.
- scaled such that terminus flux is ok (or set manually with
  ds.opts[:dhdt_scale].

- constraints: dh/dt in accum<=bdot, dh/dt in abl >=bdot
  - at least for the 1D case
- then just the usual uncertainty on btilde?
"""
function load!(ds::DataSet{DhdtData, SyntheticBenchLoader},
               bdot, terminus_flux,
               dem, glaciermask, target_proj, outline, cropbox,
               iscalving, pp::Phys)
    opts = ds.opts
    dx = step(dem.x)

    twoD_bdot = if isa(bdot.vals, Gridded) # 2D bdot
        bdot_2D = bdot.vals.v
        true
    elseif isa(bdot.vals, Gridded1d) # 1D bdot
        error("Not supported. This introduces mistakes, in particular for small glaciers.  Extrapolate a 1D bdot to 2D when loading.")
        false
    else
        error("Not supported type of bdot: $bdot, typeof = $(typeof(bdot))")
    end

    # get ELA
    if haskey(opts, :ela_elevation)
        # set by hand
        ela_elevation = opts[:ela_elevation]
    else
        if twoD_bdot
            bands, bandi = VAWTools.bin_grid(dem, 30, glaciermask, min_bands=4)
            bdot1d = VAWTools.map_onto_bands(bandi, bdot.vals.v, mean, FILL)
        else
            bdot1d = bdot.vals.v
            bands = bdot.vals.x
        end
        i = findlast(bdot1d.<=0)
        if i==0
            ela_elevation = minimum(dem.v[glaciermask])
        elseif i==length(bdot1d)
            ela_elevation = maximum(dem.v[glaciermask])
        else
            ela_elevation = bands[i]
        end
    end
    has_accumulation_area = ela_elevation < maximum(dem.v[glaciermask])
    has_ablation_area = ela_elevation > minimum(dem.v[glaciermask])

    msk_bdot = glaciermask .& (bdot_2D.!=FILL) # in case there are masked bdot values within glacier
    smb = sum(bdot_2D[msk_bdot]*dx^2)

    if has_accumulation_area
        if smb<0
            # Trial dhdt has linear dependence on elevation in ablation zone, constant in accumulation zone
            dd = dem.v-ela_elevation
            # above ELA dh/dt== small.  Do not use 0 when there is no ablation area.
            small = has_ablation_area ? 0.0 : 0.1
            for i in eachindex(dd)
                d = dd[i]
                dd[i] = d>0 ? small : d
            end
            dd[.!glaciermask] .= FILL
        else # growing glacier
            dd = dem.v-ela_elevation
        end
    else
        # Trial dhdt is flat
        dd = dem.v*0 + 1
        dd[.!glaciermask] .= FILL
    end
    # Now scale such that flux is satisfied
    if haskey(opts, :dhdt_scale)
        scale = opts[:dhdt_scale]
    else
        if twoD_bdot
            bdot_2D = bdot.vals.v
        else
            error("needs updating")
            # extrapolate to 2D
            bands, bandi = VAWTools.bin_grid(dem, -step(bdot.vals.x), glaciermask) #, min_bands=0)
            bdot_2D = VAWTools.map_back_to_2D(size(dem.v), bandi, reverse(bdot.vals.v))
            ## check that collapsing back to 1D is identity
            # @show tmp_ = VAWTools.map_onto_bands(bandi, bdot_2D)
            # @show bdot.vals.v[isnan.(tmp_)]
            # @assert all(tmp_[length.(bandi).>0] .≈ bdot.vals.v[length.(bandi).>0])
        end

        # Scale such that flux at terminus as requested
        fkind = terminus_flux[1]
        vals = terminus_flux[2].vals
        if fkind==FluxKind.abs
            # absolute flux
            target_flux = vals
        elseif fkind==FluxKind.rel_bdot_teminus # target_flux at terminus relative to flux due to expectation bdot
            bdotflux_at_terminus = sum(bdot_2D[glaciermask])*dx^2
            target_flux = bdotflux_at_terminus*vals
        elseif fkind==FluxKind.rel_bdot_max # maximum flux due to expectation bdot (== total accumulation rate)
            inds = (bdot_2D.>0) .& glaciermask
            max_bdotflux = sum(bdot_2D[inds])*dx^2
            target_flux = max_bdotflux*vals
        elseif fkind==FluxKind.rel_btilde_teminus
            error("Flux-kind: $fkind not supported with this loader $ds")
        elseif fkind==FluxKind.none && !iscalving
            # land terminating
            target_flux = 0.0
        else
            error("Unknown flux-kind: $fkind")
        end

        unscaled_dhdt = sum(dd[msk_bdot]*dx^2)
        scale = (smb-target_flux)/unscaled_dhdt
    end
    # TODO: here a offset instead of a scale could also be used.
    #       But how to determine its value?  Maybe minimize sum(abs(dhdt))?
    dd *= scale
    dd[.!glaciermask] .= FILL
    dhdt = SomeData(vals=Gridded(dem.x, dem.y, dd),
                    sigma=get(opts, :sigma, 2.0))
    return dhdt::SomeData
end
function load!(ds::DataSet{FslData, SyntheticBenchLoader}, dem, glaciermask, target_proj, outline, cropbox,
               pp::Phys)
    opts = ds.opts
    fsl = get(opts, :fsl, [0.5,0.5] ) # sliding at top and bottom
    top = maximum(dem.v[glaciermask])
    bottom = max(0,0, minimum(dem.v[glaciermask]))
    ele = linspace(bottom, top, length(fsl))
    fsl = SomeData(vals=Gridded1d(ele, collect(fsl)),
                   range_abs=(0,0.999),
                   sigma=Inf)
    return fsl::SomeData
end
function load!(ds::DataSet{TempData, SyntheticBenchLoader}, dem, glaciermask, target_proj, outline, cropbox,
              pp::Phys)
    opts = ds.opts
    temp = get!(opts, :temp, [-0.5, -0.5])
    top = maximum(dem.v[glaciermask])
    bottom = max(0,0, minimum(dem.v[glaciermask]))
    ele = get!(opts, :ele, linspace(bottom, top, length(temp)))
    temp = SomeData(vals=Gridded1d(ele, collect(temp)),
                    sigma=1,
                    range_abs=(-20,0))
    return temp::SomeData
end
function load!(ds::DataSet{TerminusFluxData, SyntheticBenchLoader}, dem, glaciermask, iscalving, pp::Phys)
    fluxkind = get(ds.opts, :fluxkind, FluxKind.rel_bdot_max)
    if iscalving
        # assume that half the balance is lost through flux
        terminus_flux = get(ds.opts, :terminus_flux, SomeData(vals=1/2,
                                                              sigma=0.1,
                                                              range_abs=(0.005,10.0)
                                                              )
                            )
    else
        # # approximately zero flux at terminus +/- 1% of max bdot flux accumulation
        # terminus_flux = get(ds.opts, :terminus_flux, SomeData(vals=0.0,
        #                                                       sigma=0.001
        #                                                       )
        #                     )
        fluxkind = FluxKind.none
        terminus_flux = SomeData(nothing)
    end

    return fluxkind, terminus_flux
end


text: largez05() = largez(5)


text:
########################################################################
#                         CONVERSION TO EXGRAPH                        #
########################################################################

make_name(id::Int) = Symbol("%$id")
make_name(op::AbstractOp) = Symbol("%$(op.id)")

unmake_name(x::Symbol) = parse(Int, string(x)[2:end])

Espresso.to_expr(op::Input) = :()
Espresso.to_expr(op::Assign) = Expr(:(=), make_name(op), make_name(op.src_id))

function Espresso.to_expr(op::Constant)
    val = op.val isa Symbol ? QuoteNode(op.val) : op.val
    return Expr(:(=), make_name(op), val)
end

function Espresso.to_expr(op::Call)
    arg_names = map(make_name, op.args)
    call = Expr(:call, op.fn, arg_names...)
    # assign = ((op.val isa Array || is_cuarray(op.val)) && op.fn == broadcast) ? :(.=) : :(=)
    assign = :(=)
    return Expr(assign, make_name(op.id), call)
end

function Espresso.to_expr(tape::Tape)
    body = Expr(:block)
    for op in tape
        push!(body.args, to_expr(op))
    end
    return body
end

# alternatives for simplegrad
function to_simple_expr(op::Call)
    arg_names = map(make_name, op.args)
    call = Expr(:call, op.fn, arg_names...)
    return Expr(:(=), make_name(op.id), call)
end

to_simple_expr(op::AbstractOp) = to_expr(op)

########################################################################
#                       CODE GENERATION                                #
########################################################################

"""
Rewrite :(Z = X * Y) into :(mul!(Z, X, Y)), but only if X and Y are arrays
"""
function rewrite_mul(tape::Tape, ex::Expr)
    @assert ex.head == :block
    new_ex = Expr(:block)
    for subex in ex.args
        if (matchingex(:(_R = $*(_X, _Y)), subex)
            && tape[unmake_name(subex.args[1])].val isa AbstractArray
            && all(tape[unmake_name(subex.args[2].args[i])].val isa AbstractArray for i=2:3))
            new_subex = rewrite(subex, :(_Z = $*(_X, _Y)), :($mul!(_Z, _X, _Y)))
            push!(new_ex.args, new_subex)
        else
            push!(new_ex.args, subex)
        end
    end
    return new_ex
end

"""
Generate initialization block of function code
"""
function generate_prologue(tape::Tape)
    ex = Expr(:body)
    for op in tape
        name = make_name(op)
        push!(ex.args, :($name = $(op).val))
    end
    return ex
end

"""
Main part of generated function code
"""
function generate_body(tape::Tape)
    # ret_var_ids = vcat(values(tape.derivs) |> collect, [tape.resultid])
    # ret_var_names = [Symbol("%$id") for id in ret_var_ids]
    # exg = to_exgraph(tape)
    # exg = Espresso.fuse_broadcasting(exg)  # TODO: this breaks test_examples.jl
    # exg = Espresso.fuse_assigned(exg; outvars=ret_var_names)
    # exg = Espresso.eliminate_common(exg)
    # ex = Espresso.to_expr_kw(exg)
    # ex = rewrite_all(ex, INPLACE_RULES; phs=Set([:X, :Y, :Z]))
    ex = to_expr(tape)
    ex = rewrite_mul(tape, ex)
    return ex
end

"""
Generate last block of function code
"""
function generate_epilogue(tape::Tape)
    ex = Expr(:body)
    for op in tape
        name = make_name(op)
        push!(ex.args, :($(op).val = $name))
    end
    return ex
end

"""
Generate function expression from the tape, binding all variables to tape's buffers
and optimizing code.
"""
function generate_function_expr(tape::Tape)
    fn_ex = :(function $(gensym("tape_fn"))() end)
    fn_ex_body = fn_ex.args[2]
    prologue = generate_prologue(tape)
    body = generate_body(tape)
    epilogue = generate_epilogue(tape)
    # abusing LineNumberNode to insert comments
    push!(fn_ex_body.args, LineNumberNode(0, :prologue))
    push!(fn_ex_body.args, prologue.args...)
    push!(fn_ex_body.args, LineNumberNode(0, :body))
    push!(fn_ex_body.args, body.args...)
    push!(fn_ex_body.args, LineNumberNode(0, :epilogue))
    push!(fn_ex_body.args, epilogue.args...)
    return fn_ex
end

"""
Generate function expression from the tape without any optimizations or binding to
tape's buffers.

`ret_grad` optional parameter controls whether function should return only
value of tape[resultid] (default) or also values of gradient nodes
"""
function generate_function_expr_unbound(tape::Tape; ret_grad=false)
    # fn_args = [Expr(:(::), make_name(inp), typeof(inp.val)) for inp in tape if isa(inp, Input)]
    fn_args = [make_name(inp) for inp in tape if isa(inp, Input)]
    fn_ex = :(function $(gensym("tape_fn"))($(fn_args...)) end)
    fn_ex_body = fn_ex.args[2]
    for op in tape
        if !isa(op, Input)
            ex = op |> to_simple_expr
            push!(fn_ex_body.args, ex)
        end
    end
    if ret_grad
        res_var = make_name(tape[tape.resultid])
        grad_vars = []
        for op in tape
            if op isa Input
                if haskey(tape.derivs, op.id)
                    push!(grad_vars, make_name(tape.derivs[op.id]))
                else
                    push!(grad_vars, nothing)
                end
            end
        end
        ret_tuple = Expr(:tuple, res_var, grad_vars...)
        push!(fn_ex_body.args, :(return $ret_tuple))
    else
        push!(fn_ex_body.args, :(return $(make_name(tape[tape.resultid]))))
    end
    return fn_ex
end

########################################################################
#                            COMPILATION                               #
########################################################################

function compile(tape::Tape; bind=true, ret_grad=false)
    fn_ex = bind ? generate_function_expr(tape) : generate_function_expr_unbound(tape; ret_grad=ret_grad)
    return Core.eval(@__MODULE__, fn_ex)
end

function compile!(tape::Tape)
    tape.compiled = compile(tape)
end


text:
export conversion,
    extend,
    extension,
    extensionsize,
    restrict,
    restriction,
    restrictionsize,
    evaluation,
    approximation

"""
```
conversion(::Type{T}, src::Dictionary, dest::Dictionary; options...)
```

Construct an operator with element type `T` that converts coefficients from the
source dictionary to coefficients of the destination dictionary.
"""
conversion

@deprecate op_eltype operatoreltype
@deprecate evaluation_operator evaluation
@deprecate interpolation_operator evaluation
@deprecate approximation_operator approximation
@deprecate transform_operator transform
@deprecate extension_operator extension
@deprecate restriction_operator restriction
@deprecate differentiation_operator differentiation
@deprecate antidifferentiation_operator antidifferentiation
@deprecate DictFun Expansion
@deprecate DictFun1d Expansion1d
@deprecate DictFun2d Expansion2d
@deprecate DictFun3d Expansion3d
@deprecate DictFun4d Expansion4d

operatoreltype(Φ::Dictionary...) = promote_type(map(coefficienttype, Φ)...)

for op in (:conversion, :extension, :restriction, :evaluation)
    # Only dictionaries are given: compute the eltype
    @eval $op(dicts::Dictionary...; options...) = $op(operatoreltype(dicts...), dicts...; options...)
    # Only grids are given: make a GridBasis dictionary
    @eval $op(::Type{T}, src::AbstractGrid, dest::AbstractGrid; options...) where {T} =
        $op(T, GridBasis{T}(src), GridBasis{T}(dest); options...)
end

############################
# Extension and restriction
############################

# Conversion between dictionaries of the same type: decide between extension and restriction.
conversion(::Type{T}, src::D, dest::D; options...) where {T,D <: Dictionary} =
    extension_restriction(T, src, dest; options...)

function extension_restriction(T, src::Dictionary, dest::Dictionary; options...)
    if dimensions(src) == dimensions(dest)
        IdentityOperator{T}(src, dest)
    elseif length(src) < length(dest)
        extension(T, src, dest; options...)
    elseif length(src) > length(dest)
        restriction(T, src, dest; options...)
    else
        # lengths are equal but dimensions are not
        error("Don't know how to convert between two `$(repr(src))` dictionaries with different structure.")
        # Perhaps we can do something elementwise
    end
end

# Transforming between dictionaries with the same type is the same as restricting
hastransform(src::D, dest::D) where {D <: Dictionary} = true

extension(::Type{T}, src::D, dest::D; options...) where {T,D} =
    error("Don't know how to extend dictionary $(name(src)) from size $(size(src)) to size $(size(dest))")

extension(::Type{T}, src::Dictionary, dest::Dictionary; options...) where {T} =
    error("Can't do extension between dictionaries of different type")

restriction(::Type{T}, src::D, dest::D; options...) where {T,D} =
        error("Don't know how to extend dictionary $(name(src)) from size $(size(src)) to size $(size(dest))")

restriction(::Type{T}, src::Dictionary, dest::Dictionary; options...) where {T} =
    error("Can't do restriction between dictionaries of different type")

"""
Return a suitable length to extend to, for example one such that the corresponding grids are nested
and function evaluations can be shared. The default is twice the length of the current set.
"""
extensionsize(s::Dictionary) = 2*length(s)

extend(Φ::Dictionary) = resize(Φ, extensionsize(Φ))

# only a single dictionary is given
extension(::Type{T}, src::Dictionary; options...) where {T} =
    extension(T, src, extend(src); options...)

"""
Return a suitable length to restrict to, for example one such that the corresponding grids are nested
and function evaluations can be shared. The default is half the length of the current set.
"""
restrictionsize(Φ::Dictionary) = length(Φ)>>1

restrict(Φ::Dictionary) = resize(Φ, restrictionsize(Φ))

restriction(::Type{T}, src::Dictionary; options...) where {T} =
    restriction(T, src, restrict(src); options...)

# For convenience with dispatch, add the grids as extra arguments when only
# GridBasis's are involved, and support AbstractGrid arguments.
extension(::Type{T}, src::GridBasis, dest::GridBasis; options...) where {T} =
    gridrestriction(T, dest, src, grid(dest), grid(src); options...)'

restriction(::Type{T}, src::GridBasis, dest::GridBasis; options...) where {T} =
    gridrestriction(T, src, dest, grid(src), grid(dest); options...)

function gridrestriction(::Type{T}, src::Dictionary, dest::Dictionary, src_grid::G, dest_grid::GridArrays.MaskedGrid{G,M,I,S}; options...) where {T,G<:AbstractGrid,M,I,S}
    @assert supergrid(dest_grid) == src_grid
    IndexRestriction{T}(src, dest, subindices(dest_grid))
end

hasextension(dg::GridBasis{T,G}) where {T,G <: GridArrays.SubGrid} = true
hasextension(dg::GridBasis{T,G}) where {T,G <: GridArrays.ProductSubGrid} = true

###############################
# Converting to grids and back
###############################

# Convert to and from grids.
# To grid: we invoke `evaluation`
# From grid: we invoke `approximation`. This defaults to inverting the corresponding
#  `evaluation` operator.

conversion(::Type{T}, src::Dictionary, dest::GridBasis; options...) where {T} =
    evaluation(T, src, dest; options...)

conversion(::Type{T}, src::GridBasis, dest::Dictionary; options...) where {T} =
    approximation(T, src, dest; options...)

# Resolve ambiguity by the above methods
conversion(::Type{T}, src::GridBasis, dest::GridBasis; options...) where {T} =
    extension_restriction(T, src, dest; options...)

evaluation(::Type{T}, src::Dictionary, dest::GridBasis; options...) where {T} =
    evaluation(T, src, dest, grid(dest); options...)


text:
export FreeVariables, NewTypeWithVar, getbinding, bind!, testbinding!

FreeVariables() = FreeVariables(Dict{Symbol,Any}())

# Working with the symbols of type variables
"""
    `push!(::FreeVariables, ::Symbol)`

Add a new type variable to stack. The variable has no value (is unbound).
"""
function Base.push!(env::FreeVariables, s::NewTypeVar)
    a = get!(env.d, s) do
        []
    end
    push!(a, nothing)
    env
end

"""
    `pop!(::FreeVariables, ::Symbol)`

Remove variable from stack of free variables.
"""
function Base.pop!(env::FreeVariables, s::NewTypeVar)
    a = getindex(env.d, s)
    ae = pop!(a)
    if isempty(a)
        delete!(env.d, s)
    end
    ae
end

function binding(env::FreeVariables, s::NewTypeVar)
    a = getindex(env.d, s)
    a[end]
end

function bind!(env::FreeVariables, s::NewTypeVar, t::Any)
    a = getindex(env.d, s)
    a[end] = t
end

"""
    `testbinding!(env::FreeVariables, s::Symbol, t::newSymbol)`

If symbol is unbound, establish binding with new `Symbol`.
If symbol is bound to other `Symbol`, verify the new binding is identical.
If symbol is bound to Type, verify the new binding is same type.
If symbol is bound to value, verify the values are identical.

usage:
```
    testbinding(env, :T, parameter) || return false
```
"""
function testbinding!(env::FreeVariables, s::NewTypeVar, t::NewType)
    # println("testbinding-type $s $t")
    # map(println, stacktrace())
    a = getindex(env.d, s)
    ae = a[end]
    r = if ae == nothing
        a[end] = t
        true
    else
        issubenv(t, ae, env)
    end
end

function testbinding!(env::FreeVariables, s::NewTypeVar, t::Any)
    # println("testbinding-any $env $s $t")
    a = getindex(env.d, s)
    ae = a[end]
    r = if ae == nothing
        a[end] = t
        true
    else
        ae == t
    end
end

"""
    `freevariables(a) -> Vector{Symbol}`

Return list of free variables in a type expression.
obsolete
"""
freevariables(a) = Symbol[]
function freevariables(a::NewUnionAll)
    x = filter(x->x!=a.var.name, freevariables(a.body))
    union(x, freevariables(a.var.ub), freevariables(a.var.lb))
end
freevariables(a::NewUnion) = union(freevariables(a.a), freevariables(a.b))
freevariables(a::NewTypeVar) = union(freevariables(a.ub), freevariables(a.lb), [a.var])
freevariables(a::NewDataType) = mapfoldl(freevariables, union, Symbol[], a.parameters)


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

using Markdown
using InteractiveUtils

# ╔═╡ 800da77a-88e9-11eb-2be8-b1eb7e3efc9f
using LinearAlgebra

# ╔═╡ 8bc33e90-88ee-11eb-221b-f7d34c0d7944
using Yao

# ╔═╡ cf7d5a10-88f3-11eb-1e8c-b521a03f976b
using Gaston; set(showable="svg")

# ╔═╡ 9e892f39-f923-4462-8dd4-1c30edff4fbc
md"""
Classical probability distributions can be written as a stochastic vector, which can be transformed to another stochastic vector by applying a stochastic matrix. In other words, the evolution of stochastic vectors can be described by a stochastic matrix.

Quantum states also evolve and their evolution is described by unitary matrices. This leads to some interesting properties in quantum computing. Unitary evolution is true for a closed system, that is, a quantum system perfectly isolated from the environment. This is not the case in the quantum computers we have today: these are open quantum systems that evolve differently due to to uncontrolled interactions with the environment. In this notebook, we take a glimpse at both types of evolution.

# Unitary evolution

A unitary matrix has the property that its conjugate transpose is its inverse. Formally, it means that a matrix $U$ is unitary if $UU^\dagger=U^\dagger U=\mathbb{1}$, where $^\dagger$ stands for conjugate transpose, and $\mathbb{1}$ is the identity matrix. A quantum computer is a machine that implements unitary operations.

As an example, we have seen the NOT operation before, which is performed by the X gate in a quantum computer. While the generic discussion on gates will only occur in a subsequent notebook, we can study the properties of the X gate. Its matrix representation is $X = \begin{bmatrix} 0 & 1\\ 1 & 0\end{bmatrix}$. Let's check if it is indeed unitary:
"""

# ╔═╡ f686db5c-88e8-11eb-272d-c9d97ade768f
contran(x) = conj(x'ᵀ)

# ╔═╡ c8090a5d-b0a7-429d-a04b-5fd863afdcb0
var_X = [0 1; 1 0]

# ╔═╡ 303eff0a-88e9-11eb-03a2-b7ba8bc0081f
md" $XX^†$"

# ╔═╡ 3c70fbca-88e9-11eb-095e-89ad0cc703fe
var_X * contran(var_X)

# ╔═╡ 4fd650fc-88e9-11eb-08e2-b9e8afbbe6af
md" $X^†X$ "

# ╔═╡ 5baa00ea-88e9-11eb-1604-11b9709e1279
contran(var_X) * var_X

# ╔═╡ f23dc244-4b63-4119-b87e-b8dda57f3a65
md"""
It looks like a legitimate unitary operation. The unitary nature ensures that the $l_2$ norm is preserved, that is, quantum states are mapped to quantum states.
"""

# ╔═╡ 89337064-88e9-11eb-2701-2d1d9ee794af
let
	zero_ket = hcat([1, 0])
Text("The norm of the state |0〉 before applying X: $(norm(zero_ket))

The norm of the state after applying X: $(norm(var_X * zero_ket))")
end

# ╔═╡ 497e5fd2-0663-48be-b7e4-23a31111903d
md"""
Furthermore, since the unitary operation is a matrix, it is linear. Measurements are also represented by matrices. These two observations imply that everything a quantum computer implements is actually linear. If we want to see some form of nonlinearity, that must involve some classical intervention.

Another consequence of the unitary operations is reversibility. Any unitary operation can be reversed. Quantum computing libraries often provide a function to reverse entire circuits. Reversing the X gate is simple: we just apply it again (its conjugate transpose is itself, therefore $X^2=\mathbb{1}$).
"""

# ╔═╡ 6bb85a60-88ec-11eb-1e3b-bf918a16f0c6
let
	q = zero_state(1)
	q = q |> chain(1, put(1=>X), put(1=>X))
	state(q)
end

# ╔═╡ 68590691-d531-457f-82b2-fea20e4efacb
md"""
which is exactly $|0\rangle$ as we would expect.

In the next notebook, you will learn about classical and quantum many-body systems and the Hamiltonian. In the notebook on adiabatic quantum computing, you will learn that a unitary operation is in fact the Schrödinger equation solved for a Hamiltonian for some duration of time. This connects the computer science way of thinking about gates and unitary operations to actual physics, but there is some learning to be done before we can make that connection. Before that, let us take another look at the interaction with the environment.
"""

# ╔═╡ a29e2c11-cae1-4436-b089-445017b23ff2
md"""
# Interaction with the environment: open systems

Actual quantum systems are seldom closed: they constantly interact with their environment in a largely uncontrolled fashion, which causes them to lose coherence. This is true for current and near-term quantum computers too.
"""

# ╔═╡ 85bffd10-88f0-11eb-2348-ab467cb4dfcd
begin
	import PlutoUI
	PlutoUI.LocalResource("../figures/blochsphere.jpg", :width => 600)
end

# ╔═╡ 7b38b2ec-88f0-11eb-2d8a-67243b44b511
md"""
This also means that their actual time evolution is not described by a unitary matrix as we would want it, but some other operator (the technical name for it is a completely positive trace-preserving map).

Quantum computing libraries often offer a variety of noise models that mimic different types of interaction, and increasing the strength of the interaction with the environment leads to faster decoherence. The timescale for decoherence is often called `` T_2 `` time. Among a couple of other parameters, `` T_2 `` time is critically important for the number of gates or the duration of the quantum computation we can perform.

A very cheap way of studying the effects of decoherence is mixing a pure state with the maximally mixed state ``\mathbb{1}/2^d``, where ``d`` is the number of qubits, with some visibility parameter in ``[0,1]``. This way we do not have to specify noise models or any other map modelling decoherence. For instance, we can mix the ``|\phi^+\rangle`` state with the maximally mixed state:
"""

# ╔═╡ bc26371a-88f1-11eb-0224-8124986b87ac
eye(x) = Float64.(Matrix(I, x, x))

# ╔═╡ e08d0c68-7a53-4303-9a8d-3b4e6cb0acc5
function mixed_state(pure_state, visibility)
    density_matrix = pure_state * contran(pure_state)
    maximally_mixed_state = eye(4)/(2^2)
    return (visibility*density_matrix + (1-visibility)*maximally_mixed_state)
end

# ╔═╡ 76592336-88f2-11eb-0ef2-1fbd3e5aa832
begin
ϕ = hcat([1, 0, 0, 1])/√2
Text("Maximum visibility is a pure state:")
end

# ╔═╡ 3298ceca-88f3-11eb-282c-53f638367451
mixed_state(ϕ, 1.0)

# ╔═╡ 3d652470-88f3-11eb-25d4-279682ede9f1
Text("The state is still entangled with visibility 0.8:")

# ╔═╡ 5f03e8c8-88f3-11eb-1041-c3e726375d98
mixed_state(ϕ, 0.8)

# ╔═╡ 6676ac1c-88f3-11eb-05bf-df507e5550b9
Text("Entanglement is lost by 0.6:")

# ╔═╡ 759c4418-88f3-11eb-1fa8-0df8c27adda0
mixed_state(ϕ, 0.6)

# ╔═╡ 7cf925fa-88f3-11eb-13a6-2d753e1475c4
Text("Barely any coherence remains by 0.2:")

# ╔═╡ 8bd2d9c4-88f3-11eb-364c-b7a76fdbed2e
mixed_state(ϕ, 0.2)

# ╔═╡ f1a12107-46e0-48b2-a58d-4ecf52a00c66
md"""
Another way to look at what happens to a quantum state in an open system is through equilibrium processes. Think of a cup of coffee: left alone, it will equilibrate with the environment, eventually reaching the temperature of the environment. This includes energy exchange. A quantum state does the same thing and the environment has a defined temperature, just like the environment of a cup of coffee.

The equilibrium state is called the thermal state. It has a very specific structure and we will revisit it, but for now, suffice to say that the energy of the samples pulled out of a thermal state follows a Boltzmann distribution. The Boltzmann -- also called Gibbs -- distribution is described as $P(E_i) = \frac {e^{-E_{i}/T}}{\sum _{j=1}^{M}{e^{-E_{j}/T}}}$, where $E_i$ is an energy, and $M$ is the total number of possible energy levels. Temperature enters the definition: the higher the temperature, the closer we are to the uniform distribution. In the infinite temperature limit, it recovers the uniform distribution. At high temperatures, all energy levels have an equal probability. In contrast, at zero temperature, the entire probability mass is concentrated on the lowest energy level, the ground state energy. To get a sense of this, let's plot the Boltzmann distribution with vastly different temperatures:
"""

# ╔═╡ a3067d19-30dd-4969-9bbf-7531acf4d7fd
begin
	temperatures = [0.5, 5, 2000]
	energies = range(0,20; length=100)

	p = let t = exp.(-energies ./ temperatures[end])
		plot([1],[1], w=:p, Axes(xrange=(0,20), yrange=(0,1.2*maximum(t)/sum(t)), xlabel="'Energy'", ylabel="'Probability'", key="box width 3")) end

	for (i,T) = enumerate(temperatures)
    	probabilities = exp.(-energies ./ T)
    	probabilities /= sum(probabilities)
    	p = plot!(energies, probabilities, lw=3, legend = "'T_$(i)'")
	end
	p
end

# ╔═╡ 3864cae0-7170-44bc-bd8f-b6118316a1b0
md"""
Here $T_1<T_2<T_3$. Notice that $T_1$ is a low temperature, and therefore it is highly peaked at low energy levels. In contrast, $T_3$ is a very high temperature and the probability distribution is almost completely flat.
"""

# ╔═╡ Cell order:
# ╟─9e892f39-f923-4462-8dd4-1c30edff4fbc
# ╠═f686db5c-88e8-11eb-272d-c9d97ade768f
# ╠═c8090a5d-b0a7-429d-a04b-5fd863afdcb0
# ╟─303eff0a-88e9-11eb-03a2-b7ba8bc0081f
# ╠═3c70fbca-88e9-11eb-095e-89ad0cc703fe
# ╟─4fd650fc-88e9-11eb-08e2-b9e8afbbe6af
# ╠═5baa00ea-88e9-11eb-1604-11b9709e1279
# ╟─f23dc244-4b63-4119-b87e-b8dda57f3a65
# ╠═800da77a-88e9-11eb-2be8-b1eb7e3efc9f
# ╠═89337064-88e9-11eb-2701-2d1d9ee794af
# ╟─497e5fd2-0663-48be-b7e4-23a31111903d
# ╠═8bc33e90-88ee-11eb-221b-f7d34c0d7944
# ╠═6bb85a60-88ec-11eb-1e3b-bf918a16f0c6
# ╟─68590691-d531-457f-82b2-fea20e4efacb
# ╟─a29e2c11-cae1-4436-b089-445017b23ff2
# ╟─85bffd10-88f0-11eb-2348-ab467cb4dfcd
# ╟─7b38b2ec-88f0-11eb-2d8a-67243b44b511
# ╠═bc26371a-88f1-11eb-0224-8124986b87ac
# ╠═e08d0c68-7a53-4303-9a8d-3b4e6cb0acc5
# ╠═76592336-88f2-11eb-0ef2-1fbd3e5aa832
# ╠═3298ceca-88f3-11eb-282c-53f638367451
# ╟─3d652470-88f3-11eb-25d4-279682ede9f1
# ╠═5f03e8c8-88f3-11eb-1041-c3e726375d98
# ╟─6676ac1c-88f3-11eb-05bf-df507e5550b9
# ╠═759c4418-88f3-11eb-1fa8-0df8c27adda0
# ╟─7cf925fa-88f3-11eb-13a6-2d753e1475c4
# ╠═8bd2d9c4-88f3-11eb-364c-b7a76fdbed2e
# ╟─f1a12107-46e0-48b2-a58d-4ecf52a00c66
# ╠═cf7d5a10-88f3-11eb-1e8c-b521a03f976b
# ╠═a3067d19-30dd-4969-9bbf-7531acf4d7fd
# ╟─3864cae0-7170-44bc-bd8f-b6118316a1b0

max_stars_repo_name: VarLad/qml-mooc

max_issues_repo_name: VarLad/qml-mooc

max_forks_repo_name: VarLad/qml-mooc