text:
# Nearest-neighbor hopping parameter
const NN_hopping = 2.8

# Graphene basis vectors
# Distance between neighboring carbon atoms
const graphene_lattice_constant = 2.46
const sublattice_shift = -1 / √(3) * graphene_lattice_constant

const graphene_d1 =
    [graphene_lattice_constant / 2, graphene_lattice_constant * √(3) / 2]

const graphene_d2 =
    [-graphene_lattice_constant / 2, graphene_lattice_constant * √(3) / 2]

#  Algebraic data type for graphene sublattices
@data Sublattice begin
    A
    B
end

"""
    GrapheneCoord(u::Int, v::Int, sublattice::Sublattice)

Lattice coordinate of a carbon atom, generated using [`graphene_A`](@ref) or
[`graphene_B`](@ref).

Each coordinate contains the sublattice index `A` or `B`, as well as
the integer coefficients of the two basis vectors
``d\\times(\\pm 1 \\hat{x} + \\sqrt{3}\\hat{y}) / 2``
(`u` for `+`, `v` for `-`), with ``d = 2.46``Å as the lattice constant.
"""
struct GrapheneCoord
    u::Int
    v::Int
    sublattice::Sublattice
end

# Define a comparison function between the coordinates to help with testing
function Base.isless(a1::GrapheneCoord, a2::GrapheneCoord)
    if a1.sublattice == A && a2.sublattice == B
        return true
    elseif a1.sublattice == B && a2.sublattice == A
        return false
    else
        if a1.v == a2.v
            return isless(a1.u, a2.u)
        else
            return isless(a1.v, a2.v)
        end
    end
end

"""
    graphene_A(u::Int, v::Int)

Create a [`GrapheneCoord`](@ref) for an atom belonging to sublattice A at the
unit cell (u, v).
"""
function graphene_A(u::Int, v::Int)
    return GrapheneCoord(u, v, A)
end

"""
    graphene_B(u::Int, v::Int)

Create a [`GrapheneCoord`](@ref) for an atom belonging to sublattice B at the
unit cell (u, v).
"""
function graphene_B(u::Int, v::Int)
    return GrapheneCoord(u, v, B)
end

"""
    graphene_neighbors(atom::GrapheneCoord)

Determine the nearest neighbors of an `atom` and return a vector of the
correspoding [`GrapheneCoord`](@ref)'s.
"""
function graphene_neighbors(atom::GrapheneCoord)
    u = atom.u
    v = atom.v
    if atom.sublattice == A
        return [
            graphene_B(u, v)
            graphene_B(u + 1, v)
            graphene_B(u, v + 1)
        ]
    elseif atom.sublattice == B
        return [
            graphene_A(u, v)
            graphene_A(u - 1, v)
            graphene_A(u, v - 1)
        ]
    else
        error("Illegal sublattice parameter")
    end
end

## Propagator
# Integrals used in computing the propagator

# When computing Ω, occasionally the integrand becomes small enough to give NaN
# The integrand helper functions help to catch these instances
@inline function Ω_Integrand(z, u, v, x::Float64)
    W = ((z / NN_hopping)^2 - 1.0) / (4.0 * cos(x)) - cos(x) |> Complex
    res =
        exp(1.0im * (u - v) * x) / cos(x) *
        ((W - √(W - 1) * √(W + 1))^abs.(u + v)) / (√(W - 1) * √(W + 1))
    return (isnan(res) ? 0.0 + 0.0im : res)
end

@inline function Ω(z, u, v)
    return ((quadgk(
        x -> 2 * Ω_Integrand(z, u, v, x) / (8.0 * π * NN_hopping^2),
        0.0,
        π,
        atol = 1e-16,
        rtol = 1e-4,
    ))[1])
end

@inline function Ωp_Integrand(z, u, v, x::Float64)
    W = ((z / NN_hopping)^2 - 1.0) / (4.0 * cos(x)) - cos(x) |> Complex
    res =
        2 *
        exp(1.0im * (u - v) * x) *
        ((W - √(W - 1) * √(W + 1))^abs.(u + v + 1)) / (√(W - 1) * √(W + 1))
    return (isnan(res) ? 0.0 + 0.0im : res)
end

@inline function Ωp(z, u, v)
    return ((quadgk(
        x -> 2 * Ωp_Integrand(z, u, v, x) / (8.0 * π * NN_hopping^2),
        0.0,
        π,
        atol = 1e-16,
        rtol = 1e-4,
    ))[1])
end

@inline function Ωn_Integrand(z, u, v, x::Float64)
    W = ((z / NN_hopping)^2 - 1.0) / (4.0 * cos(x)) - cos(x) |> Complex
    res =
        2 *
        exp(1.0im * (u - v) * x) *
        ((W - √(W - 1) * √(W + 1))^abs.(u + v - 1)) / (√(W - 1) * √(W + 1))
    return (isnan(res) ? 0.0 + 0.0im : res)
end

@inline function Ωn(z, u, v)
    return ((quadgk(
        x -> 2 * Ωn_Integrand(z, u, v, x) / (8.0 * π * NN_hopping^2),
        0.0,
        π,
        atol = 1e-16,
        rtol = 1e-4,
    ))[1])
end

# The propagator function picks out the correct element of the Ξ matrix based
# on the sublattices of the graphene coordinates.
function graphene_propagator(a_l::GrapheneCoord, a_m::GrapheneCoord, z)
    t = NN_hopping
    u = a_l.u - a_m.u
    v = a_l.v - a_m.v
    if a_l.sublattice == a_m.sublattice
        return (z * Ω(z, u, v))
    elseif ([a_l.sublattice, a_m.sublattice] == [A, B])
        return (-t * (Ω(z, u, v) + Ωp(z, u, v)))
    elseif ([a_l.sublattice, a_m.sublattice] == [B, A])
        return (-t * (Ω(z, u, v) + Ωn(z, u, v)))
    else
        error("Illegal sublattice parameter")
    end
end

"""
    crystal_to_cartesian(coord::GrapheneCoord)

Convert a [`GrapheneCoord`](@ref) to a cartesian point with lengths in Å, where
`GrapheneCoord(0, 0, A)` is at the origin.
"""
function crystal_to_cartesian(coord::GrapheneCoord)
    u = coord.u
    v = coord.v
    x = graphene_d1[1] * u + graphene_d2[1] * v
    y = graphene_d1[2] * u + graphene_d2[2] * v

    return ([x, y + (coord.sublattice == B) * sublattice_shift])
end

# Given a list of [`GrapheneCoord`](@ref), this functiohn returns a Ξ(z) propagator
# matrix. The calculation is sped up using the fact that the matrix is symmetric.
function propagator_matrix(z, Coords::Vector{GrapheneCoord})
    precomputed = Dict{Float64,ComplexF64}()
    len_coords = length(Coords)
    out = zeros(ComplexF64, len_coords, len_coords)
    for ii = 1:len_coords
        @inbounds for jj = ii:len_coords
            dist = norm(
                crystal_to_cartesian(Coords[ii]) -
                crystal_to_cartesian(Coords[jj]),
            )
            c = get(precomputed, dist, 0.0)
            if c == 0.0
                res = graphene_propagator(Coords[ii], Coords[jj], z)
                out[ii, jj] = res
                out[jj, ii] = res
                precomputed[dist] = res
            else
                out[ii, jj] = c
                out[jj, ii] = c
            end
        end
    end
    return out
end


text:
module THICK_plate_2dir_strong_examples
using FinEtools
using FinEtools.AlgoDeforLinearModule
using FinEtools.MeshUtilModule
using FinEtools.AlgoBaseModule

