text:
function get_affine_map(ds::ArchGDAL.Dataset)
    # ArchGDAL fails hard on datasets without
    # an affinemap. GDAL documents that on fail
    # a default affinemap should be returned.
    try
        global gt = ArchGDAL.getgeotransform(ds)
    catch y
        @warn y.msg
        global gt = [0.,1.,0.,0.,0.,1.]
    end
    geotransform_to_affine(SVector{6, Float64}(gt))
end

function geotransform_to_affine(gt::SVector{6, Float64})
    # See https://lists.osgeo.org/pipermail/gdal-dev/2011-July/029449.html
    # for an explanation of the geotransform format
    AffineMap(SMatrix{2,2}([gt[2] gt[3]; gt[5] gt[6]]), SVector{2}([gt[1], gt[4]]))
end
geotransform_to_affine(A::Vector{Float64}) = geotransform_to_affine(SVector{6}(A))

function affine_to_geotransform(am::AffineMap{SArray{Tuple{2,2},Float64,2,4},SArray{Tuple{2},Float64,1,2}})
    l = am.linear
    t = am.translation
    [t[1], l[1], l[3], t[2], l[2], l[4]]
end

function affine_to_geotransform(am::AffineMap{Array{Float64,2},Array{Float64,1}})
    l = am.linear
    t = am.translation
    (length(l) != 4 || length(t) != 2) || error("AffineMap has wrong dimensions.")
    [t[1], l[1], l[3], t[2], l[2], l[4]]
end

"""Check wether the AffineMap of a GeoArray contains rotations."""
function is_rotated(ga::GeoArray)
    ga.f.linear[2] != 0. || ga.f.linear[3] != 0.
end

function bbox(ga::GeoArray)
    ax, ay = ga.f(SVector(0,0))
    bx, by = ga.f(SVector(size(ga)[1:2]))
    (min_x=min(ax, bx), min_y=min(ay, by), max_x=max(ax, bx), max_y=max(ay, by))
end

function unitrange_to_affine(x::StepRangeLen, y::StepRangeLen)
    δx, δy = step(x), step(y)
    AffineMap(
        SMatrix{2,2}(δx, 0, 0, δy),
        SVector(x[1] - δx/2, y[1] - δy/2)
    )
end

function bbox_to_affine(size::Tuple{Integer, Integer}, bbox::NamedTuple{(:min_x, :min_y, :max_x, :max_y),Tuple{Float64, Float64, Float64, Float64}})
    AffineMap(
        SMatrix{2,2}((bbox.max_x - bbox.min_x) / size[1], 0, 0, (bbox.max_y - bbox.min_y)/size[2]),
        SVector(bbox.min_x, bbox.min_y)
        )
end

"""Set geotransform of `GeoArray` by specifying a bounding box.
Note that this only can result in a non-rotated or skewed `GeoArray`."""
function bbox!(ga::GeoArray, bbox::NamedTuple{(:min_x, :min_y, :max_x, :max_y),Tuple{Float64, Float64, Float64, Float64}})
    ga.f = bbox_to_affine(size(ga)[1:2], bbox)
    ga
end

"""Generate bounding boxes for GeoArray cells."""
function bboxes(ga::GeoArray)
    c = coords(ga)::Array{StaticArrays.SArray{Tuple{2},Float64,1,2},2}
    m, n = size(c)
    cellbounds = Matrix{NamedTuple}(undef, (m-1, n-1))
    for j = 1:n-1, i = 1:m-1
        v = c[i:i+1, j:j+1]
        minx, maxx = extrema(first.(v))::Tuple{Float64, Float64}
        miny, maxy = extrema(last.(v))::Tuple{Float64, Float64}
        cellbounds[i, j] = (min_x=minx, max_x=maxx, min_y=miny, max_y=maxy)
    end
    cellbounds
end

# Extend CoordinateTransformations
CoordinateTransformations.compose(ga::GeoArray, t2::AffineMap) = CoordinateTransformations.compose(ga.f, t2)
CoordinateTransformations.compose(ga::GeoArray, t2::LinearMap) = CoordinateTransformations.compose(ga.f, t2)
CoordinateTransformations.compose(ga::GeoArray, t2::Translation) = CoordinateTransformations.compose(ga.f, t2)

"""Transform an GeoArray by applying a Transformation."""
function compose!(ga::GeoArray, t2::AffineMap)
    ga.f = compose(ga, t2)
end
function compose!(ga::GeoArray, t2::LinearMap)
    ga.f = compose(ga, t2)
end
function compose!(ga::GeoArray, t2::Translation)
    ga.f = compose(ga, t2)
end

"""Function to flip GeoArray upside down to adjust to GDAL ecosystem."""
function flipud!(ga::GeoArray)
    # Flip data upside down
    ga.A = reverse(ga.A, dims=2)

    # Find new corner coordinates
    ux, uy = ga.f(SVector{2}([0, size(ga)[2]]))

    # Define y mirror and compose
    lm = LinearMap(SMatrix{2,2}([1.0 0.0; 0.0 -1.0]))
    am = compose(ga.f, lm)  # AffineMap
    translate = SVector{2}([ux, uy])
    f = AffineMap(am.linear, translate)

    ga.f = f
    ga
end


text:
# All Systems are defined in these files
include("systems/core.jl")
include("systems/camera.jl")
include("systems/mesher.jl")

abstract type AbstractRenderSystem <: System end
include("systems/rendering/prerenderer.jl")
include("systems/rendering/default.jl")
include("systems/rendering/lines.jl")
include("systems/rendering/depthpeeling.jl")
include("systems/rendering/text.jl")
include("systems/rendering/cimgui.jl")
include("systems/rendering/uniforms.jl")
include("systems/rendering/final.jl")
include("systems/rendering/uploading.jl")
include("systems/interaction.jl")
include("systems/editor.jl")

# # various simulation systems
include("systems/simulations.jl")

ext: jl


text:
module C

const parameters = [
    # Birtwistle et al., MSB (2007)
    "VmaxPY",
    "KmPY",
    "kdeg",
    "kf47",
    "Vmaxr47",
    "Kmf47",
    "Kmr47",
    "kf48",
    "Kmf48",
    "Kmr48",
    "PTEN",
    "kf49",
    "kr49",
    "Kmf49",
    "Kmr49",
    "Kmr49b",
    "kr49b",
    "kf51",
    "Vmaxr51",
    "Kmf51",
    "Kmr51",
    "Kmrb51",
    "kf52",
    "Vmaxr52",
    "Kmf52",
    "Kmr52",
    "kf54",
    "Vmaxr54",
    "Kmf54",
    "Kmr54",
    "kf55",
    "Vmaxr55",
    "Kmf55",
    "Kmr55",
    "kf38",
    "kf39",
    "kf50",
    "a98",
    "b98",
    "koff46",
    "EGF_off",
    "HRGoff_3",
    "HRGoff_4",
    "koff4",
    "koff5",
    "koff6",
    "koff7",
    "koff8",
    "koff9",
    "koff61",
    "koff62",
    "koff16",
    "koff17",
    "koff18",
    "koff19",
    "koff20",
    "koff21",
    "koff22",
    "koff23",
    "koff24",
    "koff25",
    "koff26",
    "koff27",
    "koff28",
    "koff29",
    "koff30",
    "koff31",
    "koff32",
    "koff33",
    "koff34",
    "koff35",
    "koff36",
    "koff37",
    "koff65",
    "koff66",
    "koff67",
    "koff68",
    "koff69",
    "koff70",
    "koff71",
    "koff72",
    "koff40",
    "koff41",
    "koff42",
    "koff43",
    "koff44",
    "koff45",
    "koff57",
    "koff58",
    "koff59",
    "koff60",
    "kPTP10",
    "kPTP11",
    "kPTP12",
    "kPTP13",
    "kPTP14",
    "kPTP15",
    "kPTP63",
    "kPTP64",
    "koff73",
    "koff74",
    "koff75",
    "koff76",
    "koff77",
    "koff78",
    "koff79",
    "koff80",
    "kPTP38",
    "kPTP39",
    "koff88",
    "kPTP50",
    "kf81",
    "Vmaxr81",
    "Kmf81",
    "Kmr81",
    "kf82",
    "Vmaxr82",
    "Kmf82",
    "Kmr82",
    "kf83",
    "Vmaxr83",
    "Kmf83",
    "Kmr83",
    "kf84",
    "Vmaxr84",
    "Kmf84",
    "Kmr84",
    "kf85",
    "Vmaxr85",
    "Kmf85",
    "Kmr85",
    "kcon49",
    "kon1",
    "kon86",
    "kon2",
    "kon3",
    "kon87",
    "kon4",
    "kon5",
    "kon6",
    "kon7",
    "kon8",
    "kon9",
    "kon61",
    "kon62",
    "kf10",
    "kf11",
    "kf12",
    "kf13",
    "kf14",
    "kf15",
    "kf63",
    "kf64",
    "kon16",
    "kon17",
    "kon18",
    "kon73",
    "kon19",
    "kon20",
    "kon21",
    "kon74",
    "kon22",
    "kon23",
    "kon24",
    "kon25",
    "kon75",
    "kon26",
    "kon27",
    "kon28",
    "kon29",
    "kon76",
    "kon30",
    "kon31",
    "kon32",
    "kon33",
    "kon77",
    "kon34",
    "kon35",
    "kon36",
    "kon37",
    "kon78",
    "kon65",
    "kon66",
    "kon67",
    "kon68",
    "kon79",
    "kon69",
    "kon70",
    "kon71",
    "kon72",
    "kon80",
    "kon40",
    "kon41",
    "kon42",
    "kon43",
    "kon44",
    "kon45",
    "kon88",
    "kon46",
    "kon57",
    "kon58",
    "kon59",
    "kon60",
    # Nakakuki et al., Cell (2010)
    "V1",
    "Km1",
    "V2",
    "Km2",
    "V3",
    "Km3",
    "V4",
    "Km4",
    "V5",
    "Km5",
    "V6",
    "Km6",
    "KimERK",
    "KexERK",
    "KimpERK",
    "KexpERK",
    "KimppERK",
    "KexppERK",
    "V10",
    "Km10",
    "n10",
    "p11",
    "p12",
    "p13",
    "V14",
    "Km14",
    "V15",
    "Km15",
    "p16",
    "p17",
    "KimDUSP",
    "KexDUSP",
    "KimpDUSP",
    "KexpDUSP",
    "V20",
    "Km20",
    "V21",
    "Km21",
    "p22",
    "p23",
    "V24",
    "Km24",
    "V25",
    "Km25",
    "KimRSK",
    "KexRSK",
    "V27",
    "Km27",
    "V28",
    "Km28",
    "V29",
    "Km29",
    "V30",
    "Km30",
    "V31",
    "Km31",
    "n31",
    "p32",
    "p33",
    "p34",
    "V35",
    "Km35",
    "V36",
    "Km36",
    "V37",
    "Km37",
    "p38",
    "p39",
    "KimFOS",
    "KexFOS",
    "KimpcFOS",
    "KexpcFOS",
    "V42",
    "Km42",
    "V43",
    "Km43",
    "V44",
    "Km44",
    "p45",
    "p46",
    "p47",
    "m47",
    "p48",
    "p49",
    "m49",
    "p50",
    "p51",
    "m51",
    "p52",
    "m52",
    "p53",
    "p54",
    "m54",
    "p55",
    "p56",
    "m56",
    "V57",
    "Km57",
    "n57",
    "p58",
    "p59",
    "p60",
    "p61",
    "KimF",
    "KexF",
    "p63",
    "KF31",
    "nF31",
    "Vn",
    "Vc",
    ##
    "scale_E1",
    "scale_E2",
    "scale_E3",
    "scale_E4",
    "scale_G",
    #"scale_S",
    "a_SHC1",
    "a_SHC2",
    "a_SHC3",
    "a_SHC4",
    #"scale_I",
    "a_PIK3CA",
    "a_PIK3CB",
    "a_PIK3CD",
    "a_PIK3CG",
    # ver.5
    "scale_PTEN",
    #"scale_R",
    "a_RASA1",
    "a_RASA2",
    "a_RASA3",
    #"scale_O",
    "a_SOS1",
    "a_SOS2",
    "scale_A",
    #"scale_Akt",
    "a_AKT1",
    "a_AKT2",
    #"scale_RsD",
    "a_HRAS",
    "a_KRAS",
    "a_NRAS",
    #"scale_Raf",
    "a_ARAF",
    "a_BRAF",
    "a_RAF1",
    #"scale_MEK",
    "a_MAP2K1",
    "a_MAP2K2",
    "scale_T",
    "scale_CREBn",
    #"scale_ERKc",
    "a_MAPK1",
    "a_MAPK3",
    "scale_Elk1n",
    #"scale_RSKc",
    "a_RPS6KA1",
    "a_RPS6KA2",
    "a_RPS6KA3",
    "a_RPS6KA6",
]

for (idx,name) in enumerate(parameters)
    eval(Meta.parse("const $name = $idx"))
end

const n_parameters = length(parameters)

end  # module


text:
InputFolder = './Images/VesselsSpots/';
OutputFolder = './Results/Images/VesselsSpots/Spots/';
Fill = -1;

@iA = '*.tif';

@fxm_mFilterObjSize [iA] > [M];
params.MinArea = 1;
params.MaxArea = 1;
/endf

