text:
""" Utilities for integrating GPS with InSAR maps

Links:

1. list of LLH for all gps stations: ftp://gneiss.nbmg.unr.edu/rapids/llh
Note: ^^ This file is stored in the `STATION_LLH_FILE`

2. Clickable names to explore: http://geodesy.unr.edu/NGLStationPages/GlobalStationList
3. Map of stations: http://geodesy.unr.edu/NGLStationPages/gpsnetmap/GPSNetMap.html

"""

import CSV

const GPS_BASE_URL = "http://geodesy.unr.edu/gps_timeseries/tenv3/NA12/{station}.NA12.tenv3"

const STATION_LLH_FILE = joinpath(@__DIR__, "station_llh_all.csv")
const START_YEAR = 2014  # So far I don't care about older data

"""Returns the config folder for the application.  The default behavior
is to return whatever is most appropriate for the operating system.

This is used to store gps timeseries data

the following folders could be returned:
Mac OS X:
  ``~/Library/Application Support/apertools``
Mac OS X (POSIX):
  ``~/.apertools``
Unix:
  ``~/.cache/apertools``
Unix (POSIX):
  ``~/.apertools``

force_posix: if this is set to `true` then on any POSIX system the
    folder will be stored in the home folder with a leading
    dot instead of the XDG config home or darwin's
    application support folder.

Base on: https://github.com/pallets/click/blob/master/click/utils.py#L368
"""
function get_cache_dir(force_posix=false, app_name="gps")
    if force_posix
        path = joinpath(expanduser("~/.$app_name"))
    end
    # if Sys.isapple()
    #     path = joinpth(expanduser("~/Library/Application Support"), app_name)
    # end
    path = joinpath(
        get(ENV, "XDG_CONFIG_HOME", expanduser("~/.cache")),
        app_name
    )
    !exists(path) && mkdir(path)
    return path
end

"Read in a DataFrame of gps stations with columns name, lat, lon, alt"
read_station_llas() = CSV.read(STATION_LLH_FILE)
# columns = ["name", "lat", "lon", "alt"]

"Return the (lon, lat) of a `station_name`"
function station_lonlat(station_name::AbstractString)
    df = read_station_llas()
    !(station_name in df.name)&& error("No station named $station_name found.")
    # TODO:
    # closest_names = difflib.get_close_matches(station_name, df.name, n=5)
    # error("No station named %s found. Closest: %s" % (station_name, closest_names))
    name, lat, lon, alt = df[df.name .== station_name, :][1, :]
    return lon, lat
end

station_latlon(station_name) = reverse(station_lonlat(station_name))

station_rowcol(name::AbstractString, demrsc::DemRsc) = latlon_to_rowcol(demrsc, station_latlon(name)...)
station_rowcol(name::AbstractString, img::MapImage) = latlon_to_rowcol(img.demrsc, station_latlon(name)...)

# function station_distance(station_name1, station_name2)
#     lonlat1 = station_lonlat(station_name1)
#     lonlat2 = station_lonlat(station_name2)
#     return apertools.latlon.latlon_to_dist(lonlat1[::-1], lonlat2[::-1])
# end

# function stations_within_rsc(rsc_filename=nothing, demrsc=nothing, gps_filename=nothing)
#     if demrsc is nothing
#         if rsc_filename is nothing
#             error("Need demrsc or rsc_filename")
#         end
#         demrsc = apertools.sario.load(rsc_filename)
#     end
#
#     df = read_station_llas(filename=gps_filename)
#     station_lon_lat_arr = df[["lon", "lat"]].values
#     contains_bools = [apertools.latlon.grid_contains(s, **demrsc) for s in station_lon_lat_arr]
#     return df[contains_bools][["name", "lon", "lat"]].values
# end

function download_station_data(station_name::AbstractString)
    station_name = station_name.upper()
    url = GPS_BASE_URL.format(station=station_name)
    response = requests.get(url)

    stationfile = split(url, "/")[-1]  # Just blah.NA12.tenv3
    filename = "$GPS_DIR/$stationfile"
    println("Saving to $filename")

    write(filename, response.text)
end

function load_station_enu(station::AbstractString; start_year=START_YEAR, end_year=nothing, to_cm=true)
    """Loads one gps station's ENU data since start_year until end_year
    as separate Series items

    Will scale to cm by default, and center the first data point
    to 0
    """
    station_df = load_station_data(station, to_cm=to_cm, start_year=start_year, end_year=end_year)

    start_val = station_df[["east", "north", "up"]].iloc[:10].mean()
    enu_zeroed = station_df[["east", "north", "up"]] - start_val
    dts = station_df["dt"]
    return dts, enu_zeroed
end

function load_station_data(station_name::AbstractString;
                      start_year=START_YEAR,
                      end_year=nothing,
                      download_if_missing=true,
                      to_cm=true)
    """Loads one gps station's ENU data since start_year until end_year as a dataframe

    Args:
        station_name (str): 4 Letter name of GPS station
            See http://geodesy.unr.edu/NGLStationPages/gpsnetmap/GPSNetMap.html for map
        start_year (int), default 2014, cutoff for beginning of GPS data
        end_year (int): default nothing, cut off for end of GPS data
        download_if_missing (bool): default true
    """
    station_name = station_name.upper()
    gps_data_file = os.path.join(GPS_DIR, "%s.NA12.tenv3" % station_name)
    if not os.path.exists(gps_data_file)
        logger.warning("%s does not exist.", gps_data_file)
        if download_if_missing
            logger.info("Downloading %s", station_name)
            download_station_data(station_name)
        end
    end

    df = pd.read_csv(gps_data_file, header=0, sep=r"\s+")
    clean_df = _clean_gps_df(df, start_year, end_year)
    if to_cm
        # logger.info("Converting %s GPS to cm" % station_name)
        clean_df[["east", "north", "up"]] = 100 * clean_df[["east", "north", "up"]]
    end
    return clean_df
end

function _clean_gps_df(df, start_year=nothing, end_year=nothing)
    df["dt"] = pd.to_datetime(df["YYMMMDD"], format="%y%b%d")

    if start_year
        df_ranged = df[df["dt"] >= datetime.datetime(start_year, 1, 1)]
    end
    if end_year
        df_ranged = df_ranged[df_ranged["dt"] <= datetime.datetime(end_year, 1, 1)]
    end
    df_enu = df_ranged[["dt", "__east(m)", "_north(m)", "____up(m)"]]
    df_enu = df_enu.rename(mapper=s -> replace(replace(s, "_" => ""), "(m)" => ""), axis="columns")
    df_enu.reset_index(inplace=true, drop=true)
    return df_enu
end

# function stations_within_image(image_ll=nothing, filename=nothing, mask_invalid=true, gps_filename=nothing)
#     """Given an image, find gps stations contained in area
#
#     Must have an associated .rsc file, either from being a LatlonImage (passed
#     by image_ll), or have the file in same directory as `filename`
#
#     Args
#         image_ll (LatlonImage) LatlonImage of area with data
#         filename (str) filename to load into a LatlonImage
#         mask_invalid (bool) Default true. if true, don't return stations
#             where the image value is NaN or exactly 0
#
#     Returns
#         ndarray Nx3, with columns ["name", "lon", "lat"]
#     """
#     if image_ll is nothing
#         image_ll = apertools.latlon.LatlonImage(filename=filename)
#     end
#
#     df = read_station_llas(filename=gps_filename)
#     station_lon_lat_arr = df[[:lon, :lat]].values
#     contains_bools = image_ll.contains(station_lon_lat_arr)
#     candidates = df[contains_bools][["name", "lon", "lat"]].values
#     good_stations = []
#     if mask_invalid
#         for name, lon, lat in candidates
#             val = image_ll[..., lat, lon]
#             if np.isnan(val)  # or val == 0 TODO with window 1 reference, it's 0
#                 continue
#             else
#                 good_stations.append([name, lon, lat])
#             end
#         end
#     else
#         good_stations = candidates
#     end
#     return good_stations
# end

# function save_station_points_kml(station_iter)
#     for name, lat, lon, alt in station_iter
#         apertools.kml.create_kml(
#             title=name,
#             desc="GPS station location",
#             lon_lat=(lon, lat),
#             kml_out="%s.kml" % name,
#             shape="point",
#         )
#     end
# end


text:
function build_nabla2(LF)

  Nx = LF.Nx
  Ny = LF.Ny
  Nz = LF.Nz
  Npoints = Nx*Ny*Nz

  D2jl_x = LF.LFx.D2jl
  D2jl_y = LF.LFy.D2jl
  D2jl_z = LF.LFz.D2jl

  nnzc = Nx + Ny + Nz - 2
  NNZ = nnzc*Npoints

  rowval = zeros( Int64, NNZ )
  nzval = zeros( Float64, NNZ )
  colptr = zeros( Int64, Npoints + 1 )

  rowGbl_x_orig = zeros( Int64, Nx )
  rowGbl_y_orig = zeros( Int64, Ny )
  rowGbl_z_orig = zeros( Int64, Nz )

  rowGbl_x = zeros( Int64, Nx )
  rowGbl_y = zeros( Int64, Ny )
  rowGbl_z = zeros( Int64, Nz )

  for iz = 1:Nz
    rowGbl_z_orig[iz] = iz
  end

  rowGbl_y_orig[1] = 1
  for iy = 2:Ny
    rowGbl_y_orig[iy] = rowGbl_y_orig[iy-1] + Nz
  end

  rowGbl_x_orig[1] = 1
  for ix = 2:Nx
    rowGbl_x_orig[ix] = rowGbl_x_orig[ix-1] + Ny*Nz
  end

  ip_start = 0
  ip_stop = 0

  rowval_unsrt = zeros( Int64, nnzc )
  nzval_unsrt = zeros( Float64, nnzc )

  # rowLoc = 1:Nx, 1:Ny, 1:Nz for each gblCol
  for colGbl = 1:Npoints

    colLoc_x = ceil(Int,colGbl/(Ny*Nz))

    yy = colGbl - (colLoc_x-1)*Ny*Nz
    colLoc_y = ceil( Int, yy/Nz )

    izz = ceil(Int,colGbl/Nz)
    colLoc_z = colGbl - (izz-1)*Nz

    ip = 0  # reset ip

    # diagonal element
    ip = ip + 1
    rowval_unsrt[ip] = colGbl
    nzval_unsrt[ip] = D2jl_x[colLoc_x,colLoc_x] + D2jl_y[colLoc_y,colLoc_y] +  D2jl_z[colLoc_z,colLoc_z]

    # non diagonal element
    for ir = 1:Nx
      if ir != colLoc_x
        ip = ip + 1
        rowval_unsrt[ip] = rowGbl_x_orig[ir] + colGbl - 1 - (colLoc_x-1)*Nz*Ny
        nzval_unsrt[ip] = D2jl_x[ir,colLoc_x]
      end
    end

    for ir = 1:Ny
      if ir != colLoc_y
        ip = ip + 1
        rowval_unsrt[ip] = rowGbl_y_orig[ir] + colGbl - 1 - (izz-1)*Nz + (colLoc_x-1)*Ny*Nz
        nzval_unsrt[ip] = D2jl_y[ir,colLoc_y]
      end
    end

    for ir = 1:Nz
      if ir != colLoc_z
        ip = ip + 1
        rowval_unsrt[ip] = rowGbl_z_orig[ir] + (izz-1)*Nz
        nzval_unsrt[ip] = D2jl_z[ir,colLoc_z]
      end
    end
    ip_start = (colGbl-1)*nnzc + 1
    ip_stop = ip_start + nnzc - 1

    #println( colGbl, ip_start, ip_stop )

    idx_srt = sortperm(rowval_unsrt)

    rowval[ip_start:ip_stop] = rowval_unsrt[idx_srt]
    nzval[ip_start:ip_stop] = nzval_unsrt[idx_srt]

  end

  #println("ip_stop = ", ip_stop)
  #println("NNZ = ", NNZ)

  # colptr has more predictable pattern
  colptr[1] = 1
  for i = 2:Npoints+1
    colptr[i] = colptr[i-1] + nnzc
  end

  nabla2 = SparseMatrixCSC( Npoints, Npoints, colptr, rowval, nzval )

end

ext: jl


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

JLLWrappers.@generate_main_file_header("Qt5Declarative")
JLLWrappers.@generate_main_file("Qt5Declarative", UUID("c6373c32-5b88-5913-90f5-31d7686b42da"))
end  # module Qt5Declarative_jll


text:
@testset "1D OperatorConv" begin
    modes = (16, )
    ch = 64 => 128

    m = Chain(
        Dense(2, 64),
        OperatorConv(ch, modes, FourierTransform)
    )
    @test ndims(OperatorConv(ch, modes, FourierTransform)) == 1
    @test repr(OperatorConv(ch, modes, FourierTransform)) == "OperatorConv(64 => 128, (16,), FourierTransform, permuted=false)"

    𝐱 = rand(Float32, 2, 1024, 5)
    @test size(m(𝐱)) == (128, 1024, 5)

    loss(x, y) = Flux.mse(m(x), y)
    data = [(𝐱, rand(Float32, 128, 1024, 5))]
    Flux.train!(loss, Flux.params(m), data, Flux.ADAM())