function THICK_plate_2dir_strong_MST10_conv()
    elementtag = "MST10"
    println("""
    Fiber-reinforced cantilever plate: $(elementtag)
    """)

    # This example provides three-dimensional finite element model for  the
    # transverse shear stress calculations. The problem consists of a one-, two- or
    # three-layer plate subjected to a sinusoidal distributed load, as
    # described by Pagano (1969). The resulting transverse shear and axial
    # stresses through the thickness of the plate are compared to two existing
    # analytical solutions by Pagano (1969). The first solution is derived from
    # classical laminated plate theory (CPT), while the second is an exact
    # solution from linear elasticity theory.

    for (extrap, nodevalmeth) = zip([:extrapmean, :extraptrend, :default], [:averaging, :averaging, :invdistance])
        filebase = "THICK_plate_2dir_strong_MST10_conv_$(elementtag)_$(extrap)"
        modeldatasequence = FDataDict[]
        for Refinement = [1, 2, 4, 8]

            # Orthotropic material
            E1 = 1000.0e9*phun("Pa"); E2 = 1000.0e9*phun("Pa"); E3 = 1.0e9*phun("Pa");
            G12 = 0.2e9*phun("Pa");  G13 = G12; G23 = 0.2e9*phun("Pa")
            nu12 = nu13 = nu23 =  0.25;
            # dimensions of the plate
            a = 70.0*phun("mm");
            b = 100.0*phun("mm");
            t = 50.0*phun("mm");
            # Transverse loading
            q0 = 1000.0*phun("Pa")
            # Coefficients of thermal expansion
            CTE1 = CTE2 = CTE3 = 0.0

            # Here we define the layout and the thicknesses of the layers.
            angles = vec([45.0]);
            nLayers = length(angles)
            ts = t/nLayers * ones(nLayers); # layer thicknesses

            tolerance = 0.0001*t

            # Select how fine the mesh should be
            nts= Refinement * ones(Int, nLayers);# number of elements per layer
            tnts = sum(nts)
            na, nb = 4 * tnts, 4 * tnts;

            xs = collect(linspace(0.0, a, na+1))
            ys = collect(linspace(0.0, b, nb+1))
            fens,fes = T10layeredplatex(xs, ys, ts, nts)
            println("Mesh: na, nb, nts = $na, $nb, $nts")
            println("count(fens) = $(count(fens))")

            # This is the material  model
            MR = DeforModelRed3D
            material = MatDeforElastOrtho(MR, 0.0, E1, E2, E3, nu12, nu13, nu23, G12, G13, G23, CTE1, CTE2, CTE3)

            # The material coordinate system function is defined as:
            function updatecs!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
                rotmat3!(csmatout, angles[fe_label]/180.0*pi* [0.0; 0.0; 1.0]);
            end

            # The volume integrals are evaluated using this rule
            gr = SimplexRule(3, 4)

            # We will create 3 regions, one for each of the layers
            regions = FDataDict[]
            for layer = 1:nLayers
                rls = selectelem(fens, fes, label =  layer)
                push!(regions, FDataDict("femm"=>FEMMDeforLinearMST10(MR, IntegData(subset(fes, rls), gr), CSys(3, 3, updatecs!), material)))
            end

            # The essential boundary conditions: clamped face
            lx0 = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)
            eclamped1 = FDataDict("displacement"=>  0.0, "component"=> 1, "node_list"=>lx0)
            eclamped2 = FDataDict("displacement"=>  0.0, "component"=> 2, "node_list"=>lx0)
            eclamped3 = FDataDict("displacement"=>  0.0, "component"=> 3, "node_list"=>lx0)

            # The traction boundary condition is applied at the free face opposite the clamped face.
            bfes = meshboundary(fes)
            function pfun(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
                forceout[1] = 0.0
                forceout[2] = 0.0
                forceout[3] = -q0
                return forceout
            end
            # From  the entire boundary we select those quadrilaterals that lie on the plane
            # Z = thickness
            tl = selectelem(fens, bfes, box = [a a -Inf Inf -Inf Inf], inflate=tolerance)
            Trac = FDataDict("traction_vector"=>pfun, "femm"=>FEMMBase(IntegData(subset(bfes, tl), SimplexRule(2, 3))))

            modeldata = FDataDict("fens"=>fens, "regions"=>regions, "essential_bcs"=>[eclamped1, eclamped2, eclamped3],  "traction_bcs"=> [Trac])
            modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

            # modeldata["postprocessing"] = FDataDict("file"=>filebase * "-u")
            # modeldata = AlgoDeforLinearModule.exportdeformation(modeldata)

            u = modeldata["u"]
            geom = modeldata["geom"]

            # The results of the displacement and stresses will be reported at
            # nodes located at the appropriate points.
            nbottomedge = selectnode(fens, box=[a a 0.0 0.0 0.0 0.0], inflate=tolerance)

            println("bottom edge deflection: $(mean(u.values[nbottomedge, 3], 1)/phun("mm")) [mm]")

            #  Compute  all stresses
            modeldata["postprocessing"] = FDataDict("file"=>filebase * "-s", "quantity"=>:Cauchy, "component"=>collect(1:6), "outputcsys"=>CSys(3), "nodevalmethod"=>nodevalmeth, "reportat"=>extrap)
            modeldata = AlgoDeforLinearModule.exportstress(modeldata)

            modeldata["elementsize"] = t/Refinement
            modeldata["geometricaltolerance"] = tolerance
            push!(modeldatasequence, modeldata)
        end # for refinement

        println("")
        println("Stress RMS error")
        for md = modeldatasequence
            md["targetfields"] = [e["field"] for e in md["postprocessing"]["exported"]]
        end
        elementsizes, errornorms, p = AlgoBaseModule.evalconvergencestudy(modeldatasequence)

        println("Normalized Approximate Error = $(errornorms)")

        f = log.(vec(errornorms))
        A = hcat(log.(vec(elementsizes[1:end-1])), ones(size(f)))
        p = A \ f
        println("Linear log-log fit: p = $(p)")

        csvFile = filebase * "_Stress" * ".CSV"
        savecsv(csvFile,
        elementsizes=vec(elementsizes[1:end-1]),
        elementsizes2=vec(elementsizes[1:end-1].^2),
        elementsizes3=vec(elementsizes[1:end-1].^3),
        errornorms=vec(errornorms)
        )
        println("Wrote $csvFile")

        println("")
        println("Displacement RMS error")
        for md = modeldatasequence
            md["targetfields"] = [md["u"] for r in md["regions"]]
        end
        elementsizes, errornorms, p = AlgoBaseModule.evalconvergencestudy(modeldatasequence)

        println("Normalized Approximate Error = $(errornorms)")

        f = log.(vec(errornorms))
        A = hcat(log.(vec(elementsizes[1:end-1])), ones(size(f)))
        p = A \ f
        println("Linear log-log fit: p = $(p)")

        csvFile = filebase * "_Displ" * ".CSV"
        savecsv(csvFile,
        elementsizes=vec(elementsizes[1:end-1]),
        elementsizes2=vec(elementsizes[1:end-1].^2),
        elementsizes3=vec(elementsizes[1:end-1].^3),
        errornorms=vec(errornorms)
        )
        println("Wrote $csvFile")

        # @async run(`"paraview.exe" $csvFile`)
    end # extrap
    println("Done")

end # Pagano_3lay_cyl_bend_MST10_conv

function allrun()
    println("#####################################################")
    println("# THICK_plate_2dir_strong_MST10_conv ")
    THICK_plate_2dir_strong_MST10_conv()
    return true
end # function allrun

end # module Pagano_examples

ext: jl


text:
using ApproxFunOrthogonalPolynomials, Test

c = rand(1000)
x=rand(10000)
f=Fun(Chebyshev,c)
y=f.(x)
y=f.(x)

@time y=f.(x)
println("Clenshaw large coeffs, many points: Time should be ~0.024")
# 0.012482274  with unsafe_view
# 0.024306262 with inbounds

y=f(.1)
y=f(.1)
y=f(.1)

@time y=f(.1);
println("Clenshaw large coeffs, 1 point: Time should be ~6e-6")

# @time is 8.853e-6 seconds

f=Fun(x->cos(x),20)
roots(f)
roots(f)
@time for k=1:100
    roots(f)
end
println("Small roots: Time should be ~0.015")

f=Fun(x->cos(1000x),1000)
roots(f)
roots(f)
@time roots(f)
println("Roots: Time should be ~0.08")


text:
function FHN(du,u,p,t) # FitzHugh-Nagumo oscillator
    du[1] = 30*u[1] - 10*u[1]^3 - 30*u[2]
    du[2] = 7.5*u[1]
end

function spiral(du,u,p,t) # Spiral
    true_A = [-0.1 -2.0 0; 2.0 -0.1 0; 0 0 -0.3]
    du .= true_A * (u .^ 3 + u)
end

function NDM(du,u,p,t) # Mutual Inhibition
    du[1] = 10 * (-u[1] + 1 ./ ( 1 + exp( 16 * ( u[2] - 0.5 ) ) ))
    du[2] = 10 * (-u[2] + 1 ./ ( 1 + exp( 16 * ( u[1] - 0.5 ) ) ))
end

max_stars_repo_name: timkimd/plnde

max_issues_repo_name: timkimd/plnde

max_forks_repo_name: timkimd/plnde


text:
using ApproxFun, SingularIntegralEquations, HierarchicalSingularIntegralEquations

include("HierarchicalVectorTest.jl")
include("hierarchicalsolveTest.jl")
include("FractalTest.jl")


text:
function revers(side :: HorizonSide)
    side = HorizonSide((Int(side)+2)%4)
    return side
end

function movements!(r, count, side, out_side)

    flag = true
    for i = 1:count
        if isborder(r, out_side) == false
            flag = false
            return flag
        end
        move!(r,side)
    end
end

function find_out(r)

    count = 1
    flag = true

    out_side :: HorizonSide = Nord
    side :: HorizonSide = Ost

    while flag != false

        flag = movements!(r,count, side, out_side)

        side = revers(side)

        count = count + 1

    end;

end

find_out(r)
show(r)

max_stars_repo_name: AntonDudDm/mirea-progs

max_issues_repo_name: AntonDudDm/mirea-progs

max_forks_repo_name: AntonDudDm/mirea-progs


text:
##########
## ITMIX (of all versions loaders)
#########
#

#####
# ITMIX data
#####
"Translate part of the file-name to DataKind subtypes"
const dk_dict = Dict("RES"=>ThicknessData,
                     "margin"=>OutlineData,
                     "surface"=>DEMData,
                     "mb"=>BdotData,
                     "dhdt"=>DhdtData,
                     "speed"=>IVData,
                     "azimuth"=>nothing,
                     "velocity"=>IVData)

todo = 1
"Whether a glacier is calving, is an ice cap, and its RGI 4.0 number (if there is a single one)."
const iscalving_isicecap = Dict(:Academy=>(true,true,nothing),
                                :Austfonna=>(true,true,nothing),
                                :Devon=>(true,true,"RGI40-03.02443"), # one of the major outlet glaciers in the North
                                :Hellstugubreen=>(false,false,"RGI40-08.02182"),
                                :Mocho=>(false,false,nothing), # glacier complex
                                :SouthGlacier=>(false,false,"RGI40-01.17598"),
                                :Tasman=>(false,false,"RGI40-18.02342"),
                                :Washmawapta=>(false,false,"RGI40-02.04894"),
                                :Aqqutikitsoq=>(false,false,"RGI40-05.00524"), # not sure, size does not match well
                                :Brewster=>(false,false,"RGI40-18.01130"),
                                :Elbrus=>(false,false,nothing), # glacier complex
                                :NorthGlacier=>(false,false,"RGI40-01.18240"),
                                :Starbuck=>(true,false,nothing),
                                :Unteraar=>(false,false,"RGI40-11.01328"),
                                :Columbia=>(true,false,"RGI40-01.11506"),
                                :Freya=>(false,false,"RGI40-05.19983"),
                                :Kesselwandferner=>(false,false,"RGI40-11.00787"),
                                :Urumqi=>(false,false,"RGI40-13.27054"),
                                :Synthetic1=>(false,false,nothing),
                                :Synthetic2=>(false,false,nothing),
                                :Synthetic3=>(false,true,nothing),
                                :AustreGroenfjordbreen=>(false,false,nothing),
                                :ChhotaShigri=>(false,false,nothing),
                                )

# ITMIX glacier with IV data:
const itmix_glaciers_iv = [:Synthetic1,  :Unteraar,
                           :Synthetic2, :Synthetic3,
                           :Tasman, :Austfonna, :Devon,
                           :SouthGlacier, :NorthGlacier,
                           :Brewster, :Hellstugubreen]

# of which these only have point data (at GPS locations)
const itmix_glaciers_iv_point = [:NorthGlacier, :SouthGlacier, :Hellstugubreen, :Brewster]

function LoadPara(gid::ITMIXGlacier; use_glogem=false, pl_kwargs...)
    pl_kwargs = deepcopy(pl_kwargs)
    id = getid(gid)
    data_root = gid.version==1? joinpath(LoadPara(gid=gid).data_root,"ITMIX") : joinpath(LoadPara(gid=gid).data_root,"ITMIX2")
    pl = LoadPara(; gid=gid, data_root=data_root, pl_kwargs...)
    # set the sigmas of uncertainty for the data:
    for (k,v) in Dict(ThicknessData=>10.0, # maybe change to more
                      BdotData=>0.2,  # m/a, 1m/a is what Matthias says about GloGEM
                      DhdtData=>0.2,  # m/a, e.g. 1m/a Nuth et al. 2010, Table 1, for Svalbard
                      IVData=>10.0)
        pl.dataset_opts[k][:sigma]=v
    end
    updated_pm_defaults = Dict(
        :bandsize => 30,
        :window_width_smooth=>100, # can make a noisy tau
        :cap_at_floatation=>false,
        :dist_exp => 0.5,
        #:window => 0
    )
    if iscalving_isicecap[id][2]
        merge!(updated_pm_defaults, Dict(:h_at_top => -1, :shapeF => false))
    end
    if iscalving_isicecap[id][1]
        # relative flux lost through calving
        pl.dataset_opts[TerminusFluxData][:fluxkind] =
            FluxKind.rel_btilde_teminus # compare to prior expectation of btilde integral
        #pl.dataset_opts[TerminusFluxData][:terminus_flux] = SomeData(vals=1.0, sigma=0.5, range_abs=(0.1,3.0))
        pl.dataset_opts[TerminusFluxData][:terminus_flux] =
            SomeData(vals=NaN) # require this to be set individually
    else
        # approximately zero flux at terminus +/- 1% of max bdot flux accumulation
        pl.dataset_opts[TerminusFluxData][:fluxkind] = FluxKind.rel_bdot_max
        pl.dataset_opts[TerminusFluxData][:terminus_flux] = SomeData(vals=0.0, sigma=0.01)
    end

    if id==:Austfonna
        # This makes a coarser grid as the Austfonna data is on a
        # 50m grid!
        # with 50m DEM:
        #  - 400s without cache
        #  - 35s with cache
        # with 250m DEM:
        #  - 208s without cache
        #  - 35s  with cache
        if gid.version==1
            pl.dataset_opts[DataKind][:grid_downsample_step]=5
        end
        #pl.dataset_LOADERS[DhdtData] = SyntheticBenchLoader # something is off with Austfonna's dh/dt (or bdot)

        pl.dataset_opts[TerminusFluxData][:fluxkind] = FluxKind.rel_bdot_max # compare to prior expectation of bdot integral
        # Dowdeswell at al. 2008 says this:
        pl.dataset_opts[TerminusFluxData][:terminus_flux] = SomeData(vals=1/3,
                                                                     sigma=0.4,
                                                                     range_abs=(0.05, 3.0))
    elseif id==:AustreGroenfjordbreen
        # compare to prior expectation of btilde integral because rel_bdot_max is negative
        pl.dataset_opts[TerminusFluxData][:fluxkind] = FluxKind.rel_btilde_teminus
        pl.dataset_opts[TerminusFluxData][:terminus_flux] = SomeData(vals=0.0, sigma=0.01)
    elseif id==:Columbia
        # the size of (1040, 920) make it run slow, about 0.5s per forward run
        # reduced by step=2 (4x) it's down to 0.07, in-line with the others.
        pl.dataset_opts[DataKind][:grid_downsample_step] = 2

        pl.dataset_opts[BdotData][:smb_gradient] = 0.02 # maritime
        pl.dataset_opts[TerminusFluxData][:fluxkind] = FluxKind.rel_bdot_max # compare to prior expectation of max-bdot integral
                                                                             # Note, bdot intergral is negative at terminus
        pl.dataset_opts[TerminusFluxData][:terminus_flux] = SomeData(vals=0.5,
                                                                     sigma=0.3,
                                                                     range_abs=(0.1,3.0))
        merge!(updated_pm_defaults, Dict(:h_at_top => -1)) # more ice cap like at devide
    elseif id==:Aqqutikitsoq
        if gid.version==1
            pl.dataset_opts[DataKind][:grid_downsample_step]=2
        end
    elseif id==:Hellstugubreen
        # pl.dataset_LOADERS[DhdtData] = SyntheticBenchLoader
        # pl.dataset_opts[DhdtData][:ela_elevation] = 2000
    elseif id==:Elbrus
        if gid.version==1
            pl.dataset_opts[DataKind][:grid_downsample_step]=3
        else
            # the size of (601,581) make it run slow, about 0.12s per forward run
            # reduced by step=2 (4x) it's down to 0.025, in-line with the others.
            # (now on 60m resolution vs 30m resolution)
            pl.dataset_opts[DataKind][:grid_downsample_step] = 2
        end

        pl.dataset_LOADERS[BdotData] = SyntheticBenchLoader # original only covers a small fraction of the area
                                                            # TODO: this could be turned into elevation band data:
                                                             # pl.dataset_opts[BdotData][:vals_prior] = 0.0
        # pl.dataset_LOADERS[DhdtData] = SyntheticBenchLoader # doesn't make much sense to then not use this either.
        merge!(updated_pm_defaults, Dict(:h_at_top => -1, :shapeF => false))
    elseif id==:Unteraar
        # something is amiss with dhdt or bdot data as glacier stops at confluence:
        # pl.dataset_LOADERS[DhdtData] = SyntheticBenchLoader
        # pl.dataset_opts[DhdtData][:ela_elevation] = 2900
        pl.dataset_opts[IVData][:sigma] = 4.0 # Vogel et al 2012
        pl.dataset_opts[ThicknessData][:sigma] = 27.0 # Bauder etal 2003
    elseif id==:Freya
        # uses the SyntheticBenchLoader
        pl.dataset_opts[DhdtData][:ela_elevation] = 985 # gets negative fluxes in the middle otherwise
    elseif id==:Urumqi
        pl.dataset_opts[DhdtData][:ela_elevation] = 4200
    elseif id==:Tasman
        merge!(updated_pm_defaults, Dict(:bandsize => 35))
    elseif id==:Devon
        pl.dataset_opts[TerminusFluxData][:fluxkind] = FluxKind.rel_bdot_max # compare to prior expectation of bdot integral
        # https://www.igsoc.org/annals/46/a46a288.pdf says 1/4 calving and 3/4 melt
        pl.dataset_opts[TerminusFluxData][:terminus_flux] = SomeData(vals=1/4,
                                                                     sigma=0.1,
                                                                     range_abs=(0.1,3.0))
    elseif id==:Synthetic3
        merge!(updated_pm_defaults, Dict(:h_at_top => -1))
    elseif id==:Synthetic1 || id==:Synthetic2 || id==:Synthetic3
        # there are no "measurement" errors here as no noise was put on the data.
        # Note there is still the BITEModel-error in pm.sigma_*
        pl.dataset_opts[IVData][:sigma] = 0.0
        pl.dataset_opts[ThicknessData][:sigma] = 0.0
    elseif id==:ChhotaShigri
        # flux not zero at terminus
    elseif id==:Academy
        pl.dataset_opts[TerminusFluxData][:fluxkind] = FluxKind.rel_bdot_max
        # Dowdeswell at all 2008 says this:
        pl.dataset_opts[TerminusFluxData][:terminus_flux] = SomeData(vals=1/3,
                                                                     sigma=0.3,
                                                                     range_abs=(0.1,3.0))
    elseif id==:Mocho
        # has negative mass-balance at the top!
        merge!(updated_pm_defaults, Dict(:h_at_top => -1, :shapeF => false))
    elseif id==:Starbuck
        pl.dataset_opts[TerminusFluxData][:fluxkind] = FluxKind.rel_bdot_teminus # compare to prior expectation of bdot integral
        pl.dataset_opts[TerminusFluxData][:terminus_flux] = SomeData(vals=1.0,
                                                                     sigma=0.1,
                                                                     range_abs=(0.1,3.0))
        merge!(updated_pm_defaults, Dict(:h_at_top => -1))  # more ice cap like at devide
    end
    pl.dataset_opts[DEMData][:cropbox] = :from_DEM

    # update pm but keeping any explicitly set (via kwargs) options
    merge!((a,b)->a, pl.dataset_opts[ParaData][:pm], updated_pm_defaults)

    # GloGEM
    if use_glogem && iscalving_isicecap[gid.id][3]!=nothing
        error("glogem loader disabled for ITMIX2 runs")
        pl.dataset_LOADERS[BdotData] = GloGEMLoader
        pl.dataset_LOADERS[DhdtData] = SyntheticBenchLoader
    else
        if use_glogem
            error("Not using GloGEM for $gid")
        end
    end

    # Set loaders and files:
    @unpack dataset_LOADERS, dataset_files, dataset_opts = pl
    # check what files and thus data exist
    name = getname(gid)
    dir = joinpath(pl.data_root, name)
    zip = joinpath(pl.data_root, name*".zip")
    fls = Dict()
    if isfile(zip) # zip-archive
        zipfl = ZipFile.Reader(zip)
        toread = zipfl.files
    elseif isdir(dir)
        toread = readdir(dir)
        append!(toread, joinpath.("shapefiles", readdir(joinpath(dir, "shapefiles"))))
    else
        error("Cannot read $dir or $zip")
    end
    for fl in toread
        out = parse_ITMIX_filename(fl)
        if out!=nothing
            fls[out[1]] = out[2:end]
        end
    end
    isfile(zip) && close(zipfl)
    proj_ = fls[OutlineData][end]
    args = []

    # Set the LOADERS
    for DK in datakinds()
        if !haskey(dataset_LOADERS, DK)
            if haskey(fls,DK) # ITMIX file present
                fln, nr, nameinfile, date_range, proj = fls[DK]
                dataset_LOADERS[DK] = ITMIXLoader
                dataset_files[DK] = [joinpath(dir,fln)]
                dataset_opts[DK][:proj] = proj_
            else
                if DK==ParaData
                    dataset_LOADERS[DK] = ITMIXLoader
                    if iscalving_isicecap[Symbol(name)][2]
                        dataset_opts[DK][:iv_extrapolation_scaling] = :huss
                    end
                # elseif DK==BdotData # TODO: use GloGEM by default
                #     dataset_LOADERS[DK] = GloGEMLoader
                #     dataset_files[DK] = [joinpath(dir,fln)]
                elseif DK==IVData || DK==ThicknessData
                    dataset_LOADERS[DK] = NothingLoader
                else
                    # if there is no file use the synthetic loader
                    dataset_LOADERS[DK] = SyntheticBenchLoader
                end
            end
        else
            # custom set loader & other bits, leave as is
        end
    end

    # bug in shapefile
    if id==:Tasman
        dataset_files[OutlineData] = ["data/ITMIX/Tasman/01_margin_Tasman_1986_UTM-59.csv"]
    end
    return pl
end

# Outline
function load!(ds::DataSet{OutlineData, ITMIXLoader}, gid)
    opts = ds.opts
    proj = opts[:proj]
    # fl, fid = open_file(ds.files[1])
    # outline = readdlm(fl, ';', Float64)'
    # close(fid)
    @assert length(ds.files)==1
    if endswith(ds.files[1], ".shp")
        outline, splits = load_shapefile!(ds, :outline)
    else
        outline = @open_magic readdlm(fl, ';', F)'

        if size(outline,1)==3
            # TODO
            # # process boundary ID stuff:
            # # TODO: not sure this is right yet.
            # boundary_ind = outline[3,:]
            # xy = outline[1:2,:]
            # # all boundary_ind>1 treat as outcrops
            # outline = xy[:,boundary_ind.==boundary_ind[1]]
            # outcropsx = Float64[]
            # outcropsy = Float64[]
            # splits = Int[]
            # for i=2:1000
            #     inds = boundary_ind.==i
            #     length(inds)==0 && break
            #     append!(outcropsx, xy[1,inds])
            #     append!(outcropsy, xy[2,inds])
            # end
            # outline = xy
            outline = outline[1:2,:]
        end
        if outline[:,1] !=outline[:,end]
            # close polygon
            outline = hcat(outline, outline[:,1])
        end
        splits = VAWTools.find_poly_splits(outline)
    end
    # Fix Unteraar
    if getname(gid) == "Unteraar"
        outline = outline[:, splits[end-1]:end-1]
        # remove self-intersections
        outline = outline[:,[1:1483, 1487:1667, 1671:end;]]
        splits = [1,size(outline)[2]]
        @assert outline[:,1]==outline[:,end]
    end

    iscalving,isicecap,rgi = iscalving_isicecap[getid(gid)]
    @assert isa(iscalving, Bool) "Input whether calving or not in iscalving_dict"
    return Outline(gid, outline, splits),
           iscalving::Bool, isicecap::Bool, proj::String
end

"""
opts:
- grid : if specified use a different grid from the grid of the DEM as
  computational grid.
- cropbox -- used for cropping the DEM and other 2D data.  Needs to
  contain all of the outline. `F[xmin, xmax, ymin, ymax]`, defaults
  being calculated during DEM-load

Return:
- dem::Gridded
- glaciermask::Matrix{Bool} -- true on glacier
- landmask::Matrix{Bool} -- true on solid ground (i.e. false on glaciers, ice shelves, water)
- dist_land::Matrix{F} -- distance to closest land
- dist_land_maxdist -- maximum distance to which to calculate above
- cropbox::Vector{Float64} -- box around region of interest (contains glacier + some buffer)
                              `F[xmin, xmax, ymin, ymax]`
"""
function load!(ds::DataSet{DEMData, ITMIXLoader}, outline::Outline, pp::Phys
               )::Tuple{Gridded, Matrix{Bool}, Matrix{Bool}, Matrix{F}, F, Tuple, String}
    @assert ds.opts[:cropbox]==:from_DEM "ITMIX cropbox needs to be: `:from_DEM`"
    return load!(convert(DataSet{DEMData, :ASCIIGridLoader}, ds), outline, pp)
end
function load!(ds::DataSet{IVData, ITMIXLoader}, dem, glaciermask, target_proj,
               outline, cropbox, pp::Phys)::Tuple{SomeData, Matrix{Bool}}
    if splitext(ds.files[1])[2]==".txt"
        # IV at points
        tmp = @open_magic readdlm(fl, header=true)
        speed = F.(tmp[1][:,4])
        x = F.(tmp[1][:,end-2])
        y = F.(tmp[1][:,end-1])
        iv = Traj(x,y,speed)
        ivmask = Bool.(trues(0,0))
        return SomeData(vals=iv, sigma=ds.opts[:sigma]), ivmask
    elseif splitext(ds.files[1])[2]==".asc"
        # IV in 2D
        return load!(convert(DataSet{IVData, :ASCIIGridLoader}, ds), dem, glaciermask, target_proj,
                     outline, cropbox, pp)
    else
        error("Extension of file $file not supported.  Need .txt or .asc")
    end
end
function load!(ds::DataSet{ThicknessData, ITMIXLoader}, dem, glaciermask, target_proj, outline, cropbox,
               pp::Phys)::SomeData
    @unpack files, opts, daterange = ds
    @unpack gid = outline
    sigma = opts[:sigma]
    if gid.version==1
        zb = get_cache_ram!(ds, :zb, @open_magic_fn readdlm(fl; header=true)[1])
        x,y = zb[:,1], zb[:,2]
        zb = zb[:,3]
        izs  = Interp.interpolate((dem.x, dem.y), dem.v, Interp.Gridded(Interp.Linear()) )
        zs = [izs[xx,yy] for (xx,yy) in zip(x,y)]
        thick = zs-zb
        out = Traj(x, y, thick, zs)
        split_traj!(out, 1e3) # TODO adjust split-distance
    else
        zb = get_cache_ram!(ds, :zb, @open_magic_fn readdlm(fl; header=true, skipstart=16)[1])
        lines = unique(zb[:,1])
        x,y,thick,zs = zb[:,3], zb[:,4], zb[:,end], zb[:,5]
        splits = []
        for l in lines
            push!(splits, findfirst(zb[:,1], l):findlast(zb[:,1], l))
        end
        out = Traj(x, y, thick, zs, splits)
    end
    return SomeData(vals=out, sigma=sigma)
end
function load!(ds::DataSet{BdotData, ITMIXLoader}, dem, glaciermask, target_proj, outline, cropbox,
               pp::Phys)
    bdot = load_asciigrid!(ds, :bdot, cropbox)
    # scale to m ice / a
    scale_m_ice!(bdot, pp)
    if haskey(ds.opts, :vals_prior)
        return SomeData(vals=bdot, sigma=ds.opts[:sigma], vals_prior=ds.opts[:vals_prior])
    else
        return SomeData(vals=bdot, sigma=ds.opts[:sigma])
    end
end
function load!(ds::DataSet{DhdtData, ITMIXLoader},
               bdot, terminus_flux,
               dem, glaciermask, target_proj, outline, cropbox,
               iscalving, pp::Phys)
    return load!(convert(DataSet{DhdtData, :ASCIIGridLoader}, ds),
                 dem, glaciermask, target_proj, outline, cropbox, pp)
end

function load!(ds::DataSet{FslData, ITMIXLoader}, dem, glaciermask, target_proj, outline, cropbox,
               pp::Phys)
    error("FSL data not provided in ITMIX.  Use different loader, probably syntetic")
end
function load!(ds::DataSet{TempData, ITMIXLoader}, dem, glaciermask, target_proj, outline, cropbox,
               pp::Phys)
    error("Ice-temp data not provided in ITMIX.  Use different loader, probably syntetic")
end
function load!(ds::DataSet{TerminusFluxData, ITMIXLoader}, dem, glaciermask, iscalving,
               pp::Phys)
    error("Calving flux not provided in ITMIX.  Use different loader, probably synthetic")
end

#############################
# Helper

# Each of this folders contains a combination of the following files:
#
# - 01_margin_GGGGGG_YYYY_UTMZZZ.csv
# - 02_surface_GGGGGG_YYYY_UTMZZZ.asc  #
# - 03_RES_GGGGGG.txt                  # GPR file
# - 04_mb_GGGGGG_XXXX-YYYY_UTMZZZ.asc
# - 05_dhdt_GGGGGG_XXXX-YYYY_UTMZZZ.asc
# - 06_speed_GGGGGG_XXXX-YYYY_UTMZZZ.asc
# - 07_azimuth_GGGGGG_XXXX-YYYY_UTMZZZ.asc
# - 08_velocity_GGGGGG_XXXX-YYYY_UTMZZZ.txt
#
# where
#
# GGGGGG = name of the glacier (e.g. "Brewster")
# YYYY   = year the data are referring to (e.g. "2003")
# ZZZZ   = UTM zone of the coordinate system in which the data are projected (e.g. "32") - negative values indicate the southern hemisphere.
# Files 04-08 are mean values over the period "XXXX-YYYY" (e.g. "1997-2009")
#
# The file content is as follows:
# 01: Margin of the glacier (m).
# 02: DEM of the glacier (m a.s.l.) - Note that the DEM might be cropped with the glacier margin or not
# 04: Spatially distributed average (surface) mass balance of the glacier (m w.e./a)
# 05: Spatially distributed average ice thickness change rate (m/a)
# 06: Spatially distributed average surface flow speed of the glacier (m/a)
# 07: Spatially distributed average azimuth (="direction") of the surface flow (clockwise degrees from north)
# 08: Point observations of surface flow velocity (comprising speed and azimuth)
#
# .txt-files are plain text files.
# .csv-files are comma-separated-value files.
# .asc-files are ASCII-grid files (http://en.wikipedia.org/wiki/Esri_grid)
function parse_ITMIX_filename(fln::ZipFile.ReadableFile)
    error("zip files not supported currently")
    dir,fl = splitdir(fln.name)
    if contains(dir,"shapefiles")
        return nothing
    else
        return parse_ITMIX_filename(fl)
    end
end
function parse_ITMIX_filename(fln::String)
    fln=="" && return nothing
    fl, ext = splitext(splitdir(fln)[2])
    splits = find(c->c=='_', fl)
    length(splits)<2 && return nothing
    !(ext in [".txt",".asc",".shp"]) && return nothing # drop the ".csv" outlines in favor of shp
    nr = fl[1:splits[1]-1]
    !all(isnumber, nr) && return nothing
    nr = parse(Int,nr)
    dk = fl[splits[1]+1:splits[2]-1]
    # if dk in ["margin","surface"]
    #     name = fl[splits[2]+1:splits[3]-1]
    #     DK = dk=="margin" ? Outline : DEM
    #     y = parse(Int, fl[splits[3]+1:splits[4]-1])
    #     date_range = (Date(y),Date(y))
    # else
    DK = dk_dict[dk]
    if dk=="RES"
        name = fl[splits[2]+1:end]
        date_range = (Date(),Date())
    else
        name = fl[splits[2]+1:splits[3]-1]
        sp = split(fl[splits[3]+1:splits[4]-1],'-')
        if length(sp)==1
            y1 = parse(Int,sp[1])
            y2 = y1
        else
            y1,y2 = (parse(Int,s) for s in sp)
        end
        date_range = (Date(y1),Date(y2))
    end
    if dk!="RES"
        if name=="Unteraar"
            # Unteraar is not in proper UTM projection!
            zone = fl[splits[end]+1:end]
            proj = "close to UTM zone $zone, but no cigar"
        else
            zone = fl[splits[end]+1:end]
            proj = "+proj=utm +zone=$zone +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
        end
    else
        proj = "unknown"
    end
    return DK, fln, nr, name, date_range, proj
end

"Read a agr file"
_read(name, dirs, num) = read_agr(_getfile(dirs, name, num), F, NA=FILL)

function generate_DEM_masks(gid::ITMIXGlacier, dem, outline, sealevel_cutoff=10.0)
    glaciermask, _ = maskit(dem, outline)
    # also DEM holes from glaciermask
    glaciermask[dem.v.==FILL] = false
    # remove points below sea level
    glaciermask[dem.v.<0] = false

    remove_slivers = true
    landmask = generate_landmask(glaciermask, dem, sealevel_cutoff, remove_slivers)
    return glaciermask, landmask
end

######################
# Special processing (aka hacks)
######################

function special_processing!(gid::ITMIXGlacier, outline, dem, glaciermask, landmask,
                             iv, ivmask, thickness, bdot, dhdt, fsl, temp,
                             pl::LoadPara, pp::Phys, pm::MPara, pn::Num)
    if getname(gid)=="Austfonna"
        # no easy way of figuring land vs sea out for Austfonna, set all to sea.
        landmask[:] = false
        ivmask[iv.vals.v.<=0] = false
    elseif getname(gid)=="Tasman" || getname(gid)=="Devon" || getname(gid)=="Unteraar"
        # not masked at all in the input data
        ivmask[iv.vals.v.<=0] = false
    end
    return nothing
end

###########################
## ITMIX2 Utilities
###########################
module ITMIX2
import ..BITEModel; const BM = BITEModel
import VAWTools
"Number of experiments per glacier"
const nexps = 16
using DataFrames, FileIO, CSVFiles
int2expnr(i) = i>9 ? Symbol("exp$i") : Symbol("exp0$i")

"""
    load_experiment_table(fln="data/ITMIX2/ITMIX2_experiments_key.txt")

Load ITMIX2_experiments_key.txt as a DataFrame
"""
function load_experiment_table(fln)
    tab = DataFrame(load(File(format"CSV", fln),
                         spacedelim=true,
                         header_exists=true,
                         skiplines_begin=11,
                         type_detect_rows=512,
                         #colparsers=[String,Int, fill(x->Bool(parse(Int,x)),16)..., Int]
                         ))
    tab[:glacier] = [BM.ITMIXGlacier(Symbol(gl), 2) for gl in tab[:glacier]]
    for i=1:nexps
        tab[int2expnr(i)] = Bool.(tab[int2expnr(i)])
    end
    # add exp00
    tab[:exp00] = tab[:exp01].*false
    return tab
end

"Get the priority of a glacier"
getpriority(gid) = getexps(gid)[:prio][1]
getpriority() = Dict(((gln, getexps(gln)[:prio][1]) for gln in BM.ITMIXnames))

"Get experiments for a glacier `gl`"
getexps(gid::BM.ITMIXGlacier) = DataFrames.filter(row->row[:glacier]==gid, tab)
getexps(name::Union{Symbol,String}) = getexps(BM.ITMIXGlacier(name,2))

"Get fitting radar-lines of experiment `i` for a glacier `gl`."
get_rlines(gid::BM.ITMIXGlacier, i) = i==0 ? Int[] : find(getexps(gid)[int2expnr(i)])
get_rlines(name::Union{Symbol,String}, i) =  i==0 ? Int[] : get_rlines(BM.ITMIXGlacier(name,2), i)

fln = joinpath(@__DIR__, "../../scripts/data/ITMIX2/ITMIX2_experiments_key.txt")
"The ITMIX2 experiment table"
const tab = isfile(fln) ? load_experiment_table(fln) : nothing

"The ITMIX2 glaciers sorted according to their priorities"
const priorities = tab==nothing ? nothing : [[BM.ITMIXGlacier(gln,2) for gln in keys(getpriority()) if getpriority()[gln]==i] for i=1:4]

# # iteration over the ITMIX 2 experiments
# """
#     experiments([gid])

# Return iterator for ITMIX2 experiments.

# ```
# for (gid, i, radarlines) in experiments([gid])
#    ...
# end
# ```
# """
# function experiments end

# experiments(gid) = ITMIX2Iterator(gid, getexps(gid))
# struct ITMIX2Iterator
#     gid::BM.ITMIXGlacier
# end
# Base.start(::ITMIX2Iterator) = 1
# Base.next(ii::ITMIX2Iterator, state) = ( (gid, int2expnr(state), getexp(ii.gid, state)), state+1)
# Base.done(ii::ITMIX2Iterator, state) = state==17

# experiments() = ITMIX2Iterator_All()
# struct ITMIX2Iterator_All end
# Base.start(::ITMIX2Iterator_All) = (1, BM.ITMIXv2_glaciers[1])
# function Base.next(ii::ITMIX2Iterator_All, state)
#     expnr, gid = state
#     item = (gid, int2expnr(expnr), getexp(gid, expnr))
#     state = if expnr==16 && gid!=BM.ITMIXv2_glaciers[end]
#         (1, BM.ITMIXv2_glaciers[findfirst(BM.ITMIXv2_glaciers, gid)+1])
#     else
#         (expnr+1, gid)
#     end

#     return item, state
# end
# Base.done(ii::ITMIX2Iterator_All, state) = state[1]==17

"""
    make_ITMIX2_glacier(gl, expnr::Number, rm_iv=false)
    make_ITMIX2_glacier(gl, tracks2keep, rm_iv=false)

Take the glacier with all radarlines loaded and return one which only has the
lines of that experiment loaded.  The second interface allows to specify which
tracks to keep and ditch the rest.

Also has the extra feature to make more experiments:
- exp0: no radar lines
- exp 17-32: no IV, otherwise as 1.16

Above two only make sense (in the sense that there is something to fit)
for glaciers with IV data.
"""
make_ITMIX2_glacier(gl, expnr::Number, rm_iv=false) =
    make_ITMIX2_glacier(gl, get_rlines(gl.gid, expnr), rm_iv, expnr)
function make_ITMIX2_glacier(gl, tracks, rm_iv=false, expnr=-1)
    new_h = if length(tracks)>0  # ITMIX2 experiments
        BM.SomeData(gl.h, vals=VAWTools.Traj(gl.h.vals, tracks))
    else # remove all radar-lines: name it exp00.  This is the ITMIX1 run.
        new_h = BM.SomeData(nothing)
    end
    misc = merge(gl.misc, Dict{Symbol,Any}(:expnr=>expnr, :rm_iv=>rm_iv))
    return if rm_iv
        # check that there is IV
        gl.iv isa BM.SomeData{Void} && error("Cannot remove IV from glacier $(BM.getname(gl)) as it has no IV.")
        # remove IV data
        BM.Glacier(gl, h=new_h, iv=BM.SomeData(nothing), misc=misc)
    else
        # normal run
        BM.Glacier(gl, h=new_h, misc=misc)
    end
end

##########################
## Saving
##########################
# needs to be saved on the same grid as the input data
#=
ITMIX2 - instructions for results submission

1) Results format and content
- Submitted results shall be distributed estimates of the ice thickness.
- The unit for the ice thickness is "meters".
- Results are to be submitted as Esri grids, i.e. in the same format of the input data.
- All results shall have same resolution, cellsize, extent, and coordinate system as the input grids.