/show iA > M;
/keep M > tif;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER


text:
using ..ChemistryFeaturization.AbstractType: AbstractCodec
import ..ChemistryFeaturization.encode
import ..ChemistryFeaturization.decode

export encode, decode

"""
    SimpleCodec(encode_f, decode_f)

SimpleCodec is a most simple, atomic (as in indivisible), and straightforward Codec.

A SimpleCodec is characterized by only the encoding and decoding functions it has as fields.
This means that the encoding-decoding scheme is entirely dependent only upon these functions, and the fields of the FeatureDescriptor to which this Codec is attached. No external parameters (such as number of bins, etc.) are required
by the encoding-decoding implementations.

This Codec is very useful when the encoding-decoding scheme that you want to define for a `FeatureDescriptor` is very simple.
"""
struct SimpleCodec <: AbstractCodec
    encode_f::Function
    decode_f::Function
end

encode(sc::SimpleCodec, val) = sc.encode_f(val)
decode(sc::SimpleCodec, encoded) = sc.decode_f(encoded)


text:
using UncNLPrograms
using Test, CUTEst, LinearAlgebra, Printf
import NLPModels as jo

@testset "Translation Testing" begin
    include("translation_tests.jl")
end

@testset "Interface Testing" begin
    include("interface_tests.jl")
end


text:
using Cxx
using TimeZones

for dtype in [ :Time, :Duration, :Rate ]
    cppname = string("ros::", dtype)
    cxxtval = Expr(:macrocall, Symbol("@cxxt_str"), cppname)
    cxxtref = Expr(:macrocall, Symbol("@cxxt_str"), string(cppname, "&"))
    symtval = Symbol(string("Ros", dtype))
    symtref = Symbol(string("Ros", dtype, "R"))
    symtall = Symbol(string("Ros", dtype, "U"))

    @eval begin
        global const $(symtval) = $(cxxtval)
        global const $(symtref) = $(cxxtref)
        global const $(symtall) = Union{$(cxxtval), $(cxxtref)}
        export $(symtval), $(symtref) #, $(symtall)
    end

    cxxcstr = Expr(:macrocall, Symbol("@icxx_str"), string(cppname, "((double) \$t);"))
    @eval Base.convert(::Type{$(symtval)}, t::Float64) = $(cxxcstr)
end

function Base.convert(::Type{ZonedDateTime}, rt::RosTimeU)
    dt = Dates.unix2datetime(icxx"(long)($rt.toSec());")
    astimezone(ZonedDateTime(dt, tz"UTC"), localzone())
end

Base.get(rt::RosTimeU)::Float64              = icxx"$rt.toSec();"
     set(rt::RosTimeU, t::Float64)::RosTimeU = icxx"$rt.fromSec((double) $t);"

# wrapper for ros::spin, ros::spinOnce

mutable struct RosSpin
    should_continue::Bool
    task::Any
    rate::RosRate
    RosSpin()        = new(false, nothing, (@cxx ros::Rate(1)))
    RosSpin(r::Real) = new(false, nothing, (@cxx ros::Rate(r)))
end

function spin(s::RosSpin)
    while s.should_continue
        # println("RosSpin: ", s.should_continue)
        @cxx ros::spinOnce()
        try
            sleep(.001)
        catch ex
            if isa(ex, ErrorException)
                println("got exception: ", ex.msg)
                s.should_continue = false
                break
            end
        end
        @cxx s.rate->sleep()
    end
end

function Base.start(s::RosSpin)
    s.should_continue = true
    s.task = @async spin(s)
end

function stop(s::RosSpin)
    s.should_continue = false
    #println("RosSpin: ", s.should_continue)
end

export RosSpin, stop, set

max_stars_repo_name: moyhig/CxxROS

max_issues_repo_name: moyhig/CxxROS

max_forks_repo_name: moyhig/CxxROS


text:
module Day14

using ..ReTest, ..DocStringExtensions

import ..solve, ..PuzzleInfo, ..TEST_DATA_DIRECTORY

abstract type Day14PuzzleInfo <: PuzzleInfo end

"""
Day 14 part 1 puzzle information.

$(FIELDS)
"""
struct Part1PuzzleInfo <: Day14PuzzleInfo
    """Path to input data file."""
    datapath::String
end

"""
Day 14 part 2 puzzle information.

$(FIELDS)
"""
struct Part2PuzzleInfo <: Day14PuzzleInfo
    """Path to input data file."""
    datapath::String
end

"""
    $(FUNCTIONNAME)(info)

Get number of polymerisation reaction steps to take for puzzle part described by `info`.
"""
function getnumbersteps(::Day14PuzzleInfo) end
getnumbersteps(::Part1PuzzleInfo) = 10
getnumbersteps(::Part2PuzzleInfo) = 40

"""
    $(FUNCTIONNAME)(input)

Read polymer instructions from input stream `input`.
"""
function readpolymerinstructions(input::IO)
    polymertemplate = readuntil(input, "\n\n")
    pairinsertionrules = Dict{String, Char}()
    for line in eachline(input)
        pair, insert = split(line, " -> ")
        pairinsertionrules[pair] = only(insert)
    end
    polymertemplate, pairinsertionrules
end

function solve(info::Day14PuzzleInfo)
    polymer, rules = open(readpolymerinstructions, info.datapath)
    paircounts = Dict{String, Int}()
    for i in 1:(length(polymer) - 1)
        pair = polymer[i:i+1]
        paircounts[pair] = get(paircounts, pair, 0) + 1
    end
    processpair(p) = haskey(rules, p) ? (p[1] * rules[p], rules[p] * p[2]) : (p,)
    for step in 1:getnumbersteps(info)
        newpaircounts = Dict{String, Int}()
        for pair in keys(paircounts)
            for newpair in processpair(pair)
                newpaircounts[newpair] = get(newpaircounts, newpair, 0) + paircounts[pair]
            end
        end
        paircounts = newpaircounts
    end
    elementcounts = Dict{Char, Int}()
    for pair in keys(paircounts)
        elementcounts[pair[1]] = get(elementcounts, pair[1], 0) + paircounts[pair]
    end
    elementcounts[polymer[end]] = get(elementcounts, polymer[end], 0) + 1
    maximum(values(elementcounts)) - minimum(values(elementcounts))
end

@testset begin
    @test solve(Part1PuzzleInfo(joinpath(TEST_DATA_DIRECTORY, "day14.txt"))) == 1588
    @test solve(Part2PuzzleInfo(joinpath(TEST_DATA_DIRECTORY, "day14.txt"))) == 2188189693529
end

end # module


text:
"""
    function judgeequal(potential_unit::T, topo_name) where T <: Potential_Unit
This judges wether the topology (e.g bond) is same as `potential_unit`

# Example
- `judgeequal(Potential_Charmm36_Unit(["CTL3", "CL"], [200.0, 1.522]), ["CTL3", "CL"]) == true`
- `judgeequal(Potential_Charmm36_Unit(["CTL3", "CL"], [200.0, 1.522]), ["CL", "CTL3"]) == true`
- `judgeequal(Potential_Charmm36_Unit(["CTL3", "CL"], [200.0, 1.522]), ["CTL3", "CTL3"]) == false`
"""
function judgeequal(potential_unit::T, topo_name) where T <: Potential_Unit
    if length(potential_unit.flag) != length(topo_name)
        return false
    end

    flag_foward = true
    flag_backward = true
    for (id, atom) in enumerate(potential_unit.flag)
        if atom == "X"
            continue
        elseif atom == topo_name[id]
            continue
        else
            flag_foward = false
        end
    end

    for id = 1:length(potential_unit.flag)
        atom_now = potential_unit.flag[id]
        if atom_now == "X"
            continue
        elseif atom_now == topo_name[end-id+1]
            continue
        else
            flag_backward = false
        end
    end

    return flag_foward | flag_backward
end

function strvec2intvec(raw_str_vec)
    res = [0]
    for i in raw_str_vec
        append!(res, parse(Int, i))
    end
    return res[2:end]
end


text:
using ComplexSystems
using Plots

the_walkers = walker_set(5000,-5,5)
unshift!(the_walkers,walker(0.,0.,0))

the_walkers = dla(the_walkers,0.1,-5,5)

## Plotting
the_x = [w.x for w in the_walkers]
the_y = [w.y for w in the_walkers]
p = scatter(the_x,the_y,markersize=3.,
            leg=false,ticks=nothing,
            markerstrokewidth=0.,border=false)
title!("Diffusion Limit Aggregation")
png(p,"/figs/dla")

## Minkowski dimension
min_x=-5.
max_x=5.
min_y=-5.
max_y=5.

walkers_p = map(x -> Point(x.x,x.y),the_walkers);
Δs = []
Dimension = []
for Δ ∈ 0.01:0.05:0.96
    the_grid = make_square_grid(min_x,max_x,min_y,max_y,Δ,true);
    for punto ∈ walkers_p
        map(x -> inside_box(punto,x),the_grid);
    end
    non_empty_ones = sum([!x.empty for x ∈ the_grid]);
    println(non_empty_ones)
    dimension = log(non_empty_ones)/log(1/Δ)
    push!(Δs,Δ)
    push!(Dimension,dimension)
end

mini = round(Dimension[1];digits=3)
p = scatter(Δs,Dimension,legend=false)
title!("Dimension vs box size")
yaxis!("Dimension",:log10)
xaxis!("Length of square side",:log10)

png(p,"/figs/dimension")


text:
#=
Computational kernels for pseudospectra computations.

This file is part of Pseudospectra.jl.

Julia implementation
Copyright (c) 2017-2021 Ralph A. Smith

Portions derived from EigTool:
 Copyright (c) 2002-2014, The Chancellor, Masters and Scholars
 of the University of Oxford, and the EigTool Developers. All rights reserved.
 EigTool is maintained on GitHub:  https://github.com/eigtool

SPDX-License-Identifier: BSD-3-Clause
License-Filename: LICENSES/BSD-3-Clause_Eigtool
=#

# normally hardwired, but change to get test coverage w/o huge problems
# or to get reference solutions for comparison.
mutable struct ComputeThresholds
    # FIXME: undo mutability when tests etc. are updated
    minlancs4psa::Int # use SVD for n < this
    maxstdqr4hess::Int # use HessQR for n > this (in rectangular case)
    minnev::Int # number of ew's to acquire for projection
    maxit_lancs::Int # bound on Lanczos iterations
end
const _default_thresholds = ComputeThresholds(55,200,20,99)

# FIXME: temporary alias until tests etc. are updated
const psathresholds = _default_thresholds