end

@testset "permuted 1D OperatorConv" begin
    modes = (16, )
    ch = 64 => 128

    m = Chain(
        Conv((1, ), 2=>64),
        OperatorConv(ch, modes, FourierTransform, permuted=true)
    )
    @test ndims(OperatorConv(ch, modes, FourierTransform, permuted=true)) == 1
    @test repr(OperatorConv(ch, modes, FourierTransform, permuted=true)) == "OperatorConv(64 => 128, (16,), FourierTransform, permuted=true)"

    𝐱 = rand(Float32, 2, 1024, 5)
    𝐱 = permutedims(𝐱, (2, 1, 3))
    @test size(m(𝐱)) == (1024, 128, 5)

    loss(x, y) = Flux.mse(m(x), y)
    data = [(𝐱, rand(Float32, 1024, 128, 5))]
    Flux.train!(loss, Flux.params(m), data, Flux.ADAM())
end

@testset "1D OperatorKernel" begin
    modes = (16, )
    ch = 64 => 128

    m = Chain(
        Dense(2, 64),
        OperatorKernel(ch, modes, FourierTransform)
    )
    @test repr(OperatorKernel(ch, modes, FourierTransform)) == "OperatorKernel(64 => 128, (16,), FourierTransform, σ=identity, permuted=false)"

    𝐱 = rand(Float32, 2, 1024, 5)
    @test size(m(𝐱)) == (128, 1024, 5)

    loss(x, y) = Flux.mse(m(x), y)
    data = [(𝐱, rand(Float32, 128, 1024, 5))]
    Flux.train!(loss, Flux.params(m), data, Flux.ADAM())
end

@testset "permuted 1D OperatorKernel" begin
    modes = (16, )
    ch = 64 => 128

    m = Chain(
        Conv((1, ), 2=>64),
        OperatorKernel(ch, modes, FourierTransform, permuted=true)
    )
    @test repr(OperatorKernel(ch, modes, FourierTransform, permuted=true)) == "OperatorKernel(64 => 128, (16,), FourierTransform, σ=identity, permuted=true)"

    𝐱 = rand(Float32, 2, 1024, 5)
    𝐱 = permutedims(𝐱, (2, 1, 3))
    @test size(m(𝐱)) == (1024, 128, 5)

    loss(x, y) = Flux.mse(m(x), y)
    data = [(𝐱, rand(Float32, 1024, 128, 5))]
    Flux.train!(loss, Flux.params(m), data, Flux.ADAM())
end

@testset "2D OperatorConv" begin
    modes = (16, 16)
    ch = 64 => 64

    m = Chain(
        Dense(1, 64),
        OperatorConv(ch, modes, FourierTransform)
    )
    @test ndims(OperatorConv(ch, modes, FourierTransform)) == 2

    𝐱 = rand(Float32, 1, 22, 22, 5)
    @test size(m(𝐱)) == (64, 22, 22, 5)

    loss(x, y) = Flux.mse(m(x), y)
    data = [(𝐱, rand(Float32, 64, 22, 22, 5))]
    Flux.train!(loss, Flux.params(m), data, Flux.ADAM())
end

@testset "permuted 2D OperatorConv" begin
    modes = (16, 16)
    ch = 64 => 64

    m = Chain(
        Conv((1, 1), 1=>64),
        OperatorConv(ch, modes, FourierTransform, permuted=true)
    )
    @test ndims(OperatorConv(ch, modes, FourierTransform, permuted=true)) == 2

    𝐱 = rand(Float32, 1, 22, 22, 5)
    𝐱 = permutedims(𝐱, (2, 3, 1, 4))
    @test size(m(𝐱)) == (22, 22, 64, 5)

    loss(x, y) = Flux.mse(m(x), y)
    data = [(𝐱, rand(Float32, 22, 22, 64, 5))]
    Flux.train!(loss, Flux.params(m), data, Flux.ADAM())
end

@testset "2D OperatorKernel" begin
    modes = (16, 16)
    ch = 64 => 64

    m = Chain(
        Dense(1, 64),
        OperatorKernel(ch, modes, FourierTransform)
    )

    𝐱 = rand(Float32, 1, 22, 22, 5)
    @test size(m(𝐱)) == (64, 22, 22, 5)

    loss(x, y) = Flux.mse(m(x), y)
    data = [(𝐱, rand(Float32, 64, 22, 22, 5))]
    Flux.train!(loss, Flux.params(m), data, Flux.ADAM())
end

@testset "permuted 2D OperatorKernel" begin
    modes = (16, 16)
    ch = 64 => 64

    m = Chain(
        Conv((1, 1), 1=>64),
        OperatorKernel(ch, modes, FourierTransform, permuted=true)
    )

    𝐱 = rand(Float32, 1, 22, 22, 5)
    𝐱 = permutedims(𝐱, (2, 3, 1, 4))
    @test size(m(𝐱)) == (22, 22, 64, 5)

    loss(x, y) = Flux.mse(m(x), y)
    data = [(𝐱, rand(Float32, 22, 22, 64, 5))]
    Flux.train!(loss, Flux.params(m), data, Flux.ADAM())
end

@testset "SpectralConv" begin
    modes = (16, 16)
    ch = 64 => 64
    permuted = false

    @test SpectralConv(ch, modes) isa OperatorConv
    @test SpectralConv(ch, modes).transform isa FourierTransform
end

@testset "GraphKernel" begin
    batch_size = 5
    channel = 32
    N = 10*10

    κ = Dense(2*channel, channel, relu)

    graph = grid([10, 10])
    𝐱 = rand(Float32, channel, N, batch_size)
    l = WithGraph(FeaturedGraph(graph), GraphKernel(κ, channel))
    @test repr(l.layer) == "GraphKernel(Dense(64 => 32, relu), channel=32)"
    @test size(l(𝐱)) == (channel, N, batch_size)

    g = Zygote.gradient(() -> sum(l(𝐱)), Flux.params(l))
    @test length(g.grads) == 3
end


text:
### Singletons ###
abstract type ALogicConstant <: AConstantSymbol end

# starts with capital so as not to conflict with `true`
@new_constant ALogicConstant ATrue True Statements
@new_constant ALogicConstant AFalse False Statements
# not yet computed
@new_constant ALogicConstant AUnknown Unknown Statements
# computed and the result could not be obtained within the constraints
@new_constant ALogicConstant AUnknowable Unknowable Statements

### Functions ###
abstract type ALogicFunction <: AFunctionSymbol end

@new_function(ALogicFunction, AAnd, And, (Statements, Statements), Statements, 2, ∧)
@new_function(ALogicFunction, AOr, Or, (Statements, Statements), Statements, 2, ∨)
@new_function(ALogicFunction, ANot, Not, (Statements,), Statements, 2, ~, false)
Base.show(io::IO, self::Term{<:ANot}) = print(io, "$(Meta.parse("$(data(self))($(args(self)[1]))"))")

@new_function(ALogicFunction, AImplies, Implies, (Statements, Statements), Statements, 2, ⟹)
@new_function(ALogicFunction, AEquivalent, Equivalent, (Statements, Statements), Statements, 2, ⩶)

@new_function(ALogicFunction, AElementOf, ElementOf, (Ω, Sets), Statements, 2, ∈, false)
@new_function(ALogicFunction, ASubset, Subset, (Sets, Sets), Statements, 2, ⊂)

@new_function(ALogicFunction, AExists, Exists, (Ω, Sets, Statements), Statements, 3)
Exists(x::Term{ElementOf}, y::Term) = Exists(args(x)[1], args(x)[2], y)
Base.show(io::IO, self::Term{<:AExists}) = print(io, "Exists($(args(self)[1]) ∈ $(args(self)[2]), $(args(self)[3]))")

@new_function(ALogicFunction, AForAll, ForAll, (Ω, Sets, Statements), Statements, 3)
ForAll(x::Term{ElementOf}, y::Term) = ForAll(args(x)[1], args(x)[2], y)
Base.show(io::IO, self::Term{<:AForAll}) = print(io, "ForAll($(args(self)[1]) ∈ $(args(self)[2]), $(args(self)[3]))")


text:
# using Revise

using Caesar
using KernelDensityEstimatePlotting
using KernelDensityEstimate, IncrementalInference
using DocStringExtensions, DelimitedFiles
using Gadfly, Cairo, Fontconfig
Gadfly.set_default_plot_size(35cm,25cm)

include(joinpath(@__DIR__,"slamUtils.jl"))

wstart = 415;
wlen = 171;

topdir = joinpath(ENV["HOME"],"data","kayaks")
trialstr = "20_gps_pos";
datawindow = collect(wstart:wstart+wlen-1);
datadir = joinpath(topdir,trialstr);
rangedir = joinpath(topdir, "rangeOnly_"*trialstr)

mfin = zeros(7501,2,wlen);
posData = zeros(wlen,2);
for i in 1:wlen
    dataindex = datawindow[i];
    navfile = datadir*"/nav$(dataindex).csv"
    posData[i,:] = readdlm(navfile,',',Float64,'\n')
    rangeFile = rangedir*"/range$(dataindex).txt";
    mfin[:,:,i] = readdlm(rangeFile,',',Float64,'\n')
end

dposData = deepcopy(posData)
cumulativeDrift!(dposData,[0.0;0],[0.2,0.2])

fg = initfg();

poses = [Symbol("x$i") for i in 1:wlen]
for sym in poses
  addVariable!(fg, sym, Point2)
end

priors = [1,wlen]
rtkCov = Matrix(Diagonal([0.1;0.1].^2))
#Priors
for i in priors
    pp = PriorPoint2(MvNormal(posData[i,:], rtkCov))
    addFactor!(fg, [Symbol("x$i");], pp, autoinit=true)
end

vps = 1:wlen-1
dpσ = Matrix(Diagonal([0.2;0.2].^2))
for i in vps
    dx = dposData[i+1,1] - dposData[i,1];
    dy = dposData[i+1,2] - dposData[i,2];
    p2p2 = Point2Point2(MvNormal([dx;dy],dpσ))
    addFactor!(fg, [Symbol("x$i");Symbol("x$(i+1)")], p2p2, autoinit=true)
end

beacon = :l1
addVariable!(fg, beacon, Point2 )
# manualinit!(fg,beacon,kde!(rand(MvNormal([0;0],Matrix(Diagonal([7.0;7].^2))),100)))

rangewindow = 1:3:wlen
for i in rangewindow
    sym = Symbol("x$i")
    ppR = Point2Point2Range(AliasingScalarSampler(mfin[:,1,i],exp.(mfin[:,2,i]),SNRfloor=0.8))
    addFactor!(fg, [sym;beacon], ppR, autoinit=false)
end

# writeGraphPdf(fg,viewerapp="", engine="neato", filepath="/tmp/fg.pdf")
getSolverParams(fg).drawtree = true
getSolverParams(fg).showtree = false
getSolverParams(fg).dbg = true

ensureAllInitialized!(fg)
tree, smt, hist = solveTree!(fg, maxparallel=200, recordcliqs=[:x19;:x169;:x166;:x134;:x171])

fetchCliqTaskHistoryAll!(smt, hist)
drawTree(tree, show=true)

#Plotting

plk= [];
for sym in poses[end-5:end]  #plot last 5 poses with contour
    X1 = getKDEMean(getVertKDE(fg,sym))
    push!(plk, layer(x=[X1[1];],y=[X1[2];], label=["$(sym)";], Geom.point, Geom.label, Theme(default_color=colorant"blue",point_size = 1.5pt,highlight_width = 0pt)))
    K1 = plotKDEContour(getVertKDE(fg,sym),xlbl="", ylbl="",levels=2,layers=true);
    push!(plk,K1...)
    push!(plk,Gadfly.Theme(key_position = :none));
end
for sym in poses #plotting all syms no labels
    X1 = getKDEMean(getVertKDE(fg,sym))
    push!(plk, layer(x=[X1[1];],y=[X1[2];], label=["$(sym)";], Geom.point, Theme(default_color=colorant"blue",point_size = 1.5pt,highlight_width = 0pt)))
    # K1 = plotKDEContour(getVertKDE(fg,sym),xlbl="", ylbl="",levels=2,layers=true);
    # push!(plk,K1...)
    # push!(plk,Gadfly.Theme(key_position = :none));
end
push!(plk,layer(x=posData[:,1],y=posData[:,2], Geom.path, Theme(default_color=colorant"green2")))
push!(plk,layer(x=dposData[:,1],y=dposData[:,2], Geom.path, Theme(default_color=colorant"blue2")))

# push!(plk, Coord.cartesian(xmin=-10, xmax=120, ymin=-120, ymax=-30,fixed=true)) #full plot
push!(plk, Coord.cartesian(xmin=95, xmax=105, ymin=-81, ymax=-70,fixed=true)) #crop plot to last points
plkplot = Gadfly.plot(plk...); plkplot |> PDF("/tmp/test.pdf");