2) Results nomenclature
- All results shall be contained in a single directory, called "results_[Yourname]" (e.g. "results_Farinotti")
- The directory shall be compressed (.zip) before submission. The name of the uploaded file will thus be: "results_[Yourname].zip"
- The individual files in the folder shall be named
   "thickness_[NameOfTheGlacier]_[Yourname]_exp[XX].asc"
  where
   [NameOfTheGlacier] = Name of the glacier, as defined by the input files (e.g. "NorthGlacier")
   [Yourname] = Your family name (e.g. "Farinotti")
   [XX] = Number of the experiment, as given in the "experiments_key"-file, the result is referring to (e.g. "02")
  A valid filename will thus look like, e.g., "thickness_NorthGlacier_Farinotti_exp02.asc".
- For a test case to be labelled as "considered", results for all 16 experiments will need to be submitted.
  Since there are three test cases defined as "compulsory", the submitted directory will thus contain 3*16=48 files at least.

3) Results upload
- To uplod your results, simply drag-and-drop your .zip-repository to the following address:
 https://drive.switch.ch/index.php/s/hovc5mJjSztoLVd
- When uploaded, please send a notification email to: daniel.farinotti@ethz.ch
=#

"""
    filepaths(gl::BM.Glacier, pl, expnr, rm_iv, dir)

Returns saving filenames including their full path.
The path is `joinpath(pl.data_root, dir).

Returns file names for h, h-sigma, iv, iv-sigma, and jld.
"""
function filepaths(gl::BM.Glacier, pl, expnr, rm_iv, dir)
    expnrs = expnr>9 ? Symbol("exp$expnr") : Symbol("exp0$expnr")
    glaciername = BM.getname(gl)
    fl = rm_iv ?
        joinpath(pl.data_root, dir, "thickness_$(glaciername)_Werder_$(expnrs)_rmIV") :
        joinpath(pl.data_root, dir, "thickness_$(glaciername)_Werder_$(expnrs)")
    # make the directory if it does not exist
    !isdir(pl.data_root) && error("Data-root dir does not exist: $(pl.data_root)")
    !isdir(joinpath(pl.data_root, dir)) && mkpath(joinpath(pl.data_root, dir))
    return fl*".asc", fl*"-sigma.asc", fl*"-iv.asc", fl*"-iv-sigma.asc", fl*".jld"