"""
    psa_compute(T,npts,ax,eigA,opts,S=I) -> (Z,x,y,levels,info,Tproj,eigAproj,algo)

Compute pseudospectra of a (decomposed) matrix.

Uses a modified version of the code in [^Trefethen1999].
If the matrix `T` is upper triangular (e.g. from
a Schur decomposition) the solver is much more efficient than otherwise.

# Arguments
- `T`:      input matrix, usu. from `schur()`
- `npts`:   grid will have `npts × npts` nodes
- `ax`:     axis on which to plot `[min_real, max_real, min_imag, max_imag]`
- `eigA`:   eigenvalues of the matrix, usu. also produced by `schur()`. Pass
  an empty vector if unknown.
- `S`:    2nd matrix, if this is a generalized problem arising from an
  original rectangular matrix.
- `opts`: a `Dict{Symbol,Any}` holding options. Keys used here are as follows:

| Key                 | Type   | Default | Description |
|:-----------|:---------------|:-------------------------------------|:--------|
| `:levels`  | `Vector{Real}` | auto | `log10(ϵ)` for the desired ϵ levels |
| `:recompute_levels` | `Bool` | true | automatically recompute ϵ levels? |
| `:real_matrix`      | `Bool` | `eltype(A)<:Real` | is the original matrix real? (Portrait is symmetric if so.) This is needed because `T` could be complex even if `A` was real.|
| `:proj_lev`         | `Real` |  ∞ | The proportion by which to extend the axes in all directions before projection. If negative, exclude subspace of eigenvalues smaller than inverse fraction. ∞ means no projection.|
| `:scale_equal` | `Bool` | false | force the grid to be isotropic? |
| `:threaded` | `Bool` | false | distribute computation over Julia threads?

# Notes:
- Projection is only done for square, dense matrices.  Projection for sparse
  matrices may be handled (outside this function) by a Krylov method which
  reduces the matrix to a projected Hessenberg form before invoking
  `psa_compute`.
- This function does not compute generalized pseudospectra per se. They may
  be handled by pre- and post-processing.

# Outputs:
- `Z`:         the singular values over the grid
- `x`:         the x coordinates of the grid lines
- `y`:         the y coordinates of the grid lines
- `levels`:   the levels used for the contour plot (if automatically calculated)
- `Tproj`:     the projected matrix (an alias to `T` if no projection was done)
- `eigAproj`:  eigenvalues projected onto
- `algo`: a Symbol indicating which algorithm was used
- `info`:      flag indicating where automatic level creation fails:

| info | Meaning |
|:------|:--------|
| 0 |  No error |
|-1 |  No levels in range specified (either manually, or if matrix is too normal to show levels) |
|-2 |  Matrix is so non-normal that only zero singular values were found |
|-3 |  Computation cancelled |

[^Trefethen1999]: L.N.Trefethen, "Computation of pseudospectra," Acta Numerica 8, 247-295 (1999).
"""
function psa_compute(Targ, npts::Int, ax::Vector, eigA::Vector, opts::Dict, S=I;
                     psatol = 1e-5, thresholds=_default_thresholds,
                     proglog=nothing, logger=:default,
                     ctrlflag=nothing)
    # `proglog` is a placeholder in the hope that something like ProgressMeter can
    #    be made to work in a GUI
    # `ctrlflag` allows user to force early termination.
    m,n = size(Targ)
    eigAproj = copy(eigA) # default
    if isa(S,UniformScaling)
        ms,ns = 1,1
    else
        ms,ns = size(S)
    end
    comp_opts = Dict{Symbol,Any}()
    if !haskey(opts,:recompute_levels)
        comp_opts[:recompute_levels] = false
    end
    if haskey(opts,:levels)
        levels = opts[:levels]
        if length(levels) == 1
            levels = levels * ones(Int,2)
        end
    else
        levels = -8:-1
        comp_opts[:recompute_levels] = true
    end
    all_opts = merge(comp_opts, opts)

    proj_lev = get(all_opts,:proj_lev,Inf)
    re_calc_lev = all_opts[:recompute_levels]
    verbosity = get(all_opts,:verbosity,1)
    threaded = get(all_opts,:threaded,false)
    nservers = get(all_opts,:nservers, Threads.nthreads())

    if get(all_opts,:scale_equal,false)
        y_dist = ax[4]-ax[3]
        x_dist = ax[2]-ax[1]
        if x_dist > y_dist
            x_npts = npts
            y_npts = max(5,ceil(Int,y_dist/x_dist*npts))
        else
            y_npts = npts
            x_npts = max(5,ceil(Int,x_dist/y_dist*npts))
        end
    else
        x_npts = npts
        y_npts = npts
    end
    if get(all_opts,:real_matrix,eltype(Targ)<:Real) && ax[4] > 0 && ax[3] < 0
        y, n_mirror_pts = shift_axes(ax,y_npts)
    else
        n_mirror_pts = 0
        y = collect(range(ax[3], stop=ax[4], length=y_npts))
    end
    x = collect(range(ax[1], stop=ax[2], length=x_npts))
    lx = length(x) # why??
    ly = length(y)
    Z = ones(ly,lx) .+ Inf

    if !issparse(Targ) && n==m
        Tproj, eigAproj = _maybe_project(Targ, proj_lev, ax, eigA, thresholds, verbosity)
        m = size(Tproj,1)
    else
        # sparse or rectangular
        Tproj = Targ
    end # projection branch

    # compute resolvent norms

    local progmeter

    maxit = thresholds.maxit_lancs
    warnflags = falses(2)

    if issparse(Targ)
        algo = :sparse_direct
        Tproj = Targ
        # large value used when subspace eigenproblem doesn't converg
        bigσ = 0.1*floatmax(real(eltype(Targ)))

        if proglog === nothing
            progmeter = Progress(ly,1,"Computing pseudospectra...", 20)
        end
        # reverse order so first row is likely to have a complex gridpt
        # (better timing for LU)
        for j=ly:-1:1
            # check for stop/cancel
            if (ctrlflag !== nothing) && (ctrlflag[] == 1)
                return nothing
                # Eigtool also allows for pause.
            end

            # loop over points in x-direction
            for k=1:lx
                zpt = x[k] + y[j]*im
                t0 = time()
                σ = _psa_lanczos_sparse(Targ, S, zpt, maxit, bigσ)
                Z[j,k] = 1/sqrt(σ)

            end # for k=1:lx
            if proglog === nothing
                update!(progmeter,ly-j+1)
            end
        end
    else # matrix is dense
        if proglog === nothing
            progmeter = Progress(ly,1,"Computing pseudospectra...", 20)
        end
        # Following Eigtool, we divide into batches and check for cancellation order
        # at bulk intervals.
        # Eigtool also allows for pause.
        step = _get_step_size(m,ly,real(eltype(Tproj)))
        for j=ly:-step:1
            # check for stop/cancel
            if (ctrlflag !== nothing) && (ctrlflag[] == 1)
                return nothing
            end

            last_y = max(j-step+1,1)
            q = randn(n) + randn(n)*im
            q = q / norm(q)
            t0 = time()
            Z[j:-1:last_y,:],algo,warnflags = psacore(Tproj,S,q,x,y[j:-1:last_y], m-n+1;
                                                      tol=psatol,
                                                      threaded=threaded, nt=nservers,
                                                      warned=warnflags,
                                                      thresholds=_default_thresholds)

            if proglog === nothing
                update!(progmeter, ly-j+1)
            end
        end # ly loop
    end # if sparse/dense

    # map data (and y) if accounting for symmetry
    if n_mirror_pts < 0
        # bottom half is master
        Z = vcat(Z,reverse(Z[end+n_mirror_pts+1:end,:],dims=1))
        y = vcat(y,-reverse(y[end+n_mirror_pts+1:end]))
    else
        if y[1] != 0
            Z = vcat(reverse(Z[1:n_mirror_pts,:],dims=1),Z)
            y = vcat(-reverse(y[1:n_mirror_pts]),y)
        else
            Z = vcat(reverse(Z[2:n_mirror_pts+1,:],dims=1),Z)
            y = vcat(-reverse(y[2:n_mirror_pts+1]),y)
        end
    end
    ps_tiny = 10*sqrt(floatmin(eltype(Z)))
    (verbosity > 1) && println("range of Z: ",extrema(Z))
    clamp!(Z,ps_tiny,Inf)

    err = 0
    # maybe recalc levels
    if re_calc_lev
        levels,err = recalc_levels(Z,ax)
        if err != 0
            if err == -1
                @mywarn(logger,"Range too small---no contours to plot. Refine grid or zoom out.")
            elseif err == -2
                @mywarn(logger,"Matrix too non-normal---resolvent norm is "
                * "computationally infinite within current axes. Zoom out!")
            end
            return Z,x,y,levels,err,Tproj,eigAproj,algo
        end
    else
        # check that user-supplied levels will plot something
        if ((minimum(levels) > log10(maximum(Z)))
            | (maximum(levels) < log10(minimum(Z))))
            levels, err = recalc_levels(Z,ax)
            @mywarn(logger,"No contours to plot in requested range; 'Smart' levels used.")
            return Z,x,y,levels,err,Tproj,eigAproj,algo
        end
    # check range of Z
        if minimum(levels) < log10(ps_tiny)+1
            @mywarn(logger,"Smallest level allowed by machine precision reached; "
            * "levels may be inaccurate.")
        end
    end
    return Z,x,y,levels,err,Tproj,eigAproj,algo
end

function _get_step_size(n,ly,T)
    nsmall =  precision(T) <= 53 ? 8 : 20
    if n < 100
        step = max(1,floor(Int,ly/nsmall))
    else
        step = min(ly,max(1,floor(Int,4*ly/n)))
    end
    # upstream decreases by factor of 4 if fast implementation is missing
    return step
end

# Trefethen/Wright projection scheme:
# restrict to interesting subspace by ignoring eigenvectors whose
# eigenvalues lie outside rectangle around current axes
function _maybe_project(Targ, factor, ax, eigA, thresholds, verbosity)
    m,n = size(Targ)
    axis_w = ax[2]-ax[1]
    axis_h = ax[4]-ax[3]
    if factor >= 0
        proj_w = axis_w * factor
        proj_h = axis_h * factor
    else
        proj_size = -1 / factor
    end
    np = 0
    ew_range = ax
    # iteratively extend range until 20 (or all) ews are included
    if (m > thresholds.minnev) && !isempty(eigA)
        local selection
        while np < thresholds.minnev
            if factor >= 0
                ew_range = [ew_range[1] - proj_w, ew_range[2] + proj_w,
                            ew_range[3] - proj_h, ew_range[4] + proj_h]
                selection = findall((real(eigA) .> ew_range[1])
                                    .& (real(eigA) .< ew_range[2])
                                    .& (imag(eigA) .> ew_range[3])
                                    .& (imag(eigA) .< ew_range[4]))
            else
                selection = findall(abs.(eigA) .> proj_size)
                proj_size *= (1/2)
            end
            np = length(selection)
            if factor == 0
                # restrict to ews visible in window
                break
            end
        end

    else
        np = m
    end
    # if no need to project (all ews in range)
    if m == np
        m = size(Targ,1)
        # if !opts[:no_waitbar]
        # TODO: post waitbar
        # end
        eigAproj = copy(eigA)
        Tproj = Targ # no mutation, so just dup binding
    else
        if verbosity > 1
            println("projection reduces rank $m -> $np")
        end
        m = np
        n = np
        # restrict eigenvalues and matrix
        eigAproj = eigA[selection]

        # temporarily lose triangular structure
        Tproj = copy(Matrix(Targ))
        # if we have some eigenvalues in our window
        if m>0
            # TODO: post waitbar

            # do the projection
            for i=1:m
                for k=selection[i]-1:-1:i
                    G,r = givens(conj(Tproj[k,k+1]),
                                 conj(Tproj[k,k]-Tproj[k+1,k+1]),
                                 k+1,k)
                    rmul!(Tproj,adjoint(G))
                    lmul!(G,Tproj)
                end
                # TODO: update waitbar
                # TODO: check for pause ll 291ff
                # TODO: check for stop/cancel
            end
            Tproj = UpperTriangular(triu(Tproj[1:m,1:m]))
        end
    end
    return Tproj, eigAproj
end

struct _SVD_Server{MT}
    Twt::MT
end

function (s::_SVD_Server)(c::Channel{Tuple{Int,Int}}, Twork, x, y, Z)
    m,n = size(Twork)
    try
    while true
        j,k = take!(c)
        if j < 0
            break
        end
        zpt = x[k] + y[j]*im
        copy!(s.Twt, Twork)
        s.Twt[1:m+1:end] .-= zpt
        F = svd!(s.Twt)
        Z[j,k] = minimum(F.S)
    end
    catch JE
        println("exception in svd server: "); display(JE); println()
        return false
    end
    return true
end

struct _Sq_Lancs_Server{MT,HT}
    Twt::MT
    Hwt::HT
end

function (s::_Sq_Lancs_Server)(c::Channel{Tuple{Int,Int}},
                               diaga, q0, tol, maxit, bigσ, x, y, Z)
    m,n = size(s.Twt)
    convs = 0
    fails = 0
    try
    while true
        j,k = take!(c)
        if j < 0
            break
        end
        zpt = x[k] + y[j]*im
        s.Twt[1:m+1:end] .= diaga .- zpt
        F1 = UpperTriangular(s.Twt)
        σ, conv1, fail1 = _psa_lanczos!(s.Hwt, F1, q0, tol, maxit, bigσ)
        convs += conv1
        fails += fail1
        Z[j,k] = 1/sqrt(σ)
    end
    catch JE
        println("exception in sq lancs server: "); display(JE); println()
        return false, convs, fails
    end
    return true, convs, fails
end