# igt = [17.0499;1.7832];
#
# push!(plk,layer(x=[igt[1];],y=[igt[2];], label=String["Beacon Ground Truth";],Geom.point,Geom.label(hide_overlaps=false), order=2, Theme(default_color=colorant"red",highlight_width = 0pt)));
#
# X1 = approxConv(fg,:x4l1f1,:l1,N=200)
# push!(plk,layer(x=X1[1,:],y=X1[2,:],Geom.histogram2d(xbincount=400, ybincount=400)))
# K1 = plotKDEContour(kde!(X1),xlbl="X (m)", ylbl="Y (m)",levels=4,layers=true);
# push!(plk,K1...)
# push!(plk,Gadfly.Theme(key_position = :none));

# L1 = getVal(getVariable(fg, beacon))
# L1 = rand(getVertKDE(fg,:l1),1000)
# push!(plk,layer(x=L1[1,:],y=L1[2,:],Geom.histogram2d(xbincount=300, ybincount=300)))
# K1 = plotKDEContour(getVertKDE(fg,:l1),xlbl="X (m)", ylbl="Y (m)",levels=6,layers=true);
# push!(plk,K1...)
# push!(plk,Gadfly.Theme(key_position = :none));
# push!(plk, Guide.xlabel("X (m)"),Guide.ylabel("Y (m)"))
# push!(plk, Coord.cartesian(xmin=35, xmax=45, ymin=-50, ymax=-35,fixed=true))

# Check the range distribution
bss = AliasingScalarSampler(mfin[:,1,4],exp.(mfin[:,2,4]),SNRfloor=0.8)
out = rand(bss,1000)
Gadfly.plot(x=out,Geom.histogram,Coord.cartesian(xmin=0, xmax=100, ymin=0, ymax=400,fixed=true))

## dehann dev

using RoMEPlotting

# tree = wipeBuildNewTree!(fg, maxparallel=200)
# spyCliqMat(tree, :x169)

drawTree(tree, show=true)

plotTreeProductDown(fg, tree, :x169, :x169)
plotTreeProductDown(fg, tree, :x169, :x170)

#
getCliq(tree, :x169)

# fldr = "2019-09-01T12:25:09.297"
fldr = "2019-09-08T22:16:41.83"
# fg31_au = loadDFG("/tmp/caesar/"*fldr*"/cliqSubFgs/cliq31/fg_afterupsolve", Main)
# fg31_bu = loadDFG("/tmp/caesar/"*fldr*"/cliqSubFgs/cliq31/fg_beforeupsolve", Main)

fg31bd = initfg()
loadDFG("/tmp/caesar/"*fldr*"/cliqSubFgs/cliq31/fg_beforedownsolve", Main, fg31bd)
# fg31_bd =

writeGraphPdf(fg31bd, show=true)

plotLocalProduct(fg31bd, :x169)
plotLocalProduct(fg31bd, :x170)

printCliqHistorySummary(hist, tree, :x169)

sfg31bd = buildSubgraphFromLabels(fg31bd, [:x167;:x168;:x169;:x170])

writeGraphPdf(sfg31bd, show=true)

tree31s, smt31s, hist31s = solveTree!(sfg31bd)

drawTree(tree, show=true, imgs=false)

plotKDE(sfg31bd, ls(sfg31bd), levels=1)

plotLocalProduct(sfg31bd, :x169)
plotLocalProduct(sfg31bd, :x168)

# need a new factor plot debug function
#

plotVariableGivenFactor(sfg31bd,:x168x169f1,:x169)

pts = approxConv(sfg31bd,:x168x169f1, :x169)
lie = ls(sfg31bd, :x168x169f1)
setdiff!(lie, [:x169])

fct = getFactor(sfg31bd, :x168x169f1)
@show getFactorType(fct)
getData(fct).fncargvID

fct = getFactor(fg, :x168x169f1)
getData(fct).fncargvID

ste = 12
csfg31bd = hist[31][ste][4].cliqSubFg
plotKDE(csfg31bd, [:x168;:x169;:x170], levels=1)

csmfg31bd = hist[31][ste][4].dfg
plotKDE(csmfg31bd, [:x168;:x169;:x170], levels=1)

getData(getFactor(csmfg31bd, :x168x169f1)).fncargvID
getData(getFactor(csfg31bd, :x168x169f1)).fncargvID

stuff = sandboxCliqResolveStep(hist[31], 10)

writeGraphPdf(csfg21bd, show=true)

getCliq(tree, :x170)

# what happened in the root

printCliqHistorySummary(hist, tree, :x170)

ste = 9
csfg1bd = hist[1][ste][4].cliqSubFg
plotKDE(csfg31bd, [:x168;:x169;:x170], levels=3)

plotKDE(getCliqDownMsgsAfterDownSolve(fg, getCliq(tree,:x170))[:x170], levels=3)
plotKDE(getCliqDownMsgsAfterDownSolve(fg, getCliq(tree,:x19))[:x170], levels=3)

plotTreeProductDown(fg, tree, :x169, :x169)

#


text:
################################################################################
#
#  Map/AbGrp.jl : Types for maps with domains of type GrpAbFinGen
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2015 - 2018 Tommy Hofmann, Carlo Sircana
#
################################################################################

export GrpAbFinGenMap

################################################################################
#
#  Morphisms between finitely generated abelian groups
#
################################################################################

mutable struct GrpAbFinGenMap <: Map{GrpAbFinGen, GrpAbFinGen,
                                     HeckeMap, GrpAbFinGenMap}
  header::MapHeader{GrpAbFinGen, GrpAbFinGen}

  map::fmpz_mat
  imap::fmpz_mat
  im::GrpAbFinGen  # if set
  ke::GrpAbFinGen  # if set

  function GrpAbFinGenMap(G::GrpAbFinGen)
    r = new()
    r.header = MapHeader(G, G)
    r.map = identity_matrix(FlintZZ, ngens(G))
    r.imap = identity_matrix(FlintZZ, ngens(G))
    return r
  end

  function GrpAbFinGenMap(From::GrpAbFinGen, To::GrpAbFinGen, M::fmpz_mat)
    r = new()
    r.header = MapHeader(From, To)
    r.map = M
    return r
  end

  function GrpAbFinGenMap(From::GrpAbFinGen, To::GrpAbFinGen, M::fmpz_mat, Mi::fmpz_mat)
    r = new()
    r.header = MapHeader(From, To)
    r.map = M
    r.imap = Mi
    return r
  end

  function GrpAbFinGenMap(M::T) where T <: Map{GrpAbFinGen, GrpAbFinGen}
    r = new()
    D = domain(M)
    r.header = MapHeader(D, codomain(M))
    if ngens(D) == 0
      r.map = matrix(FlintZZ, 0, ngens(codomain(M)), fmpz[])
    else
      r.map = vcat([M(D[i]).coeff for i=1:ngens(D)])
    end
    return r
  end
end

function image(f::Map(GrpAbFinGenMap), a::GrpAbFinGenElem)
  parent(a) == domain(f) || error("not in the domain")
  if !isdefined(f, :map)
    return hasimage(f, a)[2]
  end
  return GrpAbFinGenElem(codomain(f), a.coeff * f.map)
end

function image(phi::GrpAbFinGenMap, U::GrpAbFinGen, add_to_lattice::Bool = !false)
  G = domain(phi)
  fl, inj = issubgroup(U, G)
  fl || error("subgroup is not in the domain")
  return sub(codomain(phi), [phi(inj(U[i])) for i=1:ngens(U)], add_to_lattice)
end

function preimage(f::Map(GrpAbFinGenMap), a::GrpAbFinGenElem)
  if !isdefined(f, :imap)
    fl, r = haspreimage(f, a)
    fl || error("element is not in the image")
    return r
  end
  return GrpAbFinGenElem(domain(f), a.coeff * f.imap)
end

################################################################################
#
#  Morphisms from finite abelian groups into finite fields
#
################################################################################

mutable struct FiniteFieldMultGrpMap{S, T} <: Map{GrpAbFinGen, S, HeckeMap, FiniteFieldMultGrpMap{S, T}}
  header::MapHeader{GrpAbFinGen, S}
  domain::GrpAbFinGen
  codomain::S
  generator::T
  @declare_other

  function FiniteFieldMultGrpMap{S, T}(G::GrpAbFinGen, F::S, generator::T, disc_log::Function) where {S, T}
    z = new{S, T}()
    z.header = MapHeader(G, F)
    z.domain = G
    z.codomain = F
    z.header.preimage = disc_log
    z.generator = generator
    return z
  end
end

function image(f::FiniteFieldMultGrpMap, x::GrpAbFinGenElem)
  return f.generator^x[1]
end

################################################################################
#
#  Morphisms from finite abelian groups onto units of orders
#
################################################################################

mutable struct AbToNfOrdMultGrp{S, T} <: Map{GrpAbFinGen, S, SetMap, AbToNfOrdMultGrp}
  domain::GrpAbFinGen
  codomain::S
  generator::T
  @declare_other

  function AbToNfOrdMultGrp(O::NfAbsOrd, order::Int, generator::NfAbsOrdElem)
    G = abelian_group([order])
    z = new{typeof(O), typeof(generator)}()
    z.domain = G
    z.codomain = O
    z.generator = generator
    return z
  end
  function AbToNfOrdMultGrp(O::NfRelOrd, order::Int, generator::NfRelOrdElem)
    G = abelian_group([order])
    z = new{typeof(O), typeof(generator)}()
    z.domain = G
    z.codomain = O
    z.generator = generator
    return z
  end
end

function (f::AbToNfOrdMultGrp)(a::GrpAbFinGenElem)
  return f.generator^a[1]
end

mutable struct AbToNfMultGrp{S, T} <: Map{GrpAbFinGen, S, SetMap, AbToNfMultGrp}
  domain::GrpAbFinGen
  codomain::S
  generator::T
  @declare_other

  function AbToNfMultGrp(K::NumField, order::Int, generator::NumFieldElem)
    G = abelian_group(Int[order])
    z = new{typeof(K), typeof(generator)}()
    z.domain = G
    z.codomain = K
    z.generator = generator
    return z
  end

end

function (f::AbToNfMultGrp)(a::GrpAbFinGenElem)
  return f.generator^a[1]
end

################################################################################
#
#  Morphisms from finite abelian groups to orders
#
################################################################################

# S is the type of the order (the codomain) and T is the elem_type of the order
mutable struct GrpAbFinGenToAbsOrdMap{S, T} <: Map{GrpAbFinGen, S, HeckeMap, GrpAbFinGenToAbsOrdMap}
  header::MapHeader{GrpAbFinGen, S}
  generators::Vector{T}
  discrete_logarithm::Function
  modulus # this can be anything, for which powermod(::T, ::fmpz, modulus) is defined

  disc_log::GrpAbFinGenElem #Needed in the conductor computation

  function GrpAbFinGenToAbsOrdMap{S, T}(G::GrpAbFinGen, O::S, generators::Vector{T}, disc_log::Function, modulus...) where {S, T}
    @assert ngens(G) == length(generators)

    z = new{S, T}()
    modulo = false
    if length(modulus) == 1
      modulo = true
      z.modulus = modulus[1]
    end

    function _image(a::GrpAbFinGenElem)
      @assert parent(a) === G
      y = one(O)
      for i in 1:length(generators)
        a[i] == 0 && continue
        if modulo
          y *= powermod(generators[i], a[i], z.modulus)
        else
          y *= generators[i]^a[i]
        end
      end
      return y
    end

    function _preimage(a::T)
      @assert parent(a) === O
      return G(disc_log(a))
    end

    z.header = MapHeader(G, O, _image, _preimage)
    z.generators = generators
    z.discrete_logarithm = disc_log
    return z
  end

  function GrpAbFinGenToAbsOrdMap{S, T}(O::S, generators::Vector{T}, snf_structure::Vector{fmpz}, disc_log::Function, modulus...) where {S, T}
    @assert length(generators) == length(snf_structure)

    G = abelian_group(snf_structure)

    return GrpAbFinGenToAbsOrdMap{S, T}(G, O, generators, disc_log, modulus...)
  end

  function GrpAbFinGenToAbsOrdMap{S, T}(O::S, generators::Vector{T}, relation_matrix::fmpz_mat, disc_log::Function, modulus...) where {S, T}
    @assert length(generators) == nrows(relation_matrix)

    G = GrpAbFinGen(relation_matrix)

    return GrpAbFinGenToAbsOrdMap{S, T}(G, O, generators, disc_log, modulus...)
  end
end

function GrpAbFinGenToAbsOrdMap(G::GrpAbFinGen, O::S, generators::Vector{T}, disc_log::Function, modulus...) where {S, T}
  return GrpAbFinGenToAbsOrdMap{S, T}(G, O, generators, disc_log, modulus...)
end

function GrpAbFinGenToAbsOrdMap(O::S, generators::Vector{T}, snf_structure::Vector{fmpz}, disc_log::Function, modulus...) where {S, T}
  return GrpAbFinGenToAbsOrdMap{S, T}(O, generators, snf_structure, disc_log, modulus...)
end

function GrpAbFinGenToAbsOrdMap(O::S, generators::Vector{T}, relation_matrix::fmpz_mat, disc_log::Function, modulus...) where {S, T}
  return GrpAbFinGenToAbsOrdMap{S, T}(O, generators, relation_matrix, disc_log, modulus...)
end