end

"""
    save_itmix2_run(hs2d, hs2d_sigma, gl, pl, iv2d=nothing, iv2d_sigma=nothing)

Saves the thickness map according to ITMIX2 specs.
"""
function save_itmix2_run(hs2d, hs2d_sigma, gl, pl, ivs2d=nothing, ivs2d_sigma=nothing;
                         dir="results_Werder")
    expnr, rm_iv = gl.misc[:expnr], gl.misc[:rm_iv]
    fln, flns, fln_iv, fln_ivs, fln_jld = filepaths(gl, pl, expnr, rm_iv, dir)

    downsample = get(pl.dataset_opts[BM.DataKind], :grid_downsample_step, 1)
    hs2d, hs2d_sigma = copy(hs2d), copy(hs2d_sigma)
    hs2d[.!gl.glaciermask] = BM.FILL
    hs2d_sigma[.!gl.glaciermask] = BM.FILL
    h = VAWTools.Gridded(gl.dem.x, gl.dem.y, hs2d)
    hs = VAWTools.Gridded(gl.dem.x, gl.dem.y, hs2d_sigma)
    if ivs2d!=nothing
        ivs2d, ivs2d_sigma = copy(ivs2d), copy(ivs2d_sigma)
        ivs2d[.!gl.glaciermask] = BM.FILL
        ivs2d_sigma[.!gl.glaciermask] = BM.FILL
        iv = VAWTools.Gridded(gl.dem.x, gl.dem.y, ivs2d)
        ivs = VAWTools.Gridded(gl.dem.x, gl.dem.y, ivs2d_sigma)
    end
    if downsample!=1
        # up-sample output
        dem_orig = VAWTools.read_agr(pl.dataset_files[BM.DEMData][1])
        h = VAWTools.upsample(h, dem_orig.x, dem_orig.y)
        hs = VAWTools.upsample(hs, dem_orig.x, dem_orig.y)
        if ivs2d!=nothing
            iv = VAWTools.upsample(iv, dem_orig.x, dem_orig.y)
            ivs = VAWTools.upsample(ivs, dem_orig.x, dem_orig.y)
        end
    end
    VAWTools.write_agr(h, fln, NA_g=BM.FILL, NA_agr=-9999.0)
    VAWTools.write_agr(hs, flns, NA_g=BM.FILL, NA_agr=-9999.0)
    if ivs2d!=nothing
        VAWTools.write_agr(iv, fln_iv, NA_g=BM.FILL, NA_agr=-9999.0)
        VAWTools.write_agr(ivs, fln_ivs, NA_g=BM.FILL, NA_agr=-9999.0)
    end
    # TODO make sub-folders with computer-name-date
    return fln, flns, fln_iv, fln_ivs, fln_jld