"""
    psacore(T,S,q,x,y,bw;tol=1e-5,threaded=false) -> Z,algo,warninfo

Compute pseudospectra of a dense triangular matrix

# Arguments
- `T::Matrix{Number}`: long-triangular matrix whose pseudospectra to compute
- `S`: 2nd matrix from generalised pencil `zS-T`. Set to `I` if
           the problem is not generalised
- `q::Vector{Number}`: starting vector for the inverse-Lanczos iteration
           (the same vector is used to start each point in the
           grid defined by `x` and `y`). `q` **must be normalised to
           have unit length.**
- `x::Vector{Real}`: real-part grid to compute the pseudospectra over
- `y::Vector{Real}`: imaginary-part grid to compute the pseudospectra over
- `tol::Real=1e-5`:  tolerance to use to determine when to stop the
           inverse-Lanczos iteration
- `bw::Int`: lower bandwidth of the input matrix (2 for Hessenberg)
- `threaded::Bool`: whether to use multithreading

If `threaded` is `true`, computation of `Z`-values is distributed over
multiple threads. This is worthwhile if using extended precision or a
fine `Z` mesh.  If using BLAS element types, beware of oversubscription.

# Result
- `Z::Matrix{Real}`: the singular values corresponding to the grid points `x` and `y`.
- `algo::Symbol`: indicates algorithm used
- `warninfo::Vector{Bool}`: records whether warnings were issued
"""
function psacore(T, S, q0, x, y, bw; tol = 1e-5, threaded=false, nt=Threads.nthreads(),
                 warned=falses(2), thresholds=_default_thresholds)
    if isreal(T)
        Twork = T .+ complex(eltype(T))(0)
    else
        Twork = copy(T)
    end
    lx = length(x)
    ly = length(y)
    m,n = size(Twork)

    if m<n
        throw(ArgumentError("Matrix size must be m x n with m >= n"))
    end

    generalized =  !isa(S,UniformScaling)
    if generalized
        ms,ns = size(S)
        if (ms != m) || (ns != n)
            throw(ArgumentError("Dimension mismatch for S & T"))
        end
    end

    Z = zeros(ly,lx)
    diaga = diag(Twork)

    # large value used when subspace eigenproblem doesn't converg
    bigσ = 0.1*floatmax(real(eltype(T)))

    # for small matrices just use SVD
    if n < thresholds.minlancs4psa
        Twork = Matrix(Twork)
        if !generalized
            algo = :SVD
            if threaded
                svr = [_SVD_Server(similar(Twork)) for _ in 1:nt]
                c = Channel{Tuple{Int,Int}}(Inf)
                # prime the channel
                for k=1:lx
                    put!(c,(1,k))
                end
                ta = Task[]
                for it in 1:nt
                    t = Task(()->svr[it](c, Twork, x, y, Z))
                    t.sticky = false
                    schedule(t)
                    push!(ta,t)
                end
                for j=2:ly
                    for k=1:lx
                        put!(c,(j,k))
                    end
                end
                for k in 1:nt
                    put!(c,(-1,-1))
                end
                for k in 1:nt
                    r = fetch(ta[k])
                    if !isa(r, Bool) || !r
                        @warn "svd worker $k failed"
                    end
                end
            else # "serial" version
                for j=1:ly
                    for k=1:lx
                        zpt = x[k] + y[j]*im
                        Twork[1:m+1:end] .= diaga .- zpt
                        F = svd(Twork)
                        Z[j,k] = minimum(F.S)
                    end
                end
            end
        else
            algo = :SVD_gen
            for j=1:ly
                for k=1:lx
                    zpt = x[k] + y[j]*im
                    A = Twork .- zpt*S
                    F = svd!(A)
                    Z[j,k] = minimum(F.S)
                end
            end
        end
    else
        maxit = thresholds.maxit_lancs

        if (m==n) && threaded
            algo = :sq_lanc
            svr = [_Sq_Lancs_Server(copy(Twork),
                                    zeros(real(eltype(Twork)),maxit+1,maxit+1))
                   for _ in 1:nt
                   ]
            c = Channel{Tuple{Int,Int}}(Inf)
            # prime the channel
            for k=1:lx
                put!(c,(1,k))
            end
            ta = Task[]
            for it in 1:nt
                t = Task(()->svr[it](c, diaga, q0, tol, maxit, bigσ, x, y, Z))
                t.sticky = false
                schedule(t)
                push!(ta,t)
            end
            for j=2:ly
                for k=1:lx
                    put!(c,(j,k))
                end
            end
            for k in 1:nt
                put!(c,(-1,-1))
            end
            for k in 1:nt
                fetch(ta[k])
            end
        else # serial version
          H = zeros(real(eltype(Twork)),maxit+1,maxit+1)
          if m==n
            T1 = copy(Twork)
          end
          for j=1:ly
            for k=1:lx
                zpt = x[k]+y[j]*im
                if m != n
                    if !generalized
                        Twork[1:m+1:end] = diaga .- zpt
                        T1 = copy(Twork)
                    else
                        T1 = Twork - zpt*S
                    end
                    # for large rectangular Hessenberg, use HessQR algorithm
                    if (bw == 2) && (m > thresholds.maxstdqr4hess)
                        algo = :HessQR
                        for jj=1:n-1
                            # DEVNOTE: not using A_mul_B!(G,T1)
                            # because we don't want to mutate top of T1
                            # and views can be expensive
                            G,r = givens(T1[jj,jj],T1[jj+1,jj],1,2)
                            T1[jj:jj+1,jj:end] = G * T1[jj:jj+1,jj:end]
                        end
                    else
                        Qtmp,T1 = qr(T1)
                        algo = generalized ? :rect_qz : :rect_qr
                    end
                    T1 = triu(T1[1:n,1:n])
                else # square
                    algo = :sq_lanc
                    T1[1:m+1:end] .= diaga .- zpt
                end
                if !istriu(T1)
                    F1 = factorize(T1)
                else
                    F1 = UpperTriangular(T1)
                end
                σ, conv1, fail1 = _psa_lanczos!(H, F1, q0, tol, maxit, bigσ)
                if (!conv1) && (!warned[1])
                    @warn "Lanczos convergence failure(s) while computing resolvent norms"
                    warned[1] = true
                end
                if fail1 && (!warned[2])
                    @warn "Eigenvalue convergence failure(s) while computing resolvent norms"
                    warned[2] = true
                end
                Z[j,k] = 1/sqrt(σ)
            end
          end
        end
    end # svd/lanczos branch
    return Z,algo,warned
end