const GrpAbFinGenToNfAbsOrdMap = GrpAbFinGenToAbsOrdMap{NfOrd, NfOrdElem}

################################################################################
#
#  Isomorphisms from abelian groups onto multliplicative groups of residue
#  rings of orders
#
################################################################################

mutable struct GrpAbFinGenToAbsOrdQuoRingMultMap{S, T, U} <: Map{GrpAbFinGen, AbsOrdQuoRing{S, T}, HeckeMap, GrpAbFinGenToAbsOrdQuoRingMultMap}
  header::MapHeader{GrpAbFinGen, AbsOrdQuoRing{S, T}}
  generators::Vector{AbsOrdQuoRingElem{S, T, U}}
  discrete_logarithm::Function

  # Multiplicative group, tame part
  tame::Dict{T, GrpAbFinGenToAbsOrdMap{S, U}}

  # Multiplicative group, wild part
  wild::Dict{T, GrpAbFinGenToAbsOrdMap{S, U}}

  function GrpAbFinGenToAbsOrdQuoRingMultMap{S, T, U}(G::GrpAbFinGen, Q::AbsOrdQuoRing{S, T}, generators::Vector{AbsOrdQuoRingElem{S, T, U}}, disc_log::Function) where {S, T, U}
    @assert ngens(G) == length(generators)

    function _image(a::GrpAbFinGenElem)
      @assert parent(a) == G
      y = one(Q)
      for i in 1:length(generators)
        a[i] == 0 && continue
        mul!(y, y, generators[i]^a[i])
      end
      return y
    end

    function _preimage(a::AbsOrdQuoRingElem)
      @assert parent(a) === Q
      return G(disc_log(a))
    end

    z = new{S, T, U}()
    z.header = MapHeader(G, Q, _image, _preimage)
    z.generators = generators
    z.discrete_logarithm = disc_log
    return z
  end

  function GrpAbFinGenToAbsOrdQuoRingMultMap{S, T, U}(Q::AbsOrdQuoRing{S, T}, generators::Vector{AbsOrdQuoRingElem{S, T, U}}, snf_structure::Vector{fmpz}, disc_log::Function) where {S, T, U}
    @assert length(generators) == length(snf_structure)

    G = abelian_group(snf_structure)

    return GrpAbFinGenToAbsOrdQuoRingMultMap{S, T, U}(G, Q, generators, disc_log)
  end

  function GrpAbFinGenToAbsOrdQuoRingMultMap{S, T, U}(Q::AbsOrdQuoRing{S, T}, generators::Vector{AbsOrdQuoRingElem{S, T, U}}, relation_matrix::fmpz_mat, disc_log::Function) where {S, T, U}
    @assert length(generators) == nrows(relation_matrix)

    G = GrpAbFinGen(relation_matrix)

    return GrpAbFinGenToAbsOrdQuoRingMultMap{S, T, U}(G, Q, generators, disc_log)
  end

  function GrpAbFinGenToAbsOrdQuoRingMultMap{S, T, U}(G::GrpAbFinGen, Q::AbsOrdQuoRing{S, T}, exp::Function, disc_log::Function) where {S, T, U}
    z = new{S, T, U}()
    z.header = MapHeader(G, Q, exp, disc_log)
    return z
  end
end

const GrpAbFinGenToNfOrdQuoRingMultMap = GrpAbFinGenToAbsOrdQuoRingMultMap{NfOrd, NfOrdIdl, NfOrdElem}

function GrpAbFinGenToAbsOrdQuoRingMultMap(G::GrpAbFinGen, Q::AbsOrdQuoRing{S, T}, generators::Vector{AbsOrdQuoRingElem{S, T, U}}, disc_log::Function) where {S, T, U}
  return GrpAbFinGenToAbsOrdQuoRingMultMap{S, T, U}(G, Q, generators, disc_log)
end

function GrpAbFinGenToAbsOrdQuoRingMultMap(Q::AbsOrdQuoRing{S, T}, generators::Vector{AbsOrdQuoRingElem{S, T, U}}, snf_structure::Vector{fmpz}, disc_log::Function) where {S, T, U}
  return GrpAbFinGenToAbsOrdQuoRingMultMap{S, T, U}(Q, generators, snf_structure, disc_log)
end

function GrpAbFinGenToAbsOrdQuoRingMultMap(Q::AbsOrdQuoRing{S, T}, generators::Vector{AbsOrdQuoRingElem{S, T, U}}, relation_matrix::fmpz_mat, disc_log::Function) where {S, T, U}
  return GrpAbFinGenToAbsOrdQuoRingMultMap{S, T, U}(Q, generators, relation_matrix, disc_log)
end

function GrpAbFinGenToAbsOrdQuoRingMultMap(G::GrpAbFinGen, Q::AbsOrdQuoRing{S, T}, exp::Function, disc_log::Function) where {S, T}
  U = elem_type(base_ring(Q))
  return GrpAbFinGenToAbsOrdQuoRingMultMap{S, T, U}(G, Q, exp, disc_log)
end


text:
using bearingDistributions
using Interact
using Plots
using Test

const bd = bearingDistributions

@testset "bearingDistributions.jl" begin
    b1 = Bearing(atan(1,1),  .1, 0, 0);
    b2 = Bearing(atan(-1,1), .2, 0, 1);
    @test b1.x ≈ 0
    @test b1.y ≈ 0
    @test b1.θ ≈ atan(1,1)
    @test b1.σ ≈ .1
    @test b2.x ≈ 0
    @test b2.y ≈ 1
    @test b2.θ ≈ atan(-1,1)
    @test b2.σ ≈ .2
    @test bd.findIntersect(b1, b2) ≈ [.5;.5]

    b1 = Bearing(0, .1, 0, 0);
    b2 = Bearing(π, .1, 1, 0);
    @test bd.findIntersect(b1, b2) ≈ [.5, 0]

    b1 = Bearing(π, .1, 0, 0);
    b2 = Bearing(π/2, .1, 1, 0);
    @test bd.findIntersect(b1, b2) ≈ [.5, 0]

    b1 = Bearing(0, 1, 0, 1);
    b2 = Bearing(π/2, .1, 1, 0);
    @test bd.boundingBox(b1, b2) == [-.5 2.5; -.5 2.5]

    @test bd.rangeFromBox([0 2; 0 2]; length = 2) == ([0,2],[0,2])

    @test isnan(bisectionRoots((x)-> x.^2+1, -1, 1))

    @test bisectionRoots((x) -> x+1, -2, 0; xtol = 1e-9) ≈ -1

    @test bisectionRoots((x) -> x+1, -1, 0) ≈ -1
    @test bisectionRoots((x) -> x+1, -2, -1) ≈ -1

    b1 = Bearing(0, .1, 0, 1);
    b2 = Bearing(π/2, .1, 1, 0);
    len = 11;
    α = .95;
    @test typeof(plotConfidenceInterval(α,bd.intersectGrid(b1, b2; length=len)...)) <: AbstractPlot

    @test typeof(App.app([])) <: Widgets.Widget{:manipulate,Any}

    #t = App.launchApp(8001)
    #@test typeof(t) <: Task
    #@test t.state = :runnable
    #t.state = :done
    #@test t.state = :done

    #@test julia_main() == 0
end


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    AbstractBridge

Subtype of [`MathOptInterface.Bridges.AbstractBridge`](@ref) for constraint
bridges.
"""
abstract type AbstractBridge <: MOIB.AbstractBridge end

"""
    bridge_constraint(
        BT::Type{<:AbstractBridge},
        model::MOI.ModelLike,
        func::AbstractFunction,
        set::MOI.AbstractSet,
    )

Bridge the constraint `func`-in-`set` using bridge `BT` to `model` and returns
a bridge object of type `BT`. The bridge type `BT` should be a concrete type,
that is, all the type parameters of the bridge should be set. Use
[`concrete_bridge_type`](@ref) to obtain a concrete type for given function
and set types.
"""
function bridge_constraint end

"""
    MOI.get(b::AbstractBridge, ::MOI.NumberOfVariables)

The number of variables created by the bridge `b` in the model.
"""
MOI.get(::AbstractBridge, ::MOI.NumberOfVariables)::Int64 = 0

"""
    MOI.get(b::AbstractBridge, ::MOI.ListOfVariableIndices)

The list of variables created by the bridge `b` in the model.
"""
MOI.get(::AbstractBridge, ::MOI.ListOfVariableIndices) = MOI.VariableIndex[]

"""
    MOI.supports_constraint(
        BT::Type{<:AbstractBridge},
        F::Type{<:MOI.AbstractFunction},
        S::Type{<:MOI.AbstractSet},
    )::Bool

Return a `Bool` indicating whether the bridges of type `BT` support bridging
`F`-in-`S` constraints.
"""
function MOI.supports_constraint(
    ::Type{<:AbstractBridge},
    ::Type{<:MOI.AbstractFunction},
    ::Type{<:MOI.AbstractSet},
)
    return false
end

"""
    added_constrained_variable_types(
        BT::Type{<:MOI.Bridges.Constraint.AbstractBridge},
        F::Type{<:MOI.AbstractFunction},
        S::Type{<:MOI.AbstractSet},
    )

Return a list of the types of constrained variables that bridges of type `BT`
add for bridging `F`-in-`S` constraints. This falls back to
`added_constrained_variable_types(concrete_bridge_type(BT, F, S))` so bridges
should not implement this method.
"""
function MOIB.added_constrained_variable_types(
    BT::Type{<:AbstractBridge},
    F::Type{<:MOI.AbstractFunction},
    S::Type{<:MOI.AbstractSet},
)
    return MOIB.added_constrained_variable_types(concrete_bridge_type(BT, F, S))
end

"""
    added_constraint_types(
        BT::Type{<:MOI.Bridges.Constraint.AbstractBridge},
        F::Type{<:MOI.AbstractFunction},
        S::Type{<:MOI.AbstractSet},
    )

Return a list of the types of constraints that bridges of type `BT` add for
bridging `F`-in-`S` constraints. This falls back to
`added_constraint_types(concrete_bridge_type(BT, F, S))`
so bridges should not implement this method.

This function should not repeat the constraints added on variables when they
are created (this is the role of `added_constrained_variable_types`), but only
the constraints that are added with `MOI.add_constraint`.
"""
function MOIB.added_constraint_types(
    BT::Type{<:AbstractBridge},
    F::Type{<:MOI.AbstractFunction},
    S::Type{<:MOI.AbstractSet},
)
    return MOIB.added_constraint_types(concrete_bridge_type(BT, F, S))
end

"""
    concrete_bridge_type(
        BT::Type{<:AbstractBridge},
        F::Type{<:MOI.AbstractFunction},
        S::Type{<:MOI.AbstractSet}
    )::Type

Return the concrete type of the bridge supporting `F`-in-`S` constraints. This
function can only be called if `MOI.supports_constraint(BT, F, S)` is `true`.

## Examples

As a [`MathOptInterface.VariableIndex`](@ref)-in-[`MathOptInterface.Interval`](@ref)
constraint is bridged into a
[`MathOptInterface.VariableIndex`](@ref)-in-[`MathOptInterface.GreaterThan`](@ref)
and a
[`MathOptInterface.VariableIndex`](@ref)-in-[`MathOptInterface.LessThan`](@ref)
by the [`SplitIntervalBridge`](@ref):