end
end #module ITMIX2


text:
using Distributions
using Base.Test

# Core testing procedure

function test_mixture(g::UnivariateMixture, n::Int, ns::Int)
    X = zeros(n)
    for i = 1:n
        X[i] = rand(g)
    end

    K = ncomponents(g)
    pr = probs(g)
    @assert length(pr) == K

    # mean
    mu = 0.0
    for k = 1:K
        mu += pr[k] * mean(component(g, k))
    end
    @test mean(g) ≈ mu

    # evaluation of cdf
    cf = zeros(n)
    for k = 1:K
        c_k = component(g, k)
        for i = 1:n
            cf[i] += pr[k] * cdf(c_k, X[i])
        end
    end

    for i = 1:n
        @test cdf(g, X[i]) ≈ cf[i]
    end
    @test cdf.(g, X) ≈ cf

    # evaluation
    P0 = zeros(n, K)
    LP0 = zeros(n, K)
    for k = 1:K
        c_k = component(g, k)
        for i = 1:n
            x_i = X[i]
            P0[i,k] = pdf(c_k, x_i)
            LP0[i,k] = logpdf(c_k, x_i)
        end
    end

    mix_p0 = P0 * pr
    mix_lp0 = log.(mix_p0)

    for i = 1:n
        @test pdf(g, X[i])                  ≈ mix_p0[i]
        @test logpdf(g, X[i])               ≈ mix_lp0[i]
        @test componentwise_pdf(g, X[i])    ≈ vec(P0[i,:])
        @test componentwise_logpdf(g, X[i]) ≈ vec(LP0[i,:])
    end

    @test pdf.(g, X)                  ≈ mix_p0
    @test logpdf.(g, X)               ≈ mix_lp0
    @test componentwise_pdf(g, X)    ≈ P0
    @test componentwise_logpdf(g, X) ≈ LP0

    # sampling
    Xs = rand(g, ns)
    @test isa(Xs, Vector{Float64})
    @test length(Xs) == ns
    @test isapprox(mean(Xs), mean(g), atol=0.01)
end

function test_mixture(g::MultivariateMixture, n::Int, ns::Int)
    X = zeros(length(g), n)
    for i = 1:n
        X[:,i] = rand(g)
    end

    K = ncomponents(g)
    pr = probs(g)
    @assert length(pr) == K

    # mean
    mu = 0.0
    for k = 1:K
        mu += pr[k] * mean(component(g, k))
    end
    @test mean(g) ≈ mu

    # evaluation
    P0 = zeros(n, K)
    LP0 = zeros(n, K)
    for k = 1:K
        c_k = component(g, k)
        for i = 1:n
            x_i = X[:,i]
            P0[i,k] = pdf(c_k, x_i)
            LP0[i,k] = logpdf(c_k, x_i)
        end
    end

    mix_p0 = P0 * pr
    mix_lp0 = log.(mix_p0)

    for i = 1:n
        x_i = X[:,i]
        @test pdf(g, x_i)                  ≈ mix_p0[i]
        @test logpdf(g, x_i)               ≈ mix_lp0[i]
        @test componentwise_pdf(g, x_i)    ≈ vec(P0[i,:])
        @test componentwise_logpdf(g, x_i) ≈ vec(LP0[i,:])
    end

    @show g
    @show size(X)
    @show size(mix_p0)
    @test pdf(g, X)                  ≈ mix_p0
    @test logpdf(g, X)               ≈ mix_lp0
    @test componentwise_pdf(g, X)    ≈ P0
    @test componentwise_logpdf(g, X) ≈ LP0

    # sampling
    Xs = rand(g, ns)
    @test isa(Xs, Matrix{Float64})
    @test size(Xs) == (length(g), ns)
    @test isapprox(vec(mean(Xs, 2)), mean(g), atol=0.1)
    @test isapprox(cov(Xs, 2)      , cov(g) , atol=0.1)
end

function test_params(g::AbstractMixtureModel)
    C = eltype(g.components)
    pars = params(g)
    mm = MixtureModel(C, pars...)
    @test g.prior == mm.prior
    @test g.components == mm.components
end

function test_params(g::UnivariateGMM)
    pars = params(g)
    mm = UnivariateGMM(pars...)
    @test g == mm
end

# Tests

println("    testing UnivariateMixture")

g_u = MixtureModel(Normal, [(0.0, 1.0), (2.0, 1.0), (-4.0, 1.5)], [0.2, 0.5, 0.3])
@test isa(g_u, MixtureModel{Univariate, Continuous, Normal})
@test ncomponents(g_u) == 3
test_mixture(g_u, 1000, 10^6)
test_params(g_u)
@test minimum(g_u) == -Inf
@test maximum(g_u) == Inf

g_u = MixtureModel([TriangularDist(-1,2,0),TriangularDist(-.5,3,1),TriangularDist(-2,0,-1)])
@test minimum(g_u) == -2.0
@test maximum(g_u) == 3.0
@test insupport(g_u, 2.5) == true
@test insupport(g_u, 3.5) == false

g_u = UnivariateGMM([0.0, 2.0, -4.0], [1.0, 1.2, 1.5], Categorical([0.2, 0.5, 0.3]))
@test isa(g_u, UnivariateGMM)
@test ncomponents(g_u) == 3
test_mixture(g_u, 1000, 10^6)
test_params(g_u)
@test minimum(g_u) == -Inf
@test maximum(g_u) == Inf

println("    testing MultivariateMixture")
g_m = MixtureModel(
    IsoNormal[ MvNormal([0.0, 0.0], 1.0),
               MvNormal([0.2, 1.0], 1.0),
               MvNormal([-0.5, -3.0], 1.6) ],
    [0.2, 0.5, 0.3])
@test isa(g_m, MixtureModel{Multivariate, Continuous, IsoNormal})
@test length(components(g_m)) == 3
@test length(g_m) == 2
@test insupport(g_m, [0.0, 0.0]) == true
test_mixture(g_m, 1000, 10^6)
test_params(g_m)


text:
export solve_linear_static,solve_linear_eigen

function introduce_BC(K::SparseMatrixCSC{Float64},restrainedDOFs::Vector{Int})
    nDOF=size(K,1)
    mask=[(i in restrainedDOFs) ? false : true for i in 1:nDOF]
    return K[mask,mask]
end

function introduce_BC(K::Vector{Float64},restrainedDOFs::Vector{Int})
    nDOF=size(K,1)
    mask=[(i in restrainedDOFs) ? false : true for i in 1:nDOF]
    return K[mask]
end

function resolve_BC(d̄::Vector{Float64},restrainedDOFs::Vector{Int})
    rows=Array(1:length(d̄))
    cols=Array(1:length(d̄))
    vals=ones(length(d̄))
    for i in restrainedDOFs
        for j in 1:length(rows)
            if rows[j]>=i
                rows[j:end].+=1
                break
            end
        end
    end
    d̄=sparse(d̄)
    d=sparse(rows,cols,vals,length(restrainedDOFs)+length(d̄),length(d̄))*d̄
    return Array(d)
end