function _psa_lanczos!(H, F1, q0, tol, maxit, bigσ)
    m,n = size(F1)
    # q0 may be too long because of projection
    q = q0[1:n]
    qold = fill!(similar(q),zero(eltype(q)))
    β = 0.0
    σold = 0.0
    local σ
    lancz_converged = false
    H_eigs_failed = false
    for l=1:maxit
        v = (F1 \ (F1' \ q)) - β * qold
        α = real(dot(q,v)) # (q' * v)
        v .= v .- α .* q
        β = norm(v)
        copy!(qold,q)
        q .= v ./ β
        H[l+1,l] = β
        H[l,l+1] = β
        H[l,l] = α
        try
            ewp = eigvals(H[1:l,1:l])
            # eigvals may return complex eltype even if actually real
            σ = maximum(real.(ewp))
            # if !all(isreal.(ewp))
            #     @warn "psa-lancs: eigval anomaly $ewp"
            #     # Should we ask users to report this?
            # end
        catch JE
            # We want a fallback for convergence failure, but throw in
            # other cases.
            # WARNING: this is fragile, depends on library internals
            if !isa(JE, LinearAlgebra.LAPACKException)
                rethrow(JE)
            end
            H_eigs_failed = true
            σ = bigσ
            break
        end
        if (abs(σold / σ - 1) < tol || β == 0)
            lancz_converged = true
            break
        end
        σold = σ
    end
    return σ, lancz_converged, H_eigs_failed
end

function _psa_lanczos_sparse(Targ, S, zpt, maxit, bigσ)
    m,n = size(Targ)
    Tc = complex(eltype(Targ))

    F = lu(Targ - zpt*S)
    σold = 0
    qold = zeros(m)
    β = 0
    H = zeros(Tc,1,1)
    q = normalize!(randn(n) + randn(n)*im)
    w = similar(q)
    v = similar(q)
    local σ
    for l=1:maxit
        ldiv!(w,F,q)
        ldiv!(v,adjoint(F),w)
        v = v - β * qold
        α = real(dot(q,v))
        v = v - α * q
        β = norm(v)
        qold = q
        q = v * (1 / β)
        Hold = H
        H = zeros(Tc,l+1,l+1)
        copyto!(view(H,1:l,1:l),Hold)
        H[l+1,l] = β
        H[l,l+1] = β
        H[l,l] = α
        # calculate eigenvalues of H
        # if error is too big, just set a large value
        try
            ew = eigvals(H[1:l,1:l])
            σ = maximum(real.(ew))
        catch JE
            σ = bigσ
            break
        end
        if (abs(σold / σ - 1)) < 1e-3
            break
        end
        σold = σ
    end
    return σ
end


text:
using NamedDims
using Test

const testfiles = (
    "name_core.jl",
    "wrapper_array.jl",
    "functions.jl",
    "functions_dims.jl",
    "functions_math.jl",
    "broadcasting.jl",
    "tracker_compat.jl",
)

@testset "NamedDims.jl" begin
    @testset "$file" for file in testfiles
        include(file)
    end
end


text:
using LinearAlgebraicRepresentation
using Plasm
using SparseArrays
Lar = LinearAlgebraicRepresentation
Lara = Lar.Arrangement

# Data Reading
V = [
	0.0 1.5 3.0 0.0 3.0 0.0 1.5 3.0;
	0.0 0.0 0.0 1.5 1.5 3.0 3.0 3.0
	];
EV = [
	[1, 3], [4, 5],	# Orizontal
	[1, 7], [2, 8],	# Diag 1
	[2, 6], [3, 7]	# Diag 2
	];

## 1 - Input Reading
Plasm.view(Plasm.numbering(0.5)((V,[[[k] for k=1:size(V,2)], EV])));

# Planar Arrangement 1
W = convert(Lar.Points, V'); # Infering type for W = V'
copEV = Lar.coboundary_0(EV::Lar.Cells);
W1, copEV1 = Lar.planar_arrangement_1(W::Lar.Points, copEV::Lar.ChainOp);

## 2 - 2-cells fragmentation
EV1 = Lar.cop2lar(copEV1);
V1 = convert(Lar.Points, W1');
Plasm.view(Plasm.numbering(0.5)((V1,[[[k] for k=1:size(V1,2)], EV1])));

# Biconnected COmponent Evaluation
bicon_comps = Lar.Arrangement.biconnected_components(copEV1);

## 3 - Biconnected Components
hpcs = [ Plasm.lar2hpc(V1,[EV1[e] for e in comp]) for comp in bicon_comps ]
Plasm.view([
	Plasm.color(Plasm.colorkey[(k%12)==0 ? 12 : k%12])(hpcs[k])
	for k = 1 : (length(hpcs))
])

# computation of 2-cells and 2-boundary
W2, copEV2, copFE2 = Lar.planar_arrangement_2(W1, copEV1, bicon_comps)

## 4 - 3-cells identification & dangling 1-cells elimination
Plasm.view( Plasm.numbering1(0.5)((W2, copEV2, copFE2)) )

# 5 - Colorfull Representation
triangulated_faces = Lar.triangulate2D(W2, [copEV2, copFE2])
V2 = convert(Lar.Points, W2')
FVs2 = convert(Array{Lar.Cells}, triangulated_faces)
Plasm.viewcolor(V2::Lar.Points, FVs2::Array{Lar.Cells})

# polygonal face boundaries
EVs2 = Lar.FV2EVs(copEV2, copFE2)
EVs2 = convert(Array{Array{Array{Int64,1},1},1}, EVs2)
Plasm.viewcolor(V2::Lar.Points, EVs2::Array{Lar.Cells})

# 6 - Exploded Representation
model = V2,EVs2
Plasm.view(Plasm.lar_exploded(model)(1.2,1.2,1.2))


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# =================================================================================================
"""
    tess_init(
        inst::TessInst,
        languages::AbstractString = "eng",
        dataPath::AbstractString = "tessdata"
    )::Bool

Initialize the instance for the specified language(s).  Returns `false` if there was an error.

__Arguments:__

| T | Name      | Default    | Description
|:--| :-------- | :--------- | :----------
| R | inst      |            | The instance to initialize.
| O | languages | `eng`      | The language(s) to load.
| O | dataPath  | `tessdata` | The directory to look for the language files in.

__Details:__

This method can be called multiple times to reinitialize the langauges to OCR with.  Multiple
langagues can be specified by seperating them with a plus.  So if you want english and spanish you
could specify "eng+spa".  The language codes are (usually) the
[ISO 639-3](https://en.wikipedia.org/wiki/ISO_639-3) code.

Note: The language files are NOT automatically downloaded.  If you do not have them installed via
alternate means you can download them from https://github.com/tesseract-ocr/tessdata_best.

__Example:__

```jldoctest
julia> using Tesseract

julia> download_languages("eng+spa")
true

julia> instance = TessInst()
Allocated Tesseract instance.

julia> tess_init(instance, "spa")
true
```

See also: [`TessInst()`](@ref)
"""
function tess_init(
            inst::TessInst,
            languages::AbstractString = "eng",
            dataPath::AbstractString = TESS_DATA
        )::Bool
    local result = -1

    if is_valid(inst) == false
        @error "Instance has been freed."
        return false
    end

    if isempty(dataPath) && isempty(languages)
        result = ccall(
            (:TessBaseAPIInit3, TESSERACT),
            Cint,
            (Ptr{Cvoid}, Cstring, Cstring),
            inst,
            C_NULL,
            C_NULL
        )
    elseif isempty(languages)
        result = ccall(
            (:TessBaseAPIInit3, TESSERACT),
            Cint,
            (Ptr{Cvoid},Cstring,Cstring),
            inst,
            dataPath,
            C_NULL
        )
    elseif isempty(dataPath)
        result = ccall(
            (:TessBaseAPIInit3, TESSERACT),
            Cint,
            (Ptr{Cvoid},Cstring,Cstring),
            inst,
            C_NULL,
            languages
        )
    else
        result = ccall(
            (:TessBaseAPIInit3, TESSERACT),
            Cint,
            (Ptr{Cvoid},Cstring,Cstring),
            inst,
            dataPath,
            languages
        )
    end

    return result == 0
end

# =================================================================================================
"""
    tess_initialized_languages(
        inst::TessInst
    )::Union{String, Nothing}

Retrieve the last initialized language(s).  Returns `nothing` if there was an error.

__Arguments:__

| T | Name     | Default | Description
|:--| :------- | :------ | :----------
| R | inst     |         | The instance to get the languages from.

__Details:__

This method returns the language string provided in the last `tess_init()` call.

__Example:__

```jldoctest
julia> using Tesseract

julia> download_languages("eng+fra")
true

julia> instance = TessInst("eng+fra")
Allocated Tesseract instance.

julia> tess_initialized_languages(instance)
"eng+fra"
```

See also: [`tess_init`](@ref)
"""
function tess_initialized_languages(
            inst::TessInst
        )::Union{String, Nothing}

    # ---------------------------------------------------------------------------------------------
    # Make sure the Tess object is not freed.
    if is_valid(inst) == false
        @error "Instance has been freed."
        return nothing
    end

    # ---------------------------------------------------------------------------------------------
    # Call the C library.
    local retval = ccall(
        (:TessBaseAPIGetInitLanguagesAsString, TESSERACT),
        Cstring,
        (Ptr{Cvoid}, ),
        inst
    )

    # ---------------------------------------------------------------------------------------------
    # Nothing was returned.
    if retval == C_NULL
        return nothing
    end

    # ---------------------------------------------------------------------------------------------
    # Copy the string into Julia, DON'T free it.
    return unsafe_string(retval)
end

# =================================================================================================
"""
    tess_loaded_languages(
        inst::TessInst
    )::Union{Vector{String}, Nothing}

Get the the list of languages loaded into the OCR engine.  Returns `nothing` if there was an error.

__Arguments:__

| T | Name     | Default | Description
|:--| :------- | :------ | :----------
| R | inst     |         | The instance to get the languages from.

__Details:__

This method returns the language loaded into the Tesseract engine.  Some language files will load
additional languages.  Unlike `tess_initialized_languages()` this method will return all the loaded
languages not just the ones it was told to load by the client.

__Example:__

```jldoctest
julia> using Tesseract

julia> download_languages("eng+spa")
true

julia> instance = TessInst("eng+spa")
Allocated Tesseract instance.

julia> tess_loaded_languages(instance)
2-element Array{String,1}:
 "eng"
 "spa"
```

See also: [`tess_init`](@ref)
"""
function tess_loaded_languages(
            inst::TessInst
        )::Union{Vector{String}, Nothing}

    # ---------------------------------------------------------------------------------------------
    # Make sure the Tess object is not freed.
    if is_valid(inst) == false
        @error "Instance has been freed."
        return nothing
    end

    # ---------------------------------------------------------------------------------------------
    # Call the C library.
    local languages = ccall(
        (:TessBaseAPIGetLoadedLanguagesAsVector, TESSERACT),
        Ptr{Ptr{UInt8}},
        (Ptr{Cvoid}, ),
        inst
    )

    if languages == C_NULL
        return nothing
    end

    # ---------------------------------------------------------------------------------------------
    # Create the vector.
    local retval = Vector{String}()

    local pos = 1
    while unsafe_load(languages, pos) != C_NULL
        local language = unsafe_load(languages, pos)
        push!(retval, unsafe_string(language))
        pos += 1
    end

    # ---------------------------------------------------------------------------------------------
    # Free the array of strings na dreturn the vector.
    delete_array(languages)

    return retval
end

# =================================================================================================
"""
    tess_available_languages(
        inst::TessInst
    )::Union{Vector{String}, Nothing}

Get the list of available languages that can be loaded.  Returns `nothing` if there was an error.

__Arguments:__

| T | Name     | Default | Description
|:--| :------- | :------ | :----------
| R | inst     |         | The instance to query for the available languages.

__Details:__

Get the list of languages that can be used in the tess_init() function.  This method only returns
something AFTER `tess_init()`` has been called.

__Example:__

```jldoctest
julia> using Tesseract

julia> download_languages("eng+fra+spa")
true

julia> instance = TessInst("eng")
Allocated Tesseract instance.

julia> tess_available_languages(instance)
3-element Array{String,1}:
 "eng"
 "fra"
 "spa"
```

See also: [`tess_init`](@ref), [`update_languages`](@ref), [`download_languages`](@ref)
"""
function tess_available_languages(
            inst::TessInst
        )::Union{Vector{String}, Nothing}

    # ---------------------------------------------------------------------------------------------
    # Make sure the Tess object is not freed.
    if is_valid(inst) == false
        @error "Instance has been freed."
        return nothing
    end

    # ---------------------------------------------------------------------------------------------
    # Call the C library.
    local languages = ccall(
        (:TessBaseAPIGetAvailableLanguagesAsVector, TESSERACT),
        Ptr{Ptr{UInt8}},
        (Ptr{Cvoid}, ),
        inst
    )

    if languages == C_NULL
        return nothing
    end

    # ---------------------------------------------------------------------------------------------
    # Create the vector.
    local retval = Vector{String}()
    local pos    = 1

    while unsafe_load(languages, pos) != C_NULL
        local language = unsafe_load(languages, pos)
        push!(retval, unsafe_string(language))
        pos += 1
    end

    # ---------------------------------------------------------------------------------------------
    # Free the array of strings na dreturn the vector.
    delete_array(languages)

    return retval
end

# =========================================================================================
"""
    tess_params(
        inst::TessInst,
        filename::AbstractString
    )::Bool

Print out all the parameters with their values and help text to the specified file.
Returns false if there was an error.

__Arguments:__

| T | Name     | Default | Description
|:--| :------- | :------ | :----------
| R | inst     |         | The Tesseract instance to get get the parameters from.
| R | filename |         | The filename to write to.

__Details:__

For each parameter this method prints out it's name, it's value, and some descriptive text
about the variable.  Each variable is on it's own line with a tab character seperating each
value.

See also: [`tess_params_parsed`](@ref), [`tess_get_param`](@ref), [`tess_set_param`](@ref)
"""
function tess_params(
            inst::TessInst,
            filename::AbstractString
        )::Bool

    # ---------------------------------------------------------------------------------------------
    # Make sure the Tess object is not freed.
    if is_valid(inst) == false
        @error "Instance has been freed."
        return false
    end

    # ---------------------------------------------------------------------------------------------
    # Make the call.
    local retval = ccall(
        (:TessBaseAPIPrintVariablesToFile, TESSERACT),
        Cint,
        (Ptr{Cvoid}, Cstring),
        inst,
        filename
    )

    return retval == 1
end

# =========================================================================================
"""
    tess_params(
        inst::TessInst,
        stream::IO
    )::Bool

Print out all the parameters with their values and help text to the specified stream.
Returns false if there was an error.

__Arguments:__

| T | Name   | Default | Description
|:--| :----- | :------ | :----------
| R | inst   |         | The Tesseract instance to get get the parameters from.
| R | stream |         | The stream to write the files to.

__Details:__

For each parameter this method prints out it's name, it's value, and some descriptive text
about the variable.  Each variable is on it's own line with a tab character seperating each
value.

See also: [`tess_params_parsed`](@ref), [`tess_get_param`](@ref), [`tess_set_param`](@ref)
"""
function tess_params(
            inst::TessInst,
            stream::IO
        )::Bool
    local result = false

    # ---------------------------------------------------------------------------------------------
    # Make sure the Tess object is not freed.
    if is_valid(inst) == false
        @error "Instance has been freed."
        return false
    end

    # ---------------------------------------------------------------------------------------------
    # Create a temporary file to write to.
    local filename, io = mktemp(;cleanup = false)

    try
        close(io)

        if tess_params(inst, filename)
            local data = read(filename)
            write(stream, data)
            result = true
        end
    finally
        rm(filename)
    end

    return result
end

# =========================================================================================
"""
    tess_params(
        inst::TessInst
    )::Union{String, Nothing}

Print out all the parameters with their values and help text to a string.  Returns
`nothing` if there was an error.

__Arguments:__

| T | Name | Default | Description
|:--| :--- | :------ | :----------
| R | inst |         | The Tesseract instance to get the parameters from.

__Details:__

The return string will contain multiple lines, each line contains the name of a variable,
it's value, and some descriptive text about the variable.  The fields are separated by
tabs.

See also: [`tess_params_parsed`](@ref), [`tess_get_param`](@ref), [`tess_set_param`](@ref)
"""
function tess_params(
            inst::TessInst
        )::Union{String, Nothing}
    local result = nothing

    # ---------------------------------------------------------------------------------------------
    # Make sure the Tess object is not freed.
    if is_valid(inst) == false
        @error "Instance has been freed."
        return nothing
    end

    # ---------------------------------------------------------------------------------------------
    # Create a temporary file to write to.
    local filename, io = mktemp(;cleanup = false)

    try
        close(io)

        if tess_params(inst, filename)
            result = read(filename, String)
        end
    finally
        rm(filename)
    end

    return result
end

# =================================================================================================
"""
    tess_read_config(
        inst::TessInst,
        filename::AbstractString
    )::Nothing

Load configuration settings from a file into the Tesseract instance.

__Arguments:__

| T | Name     | Default | Description
|:--| :------- | :------ | :----------
| R | inst     |         | The Tesseract instance to load the settings into.
| R | filename |         | The name of the file to load the settings from.

See also: [`tess_read_debug_config`](@ref)
"""
function tess_read_config(
            inst::TessInst,
            filename::AbstractString
        )::Nothing

    # ---------------------------------------------------------------------------------------------
    # Make sure the Tess object is not freed.
    if is_valid(inst) == false
        @error "Instance has been freed."
        return nothing
    end

    # ---------------------------------------------------------------------------------------------
    # Call the C library.
    ccall(
        (:TessBaseAPIReadConfigFile, TESSERACT),
        Cvoid,
        (Ptr{Cvoid}, Cstring),
        inst,
        filename
    )
    nothing
end

# =================================================================================================
"""
    tess_read_debug_config(
        inst::TessInst,
        filename::AbstractString
    )::Nothing

Load debug configuration settings from a file into the Tesseract instance.

__Arguments:__

| T | Name     | Default | Description
|:--| :------- | :------ | :----------
| R | inst     |         | The Tesseract instance to load the settings into.
| R | filename |         | The name of the file to load the settings from.

__Details:__

Only the debug settings will be loaded, all other settings will be ignored.

See also: [`tess_read_config`](@ref)
"""
function tess_read_debug_config(
            inst::TessInst,
            filename::AbstractString
        )::Nothing

    # ---------------------------------------------------------------------------------------------
    # Make sure the Tess object is not freed.
    if is_valid(inst) == false
        @error "Instance has been freed."
        return nothing
    end

    # ---------------------------------------------------------------------------------------------
    # Call the C library.
    ccall(
        (:TessBaseAPIReadDebugConfigFile, TESSERACT),
        Cvoid,
        (Ptr{Cvoid}, Cstring),
        inst,
        filename
    )
    nothing
end


text:
# Automatically generated using Clang.jl wrap_c, version 0.0.0

using Compat

# const POLYLIB_BITS = 64
# const LONG_LONG_MAX = __LONG_LONG_MAX__
# const LONG_LONG_MIN = -LONG_LONG_MAX - 1
# const ULONG_LONG_MAX = LONG_LONG_MAX * 2 + 1
# const LINEAR_VALUE_STRING = "long int"
# const VALUE_FMT = "%ld"

# # Skipping MacroDefinition: VALUE_CONST ( val ) ( val ## L )

# const VALUE_NAN = LONG_MIN
# const VALUE_MIN = LONG_MIN + 1
# const VALUE_MAX = LONG_MAX

# # Skipping MacroDefinition: VALUE_SQRT_MIN int_to_value ( INT_MIN )
# # Skipping MacroDefinition: VALUE_SQRT_MAX int_to_value ( INT_MAX )

# const VALUE_ZERO = 0
# const VALUE_ONE = 1
# const VALUE_MONE = -1

# # Skipping MacroDefinition: VALUE_TO_LONG ( val ) ( val )
# # Skipping MacroDefinition: VALUE_TO_INT ( val ) ( ( int ) ( val ) )
# # Skipping MacroDefinition: VALUE_TO_FLOAT ( val ) ( ( float ) ( val ) )
# # Skipping MacroDefinition: VALUE_TO_DOUBLE ( val ) ( ( double ) ( val ) )
# # Skipping MacroDefinition: value_init ( val ) ( ( val ) = 0 )
# # Skipping MacroDefinition: value_assign ( v1 , v2 ) ( ( v1 ) = ( v2 ) )
# # Skipping MacroDefinition: value_set_si ( val , i ) ( ( val ) = ( Value ) ( i ) )
# # Skipping MacroDefinition: value_set_double ( val , d ) ( ( val ) = ( Value ) ( d ) )
# # Skipping MacroDefinition: value_clear ( val ) ( ( val ) = 0 )
# # Skipping MacroDefinition: value_read ( val , str ) ( sscanf ( ( str ) , VALUE_FMT , & ( val ) ) )
# # Skipping MacroDefinition: value_print ( Dst , fmt , val ) ( fprintf ( ( Dst ) , ( fmt ) , ( val ) ) )
# # Skipping MacroDefinition: value_swap ( v1 , v2 ) { Value tmp ; tmp = v2 ; v2 = v1 ; v1 = tmp ; }
# # Skipping MacroDefinition: int_to_value ( i ) ( ( Value ) ( i ) )
# # Skipping MacroDefinition: long_to_value ( l ) ( ( Value ) ( l ) )
# # Skipping MacroDefinition: float_to_value ( f ) ( ( Value ) ( f ) )
# # Skipping MacroDefinition: double_to_value ( d ) ( ( Value ) ( d ) )
# # Skipping MacroDefinition: value_eq ( v1 , v2 ) ( ( v1 ) == ( v2 ) )
# # Skipping MacroDefinition: value_ne ( v1 , v2 ) ( ( v1 ) != ( v2 ) )
# # Skipping MacroDefinition: value_gt ( v1 , v2 ) ( ( v1 ) > ( v2 ) )
# # Skipping MacroDefinition: value_ge ( v1 , v2 ) ( ( v1 ) >= ( v2 ) )
# # Skipping MacroDefinition: value_lt ( v1 , v2 ) ( ( v1 ) < ( v2 ) )
# # Skipping MacroDefinition: value_le ( v1 , v2 ) ( ( v1 ) <= ( v2 ) )
# # Skipping MacroDefinition: value_abs_eq ( v1 , v2 ) ( value_abs ( v1 ) == value_abs ( v2 ) )
# # Skipping MacroDefinition: value_abs_ne ( v1 , v2 ) ( value_abs ( v1 ) != value_abs ( v2 ) )
# # Skipping MacroDefinition: value_abs_gt ( v1 , v2 ) ( value_abs ( v1 ) > value_abs ( v2 ) )
# # Skipping MacroDefinition: value_abs_ge ( v1 , v2 ) ( value_abs ( v1 ) >= value_abs ( v2 ) )
# # Skipping MacroDefinition: value_abs_lt ( v1 , v2 ) ( value_abs ( v1 ) < value_abs ( v2 ) )
# # Skipping MacroDefinition: value_abs_le ( v1 , v2 ) ( value_abs ( v1 ) <= value_abs ( v2 ) )
# # Skipping MacroDefinition: value_sign ( v ) ( value_eq ( v , VALUE_ZERO ) ? 0 : value_lt ( v , VALUE_ZERO ) ? - 1 : 1 )
# # Skipping MacroDefinition: value_compare ( v1 , v2 ) ( value_eq ( v1 , v2 ) ? 0 : value_lt ( v1 , v2 ) ? - 1 : 1 )
# # Skipping MacroDefinition: value_plus ( v1 , v2 ) ( ( v1 ) + ( v2 ) )
# # Skipping MacroDefinition: value_div ( v1 , v2 ) ( ( v1 ) / ( v2 ) )
# # Skipping MacroDefinition: value_mod ( v1 , v2 ) ( ( v1 ) % ( v2 ) )
# # Skipping MacroDefinition: value_direct_multiply ( v1 , v2 ) ( ( v1 ) * ( v2 ) )
# # Skipping MacroDefinition: value_minus ( v1 , v2 ) ( ( v1 ) - ( v2 ) )
# # Skipping MacroDefinition: value_pdiv ( v1 , v2 ) ( DIVIDE ( ( v1 ) , ( v2 ) ) )
# # Skipping MacroDefinition: value_pmod ( v1 , v2 ) ( MODULO ( ( v1 ) , ( v2 ) ) )
# # Skipping MacroDefinition: value_min ( v1 , v2 ) ( value_le ( ( v1 ) , ( v2 ) ) ? ( v1 ) : ( v2 ) )
# # Skipping MacroDefinition: value_max ( v1 , v2 ) ( value_ge ( ( v1 ) , ( v2 ) ) ? ( v1 ) : ( v2 ) )
# # Skipping MacroDefinition: value_or ( v1 , v2 ) ( ( v1 ) | ( v2 ) )
# # Skipping MacroDefinition: value_and ( v1 , v2 ) ( ( v1 ) & ( v2 ) )
# # Skipping MacroDefinition: value_lshift ( v1 , v2 ) ( ( v1 ) << ( v2 ) )
# # Skipping MacroDefinition: value_rshift ( v1 , v2 ) ( ( v1 ) >> ( v2 ) )
# # Skipping MacroDefinition: value_addto ( ref , val1 , val2 ) ( ( ref ) = ( val1 ) + ( val2 ) )
# # Skipping MacroDefinition: value_add_int ( ref , val , vint ) ( ( ref ) = ( val ) + ( Value ) ( vint ) )
# # Skipping MacroDefinition: value_addmul ( ref , val1 , val2 ) ( ( ref ) += ( val1 ) * ( val2 ) )
# # Skipping MacroDefinition: value_increment ( ref , val ) ( ( ref ) = ( val ) + VALUE_ONE )
# # Skipping MacroDefinition: value_direct_product ( ref , val1 , val2 ) ( ( ref ) = ( val1 ) * ( val2 ) )
# # Skipping MacroDefinition: value_multiply ( ref , val1 , val2 ) ( ( ref ) = value_mult ( ( val1 ) , ( val2 ) ) )
# # Skipping MacroDefinition: value_subtract ( ref , val1 , val2 ) ( ( ref ) = ( val1 ) - ( val2 ) )
# # Skipping MacroDefinition: value_sub_int ( ref , val , vint ) ( ( ref ) = ( val ) - ( Value ) ( vint ) )
# # Skipping MacroDefinition: value_decrement ( ref , val ) ( ( ref ) = ( val ) - VALUE_ONE )
# # Skipping MacroDefinition: value_division ( ref , val1 , val2 ) ( ( ref ) = ( val1 ) / ( val2 ) )
# # Skipping MacroDefinition: value_divexact ( ref , val1 , val2 ) ( ( ref ) = ( val1 ) / ( val2 ) )
# # Skipping MacroDefinition: value_modulus ( ref , val1 , val2 ) ( ( ref ) = ( val1 ) % ( val2 ) )
# # Skipping MacroDefinition: value_pdivision ( ref , val1 , val2 ) ( ( ref ) = value_pdiv ( ( val1 ) , ( val2 ) ) )
# # Skipping MacroDefinition: value_pmodulus ( ref , val1 , val2 ) ( ( ref ) = value_pmod ( ( val1 ) , ( val2 ) ) )
# # Skipping MacroDefinition: value_oppose ( ref , val ) ( ( ref ) = value_uminus ( ( val ) ) )
# # Skipping MacroDefinition: value_absolute ( ref , val ) ( ( ref ) = value_abs ( ( val ) ) )
# # Skipping MacroDefinition: value_minimum ( ref , val1 , val2 ) ( ( ref ) = value_min ( ( val1 ) , ( val2 ) ) )
# # Skipping MacroDefinition: value_maximum ( ref , val1 , val2 ) ( ( ref ) = value_max ( ( val1 ) , ( val2 ) ) )
# # Skipping MacroDefinition: value_gcd ( ref , val1 , val2 ) Gcd ( ( val1 ) , ( val2 ) , & ( ref ) )
# # Skipping MacroDefinition: value_lcm ( ref , val1 , val2 ) Lcm3 ( ( val1 ) , ( val2 ) , & ( ref ) )
# # Skipping MacroDefinition: value_orto ( ref , val1 , val2 ) ( ( ref ) = ( val1 ) | ( val2 ) )
# # Skipping MacroDefinition: value_andto ( ref , val1 , val2 ) ( ( ref ) = ( val1 ) & ( val2 ) )
# # Skipping MacroDefinition: value_uminus ( val ) ( - ( val ) )
# # Skipping MacroDefinition: value_not ( val ) ( ~ ( val ) )
# # Skipping MacroDefinition: value_abs ( val ) ( value_posz_p ( val ) ? ( val ) : ( value_ne ( ( val ) , VALUE_NAN ) ? value_uminus ( val ) : ( THROW ( overflow_error ) , VALUE_NAN ) ) )
# # Skipping MacroDefinition: value_pos_p ( val ) value_gt ( val , VALUE_ZERO )
# # Skipping MacroDefinition: value_neg_p ( val ) value_lt ( val , VALUE_ZERO )
# # Skipping MacroDefinition: value_posz_p ( val ) value_ge ( val , VALUE_ZERO )
# # Skipping MacroDefinition: value_negz_p ( val ) value_le ( val , VALUE_ZERO )
# # Skipping MacroDefinition: value_zero_p ( val ) value_eq ( val , VALUE_ZERO )
# # Skipping MacroDefinition: value_notzero_p ( val ) value_ne ( val , VALUE_ZERO )
# # Skipping MacroDefinition: value_one_p ( val ) value_eq ( val , VALUE_ONE )
# # Skipping MacroDefinition: value_notone_p ( val ) value_ne ( val , VALUE_ONE )
# # Skipping MacroDefinition: value_mone_p ( val ) value_eq ( val , VALUE_MONE )
# # Skipping MacroDefinition: value_notmone_p ( val ) value_ne ( val , VALUE_MONE )
# # Skipping MacroDefinition: value_cmp_si ( val , n ) ( val - ( n ) )
# # Skipping MacroDefinition: value_min_p ( val ) value_eq ( val , VALUE_MIN )
# # Skipping MacroDefinition: value_max_p ( val ) value_eq ( val , VALUE_MAX )
# # Skipping MacroDefinition: value_notmin_p ( val ) value_ne ( val , VALUE_MIN )
# # Skipping MacroDefinition: value_notmax_p ( val ) value_ne ( val , VALUE_MAX )
# # Skipping MacroDefinition: EXCEPTION extern unsigned int
# # Skipping MacroDefinition: THROW ( what ) ( throw_exception ( what , __CURRENT_FUNCTION_NAME__ , __FILE__ , __LINE__ ) )
# # Skipping MacroDefinition: CATCH ( what ) if ( setjmp ( * push_exception_on_stack ( what , __CURRENT_FUNCTION_NAME__ , __FILE__ , __LINE__ ) ) )
# # Skipping MacroDefinition: UNCATCH ( what ) ( pop_exception_from_stack ( what , __CURRENT_FUNCTION_NAME__ , __FILE__ , __LINE__ ) )
# # Skipping MacroDefinition: RETHROW ( ) THROW ( the_last_just_thrown_exception )
# # Skipping MacroDefinition: value_protected_hard_idiv_multiply ( v , w , throw ) ( ( value_zero_p ( w ) || value_zero_p ( v ) ) ? VALUE_ZERO : value_lt ( value_abs ( v ) , value_div ( VALUE_MAX , value_abs ( w ) ) ) ? value_direct_multiply ( v , w ) : ( throw , VALUE_NAN ) )
# # Skipping MacroDefinition: value_protected_multiply ( v , w , throw ) value_protected_hard_idiv_multiply ( v , w , throw )
# # Skipping MacroDefinition: value_protected_mult ( v , w ) value_protected_multiply ( v , w , THROW ( overflow_error ) )
# # Skipping MacroDefinition: value_protected_product ( v , w ) v = value_protected_mult ( v , w )
# # Skipping MacroDefinition: value_mult ( v , w ) value_protected_mult ( v , w )
# # Skipping MacroDefinition: value_product ( v , w ) value_protected_product ( v , w )
# # Skipping MacroDefinition: value_substract ( ref , val1 , val2 ) ( value_subtract ( ( ref ) , ( val1 ) , ( val2 ) ) )
# # Skipping MacroDefinition: ABS ( x ) ( ( ( x ) >= 0 ) ? ( x ) : - ( x ) )
# # Skipping MacroDefinition: MIN ( x , y ) ( ( ( x ) >= ( y ) ) ? ( y ) : ( x ) )
# # Skipping MacroDefinition: MAX ( x , y ) ( ( ( x ) >= ( y ) ) ? ( x ) : ( y ) )
# # Skipping MacroDefinition: SIGN ( x ) ( ( ( x ) > 0 ) ? 1 : ( ( x ) == 0 ? 0 : - 1 ) )
# # Skipping MacroDefinition: DIVIDE ( x , y ) ( ( y ) > 0 ? POSITIVE_DIVIDE ( x , y ) : - POSITIVE_DIVIDE ( ( x ) , ( - ( y ) ) ) )
# # Skipping MacroDefinition: POSITIVE_DIVIDE ( x , y ) ( ( x ) > 0 ? ( x ) / ( y ) : - ( - ( x ) + ( y ) - 1 ) / ( y ) )
# # Skipping MacroDefinition: MODULO ( x , y ) ( ( y ) > 0 ? POSITIVE_MODULO ( x , y ) : POSITIVE_MODULO ( - x , - y ) )
# # Skipping MacroDefinition: POSITIVE_MODULO ( x , y ) ( ( x ) > 0 ? ( x ) % ( y ) : ( ( x ) % ( y ) == 0 ? 0 : ( ( y ) - ( - ( x ) ) % ( y ) ) ) )

const FIRST_PARAMETER_NAME = 'P'
const PCHAR = FIRST_PARAMETER_NAME - 1
const MAXNOOFRAYS = 200
const P_VALUE_FMT = "%4ld "
const LB_INFINITY = 1
const UB_INFINITY = 2

# # Skipping MacroDefinition: MSB ( ( unsigned ) ( ( ( unsigned ) 1 ) << ( sizeof ( int ) * 8 - 1 ) ) )
# # Skipping MacroDefinition: TOP ( ( int ) ( MSB - 1 ) )
# # Skipping MacroDefinition: NEXT ( j , b ) { if ( ! ( ( b ) >>= 1 ) ) { ( b ) = MSB ; ( j ) ++ ; } }
# # Skipping MacroDefinition: POL_HIGH_BIT ( UINT_MAX - ( UINT_MAX >> 1 ) )

# const POL_NO_DUAL = POL_HIGH_BIT | 0x0001
# const POL_INTEGER = POL_HIGH_BIT | 0x0002

# # Skipping MacroDefinition: POL_ISSET ( flags , f ) ( ( flags & f ) == f )
# # Skipping MacroDefinition: FL_INIT ( l , f ) ( l ) = ( f )
# # Skipping MacroDefinition: FL_SET ( l , f ) ( ( l ) |= ( f ) )
# # Skipping MacroDefinition: FL_CLR ( l , f ) ( ( l ) &= ~ ( f ) )
# # Skipping MacroDefinition: FL_ISSET ( l , f ) ( ( l ) & ( f ) )
# # Skipping MacroDefinition: F_INIT ( p , f ) FL_INIT ( ( p ) -> flags , f )
# # Skipping MacroDefinition: F_SET ( p , f ) FL_SET ( ( p ) -> flags , f )
# # Skipping MacroDefinition: F_CLR ( p , f ) FL_CLR ( ( p ) -> flags , f )
# # Skipping MacroDefinition: F_ISSET ( p , f ) FL_ISSET ( ( p ) -> flags , f )

# const POL_INEQUALITIES = 0x00000001
# const POL_FACETS = 0x00000002
# const POL_POINTS = 0x00000004
# const POL_VERTICES = 0x00000008
# const POL_VALID = 0x00000010

# # Skipping MacroDefinition: emptyQ ( P ) ( ( F_ISSET ( P , POL_INEQUALITIES ) && P -> NbEq > P -> Dimension ) || ( F_ISSET ( P , POL_POINTS ) && P -> NbRays == 0 ) )
# # Skipping MacroDefinition: universeQ ( P ) ( P -> Dimension == P -> NbBid )
# # Skipping MacroDefinition: FORALL_PVertex_in_ParamPolyhedron ( _V , _D , _P ) \
# #{ int _i , _ix ; unsigned _bx ; for ( _i = 0 , _ix = 0 , _bx = MSB , _V = _P -> V ; _V && ( _i < _P -> nbV ) ; _i ++ , _V = _V -> next ) { if ( _D -> F [ _ix ] & _bx ) {
# # Skipping MacroDefinition: END_FORALL_PVertex_in_ParamPolyhedron } NEXT ( _ix , _bx ) ; } \#
# }
# # Skipping MacroDefinition: FOREVER for ( ; ; )
# # Skipping MacroDefinition: POL_ENSURE_INEQUALITIES ( P ) if ( F_ISSET ( P , POL_VALID ) && ! F_ISSET ( P , POL_INEQUALITIES ) ) Polyhedron_Compute_Dual ( P ) ;
# # Skipping MacroDefinition: POL_ENSURE_POINTS ( P ) if ( F_ISSET ( P , POL_VALID ) && ! F_ISSET ( P , POL_POINTS ) ) Polyhedron_Compute_Dual ( P ) ;
# # Skipping MacroDefinition: POL_ENSURE_FACETS ( P ) if ( F_ISSET ( P , POL_VALID ) && ! F_ISSET ( P , POL_FACETS ) ) Polyhedron_Compute_Dual ( P ) ;
# # Skipping MacroDefinition: POL_ENSURE_VERTICES ( P ) if ( F_ISSET ( P , POL_VALID ) && ! F_ISSET ( P , POL_VERTICES ) ) Polyhedron_Compute_Dual ( P ) ;
# # Skipping MacroDefinition: Constraints_removeParmEqs ( a , b , c , d ) Constraints_Remove_parm_eqs ( a , b , c , d )
# # Skipping MacroDefinition: Polyhedron_removeParmEqs ( a , b , c , d , e ) Polyhedron_Remove_parm_eqs ( a , b , c , d , e )
# # Skipping MacroDefinition: show_matrix ( M ) { printf ( # M "= \n" ) ; if ( M != NULL ) { Matrix_Print ( stderr , P_VALUE_FMT , ( M ) ) ; } else { printf ( "<NULL>\n" ) ; } }
# # Skipping MacroDefinition: ensureMatrix ( M , r , c ) { if ( M == NULL ) M = Matrix_Alloc ( r , c ) ; else assert ( M -> NbRows >= r && M -> NbColumns >= c ) ; }
# # Skipping MacroDefinition: Constraints_compressLastVars ( a , b ) mpolyhedron_compress_last_vars ( a , b )
# # Skipping MacroDefinition: Constraints_eliminateFirstVars ( a , b ) mpolyhedron_eliminate_first_variables ( a , b )

typealias FILE Void
typealias Value Clong
typealias exception_callback_t Ptr{Void}

type Vector
    Size::UInt32
    p::Ptr{Value}
end

# type matrix
#     NbRows::UInt32
#     NbColumns::UInt32
#     p::Ptr{Ptr{Value}}
#     p_Init::Ptr{Value}
#     p_Init_size::Cint
# end

type Matrix
    NbRows::UInt32
    NbColumns::UInt32
    p::Ptr{Ptr{Value}}
    p_Init::Ptr{Value}
    p_Init_size::Cint
end

# type polyhedron
#     Dimension::UInt32
#     NbConstraints::UInt32
#     NbRays::UInt32
#     NbEq::UInt32
#     NbBid::UInt32
#     Constraint::Ptr{Ptr{Value}}
#     Ray::Ptr{Ptr{Value}}
#     p_Init::Ptr{Value}
#     p_Init_size::Cint
#     next::Ptr{polyhedron}
#     flags::UInt32
# end

type Polyhedron
    Dimension::UInt32
    NbConstraints::UInt32
    NbRays::UInt32
    NbEq::UInt32
    NbBid::UInt32
    Constraint::Ptr{Ptr{Value}}
    Ray::Ptr{Ptr{Value}}
    p_Init::Ptr{Value}
    p_Init_size::Cint
    next::Ptr{Polyhedron}
    flags::UInt32
end

# type interval
#     MaxN::Value
#     MaxD::Value
#     MinN::Value
#     MinD::Value
#     MaxI::Cint
#     MinI::Cint
# end

type Interval
    MaxN::Value
    MaxD::Value
    MinN::Value
    MinD::Value
    MaxI::Cint
    MinI::Cint
end

# type _Param_Vertex
#     Vertex::Ptr{Matrix}
#     Domain::Ptr{Matrix}
#     Facets::Ptr{UInt32}
#     next::Ptr{_Param_Vertex}
# end

type Param_Vertices
    Vertex::Ptr{Matrix}
    Domain::Ptr{Matrix}
    Facets::Ptr{UInt32}
    next::Ptr{Param_Vertices}
end

# type _Param_Domain
#     F::Ptr{UInt32}
#     Domain::Ptr{Polyhedron}
#     next::Ptr{_Param_Domain}
# end

type Param_Domain
    F::Ptr{UInt32}
    Domain::Ptr{Polyhedron}
    next::Ptr{Param_Domain}
end

# type _Param_Polyhedron
#     nbV::Cint
#     V::Ptr{Param_Vertices}
#     D::Ptr{Param_Domain}
#     Constraints::Ptr{Matrix}
#     Rays::Ptr{Matrix}
# end

type Param_Polyhedron
    nbV::Cint
    V::Ptr{Param_Vertices}
    D::Ptr{Param_Domain}
    Constraints::Ptr{Matrix}
    Rays::Ptr{Matrix}
end

# begin enum ANONYMOUS_1
typealias ANONYMOUS_1 UInt32
const polynomial = (UInt32)(0)
const periodic = (UInt32)(1)
const evector = (UInt32)(2)
# end enum ANONYMOUS_1

# begin enum enode_type
typealias enode_type UInt32
const polynomial = (UInt32)(0)
const periodic = (UInt32)(1)
const evector = (UInt32)(2)
# end enum enode_type

# type _evalue
#     d::Value
#     x::Void
# end

type evalue
    d::Value
    x::Void
end

immutable Array_1_evalue
    d1::evalue
end

zero(::Type{Array_1_evalue}) = begin  # /home/liu/.julia/v0.4/Clang/src/wrap_c.jl, line 266:
        Array_1_evalue(fill(zero(evalue),1)...)
    end

# type _enode
#     _type::enode_type
#     size::Cint
#     pos::Cint
#     arr::Array_1_evalue
# end

type enode
    _type::enode_type
    size::Cint
    pos::Cint
    arr::Array_1_evalue
end

# type _enumeration
#     ValidityDomain::Ptr{Polyhedron}
#     EP::evalue
#     next::Ptr{_enumeration}
# end

type Enumeration
    ValidityDomain::Ptr{Polyhedron}
    EP::evalue
    next::Ptr{Enumeration}
end

# begin enum ANONYMOUS_2
typealias ANONYMOUS_2 UInt32
const False = (UInt32)(0)
const True = (UInt32)(1)
# end enum ANONYMOUS_2

# begin enum Bool
typealias Bool UInt32
const False = (UInt32)(0)
const True = (UInt32)(1)
# end enum Bool

typealias Lattice Matrix

type LatticeUnion
    M::Ptr{Lattice}
    next::Ptr{LatticeUnion}
end

type ZPolyhedron
    Lat::Ptr{Lattice}
    P::Ptr{Polyhedron}
    next::Ptr{ZPolyhedron}
end


text:
export get_commands!

function get_commands!(model::MultiPtsDriver,scene::Union{Scene,Frame{Entity{VehicleState, BicycleModel, Int}}}, roadway::Roadway,ego_index::Int)
    veh = scene[ego_index]
    x = veh.state.posG.x
    y = veh.state.posG.y
    θ = veh.state.posG.θ
    steer=model.siji.steer
    steerdotmax=model.steerdotmax
    steermax=model.steermax
    accmax=model.accmax
    v = veh.state.v
    qStart=[x y θ steer v]';
    Pts=model.Pts
    xdata=Pts[1,:]
    ydata=Pts[2,:]
    numPts=length(xdata)

    Δt = model.Δt

    i=model.index
#    if i==numPts
#        qEnd=[xdata[i] ydata[i] 0 0 0]';
#        Y,u_solved,converged=computeTrajectory(qStart,qEnd,"xy",Δt,steermax,steerdotmax,accmax);
    if i==numPts
        qEnd=[xdata[i] ydata[i] atan((ydata[i]-ydata[i-1]),(xdata[i]-xdata[i-1])) 0 0]';
    else
        qEnd=[xdata[i] ydata[i] atan((ydata[i+1]-ydata[i]),(xdata[i+1]-xdata[i])) 0 0]';
    end
    #println("qStart : ", qStart)
    #println("qEnd : ", qEnd)

    Y1,u_solved1,converged1,cost1=computeTrajectory(veh,qStart,qEnd,"xyp",Δt,steermax,steerdotmax,accmax);

    #=
    if converged1==0
        println("xyp converge fail")
        Y2,u_solved2,converged2,cost2=computeTrajectory(qStart,qEnd,"xy",Δt,steermax,steerdotmax,accmax);
        if converged2==0
            println("xy converge fail")
            if cost1<=cost2
                println("choose xyp u")
                println("u : ",u_solved1)
                model.commands=u_solved1
            else
                println("choose xy u")
                println("u : ",u_solved2)
                model.commands=u_solved2
            end
        else
            println("xy converge success")
            println("u : ",u_solved2)
            model.commands=u_solved2
        end
    else
        println("xyp converge success")
        println("u : ",u_solved1)
        model.commands=u_solved1
    end
    =#
    model.commands=u_solved1
end

function computeTrajectory(veh::Entity,initialXYPKV::AbstractArray,targetXYPKV::AbstractArray,mode::String,Δt::Float64,kmax::Float64,kdotmax::Float64,vdotmax::Float64)
    h=Δt
    N=2
    #x=sqrt.((initialXYPKV[1]-targetXYPKV[1])^2+(initialXYPKV[2]-targetXYPKV[2])^2)
    #guess_a=2*(x-initialXYPKV[5]*h)/h^2
    #uguess = float([guess_a guess_a 0 0]')
    uguess = float([0 0 0 0]')
    Y,u_solved,converged,cost = solveTrajXYPKVLinear(veh,initialXYPKV,N,h,uguess,kmax,kdotmax,vdotmax,targetXYPKV,mode);
    #=
    if u_solved[3]*u_solved[4]<0.0 && converged==false
        println("u is ",u_solved)
        println("counter steer! average them!")
        u_solved[3]=u_solved[3]+u_solved[4]
        u_solved[4]=u_solved[3]
    end
    =#
    return Y,u_solved,converged,cost
end

function solveTrajXYPKVLinear(veh::Entity,q0::AbstractArray,N::Int,h::Float64,uguess::AbstractArray,kmax::Float64,kdotmax::Float64,vdotmax::Float64,targetXYPKV::AbstractArray,constraint_mode::String)

    w=copy(uguess)
    stepsize_multiplier = 0.5
    converged = false
    bestCost=Inf
    bestY=zeros(N+1,5);
    bestU=zeros(4,1);
    for i=1:500
        #println("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~")
        #println("w is ",w)
        u=copy(w)
        Y,dY = simXYPKVLinear(veh,q0,N,h,u,kmax,kdotmax,vdotmax)
        if constraint_mode=="xypkv"
            index=1:5
        elseif constraint_mode=="xypv"
            index=[1,2,3,5]
        elseif constraint_mode== "xyv"
            index=[1,2,5]
        elseif constraint_mode== "xyp"
            index=[1,2,3]
        elseif constraint_mode== "xy"
            index=[1,2]
        else
            index=[1,2]
        end
        f = Y[end,index] - targetXYPKV[index]

        cost=0;
        for j=1:length(index)
            if index[j]==3 || index[j]==4
                fjs=[f[j] f[j]+2*pi f[j]-2*pi];
                (abs_fj,ind_fj)=findmin(abs.(fjs))
                f[j]=fjs[ind_fj]
            end
            cost=cost+f[j]^2
        end
        dfdw = dY[index,:,end]
        #cost=maximum(abs.(f))
        dcostdw=zeros(4,1)
        for j=1:length(index)
            #println("j = ",j," f[j] is ",f[j])
            #normalize dfdw
            ####dfdw[j,:]=dfdw[j,:]/sqrt.(dfdw[j,:]'*dfdw[j,:])
            #println("dfdw is ",dfdw[j,:])
            #dcostdw=dcostdw+2*f[j]*dY[index[j],:,end]
            dcostdw=dcostdw+2*f[j]*dfdw[j,:]
        end

        if cost<bestCost
            bestCost=cost
            bestY=copy(Y)
            bestU=copy(u)
        end
        if cost<0.0001
            converged=true
            break
        end

        alpha=0.01*[vdotmax vdotmax kmax kmax]'
        #lambda=0.01*[vdotmax vdotmax kmax kmax]'
        lambda=0.001*[vdotmax vdotmax kmax kmax]'
        #cost_constrain=-log((w[1]-vdotmax)^2)-log((w[1]+vdotmax)^2)-log((w[2]-vdotmax)^2)-log((w[2]+vdotmax)^2)-log((w[3]-kmax)^2)-log((w[3]+kmax)^2)-log((w[4]-kmax)^2)-log((w[4]+kmax)^2);
        dcost_constraindw=zeros(4,1);
        for j=1:4
            if j<3 && abs.(w[j])>vdotmax*0.0
                dcost_constraindw[j]=-2/(w[j]-vdotmax)-2/(w[j]+vdotmax)
                if w[j]>(vdotmax-0.1)
                    #println("command ",j," is ",w[j])
                    dcost_constraindw[j]=-2/(-0.1)
                elseif w[j]<(-vdotmax+0.1)
                    #println("command ",j," is ",w[j])
                    dcost_constraindw[j]=-2/(0.1)
                end
            elseif j>2 && abs.(w[j])>kmax*0.0
                dcost_constraindw[j]=-2/(w[j]-kmax)-2/(w[j]+kmax)
                if w[j]>(kmax-0.1)
                    #println("command ",j," is ",w[j])
                    dcost_constraindw[j]=-2/(-0.1)
                elseif w[j]<(-kmax+0.1)
                    #println("command ",j," is ",w[j])
                    dcost_constraindw[j]=-2/(0.1)
                end
            end
        end
        #=
        println("cost at i = ",i," is ",cost)
        println("w is ",w)
        println("dcostdw is ",dcostdw)
        println("dcost_constraindw is ",lambda.*dcost_constraindw)
        =#
        gradient=dcostdw+lambda.*dcost_constraindw;
        #normalize whole gradient
        gradient=gradient/sqrt.(gradient'*gradient)
        #=
        if sqrt.(gradient'*gradient)[1]*alpha > min(vdotmax,kmax)
            println("gradient too large, is ",sqrt.(gradient'*gradient)[1]*alpha)
            gradient=gradient/sqrt.(gradient'*gradient)
        end
        =#
        #println("gradiet is : ",alpha.*gradient)
        w = w - alpha.*gradient;
    end
    Y=bestY;

    if  bestU[1]!=(bestU[1]=clamp.(bestU[1],-vdotmax,vdotmax))
        converged=false
    end
    if bestU[2]!=(bestU[2]=clamp.(bestU[2],-vdotmax,vdotmax))
        converged=false
    end
    if bestU[3]!=(bestU[3]=clamp.(bestU[3],-kmax,kmax))
        converged=false
    end
    if bestU[4]!=(bestU[4]=clamp.(bestU[4],-kmax,kmax))
        converged=false
    end
    #=
    @assert bestU[1]==(bestU[1]=clamp.(bestU[1],-vdotmax,vdotmax))

    @assert bestU[2]==(bestU[2]=clamp.(bestU[2],-vdotmax,vdotmax))

    @assert bestU[3]==(bestU[3]=clamp.(bestU[3],-kmax,kmax))

    @assert bestU[4]==(bestU[4]=clamp.(bestU[4],-kmax,kmax))
    =#
    u=bestU;
    uvec=u;
    return Y,uvec,converged,bestCost
end

function simXYPKVLinear(veh::Entity,q0::AbstractArray,N::Int,h::Float64,u::AbstractArray,kmax::Float64,kdotmax::Float64,vdotmax::Float64)
    M = length(u);
    Y = zeros(N+1,5);
    Y[1,:] = q0;
    #u=[a1 a2 k1 k2]
    # Derivative of [each state at each time step] with respect to [each decision variable]
    dY = zeros(5,M,N+1);

    # Simulate forward N timesteps
    for i=1:N
        if i==1
            steer = u[3];
            acc = u[1];
            dsteerdu=[0 0 1 0];
            daccdu=[1 0 0 0];
        else
            steer = u[4];
            acc = u[2];
            dsteerdu=[0 0 0 1];
            daccdu=[0 1 0 0];
        end

        # Max change in curvature (i.e. discrete limit on kdot)
        kdeltamax = kdotmax*h;

        # h (timestep) multiplied by f(x), where xdot=f(x)
        L = veh.def.a + veh.def.b
        l = -veh.def.b

        x=Y[i,1];
        dxdu=dY[1,:,i];
        dxdu=reshape(dxdu,1,M)
        y=Y[i,2];
        dydu=dY[2,:,i];
        dydu=reshape(dydu,1,M)
        theta=Y[i,3];
        dthetadu=dY[3,:,i];
        dthetadu=reshape(dthetadu,1,M)
        v=Y[i,5];
        dvdu=dY[5,:,i];
        dvdu=reshape(dvdu,1,M)
        acc=clamp.(acc,-vdotmax,vdotmax);
        steer=Y[i,4]+clamp.(steer-Y[i,4],-kdeltamax,kdeltamax);
        steer=clamp.(steer,-kmax,kmax);
        s=v*h+0.5*acc*h^2;
        dsdu=dvdu*h+0.5*h^2*daccdu;
        if(abs.(steer)<0.01)
         # Integrate dynamics
            theta_n=theta+steer*v*h;
            theta_n=mod(theta_n,2*pi)
            dtheta_ndu=dthetadu+dsteerdu*v*h+steer*dvdu*h;

            Y[i+1,1] = x+s*cos(theta_n);
            Y[i+1,2] = y+s*sin(theta_n);
            Y[i+1,3] = theta_n;
            Y[i+1,4] = steer;
            Y[i+1,5] = v+acc*h;

            #calculate derivative
            dY[1,:,i+1]=dxdu+dsdu*cos(theta_n)-s*sin(theta_n)*dtheta_ndu;
            dY[2,:,i+1]=dydu+dsdu*sin(theta_n)+s*cos(theta_n)*dtheta_ndu;
            dY[3,:,i+1]=dtheta_ndu;#dthetadu+dsteerdu*v*h+steer*dvdu*h;
            dY[4,:,i+1]=dsteerdu;
            dY[5,:,i+1]=dvdu+daccdu*h;
        else
            R=L/tan(steer);
            dRdu=-L*tan(steer)^(-2)*sec(steer)^2*dsteerdu;
            beta=s/R;
            dbetadu=(dsdu*R-s*dRdu)/R^2;
            xc=x-R*sin(theta)+l*cos(theta);
            dxcdu=dxdu-dRdu*sin(theta)-R*cos(theta)*dthetadu-l*sin(theta)*dthetadu;
            yc=y+R*cos(theta)+l*sin(theta);
            dycdu=dydu+dRdu*cos(theta)-R*sin(theta)*dthetadu+l*cos(theta)*dthetadu;
            theta_n=mod(theta+beta,2*pi);
            dtheta_ndu=dthetadu+dbetadu;
            # Integrate dynamics

            Y[i+1,1] = xc+R*sin(theta_n)-l*cos(theta_n);
            Y[i+1,2] = yc-R*cos(theta_n)-l*sin(theta_n);
            Y[i+1,3] = theta_n;
            Y[i+1,4] = steer;
            Y[i+1,5] = v+acc*h;

            #calculate derivative
            dY[1,:,i+1]=dxcdu+dRdu*sin(theta_n)+R*cos(theta_n)*dtheta_ndu+l*sin(theta_n)*dtheta_ndu;
            dY[2,:,i+1]=dycdu-dRdu*cos(theta_n)+R*sin(theta_n)*dtheta_ndu-l*cos(theta_n)*dtheta_ndu;
            dY[3,:,i+1]=dtheta_ndu;
            dY[4,:,i+1]=dsteerdu;
            dY[5,:,i+1]=dvdu+daccdu*h;
        end
    end
    return Y,dY
end

max_stars_repo_name: sisl/UrbanDrivingSimulation

max_issues_repo_name: sisl/UrbanDrivingSimulation

max_forks_repo_name: sisl/UrbanDrivingSimulation


text:
struct WaveletParams{TF, Tw}
    s::TF
    N::Int
    J::Int
    L::Int
    K::Int
    σ::TF

    ws::Tw

    Ψs::Vector{Matrix{TF}}

    Γ_H::Vector{NamedTuple{
        (:j, :l, :θ, :k, :j_, :l_, :θ_, :k_, :τ_),
        Tuple{Int, Int, TF, Int, Int,
        Int, TF, Int, Tuple{Int, Int}}}
        }

    ix::Vector{Tuple{Int, Int}}
    ix_shift::Vector{Tuple{Int, Int, Tuple{Int, Int}}}
    ix_0::Vector{Tuple{Int, Int}}
end

function WaveletParams(s, N, J, L, K, σ, Γ_H = default_Γ_H(J, L, K, θtype=typeof(s)))
    dx = 2s / N
    xs = range(-s; length=N, step=dx)
    ws = fftfreq(N, 1 / dx)

    g = MvNormal([0, 0], σ^2 * I)
    gaus = [pdf(g, [x, y]) for x in xs, y in xs]
    Gaus = gaus |> ifftshift |> fft |> real

    float_type = typeof(s)

    # wavelet_matrices(ws, J, L; ξ0=maximum(ws) / 2, c=1) =
    ξ0 = maximum(ws) / 2
    c = 1
    wavelet_matrices = [wavelet_matrix(ws, j, θ; ξ0, L, c) for j in 0:J - 1, θ in range(0., length=L, step=2π / L)]
    Ψs = map(Ψ -> float_type.(Ψ .* Gaus), reshape(wavelet_matrices, :))

    ix, ix_shift, ix_0, perm_ix, perm_ix_shift = get_ix_subsets_from_Γ_H(Γ_H, J, L, K)

    WaveletParams(s, N, J, L, K, σ, ws, Ψs, Γ_H, ix, ix_shift, ix_0)
end

e_θ(l, L) = (cos(π / 2 + 2π * l / L), sin(π / 2 + 2π * l / L))
τ_θ_j(l, j_, L) = round.(Int, 2^j_ .* e_θ(l, L))

function default_Γ_H(J, L, K; θtype=Float64)
    Γ_H = NamedTuple{
        (:j, :l, :θ, :k, :j_, :l_, :θ_, :k_, :τ_),
        Tuple{Int, Int, θtype, Int, Int,
        Int, θtype, Int, Tuple{Int, Int}}
    }[]
    p = Iterators.product(
        0:L - 1,
        0:L - 1,
        0:J - 1,
        0:J - 1,
        [false, true],
        0:K - 1,
        0:K - 1
    )

    for (l, l_, j, j_, isshift, k, k_) in p
        l, l_, j, j_, isshift, k, k_
        if l_ > l
            continue
        end
        if (j > j_) || (j < j_ - 2)
            continue
        end
        if k > min(k_, 1)
            continue
        end

        if isshift
            τ_ = τ_θ_j(l, j_, L)
        else
            τ_ = (0, 0)
        end

        θ = 2π * l / L
        θ_ = 2π * l_ / L
        dl = abs(l - l_)

        if j == j_
            if k == 0 && k_ > 1
                continue
            end
            if k == 1
                if k_ != 1
                    continue
                elseif dl > 2
                    continue
                end
            end
        else
            if k == 0 && k_ == 4
                continue
            end
            if k == 1
                if k_ != 2^(j_ - j)
                    continue
                elseif dl > 2
                    continue
                end
            end
        end

        push!(Γ_H, (;j, l, θ, k, j_, l_, θ_, k_, τ_))
    end
    Γ_H
end

function get_ix_from_Γ_H(p, J, L)
    ix1 = (p.j + 1) + J * (p.l) + J * L * (p.k)
    ix2 = (p.j_ + 1) + J * (p.l_) + J * L * (p.k_)
    (p.τ_ == (0, 0), ix1, ix2, p.τ_)
end

function get_ix_subsets_from_Γ_H(Γ_H, J, L, K)
    ix_all = map(p -> get_ix_from_Γ_H(p, J, L), Γ_H)

    perm_ix = findall(s -> s[1] == true, ix_all)
    ix = map(s -> s[2:3], ix_all[perm_ix])

    perm_ix_shift = findall(s -> s[1] == false, ix_all)
    ix_shift = map(s -> s[2:4], ix_all[perm_ix_shift])

    ix_0 = map(i -> (1, i), 1:J * L * K)
    ix, ix_shift, ix_0, perm_ix, perm_ix_shift
end


text:
@testset "Generalized fmap over equivalent functors" begin
  struct M{F,T,S}
    σ::F
    W::T
    b::S
  end

  @functor M

  (m::M)(x) = m.σ.(m.W * x .+ m.b)

  m = M(identity, ones(Float32, 3, 4), zeros(Float32, 3))
  x = ones(Float32, 4, 2)
  m̄, _ = gradient((m,x) -> sum(m(x)), m, x)
  m̂ = Functors.fmap(m, m̄) do x, y
    isnothing(x) && return y
    isnothing(y) && return x
    x .- 0.1f0 .* y
  end

  @test m̂.W ≈ fill(0.8f0, size(m.W))
  @test m̂.b ≈ fill(-0.2f0, size(m.b))
end