```jldoctest; setup=:(using MathOptInterface; const MOI = MathOptInterface)
MOI.Bridges.Constraint.concrete_bridge_type(
    MOI.Bridges.Constraint.SplitIntervalBridge{Float64},
    MOI.VariableIndex,
    MOI.Interval{Float64},
)

# output

MathOptInterface.Bridges.Constraint.SplitIntervalBridge{Float64, MathOptInterface.VariableIndex, MathOptInterface.Interval{Float64}, MathOptInterface.GreaterThan{Float64}, MathOptInterface.LessThan{Float64}}
```
"""
function concrete_bridge_type(
    bridge_type::Type,
    ::Type{<:MOI.AbstractFunction},
    ::Type{<:MOI.AbstractSet},
)
    return bridge_type
end

function concrete_bridge_type(
    b::MOIB.AbstractBridgeOptimizer,
    F::Type{<:MOI.AbstractFunction},
    S::Type{<:MOI.AbstractSet},
)
    return concrete_bridge_type(MOIB.bridge_type(b, F, S), F, S)
end

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
@ghdef mutable struct Team
    name::Union{String, Nothing}
    description::Union{String, Nothing}
    privacy::Union{String, Nothing}
    permission::Union{String, Nothing}
    slug::Union{String, Nothing}
    id::Union{Int, Nothing}
end

namefield(t::Team) = t.id

@api_default function members(api::GitHubAPI, team; options...)
    results, page_data = gh_get_paged_json(api, "/teams/$(name(team))/members"; options...)
    return map(Owner, results), page_data
end

@api_default function repos(api::GitHubAPI, owner::Owner, team::Team; options...)
    results, page_data = gh_get_paged_json(api, "/orgs/$(owner.login)/teams/$(team.slug)/repos"; options...)
    return map(Repo, results), page_data
end


text:
include(joinpath("..", "common", "utils.jl"))
using LinearAlgebra: Diagonal, diag

⊗ = kron

function element_operators(N, T = Float64)
    ξ, ω = lglpoints(BigFloat, N)
    D = spectralderivative(ξ)
    return (ξ = T.(ξ), ω = T.(ω), D = T.(D))
end

wave_speed(ρ, p, γ) = sqrt(γ * p ./ ρ)

function create_mesh(K, elem, z0, z1)
    ξ = elem.ξ

    # Number of pointss in element and polynomial order
    Nq = length(ξ)
    N = Nq - 1

    # Element size
    Δz = (z1 - z0) / K
    J = Δz / (ξ[Nq] - ξ[1])

    # shift ξ to go (0, 1)
    ξ01 = (ξ[1:(end - 1)] .- ξ[1]) / (ξ[Nq] - ξ[1])

    # cg to dg scatter matrix
    Q = scatter_matrix(N, K)

    # CG DOF locations
    zcg = [z0 .+ Δz * (ξ01 .+ (0:(K - 1))')[:]; z1]

    # DG DOF locations
    zdg = reshape(Q * zcg, Nq, K)

    W = Diagonal(elem.ω)
    I_KK = sparse(I, K, K)
    Wcg = Array(diag(Q' * (I_KK ⊗ (J * W)) * Q))

    return (zcg = zcg, zdg = zdg, scatter = Q, J = J, Δz = Δz, Wcg = Wcg)
end

function dg2cg_scatter(q, mesh, elem)
    Q = mesh.scatter
    ω = elem.ω
    if q isa Array
        q[:] .= Q * ((Q' * (mesh.J * ω .* q)[:]) ./ mesh.Wcg)
    elseif q isa Tuple || q isa NamedTuple
        for i in 1:length(q)
            q[i][:] .= Q * ((Q' * (mesh.J * ω .* q[i])[:]) ./ mesh.Wcg)
        end
    end
    return q
end

function legendre_polynomial(ξ, P)
    T = eltype(ξ)
    N = length(ξ) - 1
    coef = GaussQuadrature.legendre_coefs(T, max(1, P))
    v = GaussQuadrature.orthonormal_poly(ξ, coef...)[:, P + 1]
    return v
end

function balanced_rho(p, grav, mesh, elem)
    Nq = length(elem.ξ)
    N = Nq - 1

    # Compute the derivative
    ρ = -elem.D * p / (mesh.J * grav)
    ρ0 = copy(ρ)

    # make continuous by adding a higher mode on each element
    v = legendre_polynomial(elem.ξ, N)
    v ./= v[1]
    for e in 2:size(ρ, 2)
        Δρ = ρ[1, e] - ρ[Nq, e - 1]
        ρ[:, e] .-= Δρ * v
    end

    return ρ
end

"""
    element_tendency!(∂q, q, elem_operator, problem, z_elem)