function K₂(structure::Structure,d₀::Vector{Float64})::SparseMatrixCSC{Float64}
    nDOF=length(structure.nodes)*6
    restrainedDOFs=[]
    K=spzeros(nDOF,nDOF)
    for elm in values(structure.beams)
        i = elm.node1.hid
        j = elm.node2.hid
        T=sparse(elm.T)

        I=collect(1:12)
        J=[6i-5:6i;6j-5:6j]
        G=sparse(I,J,1.0,12,nDOF)

        T=sparse(elm.T)
        dᵉ=T*[d₀[i*6-5:i*6];d₀[j*6-5:j*6]]
        # Kᵉ=integrateK!(elm)
        # σ=(Kᵉ*dᵉ)[1]
        Kᵢᵉ=elm.Kᵉ
        σ=(Kᵢᵉ*dᵉ)[1]
        Kᵉ=sparse(Kᵢᵉ)+integrateKσ(elm,σ)
        A=T*G
        rDOF=findall(x->x==true,elm.release)
        if length(rDOF)!=0
            K̄ᵉ,P̄ᵉ=FEStructure.FEBeam.static_condensation(Array(Kᵉ),zeros(12),rDOF)
            K+=A'*K̄ᵉ*A
        else
            K+=A'*Kᵉ*A
        end
    end
    return K
end

function solve_linear_static(structure::Structure,loadcase,restrainedDOFs::Vector{Int})
    @info "-------------------------- 求解一阶线性工况 ---------------------------" 工况=loadcase.id
    K̄=introduce_BC(structure.K,restrainedDOFs)
    P̄=introduce_BC(loadcase.P,restrainedDOFs)
    if USE_PARDISO
        ps=Pardiso.MKLPardisoSolver()
        Pardiso.set_matrixtype!(ps,1) #实对称矩阵
        ū=Pardiso.solve(ps,K̄,P̄)
    else
        ū=Symmetric(K̄) \ P̄
    end

    u=resolve_BC(ū,restrainedDOFs)

    P=structure.K*u
    @info "------------------------------ 求解完成 ------------------------------"
    return u,P
end

function solve_linear_buckling(structure::Structure,loadcase,restrainedDOFs::Vector{Int};path=pwd())
    @info "-------------------------- 求解屈曲模态特征值 --------------------------" 工况=loadcase.id 前一步基础工况=loadcase.plc
    K=structure.K
    P=loadcase.P
    u=zero(P)
    F₀=zero(P)
    if loadcase.plc!=""
        u=read_vector(joinpath(path,".analysis"),loadcase.plc*"_u.v")
        F₀=read_vector(joinpath(path,".analysis"),loadcase.plc*"_F.v")
        Kσ=calc_Kσ(structure,u)
        K=structure.K+Kσ #初始位移及初始刚度来自前一步结果
    end
    K̄=introduce_BC(K,restrainedDOFs)

    nev=min(loadcase.nev,size(K̄,1))
    tol=loadcase.tolev
    maxiter=loadcase.maxiterev
    ω²,ϕ=eigs(K̄,M̄,nev=nev,which=:SM,tol=tol,maxiter=maxiter)

    for i in 1:size(ϕ̄,2)
        println(resolve_BC(ϕ̄[:,i],restrainedDOFs))
    end
    @info "------------------------------ 求解完成 ------------------------------"
    return ω²,ϕ
end

function solve_2nd_static(structure::Structure,loadcase,restrainedDOFs::Vector{Int};conv_tol=1e-16,steps=10,max_iter=20,path=pwd())
    @info "------------------------ 求解二阶几何非线性工况 ------------------------" 工况=loadcase.id 前一步基础工况=loadcase.plc
    K=structure.K
    P=loadcase.P
    u=zero(P)
    F₀=zero(P)
    if loadcase.plc!=""
        u=read_vector(joinpath(path,".analysis"),loadcase.plc*"_u.v")
        F₀=read_vector(joinpath(path,".analysis"),loadcase.plc*"_F.v")
        K=K₂(structure,u)
    end
    K̄=introduce_BC(K,restrainedDOFs)
    P̄=introduce_BC(P,restrainedDOFs)
    ū=introduce_BC(u,restrainedDOFs)
    F=zero(P)
    F̄=zero(ū)
    # @info "initial" F[12]

    if USE_PARDISO
        ps=Pardiso.MKLPardisoSolver()
        Pardiso.set_matrixtype!(ps,1) #实对称矩阵
    end

    for step in 1:steps #等荷载增量
        Q̄=P̄*step/steps
        iter=0
        while true
            ΔQ̄=Q̄-F̄ #当前子步不平衡力
            if USE_PARDISO
                Δū=Pardiso.solve(ps,K̄,ΔQ̄)  #线性求解子步增量位移
            else
                Δū=Symmetric(K̄)\ΔQ̄
            end
            Δu=resolve_BC(Δū,restrainedDOFs)
            ū+=Δū
            u=resolve_BC(ū,restrainedDOFs) #T.L格式，相对零位形
            K=K₂(structure,u)#应力刚化/软化
            F+=K*Δu #计算非线性内力
            K̄=introduce_BC(K,restrainedDOFs)
            F̄=introduce_BC(F,restrainedDOFs)
            iter+=1

            error=maximum(abs.(Q̄-F̄))
            if error>conv_tol
                # @printf("Iteration %d, error=%f, tol=%4.3e\n",iter,error,conv_tol)
            else
                # @info "step "*string(step) F[12]
                # @info "Converged!"
                break
            end

            if iter >= max_iter
                @error "Max iteration reached, step not converged!"
                break
            end
        end
    end
    @info "------------------------------ 求解完成 ------------------------------"
    return u,F+F₀
end


text:
# Description: startup function for solving an equation

import PDESolver.solvePDE

#=
"""
This function invokes the solver for the Elliptic equations, using the
specified input file

Inputs:
input_file: a string containing the path to an input file, or just a file
name.  If it is just a file name, it is taken to be in the
users pwd.

Outputs:
mesh: an AbstractMesh object
sbp: the SBP operator used in solving the equation
eqn: the AbstractSolutionData object used to solve the equation
opts: the options dictonary

"""
function run_elliptic(input_file::AbstractString)

  mesh, sbp, eqn, opts, pmesh = createObjects(input_file)
  solve_elliptic(mesh, sbp, eqn, opts, pmesh)

  return mesh, sbp, eqn, opts
end
=#

"""
This function creates and initializes the mesh, sbp, eqn, and opts objects

Inputs:
    opts: options dictionary

Outputs:
mesh: an AbstractMesh.  The concrete type is determined by the options
dictionary
sbp: an AbstractOperator.  The concrete type is determined by the options
dictionary
eqn: an EllipticData object
opts: the options dictionary
pmesh: mesh used for preconditioning, can be same object as mesh
"""
function createObjects(opts::Dict)

  Tdim = opts["dimensions"]
  dofpernode = 1

  sbp, mesh, pmesh, Tsol, Tres, Tmsh, mesh_time = createMeshAndOperator(opts, dofpernode)

  # TODO: input argument for dofpernode

  # create elliptic equation
  # var_type = opts["variable_type"]
  # eqn = EllipticData_{Tsol, Tres, Tdim, Tmsh, var_type}(mesh, sbp, opts)
  eqn = EllipticData_{Tsol, Tres, Tdim, Tmsh}(mesh, sbp, opts)

  # initialize physics module and populate any fields in mesh and eqn that
  # depend on the physics module
  init(mesh, sbp, eqn, opts, pmesh)

  return mesh, sbp, eqn, opts, pmesh
end

"""
Given fully initialized mesh, sbp, eqn, opts, this function solves
the Elliptic equations.  The 4 object should be obtained from createObjects().

Specifically, it applies an initial condition and invokes a nonlinear
solver according to the options dictionary.

Inputs:
mesh: an AbstractMesh
sbp: an AbstractOperator
eqn: an AbstractEllipticData
opts: the options dictionary.  This must be the options dictionary returned
by createObjects().  Changing values in the options dictionary after
calling createObjects() results in undefined behavior.
pmesh: mesh used for preconditioning, can be same object as mesh.
default value of mesh

"""
function solvePDE(mesh::AbstractMesh, sbp::AbstractOperator, eqn::AbstractEllipticData, opts::Dict, pmesh::AbstractMesh=mesh)
  #delta_t = opts["delta_t"]   # delta_t: timestep for RK

  myrank = mesh.myrank
  # flag determines whether to calculate u, dR/du, or dR/dx (1, 2, or 3)
  t_max = opts["t_max"]
  flag = opts["run_type"]
  var_type = opts["variable_type"]

  # timestepping parameters
  order = opts["order"]       # order of accuracy

  # zero some arrays out in case the user forgot
  fill!(eqn.res, 0.0)
  fill!(eqn.res_vec, 0.0)
  res_vec = eqn.res_vec
  q_vec = eqn.q_vec       # solution at current timestep

  # calculate residual of some other function for res_reltol0
  # TODO: add a boolean options here?
  Relfunc_name = opts["Relfunc_name"]
  if haskey(ICDict, Relfunc_name)
    @mpi_master println("\ncalculating residual for relative residual tolerance")
    Relfunc = ICDict[Relfunc_name]
    @mpi_master println("Relfunc = ", Relfunc)
    Relfunc(mesh, sbp, eqn, opts, q_vec)

    tmp = calcResidual(mesh, sbp, eqn, opts, evalResidual)
    res_real = real(eqn.res_vec)
    opts["res_reltol0"] = tmp
    println("res_reltol0 = ", tmp)

    saveSolutionToMesh(mesh, res_real)
    writeVisFiles(mesh, "solution_relfunc")
  end

  # populate u0 with initial condition
  @mpi_master println("\nEvaluating initial condition")
  ICfunc_name = opts["IC_name"]
  ICfunc = ICDict[ICfunc_name]
  @mpi_master println("ICfunc = ", ICfunc)
  ICfunc(mesh, sbp, eqn, opts, q_vec)

  res_vec_exact = deepcopy(q_vec)

  #writeSolutionFiles(mesh, sbp, eqn, opts, "IC")
  saveSolutionToMesh(mesh, q_vec)

  writeVisFiles(mesh, "solution_ic")

  call_nlsolver(mesh, sbp, eqn, opts, pmesh)
  postproc(mesh, sbp, eqn, opts)

  MPI.Barrier(mesh.comm)

  return mesh, sbp, eqn, opts
end  # end function
#runtest(1)

"""
This function does post processing, if requested by the input options.
Typical post processing includes calculation of errors, norms of important
quantities, writing of files. etc.

Inputs:
mesh
sbp
eqn
opts
"""
function postproc(mesh, sbp, eqn, opts)

  ##### Do postprocessing ######
  myrank = mesh.myrank

  @mpi_master println("\nDoing postprocessing")

  if haskey(opts, "exactSolution")
    l2norm, lInfnorm = calcErrorL2Norm(mesh, sbp, eqn, opts)
    @mpi_master begin
      println("L2Norm = ", l2norm)
      println("LinfNorm = ", lInfnorm)
      fname = "l2norm.dat"
      f = open(fname, "w")
      println(f, l2norm)
      close(f)
    end
  end

  if haskey(opts, "Functional")
    if haskey(opts, "exactFunctional")
      exactFunctional = opts["exactFunctional"]
    end
    functional_value = Array{typeof(eqn.q[1,1,1])}(mesh.numDofPerNode)
    eqn.functional(mesh, sbp, eqn, opts, functional_value)

    @mpi_master begin
      println("functional = ", abs(real(functional_value[1]) - exactFunctional))
      fname = "functional.dat"
      f = open(fname, "w")
      println(f, abs(real(functional_value[1]) - exactFunctional))
      close(f)
    end
  end

  return nothing
end


text:
"""
    NeighborLongitudinalResult

A structure to retrieve information about a neihbor in the longitudinal direction i.e. rear and front neighbors on the same lane.
If the neighbor index is equal to `nothing` it means there is no neighbor.

# Fields
- `ind::Union{Nothing, Int64}` index of the neighbor in the scene
- `Δs::Float64` positive distance along the lane between vehicles positions
"""
struct NeighborLongitudinalResult
    ind::Union{Nothing, Int64} # index in scene of the neighbor
    Δs::Float64 # positive distance along lane between vehicles' positions
end

"""
    VehicleTargetPoint

Defines a point on an entity that is used to measure distances.
The following target points are supported and are subtypes of VehicleTargetPoint:
- `VehicleTargetPointFront`
- `VehicleTargetPointCenter`
- `VehicleTargetPointRear`

The method `targetpoint_delta(::VehicleTargetPoint, ::Entity)` can be used to
compute the delta in the longitudinal direction to add when considering a specific target point.
"""
abstract type VehicleTargetPoint end

"""
    VehicleTargetPointFront <: VehicleTargetPoint
Set the target point at the front (and center) of the entity.
"""
struct VehicleTargetPointFront <: VehicleTargetPoint end

targetpoint_delta(::VehicleTargetPointFront, veh::Entity{S, D, I}) where {S,D<:AbstractAgentDefinition, I} = length(veh.def)/2*cos(posf(veh.state).ϕ)

"""
    VehicleTargetPointCenter <: VehicleTargetPoint
Set the target point at the center of the entity.
"""
struct VehicleTargetPointCenter <: VehicleTargetPoint end

targetpoint_delta(::VehicleTargetPointCenter, veh::Entity{S, D, I}) where {S,D<:AbstractAgentDefinition, I} = 0.0

"""
    VehicleTargetPointFront <: VehicleTargetPoint
Set the target point at the front (and center) of the entity.
"""
struct VehicleTargetPointRear <: VehicleTargetPoint end

targetpoint_delta(::VehicleTargetPointRear, veh::Entity{S, D, I}) where {S,D<:AbstractAgentDefinition, I} = -length(veh.def)/2*cos(posf(veh.state).ϕ)

const VEHICLE_TARGET_POINT_CENTER = VehicleTargetPointCenter()

"""
    find_neighbor(scene::Scene, roadway::Roawday, ego::Entity; kwargs...)

Search through lanes and road segments to find a neighbor of `ego` in the `scene`.

# Arguments

- `scene::Scene` the scene containing all the entities
- `roadway::Roadway` the road topology on which entities are driving
- `ego::Entity` the entity that we want to compute the neighbor of.

# Keyword arguments
- `lane::Union{Nothing, Lane}` the lane on which to search the neighbor, if different from the ego vehicle current lane, it uses the projection of the ego vehicle on the given lane as a reference point. If nothing, returns nothing.
- `rear::Bool = false` set to true to search for rear neighbor, search forward neighbor by default
- `max_distance::Float64 = 250.0` stop searching after this distance is reached, if the neighbor is further than max_distance, returns nothing
- `targetpoint_ego::VehicleTargetPoint` the point on the ego vehicle used for distance calculation, see `VehicleTargetPoint` for more info
- `targetpoint_neighbor::VehicleTargetPoint` the point on the neighbor vehicle used for distance calculation, see `VehicleTargetPoint` for more info
- `ids_to_ignore::Union{Nothing, Set{I}} = nothing` a list of entity ids to ignore for the search, ego is always ignored.
"""
function find_neighbor(scene::Scene, roadway::Roadway, ego::Entity{S,D,I};
                       lane::Union{Nothing, Lane} = get_lane(roadway, ego),
                       rear::Bool=false,
                       max_distance::Float64=250.0,
                       targetpoint_ego::VehicleTargetPoint = VehicleTargetPointCenter(),
                       targetpoint_neighbor::VehicleTargetPoint = VehicleTargetPointCenter(),
                       ids_to_ignore::Union{Nothing, Set{I}} = nothing) where {S,D,I}

    if lane == nothing
        return NeighborLongitudinalResult(nothing, max_distance)
    elseif get_lane(roadway, ego).tag == lane.tag
        tag_start = lane.tag
    else  # project ego on desired lane
        roadproj = proj(posg(ego.state), lane, roadway)
        tag_start = roadproj.tag
    end
    s_base = posf(ego.state).s + targetpoint_delta(targetpoint_ego, ego)

    tag_target = tag_start
    best_ind = nothing
    best_dist = max_distance

    dist_searched = 0.0
    while dist_searched < max_distance

        curr_lane = roadway[tag_target]

        for (i, veh) in enumerate(scene)

            if veh.id == ego.id
                continue
            end

            if ids_to_ignore !== nothing && veh.id ∈ ids_to_ignore
                continue
            end

            # check if veh is on thislane
            s_adjust = NaN

            if get_lane(roadway, veh).tag == curr_lane.tag
                s_adjust = 0.0

            elseif is_between_segments_hi(posf(veh.state).roadind.ind, curr_lane.curve) &&
                   is_in_entrances(roadway[tag_target], posf(veh.state).roadind.tag)

                distance_between_lanes = norm(VecE2(roadway[tag_target].curve[1].pos - roadway[posf(veh.state).roadind.tag].curve[end].pos))
                s_adjust = -(roadway[posf(veh.state).roadind.tag].curve[end].s + distance_between_lanes)

            elseif is_between_segments_lo(posf(veh.state).roadind.ind) &&
                    is_in_exits(roadway[tag_target], posf(veh.state).roadind.tag)

                distance_between_lanes = norm(VecE2(roadway[tag_target].curve[end].pos - roadway[posf(veh.state).roadind.tag].curve[1].pos))
                s_adjust = roadway[tag_target].curve[end].s + distance_between_lanes
            end

            if !isnan(s_adjust)
                s_valid = posf(veh.state).s + targetpoint_delta(targetpoint_neighbor, veh) + s_adjust
                if rear
                    dist_valid = s_base - s_valid + dist_searched
                else
                    dist_valid = s_valid - s_base + dist_searched
                end

                if dist_valid ≥ 0.0
                    s_primary = posf(veh.state).s + targetpoint_delta(targetpoint_neighbor, veh) + s_adjust
                    if rear
                        dist=  s_base - s_primary + dist_searched
                    else
                        dist = s_primary - s_base + dist_searched
                    end
                    if dist < best_dist
                        best_dist = dist
                        best_ind = i
                    end
                end
            end
        end

        # neighbor has been found exit
        if best_ind != nothing
            break
        end

        # no next lane and no neighbor found exit
        if !has_next(curr_lane) ||
            (tag_target == tag_start && dist_searched != 0.0) # exit after visiting this lane a 2nd time
            break
        end

        # go to the connected lane.
        if rear
            dist_searched += s_base
            s_base = curr_lane.curve[end].s + norm(VecE2(lane.curve[end].pos - prev_lane_point(curr_lane, roadway).pos))
            tag_target = prev_lane(lane, roadway).tag
        else
            dist_searched += (curr_lane.curve[end].s - s_base)
            s_base = -norm(VecE2(curr_lane.curve[end].pos - next_lane_point(curr_lane, roadway).pos)) # negative distance between lanes
            tag_target = next_lane(curr_lane, roadway).tag
        end
    end

    return NeighborLongitudinalResult(best_ind, best_dist)
end

"""
    FrenetRelativePosition

Contains information about the projection of a point on a lane. See `get_frenet_relative_position`.

# Fields
- `origin::RoadIndex` original roadindex used for the projection, contains the target lane ID.
- `target::RoadIndex` roadindex reached after projection
- `Δs::Float64` longitudinal distance to the original roadindex
- `t::Float64` lateral distance to the original roadindex in the frame of the target lane
- `ϕ::Float64` angle with the original roadindex in the frame of the target lane
"""
struct FrenetRelativePosition
    origin::RoadIndex
    target::RoadIndex
    Δs::Float64
    t::Float64
    ϕ::Float64
end

"""
    get_frenet_relative_position(veh_fore::Entity, veh_rear::Entity, roadway::Roadway)
return the Frenet relative position between the two vehicles. It projects the position of the first vehicle onto the lane of the second vehicle.
The result is stored as a `FrenetRelativePosition`.

Lower level:
    get_frenet_relative_position(posG::VecSE2{Float64}, roadind::RoadIndex, roadway::Roadway;
        max_distance_fore::Float64 = 250.0, # max distance to search forward [m]
        max_distance_rear::Float64 = 250.0, # max distance to search backward [m]
        improvement_threshold::Float64 = 1e-4,
    )

Project the given point to the same lane as the given RoadIndex.
This will return the projection of the point, along with the Δs along the lane from the RoadIndex.
The returned type is a `FrenetRelativePosition` object.
"""
function get_frenet_relative_position(posG::VecSE2{Float64}, roadind::RoadIndex, roadway::Roadway;
    max_distance_fore::Float64 = 250.0, # max distance to search forward [m]
    max_distance_rear::Float64 = 250.0, # max distance to search backward [m]
    improvement_threshold::Float64 = 1e-4,
    )

    # project to current lane first
    tag_start = roadind.tag
    lane_start = roadway[tag_start]
    curveproj_start = proj(posG, lane_start, roadway, move_along_curves=false).curveproj
    curvept_start = lane_start[curveproj_start.ind, roadway]
    s_base = lane_start[roadind.ind, roadway].s
    s_proj = curvept_start.s
    Δs = s_proj - s_base
    sq_dist_to_curve = normsquared(VecE2(posG - curvept_start.pos))
    retval = FrenetRelativePosition(roadind,
                RoadIndex(curveproj_start.ind, tag_start), Δs, curveproj_start.t, curveproj_start.ϕ)

    # search downstream
    if has_next(lane_start)
        dist_searched = lane_start.curve[end].s - s_base
        s_base = -norm(VecE2(lane_start.curve[end].pos - next_lane_point(lane_start, roadway).pos)) # negative distance between lanes
        tag_target = next_lane(lane_start, roadway).tag
        while dist_searched < max_distance_fore

            lane = roadway[tag_target]
            curveproj = proj(posG, lane, roadway, move_along_curves=false).curveproj
            sq_dist_to_curve2 = normsquared(VecE2(posG - lane[curveproj.ind, roadway].pos))

            if sq_dist_to_curve2 < sq_dist_to_curve - improvement_threshold

                sq_dist_to_curve = sq_dist_to_curve2
                s_proj = lane[curveproj.ind, roadway].s
                Δs = s_proj - s_base + dist_searched
                retval = FrenetRelativePosition(
                            roadind, RoadIndex(curveproj.ind, tag_target),
                            Δs, curveproj.t, curveproj.ϕ)
            end

            if !has_next(lane)
                break
            end

            dist_searched += (lane.curve[end].s - s_base)
            s_base = -norm(VecE2(lane.curve[end].pos - next_lane_point(lane, roadway).pos)) # negative distance between lanes
            tag_target = next_lane(lane, roadway).tag

            if tag_target == tag_start
                break
            end
        end
    end

    # search upstream
    if has_prev(lane_start)
        dist_searched = s_base
        tag_target = prev_lane(lane_start, roadway).tag
        s_base = roadway[tag_target].curve[end].s + norm(VecE2(lane_start.curve[1].pos - prev_lane_point(lane_start, roadway).pos)) # end of the lane

        while dist_searched < max_distance_fore

            lane = roadway[tag_target]
            curveproj = proj(posG, lane, roadway, move_along_curves=false).curveproj
            sq_dist_to_curve2 = normsquared(VecE2(posG - lane[curveproj.ind, roadway].pos))

            if sq_dist_to_curve2 < sq_dist_to_curve - improvement_threshold

                sq_dist_to_curve = sq_dist_to_curve2
                s_proj = lane[curveproj.ind, roadway].s
                dist = s_base - s_proj  + dist_searched
                Δs = -dist
                retval = FrenetRelativePosition(
                            roadind, RoadIndex(curveproj.ind, tag_target),
                            Δs, curveproj.t, curveproj.ϕ)
            end

            if !has_prev(lane)
                break
            end

            dist_searched += s_base
            s_base = lane.curve[end].s + norm(VecE2(lane.curve[1].pos - prev_lane_point(lane, roadway).pos)) # length of this lane plus crossover
            tag_target = prev_lane(lane, roadway).tag

            if tag_target == tag_start
                break
            end
        end
    end

    retval
end
get_frenet_relative_position(veh_fore::Entity{S, D, I}, veh_rear::Entity{S, D, I}, roadway::Roadway) where {S,D, I} = get_frenet_relative_position(posg(veh_fore.state), posf(veh_rear.state).roadind, roadway)

"""
    dist_to_front_neighbor(roadway::Roadway, scene::Scene, veh::Entity)
Feature function to extract the longitudinal distance to the front neighbor (in the Frenet frame).
Returns `missing` if there are no front neighbor.
"""
function dist_to_front_neighbor(roadway::Roadway, scene::Scene, veh::Entity)
    neighbor = find_neighbor(scene, roadway, veh)
    if neighbor.ind === nothing
        return missing
    else
        return neighbor.Δs
    end