Evaluate the element rhs `∂q` associated with `q` for the
`problem::AbstractProblem` using the DG `elem_operator`.
`q` and `∂q` are taken to be `NamedTuple`s with fields `(ρ, w, E)`
"""
function element_tendency!(∂q, q, O, J, z, grav, p_, Pr = I)
    ρ, w, ρe = q.ρ, q.w, q.ρe
    p = p_.(ρ, w, ρe, z)

    # ∫ ϕ J ∂ρ = ∫ (∂_ξ ϕ) ρ w
    ∂q.ρ .= (O.D' * (O.ω .* ρ .* w)) ./ (J .* O.ω)
    # ∂q.ρ .= (ρ .* (O.D' * (O.ω .* w))) ./ (J .* O.ω)
    # ∂q.ρ .-= w .* (O.D * ρ) ./ J

    #  ∫ ϕ ∂w = - ∫ ϕ (g - (∂_ξ p) / ρ̂) + ∫ (∂_ξ ϕ) w^2/2
    #  where ρ̂ has the projection operator Pr applied
    ∂q.w .= -grav .- (O.D * p) ./ (J .* Pr * ρ)

    #  ∫ ϕ ∂w = - ∫ ϕ ((g ρ̂  - ∂_ξ p) / ρ) + ∫ (∂_ξ ϕ) w^2/2
    #  where ρ̂ has the projection operator Pr applied
    # ∂q.w .= -(grav * (Pr * ρ) + (O.D * p) / J) ./ ρ
    ∂q.w .+= (O.D' * (O.ω .* w .^ 2 / 2)) ./ (J .* O.ω)

    #  ∫ ϕ ∂ρe = ∫ (∂_ξ ϕ) (ρe + p) w
    ∂q.ρe .= (O.D' * (O.ω .* (ρe + p) .* w)) ./ (J .* O.ω)

    return ∂q
end

function tendency!(
    ∂q,
    q,
    t,
    mesh,
    elem,
    grav,
    pres,
    Pr = I,
    bc = (nothing, nothing),
    forcing = nothing,
)
    # Evaluate the volume terms
    element_tendency!(∂q, q, elem, mesh.J, mesh.zdg, grav, pres, Pr)

    # Add in MMS formcing
    if !isnothing(forcing)
        ∂q.δρ .+= forcing.ρ.(mesh.zdg, t)
        ∂q.δw .+= forcing.w.(mesh.zdg, t)
        ∂q.δE .+= forcing.ρe.(mesh.zdg, t)
    end

    # Use upwind boundary treatment
    if !isnothing(bc[1])
        error("not implemented yet")
    end

    if !isnothing(bc[2])
        error("not implemented yet")
    end

    # DG -> CG projection and then scatter for storage
    dg2cg_scatter(∂q, mesh, elem)
end


text:
"""
Simple server that implements [server-sent events](https://developer.mozilla.org/en-US/docs/Web/API/Server-sent_events),
loosely following [this tutorial](https://developer.mozilla.org/en-US/docs/Web/API/Server-sent_events/Using_server-sent_events).

### Example client code (JS):
```http
<html>
<head>
    <meta charset="UTF-8">
    <title>Server-sent events demo</title>
</head>
<body>
    <h3>Fetched items:</h3>
    <ul id="list"></ul>
</body>
<script>
    const evtSource = new EventSource("http://127.0.0.1:8080/api/events")
    evtSource.onmessage = async function (event) {
        const newElement = document.createElement("li");
        const eventList = document.getElementById("list");
        if (parseFloat(event.data) > 0.5) {
            const r = await fetch("http://127.0.0.1:8080/api/getItems")
            if (r.ok) {
                const body = await r.json()
                newElement.textContent = body;
                eventList.appendChild(newElement);
            }
        }
    }
    evtSource.addEventListener("ping", function(event) {
        console.log('ping:', event.data)
    });
</script>
</html>
```

### Example client code (Julia)
```julia
using HTTP, JSON

HTTP.open("GET", "http://127.0.0.1:8080/api/events") do io
    while !eof(io)
        println(String(readavailable(io)))
    end
end
```

### Server code:
"""
using HTTP, JSON

const ROUTER = HTTP.Router()

function getItems(req::HTTP.Request)
    headers = [
        "Access-Control-Allow-Origin" => "*",
        "Access-Control-Allow-Methods" => "GET, OPTIONS"
    ]
    if HTTP.method(req) == "OPTIONS"
        return HTTP.Response(200, headers)
    end
    return HTTP.Response(200, headers; body = JSON.json(rand(2)))
end

function events(stream::HTTP.Stream)
    HTTP.setheader(stream, "Access-Control-Allow-Origin" => "*")
    HTTP.setheader(stream, "Access-Control-Allow-Methods" => "GET, OPTIONS")
    HTTP.setheader(stream, "Content-Type" => "text/event-stream")

    if HTTP.method(stream.message) == "OPTIONS"
        return nothing
    end

    HTTP.setheader(stream, "Content-Type" => "text/event-stream")
    HTTP.setheader(stream, "Cache-Control" => "no-cache")
    while true
        write(stream, "event: ping\ndata: $(round(Int, time()))\n\n")
        if rand(Bool)
            write(stream, "data: $(rand())\n\n")
        end
        sleep(1)
    end
    return nothing
end

HTTP.@register(ROUTER, "GET", "/api/getItems", getItems)
HTTP.@register(ROUTER, "/api/events", HTTP.Handlers.StreamHandlerFunction(events))

HTTP.serve(ROUTER, "127.0.0.1", 8080)


text:
@info "MTKHelpers: loading CairoMakie utils"
using .CairoMakie # syntax by Requires.jl otherwise warning
using DiffEqBase: AbstractODESolution

# pdf_figure moved to TwPrototypes

"""
    series_sol!(ax, sol::AbstractODESolution, vars; tspan=extrema(sol.t), labels=string.(vars), nt=120, kwargs...)

calls `CairoMakie.series` for a grid fo `n` points and interpolated values
from `sol`.
Currently works only with solutions created by a non-composite solver, e.g. `Tsit5`.
"""
function series_sol!(ax, sol::AbstractODESolution, vars; tspan=extrema(sol.t), labels=string.(vars), nt=120, kwargs...)
    #ts = first(tspan) .<= sol.t .<= last(tspan)
    #series!(ax, sol.t[ts], transpose(VectorOfArray(sol[vars])[ts,:]); labels, kwargs...)
    tsol = sol.t[first(tspan) .<= sol.t .<= last(tspan)]
    ts = range(first(tspan), last(tspan), length=nt)
    ts2 = sort!(vcat(tsol,ts))
    series!(ax, ts2, VectorOfArray(sol(ts2, idxs=vars).u); labels, kwargs...)
end


text:
using Test

using HalfEdges.BoundingVolumeTrees
using Base.Iterators
using AbstractTrees
using StaticArrays

const Vector3{T} = SVector{3,T}
using HalfEdges.BoundingVolumeTrees: contains

@testset "traverse" begin
  t = Node(1, Node(-2, Node(-3),Node(3,Node(-4),nothing)), Node(2,nothing,Node(4)))
  @test length(t) == 7
  @test [data(i) for i in t] == [1,-2,-3,3,-4,2,4]
  @test [i for i in t] == reduce(vcat,DepthFirst(t))
  @test drop(DepthFirst(t),2) |> first |> data == -3
  @test drop(BreadthFirst(t),2) |> first |> data == 2
  @test data.(reduce(vcat,BreadthFirst(t))) == [1, -2, 2, -3, 3, 4, -4]
  @test search(iseven∘data,t) |> data == -2

  animaltree = reduce(insert,["dog","dawg","kawt","awnt"]; init = tree("cat"))
  @test map(data, leaves(animaltree)) == ["awnt","dawg","kawt"]
  @test tree([key(n) for n in animaltree]) == animaltree

  t = tree(["cat", "awnt", "aardvark", "ant", "dog", "dag", "zat", "doog"])
  @test leaves(t, x->length(key(x)) < 4) |> first |> key == "dag"

  it = itree(["cat", "awnt", "aardvark", "ant", "dog", "dag", "zat", "doog"])
  @test map(data, leaves(it)) == map(data, leaves(t))
  @test map(data, it) == map(data, t)

  build_rosetta( (n, l, r ) ) = Node(n, build_rosetta(l), build_rosetta(r))
  build_rosetta( n::Nothing ) = nothing
  build_rosetta( n::Int ) = Node(n)
  rosetta = build_rosetta((1, (2, (4, 7, nothing), 5), (3, (6, 8, 9), nothing)))

  @test reduce( *, Traverse(rosetta, x->string(" ", data(x)), nothing, nothing)) ==
        " 1 2 4 7 5 3 6 8 9"
  @test reduce( *, Traverse(rosetta, nothing, nothing, x->string(" ", data(x)))) ==
        " 7 4 2 5 1 8 6 9 3"
  @test reduce( *, Traverse(rosetta, nothing, x->string(" ", data(x)), nothing)) ==
        " 7 4 5 2 8 9 6 3 1"

end

@testset "delete" begin
  t = tree(["cat", "dog","dawg","kawt","awnt", "doog", "zat", "aardvark"])
  @test map(key, leaves(t)) == ["aardvark", "dawg", "doog", "zat"]
  @test map(key, leaves(delete(t,"awnt"))) == ["aardvark", "dawg", "doog", "zat"]
  @test map(key, leaves(delete(t,"aardvark"))) == ["awnt", "dawg", "doog", "zat"]
  @test map(key, leaves(delete(t,"dog"))) == ["aardvark", "dawg", "zat"]
  @test map(key, leaves(delete(t,"cat"))) == ["aardvark", "doog", "zat"]
  @test map(key, leaves(delete(t,"diamond"))) == map(key, leaves(t))

  animals = ["cat", "ant", "emu", "dog", "zat", "aawk", "bat"]
  @test length(delete(delete(tree(animals),"cat"), "dog")) == length(tree(animals))-2
end

@testset "balanced" begin
  animals = ["cat", "dog","dawg","kawt","awnt", "doog", "zat", "aardvark"]
  @test tree((x->AVLData{String}(x,0)).(animals)) |> height == 3

  linear255 = tree(1:255|>collect)
  balanced255 = avltree(1:255|>collect)
  ibalanced255 = avlitree(1:255|>collect)
  @test height(linear255) == 255
  @test height(balanced255) == 7
  @test height(ibalanced255) == 7

  b15 = avltree(1:15|>collect)
  sb15 = reduce((t,i)->delete(t,i),[6];init=b15)|>left|>right
  @test sb15|>height > sb15|>left|>height

  bd = reduce((t,i)->delete(t,i),[12,6,2,13,4];init=b15)
  @test height(bd) == 3

  b31 = avltree(1:31|>collect)
  #seq = [16, 17, 15, 18, 5, 29, 20, 22, 17, 8, 29, 26, 11, 4]
  #@test length(bd) >= (31-14)
  #@test height(bd) == 4

  for i in 1:2000
    bd = reduce((acc,i)->delete(acc,rand(1:31)), 1:14; init = b31)
    @test length(bd) >= (31-14)
    @test height(bd) == 4
  end

  @test height(reduce((acc,_)->delete(acc,rand(1:255)), 1:127; init = balanced255)) == 7
  @test height(reduce((acc,i)->delete(acc,i), 1:129; init = balanced255)) == 6
end

@testset "traverse" begin
  t = Node(1, Node(-2, Node(-3),Node(3,Node(-4),nothing)), Node(2,nothing,Node(4)))
  @test length(t) == 7
  @test [data(i) for i in t] == [1,-2,-3,3,-4,2,4]
  @test [i for i in t] == reduce(vcat,DepthFirst(t))
  @test drop(DepthFirst(t),2) |> first |> data == -3
  @test drop(BreadthFirst(t),2) |> first |> data == 2
  @test data.(reduce(vcat,BreadthFirst(t))) == [1, -2, 2, -3, 3, 4, -4]
  @test search(iseven∘data,t) |> data == -2

  animaltree = reduce(insert,["dog","dawg","kawt","awnt"]; init = tree("cat"))
  @test map(data, Leaves(animaltree)) == ["awnt","dawg","kawt"]
end

@testset "AABB" begin
  a = Vector3(1.0,1,1)
  o = Vector3(0.0,0,0)
  @test overlaps(AABB(o,a),AABB(0.5*a,2.0*a)) == true
  @test overlaps(AABB(o,a),AABB(10.0*a,20.0*a)) == false
  @test volume(AABB(o,a)) == 1.0
  @test volume(inflate(AABB(o,a),1.0)) == 3.0^3
  @test contains(AABB(o,a),AABB(a*0.5,a*0.75)) == true
  @test contains(AABB(o,a),AABB(a,2.0*a)) == false
  @test contains(AABB(o,a),AABB(o,a)) == true
end

function aabbtree(extents, treebuilder)
  extents = map(x->(Vector3{Float64}(x[1]),Vector3{Float64}(x[2])), extents)
  treebuilder(map( ((n,x),i)->AABBNodeData(AABB(n,x),i), extents, 1:length(extents)))
end

boxes = [([-4,-1,-1],[-1,1,1]),([1,-1,-1],[4,1,1]),
         ([2,-1,-1],[3,1,1]), ([-10,-1,-1],[10,1,1]),
         ([11,-1,-1],[12,1,1]),([7,-1,-1],[8,1,1]),
         ([7.5,-1,-1],[7.75,1,1])]

@testset "AABBTree" begin
  abba = aabbtree(boxes, tree)
  babba = aabbtree(boxes, avltree)
  iabba = aabbtree(boxes, itree)
  biabba = aabbtree(boxes, avlitree)
  @test contains((abba |> data |> x->x.aabb), (abba |> left |> data |> x->x.aabb))
  @test contains((abba |> data |> x->x.aabb), (abba |> right |> data |> x->x.aabb))
  @test Leaves(iabba) |> collect |> x->map(data,x) == Leaves(abba) |> collect |> x->map(data,x)
  @test Leaves(biabba) |> collect |> x->map(data,x) == Leaves(babba) |> collect |> x->map(data,x)
  @test contains((babba |> key ), (abba |> left |> key ))
  @test contains((babba |> key ), (abba |> right |> key ))
  hits = query(abba, AABB(Vector3(9.0,-10.0,-10.0),Vector3(20.0,10.0,10.0)))
  @test setdiff(hits, [4,5]) |> isempty

  t = tree(AABBNodeData(AABB(Vector3(-1.0,-1.0,-1.0), Vector3(1.0,1.0,1.0)),1))
  t = insert(t,AABBNodeData(randAABB(),2))
  t = insert(t,AABBNodeData(randAABB(),3))
  t = insert(t,AABBNodeData(randAABB(),4))
  t = insert(t,AABBNodeData(randAABB(),5))
  t = insert(t,AABBNodeData(AABB(Vector3(-1.0,-1.0,-1.0), Vector3(-0.9, -0.9, -0.9)),6))
  t = insert(t,AABBNodeData(AABB(Vector3(0.9,0.9,0.9), Vector3(1.0,1.0,1.0)),7))
  getaabb(x) = x.aabb
  getaabb(x::Nothing) = x
  @test contains((t |> data |> getaabb), (t |> left |> data |> getaabb))
  @test contains((t |> data |> getaabb), (t |> left |> right |> data |> getaabb))
  @test contains((t |> left |> data |> getaabb), (t |> left |> right |> data |> getaabb))

end

@testset "leaky" begin
  @test Traverse(itree([1,3,2,4]), index) |> collect |> sort == [1,2,3,4]
  alli = Traverse(aabbtree(boxes, itree), index) |> collect |> sort
  @test alli == collect(1:length(alli))
end


text:
using HetaSimulator, Plots

p = load_platform("./sbml_00375"; source = "index.heta", type = "heta");
m = p.models[:nameless]

# res0 = sim(m; saveat = 0:0.06:3.0, observables = [:S1, :S2, :S3]) # XXX: wrong time sequence here
res0 = sim(m; tspan = (0,3), observables = [:S1, :S2, :S3])
plot(res0)
DataFrame(res0)

max_stars_repo_name: hetalang/qsp-benchmark

max_issues_repo_name: hetalang/qsp-benchmark

max_forks_repo_name: hetalang/qsp-benchmark


text:
import Base: Array, copy

export
    Array,
    AFArray,
    array,
    getBackend,
    empty,
    release!,
    dims,
    host,
    dType,
    jType,
    numdims,
    isEmpty,
    elements,
    copy,
    eval!

type AFArray{B}
    af::ArrayFire{B}
    ptr::Ptr{Void}

    # for empty
    function AFArray(af::ArrayFire{B})
        ptr = af.results.ptr
        dims = [0, 0, 0, 0]
        err = ccall(af.createHandle, Cint, (Ptr{Ptr{Void}}, Cuint, Ptr{DimT}, DType), ptr, 4, dims, f32)
        assertErr(err)
        me = new(af, ptr[])
        register!(af, me) || finalizer(me, release!)
        me
    end

    function AFArray(af::ArrayFire{B}, ptr::Ptr{Void}, wrap = true)
        me = new(af, ptr)
        if wrap
            register!(af, me) || finalizer(me, release!)
        end
        me
    end

    function AFArray{T, N}(af::ArrayFire{B}, arr::Array{T, N})
        ptr = af.results.ptr
        dims = collectDims(af.results.dims, size(arr))
        err = ccall(af.createArray, Cint, (Ptr{Ptr{Void}}, Ptr{T}, Cuint, Ptr{DimT}, DType), ptr, pointer(arr), 4, pointer(dims), asDType(T))
        assertErr(err)
        me = new(af, ptr[])
        register!(af, me) || finalizer(me, release!)
        me
    end

    function AFArray{T}(af::ArrayFire{B}, ::Type{T}, dims::DimT...)
        ptr = af.results.ptr
        dims = collectDims(af.results.dims, dims)
        err = ccall(af.createHandle, Cint, (Ptr{Ptr{Void}}, Cuint, Ptr{DimT}, DType), ptr, 4, pointer(dims), asDType(T))
        assertErr(err)
        me = new(af, ptr[])
        register!(af, me) || finalizer(me, release!)
        me
    end
end

array{B}(af::ArrayFire{B}) = AFArray{B}(af)

array{B, T, N}(af::ArrayFire{B}, arr::Array{T, N}) = AFArray{B}(af, arr)

array{B, T}(af::ArrayFire{B}, ::Type{T}, dims::DimT...) = AFArray{B}(af, T, dims...)

array{B, T}(af::ArrayFire{B}, arr::Array{T}, dims::DimT...) = array(af, reshape(arr, dimsToSize(dims...)...))

array{B}(af::ArrayFire{B}, ptr::Ptr{Void}) = AFArray{B}(af, ptr)

getBackend{B}(arr::AFArray{B}) = (verifyAccess(arr); B)

function release!(arr::AFArray)
    if (arr.ptr != C_NULL)
        release!(arr.af, arr.ptr)
        arr.ptr = C_NULL
    end
end

function release!(af, ptr)
    err = ccall(af.releaseArray, Cint, (Ptr{Void}, ), ptr)
    assertErr(err)
end

function retain!(af, ptr)
    result = af.results.ptr
    err = ccall(af.retainArray, Cint, (Ptr{Ptr{Void}}, Ptr{Void}, ), result, ptr)
    assertErr(err)
    result[]
end

verifyAccess(arr::AFArray) = arr.ptr == C_NULL && error("Cannot access to a released array.")

isEmpty(arr::AFArray) = (verifyAccess(arr); isEmpty(arr.af, arr.ptr))

function isEmpty(af::ArrayFire, ptr::Ptr{Void})
    result = af.results.bool
    err = ccall(
        af.isEmpty,
        Cint, (Ptr{Bool}, Ptr{Void}),
        result, ptr)
    assertErr(err)
    result[]
end

dims(arr::AFArray) = (verifyAccess(arr); dims(arr.af, arr.ptr))

function dims(af::ArrayFire, ptr::Ptr{Void})
    dim0 = af.results.dim0
    dim1 = af.results.dim1
    dim2 = af.results.dim2
    dim3 = af.results.dim3
    err = ccall(
        af.getDims,
        Cint, (Ptr{DimT}, Ptr{DimT}, Ptr{DimT}, Ptr{DimT}, Ptr{Void}),
        dim0, dim1, dim2, dim3, ptr)
    assertErr(err)
    [dim0[], dim1[], dim2[], dim3[]]
end

function dims(arr::AFArray, n)
    verifyAccess(arr)
    af = arr.af
    dim0 = af.results.dim0
    dim1 = af.results.dim1
    dim2 = af.results.dim2
    dim3 = af.results.dim3
    err = ccall(
        af.getDims,
        Cint, (Ptr{DimT}, Ptr{DimT}, Ptr{DimT}, Ptr{DimT}, Ptr{Void}),
        dim0, dim1, dim2, dim3, arr.ptr)
    assertErr(err)
    if n == 0
        dim0[]
    elseif n == 1
        dim1[]
    elseif n == 2
        dim2[]
    else
        dim3[]
    end
end

firstDim(arr::AFArray) = (verifyAccess(arr); dims(arr.af, arr.ptr))

function firstDim(af::ArrayFire, ptr::Ptr{Void})
    result = 1
    for dim in dims(af, ptr)
        if dim > 1
            return result
        end
        result += 1
    end
    1
end

function Base.size(arr::AFArray)
    dimsToSize(dims(arr))
end

function dType(arr::AFArray)
    verifyAccess(arr)
    dType(arr.af, arr.ptr)
end

function dType(af::ArrayFire, ptr::Ptr{Void})
    result = af.results.dType
    err = ccall(
        af.getType,
        Cint, (Ptr{DType}, Ptr{Void}),
        result, ptr)
    assertErr(err)
    result[]
end

jType(arr::AFArray) = asJType(Val{dType(arr)})

numdims(arr::AFArray) = (verifyAccess(arr); numdims(arr.af, arr.ptr))

function numdims(af::ArrayFire, ptr::Ptr{Void})
    result = af.results.dType
    err = ccall(
        af.getNumDims,
        Cint, (Ptr{UInt32}, Ptr{Void}),
        result, ptr)
    assertErr(err)
    result[]
end

elements(arr::AFArray) = (verifyAccess(arr); elements(arr.af, arr.ptr))

function elements(af::ArrayFire, ptr::Ptr{Void})
    result = af.results.dim0
    err = ccall(
        af.getElements,
        Cint, (Ptr{DimT}, Ptr{Void}),
        result, ptr)
    assertErr(err)
    result[]
end

function host(arr::AFArray)
    verifyAccess(arr)
    verifyNotEmpty(arr)
    result = Array{jType(arr)}(size(arr)...)
    _host(arr, result)
end

function host{T, N}(arr::AFArray, to::Array{T, N})
    verifyAccess(arr)
    verifyNotEmpty(arr)
    _host(arr, to)
end

function _host{T, N}(arr::AFArray, to::Array{T, N})
    err = ccall(
        arr.af.getDataPtr,
        Cint, (Ptr{T}, Ptr{Void}),
        to, arr.ptr)
    assertErr(err)
    to
end

eval!(arr::AFArray) = (verifyAccess(arr); eval!(arr.af, arr.ptr))

function eval!(af::ArrayFire, ptr::Ptr{Void})
    err = ccall(af.eval, Cint, (Ptr{Void}, ), ptr)
    assertErr(err)
end

Array(arr::AFArray) = host(arr)

copy(arr::AFArray) = arr[]

verifyNotEmpty(arr::AFArray) = isEmpty(arr) && error("Array is empty.")

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
using Zabbix
using Test

z = Zabbix.ZabbixAPI("https://www.example.com/zabbix/api_jsonrpc.php","username","password")
@test typeof(z) == ZabbixAPI


text:
function derivative_(u::PE_Unit)# Intentially changing name so this is not exported.
    if u.d_ == 0
        result_array = Array{Tuple{Float64,PE_Unit}}(undef,1)
        result_array[1] = (u.b_, PE_Unit(u.b_, u.base_, u.d_))
    else
        result_array = Array{Tuple{Float64,PE_Unit}}(undef,2)
        result_array[1] = (u.b_, PE_Unit(u.b_, u.base_, u.d_))
        result_array[2] = (u.d_, PE_Unit(u.b_, u.base_, u.d_ - 1))
    end
    return result_array
end
function derivative_(u::Array{Tuple{Float64,PE_Unit},1})# Intentially changing name so this is not exported.
    len = length(u)
    heap = Array{Tuple{Float64,PE_Unit},1}()
    for i in 1:len
        u_result = u[i]
        mult = u_result[1]
        output = derivative_(u_result[2])
        output_mults, output_units =  collect(zip(output...))
        output_mults = mult .* vcat(output_mults...)
        output_units = vcat(output_units...)
        zipped_output = collect(zip(output_mults,output_units))
        append!(heap, zipped_output)
    end
    return heap
end

"""
     derivative(f::MultivariateFunction, derivs::Dict{Symbol,Int})
This generates a function representing the derivative of function f. The derivative is
that specified by the derivs dict. So if derivs is Dict{[:x,:y] .=> [1,2]} then there
will be one derivative with respect to x and 2 with respect to y.
"""
function derivative(f::PE_Function, derivs::Dict{Symbol,Int})
    # Should always return a Sum_Of_Functions or a PE_Function.
    dims = keys(derivs)
    fdims = keys(f.units_)

    if (length(setdiff(dims,fdims)) > 0) && minimum(get.(Ref(derivs), [setdiff(dims,fdims)...],0)) < 1
        return PE_Function()
    end
    dims_dict = Dict{Symbol,Array{Tuple{Float64,PE_Unit},1}}()
    units = deepcopy(f.units_)
    for dimen in dims
        num_derivs = derivs[dimen]
        if num_derivs > 0
            if dimen in keys(units)
                unit = pop!(units, dimen)
                der = derivative_(unit)
                for i in 2:num_derivs
                    der = derivative_(der)
                end
                dims_dict[dimen] = der
            else
                dims_dict[dimen] = Array{Tuple{Float64,PE_Unit},1}([(0.0, PE_Unit())])
            end
        end
    end
    array_of_tups = [Dict(dims .=> val) for val in (collect(Iterators.product(getindex.((dims_dict,),dims)...))...,)]
    array_of_pes = PE_Function.(1.0, array_of_tups)
    remaining_dims = PE_Function(f.multiplier_, units)
    final_result = remaining_dims .* array_of_pes
    if length(final_result) == 0
        return PE_Function()
    elseif length(final_result) == 1
        return final_result[1]
    else
        return Sum_Of_Functions(final_result)
    end
end
function derivative(f::Sum_Of_Functions, derivs::Dict{Symbol,Int})
    # Should always return a Sum_Of_Functions or a PE_Function.
    deriv_funcs = derivative.(f.functions_, Ref(derivs))
    return Sum_Of_Functions(deriv_funcs)
end

function derivative(f::Piecewise_Function, derivs::Dict{Symbol,Int})
    max_order = maximum(values(derivs))
    if max_order > 0
        derivatives = derivative.(f.functions_, Ref(derivs))
        return Piecewise_Function(derivatives, f.thresholds_)
    elseif max_order == 0
        return f
    else
        error("Not possible to take a negative derivative. Try evaluating an integral instead.")
    end
end

function derivative(f::Sum_Of_Piecewise_Functions, derivs::Dict{Symbol,Int})
    return derivative(f.global_funcs_, derivs) + sum(derivative.(f.functions_, Ref(derivs)))
end

function derivative(f::Missing, derivs::Dict{Symbol,Int})
    return Missing()
end

function derivative(f::MultivariateFunction)
    if length(setdiff(underlying_dimensions(f), Set([default_symbol]))) == 0
        derivs = Dict{Symbol,Int}(default_symbol => 1)
        return derivative(f, derivs)
    else
        error("It is not possible to take the derivative without inputting a dict specifying the desired derivative unless the only variable is the default one.")
    end
end

function add_to_dict(dic::Dict{Symbol,Int}, new_symbol::Symbol)
    dd = deepcopy(dic)
    if new_symbol in keys(dd)
        dd[new_symbol] = 1 + dd[new_symbol]
        return dd
    else
        dd[new_symbol] = 1
        return dd
    end
end

"""
     all_derivatives(f::MultivariateFunction, degree::Int = 2, dimensions::Set{Symbol} = underlying_dimensions(f))
This generates a dict containing functions representing all of the deriviates of a function up to the order of degree.
"""
function all_derivatives(f::MultivariateFunction, degree::Int = 2, dimensions::Set{Symbol} = underlying_dimensions(f))
    if typeof(f) == MultivariateFunctions.Piecewise_Function
        derivs = Dict{Dict{Symbol,Int},typeof(f)}()
    else
        derivs = Dict{Dict{Symbol,Int},Sum_Of_Functions}()
    end
    derivs[Dict{Symbol,Int}()] =  f
    previous_dicts = Set(collect(keys(derivs)))
    for deg in 1:degree
        for prev in previous_dicts
            for dimension in dimensions
                proposed_dict = add_to_dict(prev,dimension)
                if !(proposed_dict in keys(derivs))
                    derivs[proposed_dict] = derivative(derivs[prev], Dict{Symbol,Int}(dimension => 1))
                end
            end
        end
        previous_dicts = setdiff(keys(derivs), Set(previous_dicts))
    end
    return derivs
end

"""
    jacobian(derivs::Union{Dict{Dict{Symbol,Int},PE_Function},Dict{Dict{Symbol,Int},Sum_Of_Functions},Dict{Dict{Symbol,Int},Piecewise_Function}}, labels::Array{Symbol,1})
    jacobian(f::MultivariateFunction, dimensions::Array{Symbol,1})
This generates an array of MultivariateFunctions representing the derivatives of a function.
This array can be evaluated with evaluate.(jacobian, Ref(coordinates))  to give a vector of the derivative values at a point.
"""
function jacobian(derivs::Union{Dict{Dict{Symbol,Int},PE_Function},Dict{Dict{Symbol,Int},Sum_Of_Functions},Dict{Dict{Symbol,Int},Piecewise_Function}}, labels::Array{Symbol,1})
    ders = Array{MultivariateFunction,1}(undef,length(labels))
    for i in 1:length(labels)
        lookup_dict = Dict{Symbol,Int}(labels[i] => 1)
        ders[i] = derivs[lookup_dict]
    end
    return ders
end

function jacobian(f::MultivariateFunction, dimensions::Array{Symbol,1})
    derivs = all_derivatives(f, 1, Set(dimensions))
    return jacobian(derivs, dimensions)
end

"""
    Hessian(f::MultivariateFunction, dimensions::Array{Symbol,1})
    Hessian(derivs::Union{Dict{Dict{Symbol,Int},PE_Function},Dict{Dict{Symbol,Int},Sum_Of_Functions},Dict{Dict{Symbol,Int},Piecewise_Function}}, labels_::Array{Symbol,1})
This represents expressions for constructing a hessian matrix for a function. It can be evaluated
to get a Hermitian matrix of the hessian at a particular location.
"""
struct Hessian
    derivs_::Union{Dict{Dict{Symbol,Int},PE_Function},Dict{Dict{Symbol,Int},Sum_Of_Functions},Dict{Dict{Symbol,Int},Piecewise_Function}}
    labels_::Array{Symbol,1}
    function Hessian(derivs::Union{Dict{Dict{Symbol,Int},PE_Function},Dict{Dict{Symbol,Int},Sum_Of_Functions},Dict{Dict{Symbol,Int},Piecewise_Function}}, labels_::Array{Symbol,1})
        new(derivs, labels_)
    end
    function Hessian(f::MultivariateFunction, dimensions::Array{Symbol,1})
        derivs = all_derivatives(f, 2, Set(dimensions))
        return Hessian(derivs, dimensions)
    end
end

"""
    evaluate(hess::Hessian, coordinates::Dict{Symbol,Float64})
This evaluates a Hessian object to create a Hermitian matrix representing the hessian
at a point.
"""
function evaluate(hess::Hessian, coordinates::Dict{Symbol,Float64})
    len = length(hess.labels_)
    second_derivs = Array{Float64,2}(undef,len,len)
    for c in 1:len
        for r in c:len
            lookup_dict = Dict{Symbol,Int}([hess.labels_[r],hess.labels_[c]] .=> [1,1])
            func = hess.derivs_[lookup_dict]
            val = evaluate(func, coordinates)
            second_derivs[r,c] = val
        end
    end
    return Hermitian(second_derivs, :L)
end

"""
    find_local_optima(func::MultivariateFunction, initial_guess::Dict{Symbol,Float64}; step_size::Float64 = 1.0, max_iters::Int = 40, convergence_tol::Float64 = 1e-10, print_reports::Bool = false)
This takes the analytical jacobian and hessian of a function and uses them to find a nearby
optima. The optima it will find are based on Newton's method. There is no way to specify whether
a minimum or a maximum is sought in Newton's method (at least the pure version of it) and thus
this function cannot selectively search for a maximum or minimum. It simply searches for a stationary point.
"""
function find_local_optima(func::MultivariateFunction, initial_guess::Dict{Symbol,Float64}; step_size::Float64 = 1.0, max_iters::Int = 40, convergence_tol::Float64 = 1e-10, print_reports::Bool = false)
    dimensions = collect(keys(initial_guess))
    dimensions_without_guess = setdiff(underlying_dimensions(func), Set(dimensions))
    if length(dimensions_without_guess) > 0
        error(string("A guess value must be input for all dimensions in the function for which a root is sought. But ", collect(dimensions_without_guess), " has/have no guess values."))
    end
    derivs = all_derivatives(func, 2, Set(dimensions))
    jacob = jacobian(derivs, dimensions)
    hess  = Hessian(derivs, dimensions)
    guess = initial_guess
    iter = 1
    while iter <= max_iters
        hessn = evaluate(hess, guess)
        inverse_hessian = inv(hessn)
        evaluated_jacobian = evaluate.(jacob, Ref(guess))
        if all(abs.(evaluated_jacobian) .< convergence_tol)
            val = evaluate(func,guess)
            if print_reports
                println(string("Converged. The converged value of f is ", val, "."))
            end
            hessian_det_sign = det(hessn) > 0
            return NamedTuple{(:coordinates, :value, :convergence, :positive_definite_hessian)}((guess, val, true, hessian_det_sign))
        end
        step_ = inverse_hessian * evaluated_jacobian
        x = get.(Ref(guess), dimensions, 0)
        x_prime = x - step_size * step_
        guess = Dict{Symbol,Float64}(dimensions .=> x_prime)
        iter = iter + 1
    end
    val = evaluate(func,guess)
    if print_reports
        println(string("Did not converge. The value of the function is ", val, "."))
    end
    return NamedTuple{(:coordinates, :value, :convergence)}((guess, val, false))
#    x - step_size * (hessian^{-1}) * jacobian.
end

"""
    uniroot(f::MultivariateFunction, initial_guess::Dict{Symbol,Float64}; step_size::Float64 = 1.0, max_iters::Int = 40, convergence_tol::Float64 = 1e-10, print_reports::Bool = false)
This takes the analytical jacobian and hessian of a function and uses them to find a nearby
root. It finds a root using Newton's method.
"""
function uniroot(f::MultivariateFunction, initial_guess::Dict{Symbol,Float64}; step_size::Float64 = 1.0, max_iters::Int = 40, convergence_tol::Float64 = 1e-10, print_reports::Bool = false)
    dimensions = collect(keys(initial_guess))
    dimensions_without_guess = setdiff(underlying_dimensions(f), Set(dimensions))
    if length(dimensions_without_guess) > 0
        error(string("A guess value must be input for all dimensions in the function for which a root is sought. But ", collect(dimensions_without_guess), " has/have no guess values."))
    end
    derivs = all_derivatives(f, 1, Set(dimensions))
    jacob = jacobian(derivs, dimensions)
    guess = initial_guess
    iter = 1
    while iter <= max_iters
        val = evaluate(f,guess)
        if abs(val) < convergence_tol
            if print_reports
                println(string("Converged. The converged value of f is ", val, "."))
            end
            return NamedTuple{(:coordinates, :value, :convergence)}((guess,val, true))
        end
        step_ = val ./ evaluate.(jacob, Ref(guess))
        if print_reports
            println(string("Iterate: ",iter,"  The value of f is ", val, " and the next step's Euclidean distance is (excluding step_size multiplier) ", sqrt(sum(step_ .^ 2))))
        end
        x = get.(Ref(guess), dimensions, 0)
        x_prime = x - step_size * step_
        guess = Dict{Symbol,Float64}(dimensions .=> x_prime)
        iter = iter + 1
    end
    val = evaluate(f,guess)
    if print_reports
        println(string("Did not converge. The value of f is ", val, "."))
    end
    return NamedTuple{(:coordinates, :value, :convergence)}((guess, val, false))
end

## Integration
function indefinite_integral(u::PE_Unit, incoming_multiplier::Float64 = 1.0)# Intentially changing name so this is not exported.
    result_array = Array{Tuple{Float64,PE_Unit}}(undef,1)
    if u.b_ ≂ 0.0
        result_array[1] = (incoming_multiplier/(u.d_+1), PE_Unit(u.b_, u.base_, u.d_+1)) # Note (u.d_+1) > 0 because u.d_ \geq 0
        return result_array
    else
        result_array[1] = (incoming_multiplier/u.b_, PE_Unit(u.b_, u.base_, u.d_))
        if u.d_ > 0
            other_pieces_multiplier = -incoming_multiplier * (u.d_ / u.b_)
            other_pieces_unit = PE_Unit(u.b_, u.base_, u.d_ - 1)
            other_pieces = indefinite_integral(other_pieces_unit, other_pieces_multiplier)
            append!(result_array, other_pieces)
        end
        return result_array
    end
end

function apply_limits(mult::Float64, indef::Array{Tuple{Float64,PE_Unit}}, left::Symbol, right::Symbol)# Intentially changing name so this is not exported.
    converted = (collect(Iterators.product(((indef,))...))...,)
    rights = [Dict{Symbol,Tuple{Float64,PE_Unit}}([right] .=> val) for val in converted]
    lefts  = [Dict{Symbol,Tuple{Float64,PE_Unit}}([left] .=> val) for val in converted]
    return Sum_Of_Functions(PE_Function.(mult, rights)) - Sum_Of_Functions(PE_Function.(mult, lefts))
end
function apply_limits(mult::Float64, indef::Array{Tuple{Float64,PE_Unit}}, left::Float64, right::Symbol)# Intentially changing name so this is not exported.
    converted = (collect(Iterators.product(((indef,))...))...,)
    rights = [Dict{Symbol,Tuple{Float64,PE_Unit}}([right] .=> val) for val in converted]
    lefts  = [Dict{Symbol,Tuple{Float64,PE_Unit}}([:left] .=> val) for val in converted]
    return Sum_Of_Functions(PE_Function.(mult, rights)) - sum(evaluate.(PE_Function.(mult, lefts), Ref(Dict{Symbol,Float64}(:left => left))))
end
function apply_limits(mult::Float64, indef::Array{Tuple{Float64,PE_Unit}}, left::Symbol, right::Float64)# Intentially changing name so this is not exported.
    converted = (collect(Iterators.product(((indef,))...))...,)
    rights = [Dict{Symbol,Tuple{Float64,PE_Unit}}([:right] .=> val) for val in converted]
    lefts  = [Dict{Symbol,Tuple{Float64,PE_Unit}}([left] .=> val) for val in converted]
    return evaluate.(Sum_Of_Functions(PE_Function.(mult, rights)), Ref(Dict{Symbol,Float64}(:right => right))   ) - Sum_Of_Functions(PE_Function.(mult, lefts))
end
function apply_limits(mult::Float64, indef::Array{Tuple{Float64,PE_Unit}}, left::Float64, right::Float64)# Intentially changing name so this is not exported.
    converted = (collect(Iterators.product(((indef,))...))...,)
    converted_again = [Dict{Symbol,Tuple{Float64,PE_Unit}}([default_symbol] .=> val) for val in converted]
    funcs = Sum_Of_Functions(PE_Function.(mult, converted_again))
    return evaluate.(funcs, right) - evaluate.(funcs, left)
end

const IntegrationLimitDict = Union{Dict{Symbol,Tuple{Union{Symbol,Float64},Union{Symbol,Float64}}},Dict{Symbol,Tuple{Symbol,Symbol}},Dict{Symbol,Tuple{Float64,Float64}},Dict{Symbol,Tuple{Symbol,Float64}},Dict{Symbol,Tuple{Float64,Symbol}}}

"""
     integral(f::PE_Function, limits::Dict{Symbol,Tuple{Union{Symbol,Float64},Union{Symbol,Float64}}})
     integral(f::Sum_Of_Functions, limits::Dict{Symbol,Tuple{Union{Symbol,Float64},Union{Symbol,Float64}}})
This gives a function representing the integral of a function, f, with limits in each dimension given by a dict. The dict should
contain a tuple for each variable in the function. The left member of the tuple contains the lower limit
and the right member the upper limite. Each can be a Float64 or a symbol. If a symbol is input then
this will get incorporated as a dimension in the MultivariateFunction created by the integral function.
"""
function integral(f::PE_Function, limits::IntegrationLimitDict)
    if length(f.units_) == 0
        volume_of_cube = 1.0
        for dimen in keys(limits)
            volume_of_cube = volume_of_cube * (limits[dimen][2] - limits[dimen][1])
        end
        return volume_of_cube * f.multiplier_
    end
    units = deepcopy(f.units_)
    result_by_dimension = Array{Union{Float64,Sum_Of_Functions,PE_Function},1}()
    for dimen in keys(limits)
        left_  = limits[dimen][1]
        right_ = limits[dimen][2]
        f_unit = Array{PE_Unit,1}(undef,1)
        if haskey(f.units_, dimen)
            f_unit = pop!(units, dimen)
        else
            f_unit = PE_Unit()
        end
        indef  = indefinite_integral(f_unit)
        result_of_applying_limits = apply_limits(1.0, indef, left_, right_)
        append!(result_by_dimension, [result_of_applying_limits])
    end
    if length(units) == 0
        return f.multiplier_ * prod(result_by_dimension)
    else
        return PE_Function(f.multiplier_, units) * prod(result_by_dimension)
    end
end

function integral(f::Sum_Of_Functions, limits::IntegrationLimitDict)
    if length(f.functions_) == 0
        return 0.0
    else
        funcs = integral.(f.functions_, Ref(limits))
        return sum(funcs)
    end
end

"""
     integral(f::Piecewise_Function, limits::Dict{Symbol,Tuple{Float64,Float64}})
     integral(f::Sum_Of_Piecewise_Functions, limits::Dict{Symbol,Tuple{Float64,Float64}})
This gives a function representing the integral of a function, f, with limits in each dimension given by a dict. The dict should
contain a tuple for each variable in the function. The left member of the tuple contains the lower limit
and the right member the upper limite. Each must be a Float64 (Support for inputting a symbol is planned but not yet implemented).
"""
function hypercubes_to_integrate(f::Piecewise_Function, lims::Dict{Symbol,Tuple{Float64,Float64}})
    limits = Dict{Symbol,Tuple{Float64,Float64}}()
    for kk in keys(f.thresholds_)
        limits[kk] = lims[kk]
    end
    ks= sort(collect(keys(limits)))
    new_dict = Dict{Symbol,Array{Float64}}()
    for dimen in ks
        lower = limits[dimen][1]
        upper = limits[dimen][2]
        if lower + eps() >= upper
            error(string("The lower limit of integration is higher than the upper for dimension ", dimen))
        end
        thresholds = f.thresholds_[dimen]
        censored_thresholds = vcat(lower, thresholds[thresholds .> lower])
        censored_thresholds2 = vcat(censored_thresholds[censored_thresholds .< upper] , upper)
        new_dict[dimen] = censored_thresholds2
    end
    lengths_minus_one = length.(get.(Ref(new_dict),ks,0)) .- 1
    index_combinations = vcat.(collect(collect(Iterators.product(range.(1,lengths_minus_one; step = 1)...))))
    indices = (collect(collect(Iterators.product(range.(1,lengths_minus_one; step = 1)...)))...,)
    indices = collect(collect.(indices))
    total_len = length(indices)
    new_cubes = Array{Dict{Symbol,Tuple{Float64,Float64}},1}(undef,total_len)
    for i in range(1, total_len; step = 1)
        new_cube = Dict{Symbol,Tuple{Float64,Float64}}()
        ind = indices[i]
        for j in 1:length(ks)
            key = ks[j]
            indd = ind[j]
            new_cube[key]  =  (new_dict[key][indd], new_dict[key][indd+1])
        end
        new_cubes[i] = new_cube
    end
    return new_cubes
end

function integral(f::Piecewise_Function, limits::Dict{Symbol,Tuple{Float64,Float64}})
    cubes_to_integrate = hypercubes_to_integrate(f, limits)
    ints = Array{Any,1}(undef,length(cubes_to_integrate))
    for i in 1:length(ints)
        piece_func = get_correct_function_from_piecewise(f, cubes_to_integrate[i])
        ints[i] = integral(piece_func, cubes_to_integrate[i])
    end
    return sum(ints)
end

function integral(f::Sum_Of_Piecewise_Functions, limits::Dict{Symbol,Tuple{Float64,Float64}})
    return integral(f.global_funcs_, limits) + sum(integral.(f.functions_, Ref(limits)))
end

function integral(f::MultivariateFunction, left_limit::Float64, right_limit::Float64)
    underlyingDims = underlying_dimensions(f)
    if length(underlyingDims) == 0
        return f.multiplier_ * (right_limit-left_limit)
    elseif underlyingDims == Set([default_symbol])
        limits_ = Dict{Symbol,Tuple{Float64,Float64}}(default_symbol => Tuple{Float64,Float64}((left_limit,right_limit)))
        return integral(f, limits_)
    else
        error("Cannot evaluate the integral of a Multivariate function without a dictionary set of coordinates unless it is a MultivariateFunction with only the default dimension being used.")
    end
end

function integral(f::MultivariateFunction, left_limit::Date, right_limit::Date)
    if underlying_dimensions(f) == Set([default_symbol])
        limits_ = Dict{Symbol,Tuple{Float64,Float64}}(default_symbol => Tuple{Float64,Float64}((years_from_global_base(left_limit),years_from_global_base(right_limit))))
        return integral(f, limits_)
    else
        error("Cannot evaluate the integral of a Multivariate function without a dictionary set of coordinates unless it is a MultivariateFunction with only the default dimension being used.")
    end
end


text:
 export empError,
        looError

# Normalized empirical error / generalization error (Blatman2009 - Dissertation)
# nemError ∈ [0,1]
function empError(Y::Vector{Float64}, Φ::Matrix{Float64}, pceCoeffs; adjusted = false)
    @assert length(Y) > 0           "Empty results vector Y"
    @assert length(Y) == size(Φ, 1) "Inconsistent number of elements"
    # TODO: Require column vectors

    # Compute PCE model response
    Y_Pce = Φ * pceCoeffs

    # Empirical variance of sampling evaluation (true model)
    n = length(Y)
    varY = n > 1 ? 1/(n-1) * sum( (Y[i] - mean(Y))^2 for i in 1:n ) : 0 # FIX: In case of one single sample this always results in a 100% fit

    # Empirical error
    empError = 1/n * sum( (Y[i] - Y_Pce[i])^2 for i in 1:n )

    if adjusted
        # TODO: adjusted empirical error (penalize larger PCE bases)
        # p = ... #bincoeff
        # n = ...
        # empError = (n-1) / (n - p - 1) * ( empError)
    end

    # Normalize with variance FIX: s.o.
    nempError = varY == 0 ? 0 : empError / varY

    return nempError
end

# Leave-one-out cross-validation error (Blatman2009 - Dissertation)
# Note: looError ∈ [0,∞]. If > 1, average modeling would be better.
function looError(Y::Vector{Float64}, Φ::Matrix{Float64}, pceCoeffs)
    @assert length(Y) > 0           "Empty results vector Y"
    @assert length(Y) == size(Φ, 1) "Inconsistent number of elements"

    # Compute PCE model response
    Y_Pce = Φ * pceCoeffs

    # h-factor for validation sets
    M = Φ * inv(Φ' * Φ) * Φ' #TODO: fix conditioning
    N = size(M, 1)
    h = ones(N) - diag(M)

    # Empirical variance of sampling evaluation (true model)
    n = length(Y)
    varY = n > 1 ? 1/(n-1) * sum( (Y[i] - mean(Y))^2 for i in 1:n ) : 0

    # compute squared error with h-factor
    loo = 1/N * sum( ( (Y[i] - Y_Pce[i]) / h[i] )^2 for i in 1:N )

    ## Debuging
    # println("M: ", M)
    # println("Y - Y_Pce: ", Y - Y_Pce)
    # println("h: ", h)
    # println("varY: ", varY)
    # println("looError: ", loo)

    # Normalize error with variance. Set error to 0, if var is 0
    looError = varY == 0 ? 0 : loo / varY

    # TODO: adjusted Error computation

    # Return loo error. Return Inf, if it is NaN. Can happen for underdetermined problems.
    return isnan(looError) ? Inf : looError
end

max_stars_repo_name: adriangrupp/sparsePolyChaos

max_issues_repo_name: adriangrupp/sparsePolyChaos

max_forks_repo_name: adriangrupp/sparsePolyChaos