end

"""
    front_neighbor_speed(roadway::Roadway, scene::Scene, veh::Entity)
Feature function to extract the velocity of the front neighbor.
Returns `missing` if there are no front neighbor.
"""
function front_neighbor_speed(roadway::Roadway, scene::Scene, veh::Entity)
    neighbor = find_neighbor(scene, roadway, veh)
    if neighbor.ind === nothing
        return missing
    else
        return vel(scene[neighbor.ind])
    end
end

"""
    time_to_collision(roadway::Roadway, scene::Scene, veh::Entity)
Feature function to extract the time to collision with the front neighbor.
Returns `missing` if there are no front neighbor.
"""
function time_to_collision(roadway::Roadway, scene::Scene, veh::Entity)
    neighbor = find_neighbor(scene, roadway, veh)
    if neighbor.ind === nothing
        return missing
    else
        len_ego = length(veh.def)
        len_oth = length(scene[neighbor.ind].def)
        Δs = neighbor.Δs - len_ego/2 - len_oth/2
        Δv = vel(scene[neighbor.ind]) - vel(veh)
        return -Δs / Δv
    end
end


text:
using Random
using MLDataUtils

"""
$(TYPEDEF)
"""
struct Dataset{F,G,H}
    xs::F
    us::G
    xsp::H
end

function rmNaN(data::Dataset)
    goodidx1 = Set([a[2] for a in findall(sum((!isnan).(data.xs), dims=1) .== size(data.xs, 1))])
    goodidx2 = Set([a[2] for a in findall(sum((!isnan).(data.xsp), dims=1) .== size(data.xsp, 1))])
    allgood = collect(intersect(goodidx1, goodidx2))
    newxs = data.xs[:, allgood]
    newus = data.us[:, allgood]
    newxsp = data.xsp[:, allgood]
    @assert !any(isnan, newxs)
    @assert !any(isnan, newus)
    @assert !any(isnan, newxsp)
    newd = Dataset(newxs, newus, newxsp)
    oldlen = length(data)
    newlen = length(newd)
    @info "Previously, $oldlen. Now, $newlen"
end

"""
Adds noise to position and velocity measuresments in a dataset

$(SIGNATURES)
"""
function addnoise(data::Dataset, std::Float64=0.001)
    Dataset(data.xs .+ std .* randn(size(data.xs)), data.us, data.xsp .+ std .* randn(size(data.xsp)))
end

function MLDataUtils.splitobs(data::Dataset, at::AbstractFloat=0.7)
    train, test = MLDataUtils.splitobs((data.xs, data.us, data.xsp))
    return Dataset(train...), Dataset(test...)
end

function Random.shuffle(r::AbstractRNG, d::Dataset)
    newbid = randperm(r, length(d))
    return Dataset(d.xs[:, newbid], d.us[:, newbid], d.xsp[:, newbid])
end

Base.length(d::Dataset) = size(d.xs)[end]

function Base.hcat(ds::Dataset...)
    newxs = reduce(hcat, [d.xs for d in ds])
    newus = reduce(hcat, [d.us for d in ds])
    newxsp = reduce(hcat, [d.xsp for d in ds])
    return Dataset(newxs, newus, newxsp)
end

function MLDataUtils.eachbatch(data::Dataset, size::Int)
    return MLDataUtils.eachbatch((data.xs, data.us, data.xsp), size=size)
end

function MLDataUtils.shuffleobs(data::Dataset)
    return MLDataUtils.shuffleobs((data.xs, data.us, data.xsp))
end


text:
# Derived from https://github.com/zachallaun/FactCheck.jl/blob/master/src/FactCheck.jl, hacked for CodeWars

module Test

export @fact,
       @fact_throws,
       facts,
       context,
       getstats,
       exitstatus,
       # assertion helpers
       not,
       truthy,
       falsey,
       falsy,
       anything,
       irrelevant,
       exactly,
       roughly,
       @runtest

allresults = {}

# HACK: get the current line number
#
# This only works inside of a function body:
#
#     julia> hmm = function()
#                2
#                3
#                getline()
#            end
#
#     julia> hmm()
#     4
#
function getline()
    bt = backtrace()
    issecond = false
    for frame in bt
        lookup = ccall(:jl_lookup_code_address, Any, (Ptr{Void}, Int32), frame, 0)
        if lookup != ()
            if issecond
                return lookup[3]
            else
                issecond = true
            end
        end
    end
end

# Represents the result of a test. The `meta` dictionary is used to retain
# information about the test, such as its file, line number, description, etc.
#
abstract Result
type Success <: Result
    expr::Expr
    val
    meta::Dict
end
type Failure <: Result
    expr::Expr
    val
    meta::Dict
end
type Error <: Result
    expr::Expr
    err::Exception
    backtrace
    meta::Dict
end

LF = "<:LF:>\n"

function format_line(s...)
    context = isempty(contexts) ? () : ("<IT::>$(contexts[end])",LF)
    apply(string,tuple(context...,s...,LF))
end

# A TestSuite collects the results of a series of tests, as well as some
# information about the tests such as their file and description.
#
type TestSuite
    filename
    desc
    successes::Array{Success}
    failures::Array{Failure}
    errors::Array{Error}
end
function TestSuite(filename, desc)
    TestSuite(filename, desc, Success[], Failure[], Error[])
end

pluralize(s::String, n::Number) = n == 1 ? s : string(s, "s")

# Implementing Base.show(io::IO, t::SomeType) gives you control over the
# printed representation of that type. For example:
#
#     type Foo
#     a
#     end
#
#     show(io::IO, f::Foo) = print("Foo: a=$(repr(f.a))")
#
#     print(Foo("attr"))
#     # prints Foo: a="attr"
#
import Base.show

function mistake(ex::Expr, actual)
    x,y = ex.args
    "$(repr(x)) - expected: $(repr(y)) actual: $(actual)$(LF)"
end

function show(io::IO, f::Failure)
    print(io, format_line("<FAILURE::>", mistake(f.expr, f.val)))
end

import Base.error
function show(io::IO, e::Error)
    bt = replace(strip(sprint(io->Base.show_backtrace(io, e.backtrace)),'\n'), "\n", LF)
    print(io, format_line("<ERROR::>", mistake(e.expr, e.err), bt))
end

function show(io::IO, s::Success)
    print(io, format_line("<PASSED::>Test Passed"))
end

function format_suite(suite::TestSuite)
    suite.desc != nothing ? "<DESCRIBE::>$(suite.desc)$(LF)" : ""
end

# FactCheck core functions and macros
# ========================================

# The last handler function found in `handlers` will be passed test results.
# This means the default handler set up by FactCheck could be overridden with
# `push!(FactCheck.handlers, my_custom_handler)`.
#
const handlers = Function[]

# A list of test contexts. `contexts[end]` should be the inner-most context.
#
const contexts = String[]

# `do_fact` constructs a Success, Failure, or Error depending on the outcome
# of a test and passes it off to the active test handler (`FactCheck.handlers[end]`).
#
# `thunk` should be a parameterless boolean function representing a test.
# `factex` should be the Expr from which `thunk` was constructed.
# `meta` should contain meta information about the test.
#
function do_fact(thunk::Function, factex::Expr, meta::Dict)
    result = try
        res, val = thunk()
        res ? Success(factex, val, meta) : Failure(factex, val, meta)
    catch err
        Error(factex, err, catch_backtrace(), meta)
    end

    !isempty(handlers) && handlers[end](result)
    push!(allresults, result)
    result
end

# Constructs a boolean expression from a given expression `ex` that, when
# evaluated, returns true if `ex` throws an error and false if `ex` does not.
#
throws_pred(ex) = quote
    try
        $(esc(ex))
        (false, "no error")
    catch e
        (true, "error")
    end
end

# Constructs a boolean expression from two values that works differently
# depending on what `assertion` evaluates to.
#
# If `assertion` evaluates to a function, the result of the expression will be
# `assertion(ex)`. Otherwise, the result of the expression will be
# `assertion == ex`.
#
function fact_pred(ex, assertion)
    quote
        pred = function(t)
            e = $(esc(assertion))
            isa(e, Function) ? (e(t), t) : (e == t, t)
        end
        pred($(esc(ex)))
    end
end

# `@fact` rewrites assertions and generates calls to `do_fact`, which
# is responsible for actually running the test.
#
#     macroexpand(:(@fact 1 => 1))
#     #=> do_fact( () -> 1 == 1, :(1 => 1), ...)
#
macro fact(factex::Expr)
    if factex.head == :(=>)
        :(do_fact(() -> $(fact_pred(factex.args...)),
                  $(Expr(:quote, factex)),
                  {"line" => getline()}))
    else
        error("@fact doesn't support expression: $factex")
    end
end

macro fact_throws(factex::Expr)
    :(do_fact(() -> $(throws_pred(factex)),
              $(Expr(:quote, factex)),
              {"line" => getline()}))
end

# Constructs a function that handles Successes, Failures, and Errors,
# pushing them into a given TestSuite and printing Failures and Errors
# as they arrive.

function make_handler(suite::TestSuite)
    function delayed_handler(r::Success)
        push!(suite.successes, r)
        print(r)
    end
    function delayed_handler(r::Failure)
        push!(suite.failures, r)
        print(r)
    end
    function delayed_handler(r::Error)
        push!(suite.errors, r)
        print(r)
    end
    delayed_handler
end

# Executes a battery of tests in some descriptive context.
#
function context(f::Function, desc)
    push!(contexts, desc)
    f()
    pop!(contexts)
end
context(f::Function) = f()

# `facts` creates test scope. It is responsible for setting up a testing
# environment, which means constructing a `TestSuite`, generating and
# registering test handlers, and reporting results.
#
# `f` should be a function containing `@fact` invocations.
#
facts(f::Function) = facts(f, nothing)
function facts(f::Function, desc)
    suite = TestSuite(nothing, desc)
    test_handler = make_handler(suite)
    push!(handlers, test_handler)

    print(format_suite(suite))

    f()

    pop!(handlers)
end

# `getstats` return a dictionary with a summary over all tests run

function getstats()
    s = 0
    f = 0
    e = 0
    ns = 0
    for r in allresults
        if isa(r, Success)
            s += 1
        elseif isa(r, Failure)
            f += 1
            ns += 1
        elseif isa(r, Error)
            e += 1
            ns += 1
        end
    end
    assert(s+f+e == length(allresults) == s+ns)
    {"nSuccesses" => s, "nFailures" => f, "nErrors" => e, "nNonSuccessful" => ns}
end

exitstatus() = exit(getstats()["nNonSuccessful"])

# Assertion helpers
# =================

# Logical not for values and functions.
not(x) = isa(x, Function) ? (y) -> !x(y) : (y) -> x != y

# Truthiness is defined as not `nothing` or `false` (which is 0).
# Falsiness is its opposite.
#
truthy(x) = nothing != x != false
falsey = falsy = not(truthy)

irrelevant = anything(x) = true

# Can be used to test object/function equality:
#
#     @fact iseven => exactly(iseven)
#
exactly(x) = (y) -> is(x, y)

# Useful for comparing floating point numbers:
#
#     @fact 4.99999 => roughly(5)
#

roughly(n::Number; kvtols...) = i::Number -> isapprox(i,n; kvtols...)

roughly(X::AbstractArray; kvtols...) = Y::AbstractArray -> begin
    if size(X) != size(Y)
        return false
    end

    for i in 1:length(X)
        if !isapprox(X[i], Y[i]; kvtols...)
            return false
        end
    end
    return true
end

macro runtest(pkg, files...)
  for f in files
    include(Pkg.dir("$pkg/test/$f.jl"))
  end
end

end # module FactCheck

max_stars_repo_name: xcthulhu/codewars-runner

max_issues_repo_name: xcthulhu/codewars-runner

max_forks_repo_name: xcthulhu/codewars-runner


text:
#Jupyter Notebook display
function show(io::IO, ::MIME"text/html", ec::EChart)

    divid = "Echart" * randstring(10)
    option = json(ec.options)
    width = ec.width
    height = ec.height

        display("text/html", """

              """)
end

function html(ec::EChart;public=false)

    divid = "Echart" * randstring(10)
    option = json(ec.options)
    width = ec.width
    height = ec.height

    public==true ? src="""""" : src=""

        return """

                $src

              """
end