text:
################################################################################
# Interface
# A Measurement must inherit from AbstractMeasurement and implement the
# following functions.

abstract type AbstractMeasurement end

###################################
# You must implement
###################################

"""
    measure!(measurement, mc, model, sweep_index)

Performs a measurement during the measurement phase.
"""
function measure!(m::AbstractMeasurement, mc, model, i)
    throw(MethodError(measure!, (m, mc, model)))
end

###################################
# You may implement
###################################

"""
    prepare!(measurement, mc, model)

This method is called between the thermalization and the measurement phase.
"""
prepare!(m::AbstractMeasurement, mc, model) = nothing

"""
    finish!(measurement, mc, model)

Finish a measurement. This method is called after the measurement phase.
"""
finish!(m::AbstractMeasurement, mc, model) = nothing

"""
    observables(measurement)

Returns a dictionary of observables.

The default implementation searches for all fields `<: AbstractObservable` and
builds a dictionary pairs `name(obs) => obs`.
"""
function observables(m::AbstractMeasurement)
    os = [getfield(m, n) for n in obs_fieldnames_from_obj(m)]
    Dict{String, AbstractObservable}(MonteCarloObservable.name(o) => o for o in os)
end

"""
    _save(filename, measurement, entryname)

Saves a measurement to a jld-file `filename` in group `entryname`.

The default implementation saves the full measurement object. This requires the
default constructor to be available, i.e. construction by passing fields. The
group `entryname` follows the structure given by `measurements()`, e.g.
`ME/config` for a `ConfigurationMeasurement` with tag `:config` in the
measurement phase.

See also [`save_measurements`](@ref), [`measurements`](@ref), [`_load`](@ref)
"""
function _save(file::JLDFile, m::AbstractMeasurement, entryname::String)
    # NOTE: `VERSION` and `type` are necessary
    write(file, entryname * "/VERSION", 0)
    write(file, entryname * "/tag", "Generic")
    write(file, entryname * "/data", m)
    nothing
end

# Statistics forwarded from MonteCarloObservable/BinningAnalysis
# Generates functions
#   mean(measurement)       - returns the mean/expectation value of a measurement
#   var(measurement)        - returns the variance of a measurement
#   std_error(measurement)  - returns the standard error of the mean of a measurement
#   tau(measurement)        - return the autocorrelation time of a measurement
for (func, name) in zip(
        (:mean, :var, :std_error, :tau),
        ("mean", "variance", "standard error", "autocorrelation time")
    )
    docstring = """
        $func(measurement)

    Returns the $name of a given `measurement`.

    The default implementation searches for all fields `<: AbstractObservable`
    and returns `$func(x)` for each field `x`. If there are multiple fields
    with type `<: AbstractObservable` a dictionary will be return.
    """
    @eval begin
        @doc $docstring $func
        function MonteCarloObservable.$(func)(m::AbstractMeasurement)
            fn = obs_fieldnames_from_obj(m)
            os = [getfield(m, n) for n in fn]
            if isempty(os)
                throw(error("Did not find any observables in $m."))
            elseif length(os) == 1
                return $(func)(os[1])
            else
                return Dict(MonteCarloObservable.name(o) => $(func)(o) for o in os)
            end
        end
    end
end

for func in (:length, :isempty, :empty!)
    @eval begin
        function Base.$(func)(m::AbstractMeasurement)
            fn = obs_fieldnames_from_obj(m)
            os = [getfield(m, n) for n in fn]
            if isempty(os)
                throw(error("Did not find any observables in $m."))
            elseif length(os) == 1
                return $(func)(os[1])
            else
                return Dict(MonteCarloObservable.name(o) => $(func)(o) for o in os)
            end
        end
    end
end

################################################################################
# A new model may implement the following for convenience

"""
    default_measurements(mc, model)

Return a dictionary of default measurements for a given Monte Carlo flavour and
model. If there is no implementation given for the specific Monte Carlo flavour
an empty dictionary will be returned.
"""
default_measurements(mc, model) = Dict{Symbol, AbstractMeasurement}()

################################################################################
# mc based, default measurements

"""
    ConfigurationMeasurement(mc, model[, rate=1])

Measures configurations of the given Monte Carlo flavour and model. The rate of
measurements can be reduced with `rate`. (e.g. `rate=10` means 1 measurement per
10 sweeps)
"""
struct ConfigurationMeasurement{OT <: Observable} <: AbstractMeasurement
    obs::OT
    rate::Int64
    function ConfigurationMeasurement{OT}(obs, rate) where {OT <: Observable}
        @warn(
            "`ConfigurationMeasurement` is deprecated. Configurations can now " *
            "be collected in `mc.configs` using an `<: AbstractConfiguartionAccumulator`"
        )
        new{OT}(obs, rate)
    end
end
function ConfigurationMeasurement(mc::MonteCarloFlavor, model::Model, rate=1)
    o = Observable(typeof(mc.conf), "Configurations")
    ConfigurationMeasurement{typeof(o)}(o, rate)
end
@bm function measure!(m::ConfigurationMeasurement, mc, model, i::Int64)
    (i % m.rate == 0) && push!(m.obs, conf(mc))
    nothing
end

################################################################################
# called by simulation

for function_name in (:prepare!, :finish!)
    @eval begin
        function $(function_name)(
                measurements::Dict{Symbol, AbstractMeasurement},
                mc, model
            )
            for (k, m) in measurements
                $(function_name)(m, mc, model)
            end
            nothing
        end
    end
end

function measure!(
        measurements::Dict{Symbol, AbstractMeasurement},
        mc, model, sweep_index::Int64
    )
    for (k, m) in measurements
        measure!(m, mc, model, sweep_index)
    end
    nothing
end

################################################################################
# other convenience functions

# Searches `obj` for fields which <: AbstractObservables
# if get_all_names = false, return array of fieldnames <: AbstractObservable
# if get_all_names = true, also return other fieldnames
function obs_fieldnames_from_obj(obj, get_all_names=false)
    ObsTypes = Union{AbstractObservable, LogBinner}
    fnames = fieldnames(typeof(obj))
    obs_names = [
        s for s in fnames if getfield(obj, s) isa ObsTypes
    ]
    if get_all_names
        other_names = [
            s for s in fnames if !(getfield(obj, s) isa ObsTypes)
        ]
        return obs_names, other_names
    else
        return obs_names
    end
end

# printing
# function Base.show(io::IO, m::AbstractMeasurement)
#     #  no parametrization -v    v- no MonteCarlo.
#     typename = typeof(m).name.name
#     observables, other = obs_fieldnames_from_obj(m, true)
#     println(io, typename)
#     for obs_fieldname in observables
#         o = getfield(m, obs_fieldname)
#         oname = MonteCarloObservable.name(o)
#         otypename = typeof(o).name.name
#         println(io, "\t", obs_fieldname, "::", otypename, "\t → \"", oname, "\"")
#     end
#     for fieldname in other
#         println(io, "\t", fieldname, "::", typeof(getfield(m, fieldname)))
#     end
#     nothing
# end

function Base.show(io::IO, ::MIME"text/plain", m::AbstractMeasurement)
    #small
    #  no parametrization -v    v- no MonteCarlo.
    typename = typeof(m).name.name
    temp = obs_fieldnames_from_obj(m)
    observable_names = map(temp) do obs_fieldname
        MonteCarloObservable.name(getfield(m, obs_fieldname))
    end
    print(io, typename, "(\"", join(observable_names, "\", \""), "\")")
    nothing
end

"""
    measurements(mc, stage = :ME)

Returns a nested dictionary of measurements used in a given Monte Carlo
simulation `mc`.
By default, only measurements from the measurement stage (`:ME`) are returned.
To get measurements from the thermalization stage, use `stage = :TH`. To get all
measurements, use `stage = :all`.

```julia
julia> julia> measurements(dqmc)
Dict{Symbol,MonteCarlo.AbstractMeasurement} with 3 entries:
  :conf              => ConfigurationMeasurement…
  :Greens            => GreensMeasurement…
  :BosonEnergy       => BosonEnergyMeasurement…
```
"""
function measurements(mc::MonteCarloFlavor, stage = :ME)
    measurement_stage = (:ME, :me, :Measurement, :measurement)
    thermalization_stage = (:TH, :th, :Thermalization, :thermalization, :thermalization_measurements)
    all_stages = (:all, :ALL)

    if stage in measurement_stage
        return mc.measurements
    elseif stage in thermalization_stage
        return mc.thermalization_measurements
    elseif stage in all_stages
        return Dict(
            :TH => mc.thermalization_measurements,
            :ME => mc.measurements
        )
    else
        throw(error(
            "The given stage `$stage` was not recognized. Try one of (" *
            join(string.([
                measurement_stage...,
                thermalization_stage...,
                all_stages...
            ]), ", ") * ")"
        ))
    end
end

"""
    observables(mc::MonteCarloFlavor, stage = :ME)

Returns a nested dictionary of all observables used during a given `stage` of
the given Monte Carlo simulation `mc`. By default the resulting dictionary will
contain measurements from the measurement (`:ME`) stage.

The thermalization stage can be accessed with `stage = :TH` and both stages can
be retrieved together using `stage = :all`.

The dictionary is generated by calling `observables(measurement)` for each
`measurement` in the given `stage` of the Monte Carlo simulation `mc`. By
default this will generate a nested dictionary indexed as
`observables(mc)[measurement_name][observable_name]`.

```julia
julia> observables(dqmc)
Dict{Symbol,Dict{String,MonteCarloObservable.AbstractObservable}} with 3 entries:
  :BosonEnergy       => Dict{String,MonteCarloObservable.AbstractObservable}("Bosonic Energy"=>LightObservable{Float64,20}())
  :conf              => Dict{String,MonteCarloObservable.AbstractObservable}("Configurations"=>Array{Int8,2} Observable…
  :Greens            => Dict{String,MonteCarloObservable.AbstractObservable}("Equal-times Green's function"=>LightObservable{Array{Complex{Float64},2},20}())
```
"""
function observables(mc::MonteCarloFlavor, stage = :ME)
    measurement_stage = (:ME, :me, :Measurement, :measurement)
    thermalization_stage = (:TH, :th, :Thermalization, :thermalization, :thermalization_measurements)
    all_stages = (:all, :ALL)

    if stage in measurement_stage
        me_obs = Dict{Symbol, Dict{String, AbstractObservable}}(
            k => observables(m) for (k, m) in mc.measurements
        )
        return me_obs

    elseif stage in thermalization_stage
        th_obs = Dict{Symbol, Dict{String, AbstractObservable}}(
            k => observables(m) for (k, m) in mc.thermalization_measurements
        )
        return th_obs

    elseif stage in all_stages
        th_obs = Dict{Symbol, Dict{String, AbstractObservable}}(
            k => observables(m) for (k, m) in mc.thermalization_measurements
        )
        me_obs = Dict{Symbol, Dict{String, AbstractObservable}}(
            k => observables(m) for (k, m) in mc.measurements
        )
        return Dict(:TH => th_obs, :ME => me_obs)

    else
        throw(error(
            "The given stage `$stage` was not recognized. Try one of (" *
            join(string.([
                measurement_stage...,
                thermalization_stage...,
                all_stages...
            ]), ", ") * ")"
        ))
    end
end

# For `mc[:M]` instead of `measurements(mc)[:M]`. How convenient!
Base.getindex(mc::MonteCarloFlavor, k) = getindex(measurements(mc), k)
# Allow `mc[:M] = MagnetizationMeasurement(mc, model)` to add measurements
Base.setindex!(mc::MonteCarloFlavor, v, k) = setindex!(measurements(mc), v, k)
# Allow `keys(dqmc)` to get measurement keys
Base.keys(mc::MonteCarloFlavor) = keys(measurements(mc))
Base.haskey(mc::MonteCarloFlavor, key) = haskey(measurements(mc), key)

"""
    push!(mc, tag::Symbol => MT::Type{<:AbstractMeasurement}[, stage=:ME; paassthrough...])

Adds a new pair `tag => MT(mc, model; passthorugh...)`, where `MT` is a type
`<: AbstractMeasurement`, to either the thermalization or measurement `stage`
(`:TH` or `:ME`) of the simulation `mc`.

Examples:
```
push!(mc, :conf => ConfigurationMeasurement)
push!(mc, :conf => ConfigurationMeasurement, stage=:ME)
```

See also: [`unsafe_push!`](@ref)
"""
function Base.push!(mc::MonteCarloFlavor, p::Pair{Symbol, T}, stage=:ME; passthrough...) where T
    tag, MT = p
    p[2] <: AbstractMeasurement || throw(ErrorException(
        "The given `tag => MT` pair should be of type " *
        "`Pair{Symbol, Type(<: AbstractMeasurement)}`, but is " *
        "`Pair{Symbol, Type{$(p[2])}}`."
    ))
    unsafe_push!(mc, tag => MT(mc, mc.model; passthrough...), stage)
end

"""
    unsafe_push!(mc, tag::Symbol => m::AbstractMeasurement[, stage=:ME])

Adds a pair `tag => m` to either the thermalization or measurement stage (`:TH`
or `:ME`) of the given simulation `mc`.

Note that this function is unsafe as it does not test whether `m` is a valid
measurement for the given simulation.

Examples:
```
push!(mc, :conf => ConfigurationMeasurement(mc, model))
push!(mc, :conf => ConfigurationMeasurement(mc, model), stage=:ME)
```

See also: [`MonteCarlo.push!`](@ref)
"""
function unsafe_push!(mc::MonteCarloFlavor, p::Pair{Symbol, <:AbstractMeasurement}, stage=:ME)
    measurement_stage = (:ME, :me, :Measurement, :measurement)
    thermalization_stage = (:TH, :th, :Thermalization, :thermalization, :thermalization_measurements)
    all_stages = (:all, :ALL)

    if stage in measurement_stage
        push!(mc.measurements, p)
    elseif stage in thermalization_stage
        push!(mc.thermalization_measurements, p)
    elseif stage in all_stages
        push!(mc.thermalization_measurements, p)
        push!(mc.measurements, p)
    else
        throw(error(
            "The given stage `$stage` was not recognized. Try one of (" *
            join(string.([
                measurement_stage...,
                thermalization_stage...,
                all_stages...
            ]), ", ") * ")"
        ))
    end
end

"""
    delete!(mc, key[, stage=:MC])
    delete!(mc, MT::Type{<:AbstractMeasurement}[, stage=:MC])

Deletes a measurement from the given Monte Carlo simulation by key or by type.
When deleting by type, multiple measurements can be targeted using inheritance.
For example, `delete!(mc, IsingMeasurement)` wil delete all
`IsingEnergyMeasurement` and `IsingMagnmetizationMeasurement` objects.

Examples:
```
delete!(mc, :conf)
delete!(mc, ConfigurationMeasurement)
delete!(mc, ConfigurationMeasurement, stage=:ME)
```
"""
function Base.delete!(mc::MonteCarloFlavor, key::Symbol, stage=:ME)
    measurement_stage = (:ME, :me, :Measurement, :measurement)
    thermalization_stage = (:TH, :th, :Thermalization, :thermalization, :thermalization_measurements)

    if stage in measurement_stage
        delete!(mc.measurements, key)
    elseif stage in thermalization_stage
        delete!(mc.thermalization_measurements, key)
    else
        throw(error(
            "The given stage `$stage` was not recognized. Try one of (" *
            join(string.([
                measurement_stage...,
                thermalization_stage...
            ]), ", ") * ")"
        ))
    end
end

function Base.delete!(mc::MonteCarloFlavor, MT::Type{<: AbstractMeasurement}, stage=:ME)
    measurement_stage = (:ME, :me, :Measurement, :measurement)
    thermalization_stage = (:TH, :th, :Thermalization, :thermalization, :thermalization_measurements)

    if stage in measurement_stage
        ks = collect(keys(mc.measurements))
        for k in ks
            if typeof(mc.measurements[k]) <: MT
                delete!(mc.measurements, k)
            end
        end
        mc.measurements
    elseif stage in thermalization_stage
        ks = collect(keys(mc.thermalization_measurements))
        for k in ks
            if typeof(mc.thermalization_measurements[k]) <: MT
                delete!(mc.thermalization_measurements, k)
            end
        end
        mc.thermalization_measurements
    else
        throw(error(
            "The given stage `$stage` was not recognized. Try one of (" *
            join(string.([
                measurement_stage...,
                thermalization_stage...
            ]), ", ") * ")"
        ))
    end
end

################################################################################
### FileIO
################################################################################

# This is just for dispath
struct Measurements end

"""
    save_measurements!(mc, filename[, entryname=""; overwrite=false, rename=true])

Saves all measurements to `filename`.

If `overwrite = true` the file will
be overwritten if it already exists. If `rename = true` random characters
will be added to the filename until it becomes unique.
"""
function save_measurements(
        filename::String, mc::MonteCarloFlavor, entryname::String="";
        backend = endswith(filename, "jld2") ? JLD2 : JLD,
        overwrite = false, rename = true
    )
    isfile(filename) && !overwrite && !rename && throw(ErrorException(
        "Cannot save because \"$filename\" already exists. Consider setting " *
        "`reanme = true` to adjust the filename or `overwrite = true`" *
        " to overwrite the file."
    ))
    if isfile(filename) && !overwrite && rename
        while isfile(filename)
            # those map to 0-9, A-Z, a-z
            x = rand([(48:57)..., (65:90)..., (97:122)...])
            s = string(Char(x))
            filename = filename[1:end-4] * s * ".jld"
        end
    end

    mode = isfile(filename) ? "r+" : "w"
    file = backend.jldopen(filename, mode)
    save_measurements(file, mc, entryname)
    close(file)
    filename
end
function save_measurements(file::JLDFile, mc::MonteCarloFlavor, entryname::String="")
    !isempty(entryname) && !endswith(entryname, "/") && (entryname *= "/")
    write(file, entryname * "VERSION", 1)
    write(file, entryname * "tag", "Measurements")
    measurement_dict = measurements(mc, :all)
    for (k0, v0) in measurement_dict # :TH or :ME
        for (k1, meas) in v0 # Measurement name
            _entryname = entryname * "$k0/$k1"
            _save(file, meas, _entryname)
        end
    end
end

to_tag(::Type{<: Measurements}) = Val(:Measurements)
function _load(data, ::Val{:Measurements})
    if !(data["VERSION"] == 1)
        throw(ErrorException("Failed to load measurements version $(data["VERSION"])"))
    end
    !haskey(data, "ME") && @warn "No measurement stage found (key \"ME\" missing)"
    !haskey(data, "TH") && @debug "No thermalization stage found (key \"TH\" missing)"

    Dict{Symbol, Dict{Symbol, AbstractMeasurement}}(
        :TH => if haskey(data, "TH")
            Dict{Symbol, AbstractMeasurement}(
                Symbol(k) => _load(data["TH"][k], to_tag(data["TH"][k]))
                for k in keys(data["TH"])
            )
        else
            Dict{Symbol, AbstractMeasurement}()
        end,
        :ME => if haskey(data, "ME")
            Dict{Symbol, AbstractMeasurement}(
                Symbol(k) => _load(data["ME"][k], to_tag(data["ME"][k]))
                for k in keys(data["ME"])
            )
        else
            Dict{Symbol, AbstractMeasurement}()
        end
    )
end


text:
using MOCNeutronTransport
using StaticArrays
@testset "Point_2D" begin
    @testset "$F" for F in [Float32, Float64, BigFloat]
        @testset "Operators" begin
            p₁ = Point_2D{F}(1, 2)
            p₂ = Point_2D{F}(2, 4)

            # Point_2D isapprox
            p = Point_2D{F}(1, 2 - 10*eps(F))
            @test F(2) ≈ 2 - 10*eps(F)
            @test p ≈ p₁

            # Point_2D addition
            p = p₁ + p₂
            @test p == [3, 6]

            # Point_2D subtraction
            p = p₁ - p₂
            @test p == [-1, -2]

            # Cross product
            p₁ = Point_2D{F}(2, 3)
            p₂ = Point_2D{F}(5, 6)
            v = p₁ × p₂
            @test v ≈ -3

            # Dot product
            @test p₁ ⋅ p₂ ≈ 10 + 18

            # Number addition
            p₁ = Point_2D{F}(1, 2)
            p₂ = Point_2D{F}(2, 4)
            p = p₁ .+ 1
            @test p  ≈ [2, 3]

            # Subtraction
            p = p₁ .- 1
            @test p == [0, 1]

            # Multiplication
            p = 4*p₁
            @test p == [4, 8]

            # Division
            p = p₁/4
            @test p == [1//4, 1//2]

            # Unary -
            p = -p₁
            @test p == [-1, -2]

            # Matrix multiplcation
            A = SMatrix{2, 2, F, 4}(1, 2, 3, 4)
            p = Point_2D{F}(1, 2)
            q = A*p
            @test q[1] ≈ 7
            @test q[2] ≈ 10
            @test typeof(q) == Point_2D{F}
        end

        @testset "Methods" begin
            p₁ = Point_2D{F}(1, 2)
            p₂ = Point_2D{F}(2, 4)

            # distance
            @test distance(p₁, p₂) ≈ sqrt(5)

            # distance²
            @test distance²(p₁, p₂) ≈ 5

            # norm
            @test norm(p₁) ≈ sqrt(5)
            @test norm(p₂) ≈ sqrt(4 + 16)

            # norm²
            @test norm²(p₁) ≈ 5
            @test norm²(p₂) ≈ 4 + 16

            # midpoint
            mp = midpoint(p₁, p₂)
            @test mp[1] ≈ 3//2
            @test mp[2] ≈ 3

            # sort_points
            p₁ = Point_2D{F}(1, 0)
            p₂ = Point_2D{F}(2, 0)
            p₃ = Point_2D{F}(3, 0)
            points = [p₃, p₁, p₂]
            points_sorted = sortpoints(Point_2D{F}(0, 0), points)
            @test points_sorted[1] == p₁
            @test points_sorted[2] == p₂
            @test points_sorted[3] == p₃
            # mutating
            sortpoints!(Point_2D{F}(0, 0), points)
            @test points[1] == p₁
            @test points[2] == p₂
            @test points[3] == p₃
        end
    end
end

ext: jl

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
export iauGst94
"""
Greenwich apparent sidereal time (consistent with IAU 1982/94
resolutions).

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   uta,utb    double    UT1 as a 2-part Julian Date (Notes 1,2)

Returned (function value):
               double    Greenwich apparent sidereal time (radians)

Notes:

   1. The UT1 date uta+utb is a Julian Date, apportioned in any
   convenient way between the argument pair.  For example,
   JD=2450123.7 could be expressed in any of these ways, among
   others:

            uta            utb

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in cases
   where the loss of several decimal digits of resolution is
   acceptable.  The J2000 and MJD methods are good compromises
   between resolution and convenience.  For UT, the date & time
   method is best matched to the algorithm that is used by the Earth
   Rotation Angle function, called internally:  maximum precision is
   delivered when the uta argument is for 0hrs UT1 on the day in
   question and the utb argument lies in the range 0 to 1, or vice
   versa.

   2. The result is compatible with the IAU 1982 and 1994 resolutions,
   except that accuracy has been compromised for the sake of
   convenience in that UT is used instead of TDB (or TT) to compute
   the equation of the equinoxes.

   3. This GAST must be used only in conjunction with contemporaneous
   IAU standards such as 1976 precession, 1980 obliquity and 1982
   nutation.  It is not compatible with the IAU 2000 resolutions.

   4. The result is returned in the range 0 to 2pi.

Called:
   iauGmst82    Greenwich mean sidereal time, IAU 1982
   iauEqeq94    equation of the equinoxes, IAU 1994
   iauAnp       normalize angle into range 0 to 2pi

References:

   Explanatory Supplement to the Astronomical Almanac,
   P. Kenneth Seidelmann (ed), University Science Books (1992)

   IAU Resolution C7, Recommendation 3 (1994)

This revision:  2008 May 16

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# double iauGst94(double uta, double utb)

function iauGst94(uta::Real, utb::Real)
   return ccall((:iauGst94, libsofa_c), Cdouble,
               (Cdouble, Cdouble),
               convert(Float64, uta), convert(Float64, utb))
end


text:
function differentiate_node!(program::Program, derived_values::Dict{Int,Int}, call_node::CallNode)
    args = call_node.args

    derived_value = nothing

    for arg_i in 1:length(args)
        arg = args[arg_i]

        if !(arg isa TrackedValue)
            continue
        end

        arg::TrackedValue

        if !haskey(derived_values, arg.id)
            continue
        end

        arg_derived = derived_values[arg.id]

        f_der_lin = diff_rule(call_node.f, Val(arg_i), args...)

        # TODO: handle case that not TrackedValue
        return_value = f_der_lin(TrackedValue(arg_derived, program))::TrackedValue

        if derived_value === nothing
            derived_value = return_value
        else
            derived_value = derived_value::TrackedValue + return_value
        end
    end

    derived_value::TrackedValue
    return derived_value.id
end

function differentiate_node!(program::Program, derived_values::Dict{Int,Int}, tuple_node::TupleNode)
    args = tuple_node.args

    derived_args = Any[]

    for arg in tuple_node.args
        if !(arg isa TrackedValue)
            push!(derived_args, zero(arg))
            continue
        end

        arg::TrackedValue

        if !haskey(derived_values, arg.id)
            push!(derived_args, add_zero!(program, arg.id))
            continue
        end

        push!(derived_args, TrackedValue(derived_values[arg.id], program))
    end

    derived_value = add_tuple!(program, derived_args)
    return derived_value.id
end

function differentiate!(program::Program)
    # TODO: remember when program is in a compact state?
    compact!(program)

    executable_order = topological_sort(program.graph)

    # we only consider the nodes that depend on x
    intersect!(executable_order, descendants(program.graph, Set(program.args[1][1])))

    push!(program.nodes, ArgumentNode())
    y_node_id = create_node!(program.graph)

    y_value = TrackedValue(y_node_id, program)

    push!(program.args, Int[y_node_id])

    # denote by x the first argument
    # denote by y the argument to the linearized map
    # denote by f_i the map x -> si. then derived_values[i] represents f_i'(x)[y]
    derived_values = Dict{Int,Int}()
    derived_values[program.args[1][1]] = y_node_id

    for node_id in executable_order
        node = program.nodes[node_id]
        derived_values[node_id] = differentiate_node!(program, derived_values, node)
    end

    if !haskey(derived_values, program.return_node)
        # return node does not depend on x
        derived_values[program.return_node] = add_zero!(program, program.return_node).id
    end

    program.return_node = derived_values[program.return_node]

    compact!(program)

    return program
end

differentiate(program::Program) = differentiate!(deepcopy(program))


text:
# Reference journal article:
# Nikkhoo, M., Walter, T. R., Lundgren, P. R., Prats-Iraola, P. (2016):
# Compound dislocation models (CDMs) for volcano deformation analyses.
# Submitted to Geophysical Journal International
#
# Copyright (c) 2016 Mehdi Nikkhoo
#
# 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.
#
# I appreciate any comments or bug reports.
#
# Mehdi Nikkhoo
# created: 2013.1.25
# Last modified: 2016.10.18
#
# Section 2.1, Physics of Earthquakes and Volcanoes
# Department 2, Geophysics
# Helmholtz Centre Potsdam
# German Research Centre for Geosciences (GFZ)
#
# email:
# mehdi.nikkhoo@gfz-potsdam.de
# mehdi.nikkhoo@gmail.com
#
# website:
# http://www.volcanodeformation.com
#
# Translated by Pengcheng Shi (shipengcheng1230@gmail.com) 12/2018
# The broadcast here causes tremendous amount of compiliation time, needs rewriting!

# export disp_rect_fs, strain_rect_fs, stress_rect_fs

# using LinearAlgebra

# function disp_rect_fs(x::V, y::V, z::V, x₀::T, y₀::T, depth::T, L::T, W::T, plunge::T, dip::T, strike::T, rake::T, slip::T, opening::T, nu::T, ref::S) where {V<:AbstractVector, T, S}
#     @assert size(x) == size(y) == size(z) "Unmatched length of observational coordinates."

#     bx = opening
#     by = slip * cosd(rake)
#     bz = slip * sind(rake)
#     Rz1 = [cosd(plunge) sind(plunge) zero(T); -sind(plunge) cosd(plunge) zero(T); zero(T) zero(T) one(T)]
#     Ry = [cosd(dip) zero(T) sind(dip); zero(T) one(T) zero(T); -sind(dip) zero(T) cosd(dip)]
#     Rz2 = [cosd(strike) sind(strike) zero(T); -sind(strike) cosd(strike) zero(T); zero(T) zero(T) one(T)]
#     Rt = Rz2 * Ry * Rz1

#     Pt1 = [-W/2, L/2, zero(T)]
#     Pt2 = [-W/2, -L/2, zero(T)]
#     Pt3 = [W/2, -L/2, zero(T)]
#     Pt4 = [W/2, L/2, zero(T)]

#     Ptr = _ref_pts(ref, L, W)
#     Pr = [x₀, y₀, -depth] - Rt * Ptr
#     P1 = Rt * Pt1 + Pr
#     P2 = Rt * Pt2 + Pr
#     P3 = Rt * Pt3 + Pr
#     P4 = Rt * Pt4 + Pr

#     ez = [zero(T), zero(T), one(T)]
#     Vnorm = Rt * ez
#     Vstrike = [sind(strike), cosd(strike), zero(T)]
#     Vdip = Vnorm × Vstrike

#     Pm = (P1 + P2 + P3 + P4) / 4
#     p1, p2, p3, p4 = [zeros(T, 3) for _ in 1: 4]
#     At = hcat(Vnorm, Vstrike, Vdip)'

#     x, y, z = coord_trans(x .- Pm[1], y .- Pm[2], z .- Pm[3], At)
#     p1[1], p1[2], p1[3] = coord_trans(P1[1] - Pm[1], P1[2] - Pm[2], P1[3] - Pm[3], At)
#     p2[1], p2[2], p2[3] = coord_trans(P2[1] - Pm[1], P2[2] - Pm[2], P2[3] - Pm[3], At)
#     p3[1], p3[2], p3[3] = coord_trans(P3[1] - Pm[1], P3[2] - Pm[2], P3[3] - Pm[3], At)
#     p4[1], p4[2], p4[3] = coord_trans(P4[1] - Pm[1], P4[2] - Pm[2], P4[3] - Pm[3], At)

#     e12 = normalize(p2 - p1)
#     e23 = normalize(p3 - p2)
#     e34 = normalize(p4 - p3)
#     e14 = normalize(p4 - p1)

#     Rectmode = @views rectmodefInder(y, z, x, p1[2:3], p2[2:3], p3[2:3], p4[2:3])

#     # use list comprehension causes being unable infer type (until Julia v1.2), see one example mentioned in the link below
#     # https://discourse.julialang.org/t/comprehension-type-inference/4324
#     # u, v, w = [zeros(T, length(x)) for _ ∈ 1: 3]
#     u = zeros(T, length(x)); v = zeros(T, length(x)); w = zeros(T, length(x))

#     casepLog = Rectmode .== 1
#     casenLog = Rectmode .== -1
#     casezLog = Rectmode .== 0

#     if any(casepLog)
#         xp, yp, zp = x[casepLog], y[casepLog], z[casepLog]
#         utmp, vtmp, wtmp = [Vector{T}(undef, length(xp)) for _ ∈ 1: 3]

#         utmp, vtmp, wtmp = RDSetupD(xp, yp, zp, bx, by, bz, nu, p1, e14)
#         u[casepLog] .+= utmp; v[casepLog] .+= vtmp; w[casepLog] .+= wtmp

#         utmp, vtmp, wtmp = RDSetupD(xp, yp, zp, bx, by, bz, nu, p2, -e12)
#         u[casepLog] .+= utmp; v[casepLog] .+= vtmp; w[casepLog] .+= wtmp

#         utmp, vtmp, wtmp = RDSetupD(xp, yp, zp, bx, by, bz, nu, p3, -e23)
#         u[casepLog] .+= utmp; v[casepLog] .+= vtmp; w[casepLog] .+= wtmp

#         utmp, vtmp, wtmp = RDSetupD(xp, yp, zp, bx, by, bz, nu, p4, -e34)
#         u[casepLog] .+= utmp; v[casepLog] .+= vtmp; w[casepLog] .+= wtmp
#     end
#     if any(casenLog)
#         xn, yn, zn = x[casenLog], y[casenLog], z[casenLog]
#         utmp, vtmp, wtmp = [Vector{T}(undef, length(xn)) for _ ∈ 1: 3]

#         utmp, vtmp, wtmp = RDSetupD(xn, yn, zn, bx, by, bz, nu, p1, -e14)
#         u[casenLog] .+= utmp; v[casenLog] .+= vtmp; w[casenLog] .+= wtmp

#         utmp, vtmp, wtmp = RDSetupD(xn, yn, zn, bx, by, bz, nu, p2, e12)
#         u[casenLog] .+= utmp; v[casenLog] .+= vtmp; w[casenLog] .+= wtmp

#         utmp, vtmp, wtmp = RDSetupD(xn, yn, zn, bx, by, bz, nu, p3, e23)
#         u[casenLog] .+= utmp; v[casenLog] .+= vtmp; w[casenLog] .+= wtmp

#         utmp, vtmp, wtmp = RDSetupD(xn, yn, zn, bx, by, bz, nu, p4, e34)
#         u[casenLog] .+= utmp; v[casenLog] .+= vtmp; w[casenLog] .+= wtmp
#     end
#     if any(casezLog)
#         u[casezLog] .= NaN; v[casezLog] .= NaN; w[casezLog] .= NaN
#     end

#     Fi = BurgersFuncRD(x, y, z, p1, p2, p3, p4)
#     @. u += Fi * bx
#     @. v += Fi * by
#     @. w += Fi * bz
#     u, v, w = coord_trans(u, v, w, At')
#     return u, v, w
# end

# function rd_disp_hs()
#     #TODO
# end

# function strain_rect_fs(x::V, y::V, z::V, x₀::T, y₀::T, depth::T, L::T, W::T, plunge::T, dip::T, strike::T, rake::T, slip::T, opening::T, λ::T, μ::T, ref::S) where {V<:AbstractVector, T, S}
#     nu = λ / (λ + μ) / 2

#     # lots of repeated codes here
#     bx = opening
#     by = slip * cosd(rake)
#     bz = slip * sind(rake)
#     Rz1 = [cosd(plunge) sind(plunge) zero(T); -sind(plunge) cosd(plunge) zero(T); zero(T) zero(T) one(T)]
#     Ry = [cosd(dip) zero(T) sind(dip); zero(T) one(T) zero(T); -sind(dip) zero(T) cosd(dip)]
#     Rz2 = [cosd(strike) sind(strike) zero(T); -sind(strike) cosd(strike) zero(T); zero(T) zero(T) one(T)]
#     Rt = Rz2 * Ry * Rz1

#     Pt1 = [-W/2, L/2, zero(T)]
#     Pt2 = [-W/2, -L/2, zero(T)]
#     Pt3 = [W/2, -L/2, zero(T)]
#     Pt4 = [W/2, L/2, zero(T)]

#     Ptr = _ref_pts(ref, L, W)
#     Pr = [x₀, y₀, -depth] - Rt * Ptr
#     P1 = Rt * Pt1 + Pr
#     P2 = Rt * Pt2 + Pr
#     P3 = Rt * Pt3 + Pr
#     P4 = Rt * Pt4 + Pr

#     ez = [zero(T), zero(T), one(T)]
#     Vnorm = Rt * ez
#     Vstrike = [sind(strike), cosd(strike), zero(T)]
#     Vdip = Vnorm × Vstrike

#     Pm = (P1 + P2 + P3 + P4) / 4
#     p1, p2, p3, p4 = [zeros(T, 3) for _ in 1: 4]
#     At = hcat(Vnorm, Vstrike, Vdip)'

#     x, y, z = coord_trans(x .- Pm[1], y .- Pm[2], z .- Pm[3], At)
#     p1[1], p1[2], p1[3] = coord_trans(P1[1] - Pm[1], P1[2] - Pm[2], P1[3] - Pm[3], At)
#     p2[1], p2[2], p2[3] = coord_trans(P2[1] - Pm[1], P2[2] - Pm[2], P2[3] - Pm[3], At)
#     p3[1], p3[2], p3[3] = coord_trans(P3[1] - Pm[1], P3[2] - Pm[2], P3[3] - Pm[3], At)
#     p4[1], p4[2], p4[3] = coord_trans(P4[1] - Pm[1], P4[2] - Pm[2], P4[3] - Pm[3], At)

#     e12 = normalize(p2 - p1)
#     e23 = normalize(p3 - p2)
#     e34 = normalize(p4 - p3)
#     e14 = normalize(p4 - p1)

#     Rectmode = @views rectmodefInder(y, z, x, p1[2:3], p2[2:3], p3[2:3], p4[2:3])

#     casepLog = Rectmode .== 1
#     casenLog = Rectmode .== -1
#     casezLog = Rectmode .== 0
#     lenx = length(x)

#     # use list comprehension causes being unable infer type (until Julia v1.2), see one example mentioned in the link below
#     # https://discourse.julialang.org/t/comprehension-type-inference/4324
#     # exx, eyy, ezz, exy, exz, eyz = [zeros(T, length(x)) for _ in 1: 6]
#     exx, eyy, ezz, exy, exz, eyz = zeros(T, lenx), zeros(T, lenx), zeros(T, lenx), zeros(T, lenx), zeros(T, lenx), zeros(T, lenx)

#     if any(casepLog)
#         xp, yp, zp = x[casepLog], y[casepLog], z[casepLog]
#         exxt, eyyt, ezzt, exyt, exzt, eyzt = [Vector{T}(undef, length(xp)) for _ ∈ 1: 6]

#         exxt, eyyt, ezzt, exyt, exzt, eyzt = RDSetupS(xp, yp, zp, bx, by, bz, nu, p1, e14)
#         exx[casepLog] .+= exxt; eyy[casepLog] .+= eyyt; ezz[casepLog] .+= ezzt; exy[casepLog] .+= exyt; exz[casepLog] .+= exzt; eyz[casepLog] .+= eyzt

#         exxt, eyyt, ezzt, exyt, exzt, eyzt = RDSetupS(xp, yp, zp, bx, by, bz, nu, p2, -e12)
#         exx[casepLog] .+= exxt; eyy[casepLog] .+= eyyt; ezz[casepLog] .+= ezzt; exy[casepLog] .+= exyt; exz[casepLog] .+= exzt; eyz[casepLog] .+= eyzt

#         exxt, eyyt, ezzt, exyt, exzt, eyzt = RDSetupS(xp, yp, zp, bx, by, bz, nu, p3, -e23)
#         exx[casepLog] .+= exxt; eyy[casepLog] .+= eyyt; ezz[casepLog] .+= ezzt; exy[casepLog] .+= exyt; exz[casepLog] .+= exzt; eyz[casepLog] .+= eyzt

#         exxt, eyyt, ezzt, exyt, exzt, eyzt = RDSetupS(xp, yp, zp, bx, by, bz, nu, p4, -e34)
#         exx[casepLog] .+= exxt; eyy[casepLog] .+= eyyt; ezz[casepLog] .+= ezzt; exy[casepLog] .+= exyt; exz[casepLog] .+= exzt; eyz[casepLog] .+= eyzt
#     end
#     if any(casenLog)
#         xn, yn, zn = x[casenLog], y[casenLog], z[casenLog]
#         exxt, eyyt, ezzt, exyt, exzt, eyzt = [Vector{T}(undef, length(xn)) for _ ∈ 1: 6]

#         exxt, eyyt, ezzt, exyt, exzt, eyzt = RDSetupS(xn, yn, zn, bx, by, bz, nu, p1, -e14)
#         exx[casenLog] .+= exxt; eyy[casenLog] .+= eyyt; ezz[casenLog] .+= ezzt; exy[casenLog] .+= exyt; exz[casenLog] .+= exzt; eyz[casenLog] .+= eyzt

#         exxt, eyyt, ezzt, exyt, exzt, eyzt = RDSetupS(xn, yn, zn, bx, by, bz, nu, p2, e12)
#         exx[casenLog] .+= exxt; eyy[casenLog] .+= eyyt; ezz[casenLog] .+= ezzt; exy[casenLog] .+= exyt; exz[casenLog] .+= exzt; eyz[casenLog] .+= eyzt

#         exxt, eyyt, ezzt, exyt, exzt, eyzt = RDSetupS(xn, yn, zn, bx, by, bz, nu, p3, e23)
#         exx[casenLog] .+= exxt; eyy[casenLog] .+= eyyt; ezz[casenLog] .+= ezzt; exy[casenLog] .+= exyt; exz[casenLog] .+= exzt; eyz[casenLog] .+= eyzt

#         exxt, eyyt, ezzt, exyt, exzt, eyzt = RDSetupS(xn, yn, zn, bx, by, bz, nu, p4, e34)
#         exx[casenLog] .+= exxt; eyy[casenLog] .+= eyyt; ezz[casenLog] .+= ezzt; exy[casenLog] .+= exyt; exz[casenLog] .+= exzt; eyz[casenLog] .+= eyzt
#     end
#     if any(casezLog)
#         exx[casezLog] .= NaN; eyy[casezLog] .= NaN; ezz[casezLog] .= NaN; exy[casezLog] .= NaN; exz[casezLog] .= NaN; eyz[casezLog] .= NaN
#     end
#     exx, eyy, ezz, exy, exz, eyz = TensTrans(exx, eyy, ezz, exy, exz, eyz, At')
# end

# function rd_strain_hs()
#     #TODO
# end

# function stress_rect_fs(x::V, y::V, z::V, x₀::T, y₀::T, depth::T, L::T, W::T, plunge::T, dip::T, strike::T, rake::T, slip::T, opening::T, λ::T, μ::T, ref::S) where {V, T, S}
#     exx, eyy, ezz, exy, exz, eyz = strain_rect_fs(x, y, z, x₀, y₀, depth, L, W, plunge, dip, strike, rake, slip, opening, λ, μ, ref)
#     ekk = exx + eyy + ezz
#     Sxx = @. 2μ * exx + λ * ekk
#     Syy = @. 2μ * eyy + λ * ekk
#     Szz = @. 2μ * ezz + λ * ekk
#     Sxy = @. 2μ * exy
#     Sxz = @. 2μ * exz
#     Syz = @. 2μ * eyz
#     return Sxx, Syy, Szz, Sxy, Sxz, Syz
# end

# function rd_stress_hs()
#     #TODO
# end

# @inline _ref_pts(::Val{:pc}, L::T, W::T) where T = zeros(T, 3)
# @inline _ref_pts(::Val{:p1}, L::T, W::T) where T = [-W/2, L/2, zero(T)]
# @inline _ref_pts(::Val{:p2}, L::T, W::T) where T = [-W/2, -L/2, zero(T)]
# @inline _ref_pts(::Val{:p3}, L::T, W::T) where T = [W/2, -L/2, zero(T)]
# @inline _ref_pts(::Val{:p4}, L::T, W::T) where T = [W/2, L/2, zero(T)]
# @inline _ref_pts(::Val{:mp12}, L::T, W::T) where T = [-W/2, zero(T), zero(T)]
# @inline _ref_pts(::Val{:mp23}, L::T, W::T) where T = [zero(T), -L/2, zero(T)]
# @inline _ref_pts(::Val{:mp34}, L::T, W::T) where T = [W/2, zero(T), zero(T)]
# @inline _ref_pts(::Val{:mp41}, L::T, W::T) where T = [zero(T), L/2, zero(T)]
# @inline _ref_pts(::Val{:mp21}, L::T, W::T) where T = _ref_pts(Val(:mp12), L, W)
# @inline _ref_pts(::Val{:mp32}, L::T, W::T) where T = _ref_pts(Val(:mp23), L, W)
# @inline _ref_pts(::Val{:mp43}, L::T, W::T) where T = _ref_pts(Val(:mp34), L, W)
# @inline _ref_pts(::Val{:mp14}, L::T, W::T) where T = _ref_pts(Val(:mp14), L, W)

# @inline function rectmodefInder(x::T, y::T, z::T, p1::U, p2::U, p3::U, p4::U) where {T, U}
#     pm = (p1 + p2 + p3 + p4) / 4
#     e21 = normalize(p1 - p2)
#     e41 = normalize(p1 - p4)
#     A = hcat(e21, e41)
#     rectmode = ones(Int, length(x))
#     r = A' * hcat(p1 - pm, p2 - pm, p3 - pm, p4 - pm)
#     P1, P2, P3, P4 = @views r[:,1], r[:,2], r[:,3], r[:,4]

#     @inbounds @simd for i ∈ eachindex(rectmode)
#         _x, _y = x[i] - pm[1], y[i] - pm[2]
#         r1 = A[1,1] * _x + A[2,1] * _y
#         r2 = A[1,2] * _x + A[2,2] * _y
#         _x, _y = r1, r2
#         if (_x ≥ 0 && _y ≥ 0 && ((_y - P1[2]) < (_x - P1[1]))) ||
#             (_x ≤ 0 && _y ≥ 0 && ((_y - P2[2]) > -(_x - P2[1]))) ||
#             (_x ≤ 0 && _y ≤ 0 && ((_y - P3[2]) > (_x - P3[1]))) ||
#             (_x ≥ 0 && _y ≤ 0 && ((_y - P4[2]) < -(_x - P4[1]))) ||
#             (_x < P1[1] && _x > P3[1] && _y < P1[2] && _y > P3[2])
#             rectmode[i] = -1
#         elseif ((_x == P1[1] || _x == P3[1]) && _y ≤ P1[2] && _y ≥ P3[2] && z[i] == 0) ||
#             ((_y == P1[2] || _y == P3[2]) && _x ≤ P1[1] && _x ≥ P3[1] && z[i] == 0)
#             rectmode[i] = 0
#         end
#     end
#     rectmode
# end

# @inline function RDSetupD(x::U, y::U, z::U, bx::T, by::T, bz::T, nu::T, RDVertex::V, SideVec::V) where {T, U, V}
#     y1 = @. SideVec[3] * (y - RDVertex[2]) - SideVec[2] * (z - RDVertex[3])
#     z1 = @. SideVec[2] * (y - RDVertex[2]) + SideVec[3] * (z - RDVertex[3])
#     by1 = @. SideVec[3] * by - SideVec[2] * bz
#     bz1 = @. SideVec[2] * by + SideVec[3] * bz

#     u, v0, w0 = AngDisDisp(x, y1, z1, -π/2, bx, by1, bz1, nu)
#     v = @. SideVec[3] * v0 + SideVec[2] * w0
#     w = @. -SideVec[2] * v0 + SideVec[3] * w0
#     return u, v, w
# end

# @inline function BurgersFuncRD(x::T, y::T, z::T, p1::U, p2::U, p3::U, p4::U) where {T, U}
#     Fi = zeros(size(x))
#     Ind = @. (abs(y) ≤ abs(p1[2] - p2[2])/2) & (abs(z) ≤ abs(p1[3] - p4[3])/2)

#     if any(Ind)
#         xI, yI, zI = x[Ind], y[Ind], z[Ind]

#         FiD1 = @. √(xI ^ 2 + (yI - p1[2]) ^ 2 + (zI - p1[3]) ^ 2) - (zI - p1[3]) - (yI - p1[2])
#         FiD2 = @. √(xI ^ 2 + (yI - p2[2]) ^ 2 + (zI - p2[3]) ^ 2) - (zI - p2[3]) - (yI - p2[2])
#         FiD3 = @. √(xI ^ 2 + (yI - p3[2]) ^ 2 + (zI - p3[3]) ^ 2) - (zI - p3[3]) - (yI - p3[2])
#         FiD4 = @. √(xI ^ 2 + (yI - p4[2]) ^ 2 + (zI - p4[3]) ^ 2) - (zI - p4[3]) - (yI - p4[2])

#         FiNt = @. xI * (FiD1 + FiD3) * (FiD2 * FiD4 - xI ^ 2) - xI * (FiD2 + FiD4) * (FiD1 * FiD3 - xI ^ 2)
#         FiDt = @. (FiD1 * FiD3 - xI ^ 2) * (FiD2 * FiD4 - xI ^ 2) + xI ^ 2 * (FiD1 + FiD3) * (FiD2 + FiD4)
#         Fi[Ind] = @. atan(FiNt, FiDt) / 2 / π
#     elseif !all(Ind)
#         bool = .!Ind
#         xO, yO, zO = x[bool], y[bool], z[bool]
#         a = cat(-xO, p1[2] .- yO, p1[3] .- zO; dims=2)
#         b = cat(-xO, p2[2] .- yO, p2[3] .- zO; dims=2)
#         c = cat(-xO, p3[2] .- yO, p3[3] .- zO; dims=2)
#         d = cat(-xO, p4[2] .- yO, p4[3] .- zO; dims=2)
#         na, nb, nc, nd = vec(.√(sum(abs2, a; dims=2))), vec(.√(sum(abs2, b; dims=2))), vec(.√(sum(abs2, c; dims=2))), vec(.√(sum(abs2, d; dims=2)))
#         @views begin
#             FiN1t = @. a[:,1] * (b[:,2] * c[:,3] - b[:,3] * c[:,2]) - a[:,2] * (b[:,1] * c[:,3]-b[:,3] * c[:,1]) + a[:,3] * (b[:,1] * c[:,2]-b[:,2] * c[:,1])
#             FiD1t = @. na * nb * nc + $vec($sum(a * b; dims=2)) * nc + $vec($sum(a * c; dims=2)) * nb + $vec($sum(b * c; dims=2)) * na
#             FiN2t = @. a[:,1] * (c[:,2] * d[:,3] - c[:,3] * d[:,2]) - a[:,2] * (c[:,1] * d[:,3] - c[:,3] * d[:,1]) + a[:,3] * (c[:,1] * d[:,2]-c[:,2] * d[:,1])
#             FiD2t = @. na * nc * nd + $vec($sum(a * c; dims=2)) * nd + $vec($sum(a * d; dims=2)) * nc + $vec($sum(c * d; dims=2)) * na
#         end
#         Fi[bool] = @. -2 * atan(FiN1t * FiD2t + FiN2t * FiD1t, FiD1t * FiD2t - FiN1t * FiN2t) / 4 / π
#     end
#     return Fi
# end

# @inline function RDSetupS(x::T, y::T, z::T, bx::U, by::U, bz::U, nu::U, RDVertex::V, SideVec::V) where {T, U, V}
#     y1 = @. SideVec[3] * (y - RDVertex[2]) - SideVec[2] * (z - RDVertex[3])
#     z1 = @. SideVec[2] * (y - RDVertex[2]) + SideVec[3] * (z - RDVertex[3])
#     by1 = @. SideVec[3] * by - SideVec[2] * bz
#     bz1 = @. SideVec[2] * by + SideVec[3] * bz

#     exx, eyy, ezz, exy, exz, eyz = AngDisStrain(x, y1, z1, -pi/2, bx, by1, bz1, nu)

#     B = [1 0 0; 0 SideVec[3] SideVec[2]; 0 -SideVec[2] SideVec[3]]
#     exx, eyy, ezz, exy, exz, eyz = TensTrans(exx, eyy, ezz, exy, exz, eyz, B)
#     return exx, eyy, ezz, exy, exz, eyz
# end

# function disp_rect_fs(x::V, y::V, z::V, x₀::T, y₀::T, depth::T, L::T, W::T, plunge::T, dip::T, strike::T, rake::T, slip::T, opening::T, nu::T, ref::S) where {V<:Number, T, S}
#     u, v, w = disp_rect_fs([x], [y], [z], x₀, y₀, depth, L, W, plunge, dip, strike, rake, slip, opening, nu, ref)
#     u[1], v[1], w[1]
# end

# function strain_rect_fs(x::V, y::V, z::V, x₀::T, y₀::T, depth::T, L::T, W::T, plunge::T, dip::T, strike::T, rake::T, slip::T, opening::T, λ::T, μ::T, ref::S) where {V<:Number, T, S}
#     out = strain_rect_fs([x], [y], [z], x₀, y₀, depth, L, W, plunge, dip, strike, rake, slip, opening, λ, μ, ref)
#     [q[1] for q in out]
# end


text:
# Julia wrapper for header: /usr/include/scip/scip_relax.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeRelax(scip, name, desc, priority, freq, relaxcopy, relaxfree, relaxinit, relaxexit, relaxinitsol, relaxexitsol, relaxexec, relaxdata)
    ccall((:SCIPincludeRelax, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Cstring, Cint, Cint, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{SCIP_RELAXDATA}), scip, name, desc, priority, freq, relaxcopy, relaxfree, relaxinit, relaxexit, relaxinitsol, relaxexitsol, relaxexec, relaxdata)
end

function SCIPincludeRelaxBasic(scip, relaxptr, name, desc, priority, freq, relaxexec, relaxdata)
    ccall((:SCIPincludeRelaxBasic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_RELAX}}, Cstring, Cstring, Cint, Cint, Ptr{Cvoid}, Ptr{SCIP_RELAXDATA}), scip, relaxptr, name, desc, priority, freq, relaxexec, relaxdata)
end

function SCIPsetRelaxCopy(scip, relax, relaxcopy)
    ccall((:SCIPsetRelaxCopy, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_RELAX}, Ptr{Cvoid}), scip, relax, relaxcopy)
end

function SCIPsetRelaxFree(scip, relax, relaxfree)
    ccall((:SCIPsetRelaxFree, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_RELAX}, Ptr{Cvoid}), scip, relax, relaxfree)
end

function SCIPsetRelaxInit(scip, relax, relaxinit)
    ccall((:SCIPsetRelaxInit, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_RELAX}, Ptr{Cvoid}), scip, relax, relaxinit)
end

function SCIPsetRelaxExit(scip, relax, relaxexit)
    ccall((:SCIPsetRelaxExit, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_RELAX}, Ptr{Cvoid}), scip, relax, relaxexit)
end

function SCIPsetRelaxInitsol(scip, relax, relaxinitsol)
    ccall((:SCIPsetRelaxInitsol, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_RELAX}, Ptr{Cvoid}), scip, relax, relaxinitsol)
end

function SCIPsetRelaxExitsol(scip, relax, relaxexitsol)
    ccall((:SCIPsetRelaxExitsol, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_RELAX}, Ptr{Cvoid}), scip, relax, relaxexitsol)
end

function SCIPfindRelax(scip, name)
    ccall((:SCIPfindRelax, libscip), Ptr{SCIP_RELAX}, (Ptr{SCIP_}, Cstring), scip, name)
end

function SCIPgetRelaxs(scip)
    ccall((:SCIPgetRelaxs, libscip), Ptr{Ptr{SCIP_RELAX}}, (Ptr{SCIP_},), scip)
end

function SCIPgetNRelaxs(scip)
    ccall((:SCIPgetNRelaxs, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPsetRelaxPriority(scip, relax, priority)
    ccall((:SCIPsetRelaxPriority, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_RELAX}, Cint), scip, relax, priority)
end


text:
using CellListMap, StaticArrays

# System properties
n = 100_000
N = n
sides = [250,250,250]
cutoff = 10

# Initialize linked lists and box structures
box = Box(sides,cutoff)
# masses
const mass = rand(N)

# Particle positions
x = [ box.sides .* rand(SVector{3,Float64}) for i in 1:n ]

cl = CellList(x,box)

# Function to be evalulated for each pair
function potential(x,y,i,j,d2,u,mass)
  d = sqrt(d2)
  u = u - 9.8*mass[i]*mass[j]/d
  return u
end

# Run pairwise computation
u = map_pairwise!((x,y,i,j,d2,u) -> potential(x,y,i,j,d2,u,mass),0.0,box,cl)


text:
sourcedir = joinpath(dirname(@__DIR__), "cmt")

builddir = joinpath(@__DIR__, "build")
if isdir(builddir)
    rm(builddir; recursive=true)
end
mkdir(builddir)

run(`cmake -S $sourcedir -B $builddir -DCMAKE_LIBRARY_OUTPUT_DIRECTORY=$(@__DIR__) -DCMAKE_BUILD_TYPE=Debug`)
run(`make -C $builddir -j$(Sys.CPU_THREADS) VERBOSE=1`)


text:
# Tests for conditional HINT layer (Kruse et al, 2020)
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

using InvertibleNetworks, LinearAlgebra, Test

#######################################################################################################################
# Test invertibility

# Use affine or additive coupling ConditionalLayerSLIM
type = "affine" # "affine", "additive" or "learned"

# X dimensions
nx1 = 32
nx2 = 32
nx_channel = 32
nx_hidden = 64
batchsize = 2

# Y dimensions
ny1 = 128
ny2 = 64
ny_channel = 1
ny_hidden = 32

# Linear operator
A = randn(Float32, ny1*ny2, nx1*nx2)

# Input image
X = glorot_uniform(nx1, nx2, nx_channel, batchsize)
Y = glorot_uniform(ny1, ny2, ny_channel, batchsize)

# Conditional HINT layer
CI = ConditionalLayerSLIM(nx1, nx2, nx_channel, nx_hidden, ny1, ny2, ny_channel, ny_hidden, batchsize; type=type)

# Forward/inverse
Zx, Zy, logdet = CI.forward(X, Y, A)
X_, Y_ = CI.inverse(Zx, Zy, A)

@test isapprox(norm(X - X_)/norm(X), 0f0; atol=1f-3)
@test isapprox(norm(Y - Y_)/norm(Y), 0f0; atol=1f-3)

# Forward/backward
Zx, Zy, logdet = CI.forward(X, Y, A)
X_, Y_ = CI.backward(0f0*Zx, 0f0*Zy, Zx, Zy, A)[3:4]

@test isapprox(norm(X - X_)/norm(X), 0f0; atol=1f-3)
@test isapprox(norm(Y - Y_)/norm(Y), 0f0; atol=1f-3)

# Forward/inverse Y-lane only
Zy = CI.forward_Y(Y)
Y_ = CI.inverse_Y(Zy)

@test isapprox(norm(Y - Y_)/norm(Y), 0f0; atol=1f-6)

#######################################################################################################################
# Gradient test

# Input image
X = randn(Float32, nx1, nx2, nx_channel, batchsize)
X0 = randn(Float32, nx1, nx2, nx_channel, batchsize)
dX = X - X0

# Input data
Y = randn(Float32, ny1, ny2, ny_channel, batchsize)
Y0 = randn(Float32, ny1, ny2, ny_channel, batchsize)
dY = Y - Y0

function loss(CI, X, Y, A)
    Zx, Zy, logdet = CI.forward(X, Y, A)
    f = -log_likelihood(Zx) -log_likelihood(Zy) - logdet
    ΔZx = -∇log_likelihood(Zx)
    ΔZy = -∇log_likelihood(Zy)
    ΔX, ΔY = CI.backward(ΔZx, ΔZy, Zx, Zy, A)[1:2]
    return f, ΔX, ΔY, CI.CL_X.CL[1].RB.W1.grad, CI.C_X.v1.grad
end

# Gradient test for input X, Y
CI = ConditionalLayerSLIM(nx1, nx2, nx_channel, nx_hidden, ny1, ny2, ny_channel, ny_hidden, batchsize; type=type)
f0, gX, gY = loss(CI, X0, Y0, A)[1:3]

maxiter = 5
h = 0.25f0
err1 = zeros(Float32, maxiter)
err2 = zeros(Float32, maxiter)

print("Gradient test ΔX\n")
for j=1:maxiter
    f = loss(CI, X0 + h*dX, Y0 + h*dY, A)[1]
    err1[j] = abs(f - f0)
    err2[j] = abs(f - f0 - h*dot(dX, gX) - h*dot(dY, gY))
    print(err1[j], "; ", err2[j], "\n")
    global h = h/2f0
end

@test isapprox(err1[end] / (err1[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err2[end] / (err2[1]/4^(maxiter-1)), 1f0; atol=1f1)

# Test for weights
X = randn(Float32, nx1, nx2, nx_channel, batchsize)
Y = randn(Float32, ny1, ny2, ny_channel, batchsize)
CI = ConditionalLayerSLIM(nx1, nx2, nx_channel, nx_hidden, ny1, ny2, ny_channel, ny_hidden, batchsize; type=type)
CI0 = ConditionalLayerSLIM(nx1, nx2, nx_channel, nx_hidden, ny1, ny2, ny_channel, ny_hidden, batchsize; type=type)

# Make weights larger (otherweise too close to zero after initialization)
CI.CL_X.CL[1].RB.W1.data .*= 4; CI0.CL_X.CL[1].RB.W1.data .*= 4
CI.C_X.v1.data .*= 4; CI0.C_X.v1.data .*= 4
CIini = deepcopy(CI0)

dW = CI.CL_X.CL[1].RB.W1.data - CI0.CL_X.CL[1].RB.W1.data
dv = CI.C_X.v1.data - CI0.C_X.v1.data

f0, gW, gv = loss(CI0, X, Y, A)[[1,4,5]]

maxiter = 5
h = 0.25f0
err3 = zeros(Float32, maxiter)
err4 = zeros(Float32, maxiter)

print("\nGradient test weights\n")
for j=1:maxiter
    CI0.CL_X.CL[1].RB.W1.data = CIini.CL_X.CL[1].RB.W1.data + h*dW
    CI0.C_X.v1.data = CIini.C_X.v1.data + h*dv
    f = loss(CI0, X, Y, A)[1]
    err3[j] = abs(f - f0)
    err4[j] = abs(f - f0 - h*dot(gW, dW) - h*dot(gv, dv))
    print(err3[j], "; ", err4[j], "\n")
    global h = h/2f0
end

@test isapprox(err3[end] / (err3[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err4[end] / (err4[1]/4^(maxiter-1)), 1f0; atol=1f1)


text:
module MaxHelpingHand
# *waves at Communal Helper*
end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
using Printf, LinearAlgebra
import FixedPointNumbers: N0f8
import Images
import CSV, DataFrames
import Cairo

toRadians(α::Float64) = α / 180.0 * π
toDegrees(α::Float64) = α * 180.0 / π

#===========================================
  ____            _   _ _   _ _
 / ___| ___  ___ | | | | |_(_) |___
| |  _ / _ \/ _ \| | | | __| | / __|
| |_| |  __/ (_) | |_| | |_| | \__ \
 \____|\___|\___/ \___/ \__|_|_|___/

=============================================#

struct PositionLLH
    lat::Float64
    lon::Float64
    height::Float64
end

struct PositionXYZ
    x::Float64
    y::Float64
    z::Float64
end

abstract type Ellipsoid
end

struct Wgs84 <: Ellipsoid
    a::Float64
    b::Float64
    f::Float64
    e::Float64
    e²::Float64
    function Wgs84()
        a  = 6378137.0
        f  = 1.0/298.257223563
        b  = a*(1.0-f)
        e² = 1.0 - b*b/(a*a) # or f * (2.0-f)
        e  = sqrt(e²)
        new(a,b,f,e,e²)
    end
end

struct SphericalEarth <: Ellipsoid
    r::Float64
    SphericalEarth() = new(6378137.0)
end

function llh_to_xyz(ellipsoid::Wgs84, llh::PositionLLH)::PositionXYZ
    Φ    = toRadians(llh.lat)
    λ    = toRadians(llh.lon)
    h    = llh.height
    e²   = ellipsoid.e²
    sinΦ = sin(Φ)
    cosΦ = cos(Φ)
    sinλ = sin(λ)
    cosλ = cos(λ)
    v = ellipsoid.a / sqrt(1.0 - e² * sinΦ * sinΦ)
    x = (v+h) * cosΦ * cosλ
    y = (v+h) * cosΦ * sinλ
    z = ((1.0 - e²) * v + h) * sinΦ;
    return PositionXYZ(x,y,z)
end

# TODO: make it work for south and west
function xyz_to_llh(ellipsoid::Wgs84, xyz::PositionXYZ)::PositionLLH
    e²       = ellipsoid.e²
    p        = sqrt(xyz.x*xyz.x + xyz.y*xyz.y)
    lon      = atan(xyz.y / xyz.x)
    lat_init = atan(xyz.z/(p*(1.0 - e²)))
    v        = ellipsoid.a / sqrt(1.0-e²*sin(lat_init)*sin(lat_init))
    lat      = atan((xyz.z + e²*v*sin(lat_init))/p)
    height   = (p/cos(lat))-v
    return PositionLLH(toDegrees(lat), toDegrees(lon), height)
end

function llh_to_xyz(ellipsioid::SphericalEarth, llh::PositionLLH)::PositionXYZ
    Φ    = toRadians(llh.lat)
    λ    = toRadians(llh.lon)
    h    = llh.height
    cosΦ = cos(Φ)
    v    = ellipsioid.r + h
    x    = v * cosΦ * cos(λ)
    y    = v * cosΦ * sin(λ)
    z    = v * sin(Φ)
    return PositionXYZ(x,y,z)
end

function xyz_to_llh(ellipsoid::SphericalEarth, xyz::PositionXYZ)::PositionLLH
    v      = sqrt(xyz.x*xyz.x + xyz.y*xyz.y + xyz.z*xyz.z)
    height = v - ellipsoid.r
    lon    = atan(xyz.y, xyz.x)
    lat    = asin(xyz.z/v)
    return PositionLLH(toDegrees(lat), toDegrees(lon), height)
end

# https://www.movable-type.co.uk/scripts/latlong.html
# for SphericalEarth, but best source found
function bearing(p1::PositionLLH, p2::PositionLLH)::Float64
    φ1 = toRadians(p1.lat)
    φ2 = toRadians(p2.lat)
    Δλ = toRadians(p2.lon-p1.lon)
    y  = sin(Δλ)*cos(φ2)
    x  = cos(φ1)*sin(φ2) - sin(φ1)*cos(φ2)*cos(Δλ)
    θ  = atan(y, x)
    return mod(toDegrees(θ) + 360.0, 360.0)
end

xyz_to_vector(xyz::PositionXYZ)::Vector{Float64} = [xyz.x; xyz.y; xyz.z]

# GeoUtils test
#testXYZ = llh_to_xyz(wgs84, PositionLLH(49.0, 16.0, 225.0))
#println(testXYZ)
#testLLH= xyz_to_llh(wgs84, testXYZ)
#println(testLLH)

#================================================================
 _   _      _       _     _   __  __
| | | | ___(_) __ _| |__ | |_|  \/  | __ _ _ __
| |_| |/ _ \ |/ _` | '_ \| __| |\/| |/ _` | '_ \
|  _  |  __/ | (_| | | | | |_| |  | | (_| | |_) |
|_| |_|\___|_|\__, |_| |_|\__|_|  |_|\__,_| .__/
              |___/                       |_|
===============================================================#

struct LatLonRange
    minLat::Int16
    minLon::Int16
    maxLat::Int16
    maxLon::Int16
end

function loadData(range::LatLonRange, tileDir)
    getHgtFileName(lat, lon) = @sprintf "N%02dE%03d.hgt" lat lon
    getHgtFilePath(lat, lon, tileDir) = @sprintf "%s/%s" tileDir getHgtFileName(lat, lon)

    nTilesHoriz = range.maxLon - range.minLon + 1
    nTilesVert  = range.maxLat - range.minLat + 1
    nTilesTotal = nTilesHoriz * nTilesVert
    dataWidth   = nTilesHoriz*1200+1
    dataHeight  = nTilesVert*1200+1

    @printf("Requesting data for area %d°N %d°E - %d°N %d°E ... (aprox. %3.0fx%3.0f km).\n",
        range.minLat, range.minLon,
        range.maxLat, range.maxLon,
        (nTilesHoriz)*111.1*cos(range.minLat/180.0*π),
        (nTilesVert)*111.1
        )
    @printf("I will read %dx%d=%d tiles, heightmap size is %dx%d (%d MB).\n",
        nTilesHoriz, nTilesVert, nTilesTotal,
        dataWidth, dataHeight, dataWidth*dataHeight*2/1000000
        )

    # Note: array is indexed by row, column and starting index is 1
    data = Array{UInt16}(undef, dataHeight, dataWidth)
    progress = 0
    progressLock = Threads.SpinLock()
    Threads.@threads for i in 0:(nTilesTotal-1)
        lat = range.minLat + div(i, nTilesHoriz)
        lon = range.minLon + mod(i, nTilesHoriz)
        # Print progress
        lock(progressLock) do
            progress = progress + 1
            @printf("Loading tile %03d/%03d lat=%02d, lon=%02d\n", progress, nTilesTotal, lat, lon)
        end
        # Load tile
        tile = Array{Int16}(undef, 1201, 1201)
        path = getHgtFilePath(lat, lon, tileDir)
        io = open(path, "r")
        read!(io, tile)
        close(io)
        # Fix endianity, clamp, transpose and copy to data
        tile .= ntoh.(tile)
        clamp!(tile, 0, 6000)
        tile = transpose(tile)
        rOffset = (range.maxLat-lat) * 1200
        cOffset = (lon-range.minLon) * 1200
        for r in 1:1201
            @simd for c in 1:1201
                data[rOffset+r, cOffset+c] = tile[r, c]
            end
        end
    end
    return data
end

function saveHeightMap(data)
   maxValue = maximum(data)
   norm = Images.Gray.(data/maxValue)
   println("Saving heightmap-gray.png")
   save("heightmap-gray.png", norm)
end

# getHeight(srtmRange, heightMap, 49.142158, 16.627978) -> 192 Svratka, Svitava
function getHeight(range::LatLonRange, data::Matrix{UInt16}, lat::Float64, lon::Float64)
    r = Int64(trunc((range.maxLat+1 - lat)*1200))
    c = Int64(trunc((lon-range.minLon)*1200))
    r = clamp(r, 1, size(data)[1])
    c = clamp(c, 1, size(data)[2])
    return data[r,c]
end

makeEarthCurve(radius, distMax, distStep) = [sqrt(radius*radius-x*x)-radius for x=range(0, distMax, step=distStep)]

struct ViewPort
    ellipsoid::Ellipsoid
    eye::PositionLLH

    angleMin::Float64
    angleMax::Float64
    vertAngleMin::Float64
    vertAngleMax::Float64
    angleStep::Float64

    distMax::Float64
    distStep::Float64

    refractionCoef::Float64

    vUp::Vector{Float64}
    vNorth::Vector{Float64}
    vEast::Vector{Float64}

    outWidth::Int
    outHeight::Int

    function ViewPort(ellipsoid::Ellipsoid, eye::PositionLLH, azimuthMinR, azimuthMaxR, elevationMinR, elevationMaxR, angularStepR, distMaxM, refractionCoef)
        distStep = 50.0
        pRef   = xyz_to_vector(llh_to_xyz(ellipsoid, eye))
        vZ     = [0.0; 0.0; 1.0]
        vUp    = normalize(pRef)
        vEast  = normalize(cross(-vUp,vZ))
        vNorth = normalize(cross(vEast,-vUp))

        azimuthDiff = abs(azimuthMinR-azimuthMaxR)
        if (azimuthDiff > 2.0*π)
            throw(DomainError(azimuthDiff, "Azimuth difference should not exceed 2π"))
        end
        if (azimuthMinR > azimuthMaxR)
            azimuthMinR = azimuthMinR - 2.0*π
        end
        if (azimuthMinR < 0.0 && azimuthMaxR < 0.0)
            azimuthMinR = azimuthMinR + 2.0*π
            azimuthMaxR = azimuthMaxR + 2.0*π
        end
        xMax        = Int64(trunc( (azimuthMaxR-azimuthMinR)/angularStepR ))+1
        outWidth    = xMax+1
        outHeight   = size(range(elevationMinR, elevationMaxR, step=angularStepR))[1]+1
        new(ellipsoid, eye,
            azimuthMinR, azimuthMaxR, elevationMinR, elevationMaxR, angularStepR,
            distMaxM, distStep,
            refractionCoef,
            vUp, vNorth, vEast,
            outWidth, outHeight)
    end
end

function eyeVec(vp::ViewPort)::Vector{Float64}
    eyeXYZ = llh_to_xyz(vp.ellipsoid, vp.eye)
    return [eyeXYZ.x; eyeXYZ.y; eyeXYZ.z]
end

function makeDistMap(vp::ViewPort, latLonRange::LatLonRange, heightMap::Matrix{UInt16})::Matrix{UInt16}
    pRef = eyeVec(vp)
    earthRadius = sqrt(dot(pRef, pRef)) * vp.refractionCoef
    earthCurve  = makeEarthCurve(earthRadius, vp.distMax, vp.distStep)
    @printf("Earth radius is %6.1f km (refraction x%4.2f)\n", earthRadius/vp.refractionCoef/1000.0, vp.refractionCoef)
    @printf("Output size is %d x %d pixels\n", vp.outWidth, vp.outHeight)
    @printf("Output resolution is %f mrad per pixel or %f pixels per degree\n", vp.angleStep * 1000.0, 1.0/toDegrees(vp.angleStep))
    output      = zeros(UInt16, vp.outHeight, vp.outWidth)
    distances   = range(0.0, vp.distMax, step=vp.distStep)
    Threads.@threads for x in 0:(vp.outWidth-1)
        azimuth = vp.angleMin + x*vp.angleStep
        cosAz   = cos(azimuth)
        sinAz   = sin(azimuth)

        vertAngle     = vp.vertAngleMin
        h0            = vp.eye.height
        rayCastHeight = h0
        index         = 1
        direction     = vp.vNorth*cosAz + vp.vEast*sinAz
        point = [0.0,0.0,0.0]
        for dist in distances
            #point = pRef + dist * direction;
            #point = [pRef[1]+dist*direction[1], pRef[2]+dist*direction[2], pRef[3]+dist*direction[3]]
            point[1] = pRef[1]+dist*direction[1]
            point[2] = pRef[2]+dist*direction[2]
            point[3] = pRef[3]+dist*direction[3]
            llh = xyz_to_llh(vp.ellipsoid, PositionXYZ(point[1], point[2], point[3]))
            rayCastHeight = h0 + sin(vertAngle) * dist
            terrainHeight = earthCurve[index] + getHeight(latLonRange, heightMap, llh.lat, llh.lon)
            if terrainHeight > rayCastHeight
                newVertAngle = atan((terrainHeight-h0)/dist)
                yTop = Int64(trunc( (vp.vertAngleMax-newVertAngle)/vp.angleStep ))
                yBot = Int64(trunc( (vp.vertAngleMax-   vertAngle)/vp.angleStep ))
                v = UInt16(trunc( dist / vp.distStep ))
                for y in yTop:yBot
                    output[y, x+1] = v
                end
                vertAngle = newVertAngle
            end
            index = index + 1
        end
    end
    return output
end

function extractOutlines(distMap::Matrix{UInt16})
    nRows = size(distMap)[1]
    nCols = size(distMap)[2]
    #ouput = Array{Int16}(undef, nRows, nCols)
    output = zeros(N0f8, nRows, nCols)
    Threads.@threads for row in 2:nRows
        for col in 1:nCols
            diff::UInt8 = 255-clamp(abs(reinterpret(Int16,distMap[row-1,col]) - reinterpret(Int16,distMap[row,col])), 0, 255)
            output[row, col] = reinterpret(N0f8, diff)
        end
    end
    return output
end

function testPixel(distMap::Matrix{UInt16}, x::UInt64, y::UInt64, radius::Int64, value::UInt16, valueTolerance::UInt16)::Bool
    if (x < radius+1) || (y < radius+1) || (x+radius> size(distMap)[2]) || (y+radius> size(distMap)[1])
        return false
    end
    for row in (y-radius):(y+radius)
        for col in (x-radius):(x+radius)
            mapValue = distMap[row, col]
            if (mapValue > value) && ((mapValue - value) ≤ valueTolerance)
                return true
            end
            if (mapValue < value) && ((value - mapValue) ≤ valueTolerance)
                return true
            end
        end
    end
    return false
end

function drawSummits(vp::ViewPort, distMap::Matrix{UInt16})
    println("Loading data")
    dfFiltered = DataFrames.DataFrame(Summit = String[], Elevation = Float64[], Distance=Float64[], X=UInt64[], Y=UInt64[])

    # TODO: options
    #hillsCZ = CSV.File("data-cz-prom100.tsv") |> DataFrames.DataFrame
    #hillsSK = CSV.File("data-sk-prom200.tsv") |> DataFrames.DataFrame
    #hills = vcat(hillsCZ, hillsSK)
    hills =  CSV.File("osm-cz-sk.tsv") |> DataFrames.DataFrame
    # convert to ours azimuth, angle above horizon and distance - project into
    hill_to_xyz(ellipsoid::Ellipsoid, dfRow)::PositionXYZ = llh_to_xyz(ellipsoid, PositionLLH(dfRow["Latitude"], dfRow["Longitude"], dfRow["Elevation"]))
    # difference between true and seen earth curvature
    elevationDropAtDistance(distance::Float64, radius::Float64)::Float64 = sqrt(radius*radius-distance*distance)-radius
    elevationDropCompensation(distance::Float64, radius::Float64, refractionCoef::Float64)::Float64 = elevationDropAtDistance(distance, radius*refractionCoef) - elevationDropAtDistance(distance, radius)
    # TODO: make distance function
    mLocalToWorld = hcat(vp.vEast,vp.vNorth,vp.vUp)
    mWorldToLocal = inv(mLocalToWorld)
    #ground = PositionLLH(vp.eye.lat, vp.eye.lon, 0.0)
    pRef = xyz_to_vector(llh_to_xyz(vp.ellipsoid, vp.eye))
    earthRadius = sqrt(dot(pRef, pRef))
    for hill in eachrow(hills)
        hill_world = hill_to_xyz(vp.ellipsoid, hill)
        hill_local_xyz = mWorldToLocal * [hill_world.x; hill_world.y; hill_world.z]
#        hill_local_xyz[3] = hill_local_xyz[3] - earthRadius  <---- this gives weird altitutude, let's restore it in the same way as raytracer
        hill_local_xyz[3] = 0.0
        distance = sqrt(dot(hill_local_xyz, hill_local_xyz))
        if distance > vp.distMax
            continue
        end
        hill_local_xyz[3] = hill["Elevation"] + elevationDropAtDistance(distance, earthRadius * vp.refractionCoef) - vp.eye.height
        azimuth = bearing(vp.eye, PositionLLH(hill["Latitude"], hill["Longitude"], 0.0))
        if (azimuth < toDegrees(vp.angleMin) || azimuth > toDegrees(vp.angleMax)) # FIXME: test for weird angles
            continue
        end
        @printf("%20s is possibly visible at azimuth %6.2f, distance %6.2f km", hill["Summit"], azimuth, distance/1000.0)
        elevationAngle=atan(hill_local_xyz[3], distance)
#        @printf("              hill=%f+%f, dist=%f\n", hill_local_xyz[3], elevationDropCompensation(distance, earthRadius, vp.refractionCoef), distance)
#        @printf("              , pixel.x,y=%5.0f,%5.0f\n", (toRadians(azimuth)-vp.angleMin)/vp.angleStep , (vp.vertAngleMax-elevationAngle)/vp.angleStep )
        testX::UInt64 = round((toRadians(azimuth)-vp.angleMin)/vp.angleStep)
        testY::UInt64 = round((vp.vertAngleMax-elevationAngle)/vp.angleStep)
        visible::Bool = testPixel(distMap, testX, testY, 4, UInt16(trunc(distance/vp.distStep)), UInt16(5))
        @printf(", visible=%s\n", visible)
        if visible
            push!(dfFiltered, (hill["Summit"], hill["Elevation"], distance, testX, testY))
        end
    end
    println("Drawing ...")
    # Initialize
    Cairo_set_line_color(ctx::Cairo.CairoContext) = Cairo.set_source_rgb(ctx, 131/255, 148/255, 150/255)
    Cairo_set_text_color(ctx::Cairo.CairoContext) = Cairo.set_source_rgb(ctx,  38/255, 139/255, 210/255)
    function Cairo_line(ctx::Cairo.CairoContext, x1::Core.Real, y1::Core.Real, x2::Core.Real, y2::Core.Real)
        Cairo.move_to(ctx, x1, y1)
        Cairo.line_to(ctx, x2, y2)
        Cairo.stroke(ctx)
    end
    surf = Cairo.CairoARGBSurface(vp.outWidth, vp.outHeight)
    ctx  = Cairo.CairoContext(surf)
    # Background (previous image)
    bgSurf = Cairo.read_from_png("outlines.png")
    Cairo.set_source_surface(ctx, bgSurf, 0.0, 0.0)
    Cairo.paint(ctx)
    # Annotations
    Cairo.select_font_face(ctx, "Fira Sans", Cairo.FONT_SLANT_NORMAL, Cairo.FONT_WEIGHT_NORMAL)
    Cairo.set_font_size(ctx, 18.0)
    Cairo.set_line_width(ctx, 1.0)
    for poi in eachrow(dfFiltered)
        Cairo_set_line_color(ctx)
        Cairo_line(ctx, poi["X"]+0.5, poi["Y"], poi["X"]+0.5, 300.0)
        Cairo_set_text_color(ctx)
        Cairo.move_to(ctx, poi["X"]+5, 300.0-5.0)
        Cairo.save(ctx)
        Cairo.rotate(ctx, toRadians(-45.0))
        Cairo.show_text(ctx, poi["Summit"])
        Cairo.show_text(ctx, @sprintf(" (%3.0f km)", poi["Distance"]/1000.0))
        Cairo.restore(ctx)
    end
    # Azimuth ticks
    azMinD::Int = Int(ceil(toDegrees(vp.angleMin)))
    azMaxD::Int = Int(floor(toDegrees(vp.angleMax)))
    for az in azMinD:azMaxD
        x = round((toRadians(Float64(az))-vp.angleMin)/vp.angleStep)+0.5
        Cairo_set_line_color(ctx)
        Cairo_line(ctx, x, 38, x, 42)
        Cairo_line(ctx, x, 63, x, 68)
        Cairo_set_text_color(ctx)
        ext = Cairo.text_extents(ctx, @sprintf("%d", az))
        Cairo.move_to(ctx, x - ext[3]/2, 58)
        Cairo.show_text(ctx, @sprintf("%d °", az))
    end
    # Horizon line
    horizY = round(vp.vertAngleMax/vp.angleStep)+0.5
    Cairo.set_source_rgb(ctx, 238/255, 232/255, 213/255)
    Cairo_line(ctx, 0.0, horizY, vp.outWidth, horizY)
    # Finalize
    Cairo.write_to_png(surf, "outline-with-annotations.png" )
end

# Info (from https://www.udeuschle.de/panoramas/makepanoramas_en.htm)
# Lat: 50.08309 Lon 17.23094 Alt(auto+10m): 1500+10
# View direction: 112.5, extension 45 left: 90, right: 135, resolution 20pix/deg
# Tilt, range, vert. exaggeration 1.2

# TODO: determine lat/long range automatically
function main()
    tileDir     = "d:/_disk_d_old/devel-python/panorama/data_srtm"
	#tileDir     = "data_srtm"
    latLonRange = LatLonRange(47, 15, 50, 21)
    eye         = PositionLLH(50.08309, 17.23094, 1510)

    heightMap = loadData(latLonRange, tileDir)
    #saveHeightMap(data)
    ellipsoid = SphericalEarth()
    #ellipsoid = Wgs84()
    vp = ViewPort(ellipsoid, eye, toRadians(90.0), toRadians(135.0), -0.0560, 0.0339, 0.0001, 250.0e3, 1.18)
    distMap   = makeDistMap(vp, latLonRange, heightMap)

    minValue = minimum(distMap)
    maxValue = maximum(distMap)
    @printf("min=%d max=%d\n", minValue, maxValue)
    println("Saving distmap-gray.png")
    Images.save("distmap-gray.png", Images.Gray.(distMap/maxValue))

    println("Extracting outlines ...")
    outlines=extractOutlines(distMap)
    println("Saving outlines.png")
    Images.save("outlines.png", outlines)

    println("Creating annotations")
    drawSummits(vp, distMap)
    println("All done")

    # This can be fun: https://wiki.flightgear.org/Atmospheric_light_scattering
    # http://www.science-and-fiction.org/rendering/als.html
end

main()

max_stars_repo_name: pavel-perina/panorama-jl

max_issues_repo_name: pavel-perina/panorama-jl

max_forks_repo_name: pavel-perina/panorama-jl


text:
# -*- coding: utf-8 -*-
using Test
using Polynomials
using DifferentialEquations
import DifferentialEquations: solve

# f = @ode_def LotkaVolterra begin
#   dx = a*x - b*x*y
#   dy = -c*y + d*x*y
# end a b c d

"""   hookeslaw

the ODE representation of a spring. The solutions are periodic functions.
"""
function hookeslaw(du, u, p, t)
    # u[1] = p
    # u[2] = v
    du[1] = u[2]
    du[2] = -p[1] * u[1]
end

"""    SpringModel

represents the second order linear ODE goverened by [`hookeslaw`](@ref).

"""
mutable struct SpringModel{F,D,U}
    frequency::F
    domain::D
    initial_position::U
end

parameters(spm::SpringModel) = spm.frequency
domain(spm::SpringModel) = spm.domain
flux(spm::SpringModel) = hookeslaw
initial_conditions(spm::SpringModel) = spm.initial_position
odeproblem(spm::SpringModel) = ODEProblem(flux(spm), initial_conditions(spm), domain(spm), parameters(spm))
solve(spm::SpringModel, alg=Vern7()) = solve(odeproblem(spm), alg)
accelerate!(spm::SpringModel, factor::Number) = begin
    spm.frequency = factor*spm.frequency
end

struct SIRParams{T,U}
    β::T
    γ::U
end

struct SIRSimulation{V, T, P}
    initial_populations::V
    tspan::T
    params::P
end

parameters(sir::SIRSimulation) = sir.params
domain(sir::SIRSimulation) = sir.tspan
initial_conditions(sir::SIRSimulation) = sir.initial_populations

function flux(sir::SIRSimulation)
    function sirflux(du, u, p, t)
        # println("sirflux call")
        # @show t
        β, γ = p.β, p.γ
        # println(β, γ)
        N = sum(u)
        # println(N)
        infections = β * (u[1]*u[2]/N)
        # println(infections)
        recoveries = γ * (u[2]/N)
        # println(recoveries)
        du[1] = -infections
        du[2] =  infections - recoveries
        du[3] =  recoveries
        # println(du)
        # println(u)
        return du
    end
    return sirflux
end

function odeproblem(sir::SIRSimulation)
    problem = ODEProblem(flux(sir), initial_conditions(sir), domain(sir), parameters(sir))
    return problem
end

solve(sir::SIRSimulation) = solve(odeproblem(sir), alg=Vern9())
function RealSIR()
    # β = NumParameter(:β, u"person/s", TransitionRate)
    # γ = NumParameter(:γ, u"person/s", TransitionRate)
    # @show β, γ

    # S = NumVariable(:S, u"person", Amount)
    # I = NumVariable(:I, u"person", Amount)
    # R = NumVariable(:R, u"person", Amount)

    # sir = SIR([ODE(:(dS/dt -> -βSI/N),
    #             [β], [S]),
    #             ODE(:(dI/dt -> βSI/N -γI),
    #                 [β, γ], [S, I]),
    #         ODE(:(dR/dt -> γI), [γ], [I])],
    #         NumVariable.([:N,:S,:I,:R], u"person", Amount),
    #         NumVariable(:t, u"s", Amount))

    return odeproblem(SIRSimulation([100,1,0.0],
                      (0,100.0),
                       SIRParams(0.50,5.0)))
end
# test case where nothing happens because of *0.0
    # return odeproblem(SIRSimulation([100,0,0.0],
    #                   (0,10.0),
    #                    SIRParams(1.0,0.2)))

# this works at this point the next step is to run it with units
# β in person^-2/day
# γ in person/day

# then use the number class to encode more information

springmodel = SpringModel([1.0], (0,4π), [1.0, 0.0])
solslow = solve(odeproblem(springmodel))
@test abs(solslow(π/2)[1]) < 1e-6
@test abs(solslow(3π/2)[1]) < 1e-4

accelerate!(springmodel, 4.0)
solfast = solve(odeproblem(springmodel), atol=1e-8)
@test abs(solfast(π/2)[1]+1) < 1e-5
@test_skip abs(solfast(3π/2)[1]+1) < 1e-4
@test abs(solfast(π/4)[1]) < 1e-5
@test abs(solfast(3π/4)[1]) < 1e-4


text:
using Eirene
using Plots
include("clique_top_Julia/CliqueTop.jl")
include("VideoProcessing.jl")
include("MatrixToolbox.jl")
include("Settings.jl")

function testing_pariwise_corr()
    do_clique_top = test_params["do_clique_top"]
    do_eirene =     test_params["do_eirene"]
    save_figures = test_params["save_figures"]
    plot_betti_figrues = test_params["plot_betti_figrues"]
    plot_vectorized_video = test_params["plot_vectorized_video"]
    size_limiter = test_params["size_limiter"]
    ind_distrib = test_params["ind_distrib"]
    videos_set = test_params["videos_set"]
    tau_max_set = test_params["tau_max_set"]
    points_per_dim_set = test_params["points_per_dim_set"]
    shifts_set = test_params["shifts_set"]
    patch_params = test_params["patch_params"]
    video_path = test_params["video_path"]
    results_path = test_params["results_path"]
    videos = test_params["videos_names"]

    do_local_corr = false
    do_local_grad = false

    if ind_distrib == "local_corr"
        shift_set = test_params["shift_set"]
        sub_img_size_set = [9]
        do_local_corr = true
        do_local_grad = false
        @debug "Doing local correlation" do_local_corr
    elseif ind_distrib == "local_grad"
        shift_set = [1]
        sub_img_size_set = test_params["sub_img_size_set"]
        do_local_corr = false
        do_local_grad = true
        @debug "Doing local gradient" do_local_grad
    else
        shift_set = [1]
        sub_img_size_set = [9]
        do_local_corr = false
        do_local_grad = false
    end
    @info "Using following distribution: " test_params["ind_distrib"]

    @debug "All videos are: " videos_names
    @debug "Video set is : " videos_set
    for video in videos_set
        choice = videos_names[video]
        @info "Selected video: " choice
        @debug "Path and choice is:" video_path*choice

        video_array = get_video_array_from_file(video_path*choice)
        @info "Array extracted."

        video_dimensions = get_video_dimension(video_array)
        for points_per_dim in points_per_dim_set
            for shift in shift_set, sub_img_size in sub_img_size_set
                if do_local_corr
                    centers = get_local_centers(points_per_dim, video_dimensions, shift, sub_img_size)

                    extracted_pixels_matrix = get_subimg_correlations(video_array, centers, sub_img_size, shift)
                elseif do_local_grad
                    centers = get_local_centers(points_per_dim, video_dimensions, shift, sub_img_size)

                    extracted_pixels_matrix = get_local_gradients(video_array, centers, sub_img_size)
                else
                    indicies_set = get_video_mask(points_per_dim, video_dimensions,  distribution=ind_distrib, patch_params)
                    extracted_pixels_matrix = extract_pixels_from_video(video_array, indicies_set, video_dimensions)
                end
                @info "Pixels extracted."

                vectorized_video = vectorize_video(extracted_pixels_matrix)
                @info "Video is vectorized, proceeding to Pairwise correlation."

                for tau in tau_max_set
                    ## Compute pairwise correlation
                    C_ij = get_pairwise_correlation_matrix(vectorized_video, tau)

                    # set the diagonal to zero
                    for diag_elem in 1:size(C_ij,1)
                        C_ij[diag_elem,diag_elem] = 0
                    end
                    @info "Pairwise correlation finished, proceeding to persistance homology."

                    # Compute persistance homology with CliqueTopJulia
                    size_limiter = test_params["size_limiter"]
                    @debug "using size limiter = " size_limiter

                    if size_limiter > size(C_ij,1)
                        @warn "Used size limiter is larger than matrix dimension: " size_limiter size(C_ij,1)
                        @warn "Using maximal size instead"
                        size_limiter = size(C_ij,1)
                    end

                    @debug "do_clique_top: " do_clique_top
                    @debug "test_params['do_clique_top']: " test_params["do_clique_top"]
                    if do_clique_top
                        @debug pwd()
                        @time c_ij_betti_num, edge_density, persistence_intervals, unbounded_intervals = compute_clique_topology(C_ij[1:size_limiter, 1:size_limiter], edgeDensity = 0.6)
                    end

                    @debug "do_eirene: " do_eirene
                    if do_eirene
                        C = eirene(C_ij[1:size_limiter, 1:size_limiter],maxdim=3,model="vr")
                    end

                    # ---------------------------------------------------------
                    # Plot results
                    @debug "Proceeding to plotting."
                    if plot_vectorized_video
                        vector_plot_ref = heatmap(vectorized_video, color=:grays)
                        if save_figures
                            name = split(choice, ".")[1]
                            name = "vec_" * name * "_sz$(size_limiter)_p$(points_per_dim)_tau$(tau).png"
                            savefig(vector_plot_ref, name)
                        end #save vec
                    end #plot vec

                    if plot_betti_figrues && do_clique_top
                        betti_numbers = c_ij_betti_num
                        title = "Betti curves for pairwise corr. matrix"
                        p1 = plot_betti_numbers(c_ij_betti_num, edge_density, title);

                        heat_map1 = heatmap(C_ij,  color=:lightrainbow, title="Pariwise Correlation matrix, number of points: $(points_per_dim)");

                        betti_plot_clq_ref = plot(p1, heat_map1, layout = (2,1))

                        if save_figures
                            saving_figures(betti_plot_clq_ref, results_cliq_path, choice, points_per_dim, tau, size_limiter)
                        end#save fig
                    end #plot cliq

                    if plot_betti_figrues && do_eirene
                        p1, heat_map1 = plot_eirene_betti_curves(C, C_ij)
                        betti_plot_ei_ref = plot(p1, heat_map1, layout = (2,1))

                        if save_figures
                            saving_figures(betti_plot_ei_ref, results_cliq_path, choice, points_per_dim, tau, size_limiter)
                        end#save fig
                    end #plot eirene
                end #for tau
            end #for shift
        end #for points_per_dim
    end #for video set
    @info "Finished testing"
end #func


text:
# Modified from Flux.jl to accept the td error instead of two inputs
function Flux.huber_loss(y; δ=eltype(y)(1))
    abs_error = abs.(y)
    temp = abs_error .< δ
    x = eltype(y)(0.5)
    hub_loss = sum(((abs_error.^2) .* temp) .* x .+ δ*(abs_error .- x*δ) .* (1 .- temp)) * 1 // length(y)
end


text:
matchesprefix(c::AString, pre::AString) = isempty(pre) || lowercase(c[1]) == lowercase(pre[1])
matchesprefix(c::Dict, pre::AString) = matchesprefix(c[:text], pre)
matchesprefix(c, ::Void) = true

handle("completions") do data
  @destruct [path || nothing, mod || "Main", line, column, force] = data
  withpath(path) do
    pre = CodeTools.prefix(line[1:column-1])
    pre = isempty(pre) ? nothing : pre[end]
    cs = CodeTools.completions(line[1:column-1], getthing(mod), default = false)
    cs == nothing && pre == nothing && !force && (cs = [])
    d(:completions => cs,
      :prefix      => pre,
      :mod         => mod)
  end
end

handle("cacheCompletions") do mod
  CodeTools.completions(getthing(mod))
end


text:
using Distributions
using ConjugatePriors: NormalInverseWishart, posterior_canon
using Particles: marginal_lhood, marginal_log_lhood

ConjugatePriors.NormalInverseWishart(nix2::NormalInverseChisq) =
    NormalInverseWishart([nix2.μ], nix2.κ, nix2.ν*reshape([nix2.σ2], 1, 1), nix2.ν)

@testset "NormalInverseWishart" begin

    @testset "Posterior predictive" begin
        d = NormalInverseWishart([0., 1.], 10., [4. 2.; 2. 3.], 11.)
        d_pp_samps = zeros(2, 1_000_000)
        for i in 1:1_000_000
            d_pp_samps[:,i] = rand(MvNormal(rand(d)...))
        end
        pp_samps = rand(posterior_predictive(d), 1_000_000)

        @test isapprox(cov(d_pp_samps'), cov(pp_samps'), atol=0.01)
        @test isapprox(mean(d_pp_samps, dims=2), mean(pp_samps, dims=2), atol=0.01)
    end

    @testset "Marginal likelihood" begin
        Random.seed!(1001)

        prior = NormalInverseWishart([0., 1.], 10., Matrix(1.0I, 2,2), 50.)
        x = rand(MvNormal(rand(prior)...), 10)
        ss = suffstats(MvNormal, x)
        # this is trivially true at the moment since one just calls the other
        @test marginal_log_lhood(prior, ss) ≈ log(marginal_lhood(prior, ss))
        xvecs = [x[:,i] for i in 1:size(x,2)]
        logpdfs = [sum(logpdf.(Ref(MvNormal(prior_samp...)), xvecs))
                   for prior_samp
                   in (rand(prior) for _ in 1:1000)]
        @test isapprox(marginal_log_lhood(prior, ss), log(mean(exp.(logpdfs))), rtol=0.1)
        # sanity check: lhood under samples from prior is more variable than
        # error from analytical expression
        @test std(logpdfs) > abs(marginal_log_lhood(prior, ss) - log(mean(exp.(logpdfs))))
    end

    # I think the parameter space is too high-dimensional to have any hope of
    # testing reliably with sampling like this.  Next best is to check against
    # the 1-D case which is equivalent to the Normal-Inverse Chi-squared
    @testset "Equivalence with Normal-Inverse Chi-squared" begin
        nix2 = NormalInverseChisq(1., 2., 3., 4.)
        niw = NormalInverseWishart(nix2)

        x = rand(Normal(3, 2), 100)

        ss = suffstats(Normal, x)
        ss_mv = suffstats(MvNormal, reshape(x, 1, :))

        post_nix2 = posterior_canon(nix2, ss)
        post_niw = posterior_canon(niw, ss_mv)

        @test all(post_nix2.μ .≈ post_niw.mu)
        @test post_nix2.κ ≈ post_niw.kappa
        @test post_nix2.ν ≈ post_niw.nu
        @test all(post_nix2.σ2 .≈ Matrix(post_niw.Lamchol)[1] ./ post_niw.nu)

        μ, σ2 = rand(post_nix2)
        @test logpdf(post_nix2, μ, σ2) ≈ logpdf(post_niw, [μ], reshape([σ2], 1, 1))
        @test logpdf(nix2, μ, σ2) ≈ logpdf(niw, [μ], reshape([σ2], 1, 1))

        @test marginal_log_lhood(nix2, ss) ≈ marginal_log_lhood(niw, ss_mv)

    end
end


text:
function shownode(io::IO, node, depth::Int, item_prefix::String, prefix::String)
    buf = PipeBuffer()
    printnode(buf, node)
    for (i, line) in enumerate(eachline(buf, keep=true))
        if i == 1
            print(io, item_prefix*line)
        else
            print(io, prefix*line)
        end
    end
    if depth <= 0
        println(io, " ($(n_children(node)) children)")
    else
        println(io)
        n = n_children(node)
        for (i, c) in enumerate(children(node))
            if i < n
                shownode(io, c, depth-1, prefix*"├──", prefix*"│  ")
            else # last one
                shownode(io, c, depth-1, prefix*"└──", prefix*"   ")
            end
        end
    end
end


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

"""
    AdaMax(; kwargs...)

This is a variant of of the Adam algorithm based on the infinity norm.
See [1] for further description.

### Arguments
* `η`: default `0.002`, learning rate.
* `β1`: default `0.9`, exponential decay rate for the first moment estimates.
* `β2`: default `0.999`, exponential decay rate for the weighted
  infinity norm estimates.
* `ϵ`: default `1e-8`, small value added for numerical stability.
* `clip`: default `0`, gradient clipping.
  If positive, will clip the gradient into the range `[-clip, clip]`.
* `scale`: default `0`, gradient rescaling.
  If != 0, multiply the gradient with `scale` before updating.
  Often choose to be `1.0 / batch_size`.
  If leave it default, high-level API like `fit!` will set it to
  `1.0 / batch_size`, since `fit!` knows the `batch_size`.
* `λ`: default `0.00001`, weight decay is equivalent
  to adding a global l2 regularizer for all the parameters.

### References
1. Kingma, Diederik, and Jimmy Ba (2014):
   Adam: A Method for Stochastic Optimization. Section 7.
   [http://arxiv.org/abs/1412.6980]
   (http://arxiv.org/abs/1412.6980).
"""
AdaMax

@defstruct AdaMax <: AbstractOptimizer (
  (η      :: Real = 0.002, η > 0),
  (β1     :: Real = 0.9,   0 <= β1 < 1),
  (β2     :: Real = 0.999, 0 <= β2 < 1),
  (ϵ      :: Real = 1e-8,  ϵ > 0),
  (clip   :: Real = 0,     clip >= 0),
   scale  :: Real = 0,
  (λ      :: Real = 1e-5,  λ >= 0),
  η_sched :: Any  = initlrsched(η)
)

mutable struct AdaMaxState
  mₜ  :: NDArray
  uₜ  :: NDArray
  β1ᵗ :: Float64
end

create_state(ada::AdaMax, ::Int, W::NDArray) =
  AdaMaxState(zeros(size(W), context(W)),
              zeros(size(W), context(W)),
              ada.β1)

function update!(ada::AdaMax, ::Int, W::NDArray, ∇::NDArray, s::AdaMaxState)
  η = get(ada.η_sched)
  β1 = ada.β1
  β2 = ada.β2
  ϵ = ada.ϵ

  normgrad!(ada, W, ∇)

  s.mₜ = β1 * s.mₜ .+ (1 - β1) .* ∇
  s.uₜ = _maximum(β2 * s.uₜ, abs(∇))  # FIXME abs dot-call

  @inplace W .+= -η / (1 - s.β1ᵗ) * s.mₜ ./ (s.uₜ + ϵ)

  s.β1ᵗ *= ada.β1
end

max_stars_repo_name: Vikas-kum/incubator-mxnet

max_issues_repo_name: Vikas-kum/incubator-mxnet

max_forks_repo_name: Vikas-kum/incubator-mxnet


text:
module Utilities
using ..Concepts

import Distributions
import LinearAlgebra

# export ErrorMatrix,
#     relative_error,
#     total_error

include("./Misc.jl")
include("./Diagnostics.jl")
include("./ExponentialFamily.jl")
include("./FastEigen.jl")
include("./RandomMatrices.jl")
include("./Indexing.jl")
include("./Sampling.jl")
include("./BatchUtils.jl")
include("./PrettyPrinter.jl")

# include("./TestModule.jl")

end


text:
using PyCall
using HDF5
using ImpulseResponse

@pyimport tkinter as tk

speakermix2mat(x) = permutedims([parse(Int,i) for i in split(x[2:end-1])])
microphonemix2mat(x) = [parse(Int,i) for i in split(x[2:end-1],',')]

f = Any
m = Any
ϕ = Any
xds = Any
h = Any

function measureclicked()
end
function calculateclicked()
end
function validateclicked()
end
function saveclicked()
end

@pydef mutable struct SampleApp <: tk.Tk
    __init__(self, args...; kwargs...) = begin
        tk.Tk[:__init__](self, args...; kwargs...)

        # label frame parameters
        self[:lf0] = tk.LabelFrame(text="parameters")
        self[:lf0][:pack](fill="both", expand="yes")

        self[:lab00] = tk.Label(self[:lf0], text="sample rate")
        self[:lab00][:pack](side="left")
        self[:ent00v] = tk.StringVar(value="48000")
        self[:ent00] = tk.Entry(self[:lf0], width=10, textvariable=self[:ent00v])
        self[:ent00][:pack](side="left")

        self[:lab01] = tk.Label(self[:lf0], text="nfft/2")
        self[:lab01][:pack](side="left")
        self[:ent01v] = tk.StringVar(value="1024")
        self[:ent01] = tk.Entry(self[:lf0], width=10, textvariable=self[:ent01v])
        self[:ent01][:pack](side="left")

        self[:lab02] = tk.Label(self[:lf0], text="attenuate(dB)")
        self[:lab02][:pack](side="left")
        self[:ent02v] = tk.StringVar(value="-20")
        self[:ent02] = tk.Entry(self[:lf0], width=10, textvariable=self[:ent02v])
        self[:ent02][:pack](side="left")

        self[:lab03] = tk.Label(self[:lf0], text="speaker mix")
        self[:lab03][:pack](side="left")
        self[:ent03v] = tk.StringVar(value="[0 0 0 0 1 0 0 0]")
        self[:ent03] = tk.Entry(self[:lf0], width=20, textvariable=self[:ent03v])
        self[:ent03][:pack](side="left")

        self[:lab04] = tk.Label(self[:lf0], text="microphone mix")
        self[:lab04][:pack](side="left")
        self[:ent04v] = tk.StringVar(value="[1,0,0,0,0,0,0,0]")
        self[:ent04] = tk.Entry(self[:lf0], width=20, textvariable=self[:ent04v])
        self[:ent04][:pack](side="left")

        # label farme impulse response
        self[:lf1] = tk.LabelFrame(text="impulse response of system WITHOUT equalization")
        self[:lf1][:pack](fill="both", expand="yes")
        self[:lf2] = tk.LabelFrame(text="minimal phase fir equalization filter settings")
        self[:lf2][:pack](fill="both", expand="yes")
        self[:lf3] = tk.LabelFrame(text="impulse response of system WITH equalization")
        self[:lf3][:pack](fill="both", expand="yes")

        self[:lab1] = tk.Label(self[:lf1], text="start(Hz)")
        self[:lab1][:pack](side="left")
        self[:ent1v] = tk.StringVar(value="22")
        self[:ent1] = tk.Entry(self[:lf1], width=10, textvariable=self[:ent1v])
        self[:ent1][:pack](side="left")

        self[:lab2] = tk.Label(self[:lf1], text="stop(Hz)")
        self[:lab2][:pack](side="left")
        self[:ent2v] = tk.StringVar(value="22000")
        self[:ent2] = tk.Entry(self[:lf1], width=10, textvariable=self[:ent2v])
        self[:ent2][:pack](side="left")

        self[:btn1] = tk.Button(self[:lf1], text="measure", command=measureclicked)
        self[:btn1][:pack](side="left")

        # label frame filter design
        self[:lab3] = tk.Label(self[:lf2], text="start(Hz)")
        self[:lab3][:pack](side="left")
        self[:ent3v] = tk.StringVar(value="142")
        self[:ent3] = tk.Entry(self[:lf2], width=10, textvariable=self[:ent3v])
        self[:ent3][:pack](side="left")

        self[:lab4] = tk.Label(self[:lf2], text="anchor(Hz)")
        self[:lab4][:pack](side="left")
        self[:ent4v] = tk.StringVar(value="214")
        self[:ent4] = tk.Entry(self[:lf2], width=10, textvariable=self[:ent4v])
        self[:ent4][:pack](side="left")

        self[:lab5] = tk.Label(self[:lf2], text="stop(Hz)")
        self[:lab5][:pack](side="left")
        self[:ent5v] = tk.StringVar(value="13564")
        self[:ent5] = tk.Entry(self[:lf2], width=10, textvariable=self[:ent5v])
        self[:ent5][:pack](side="left")

        self[:btn2] = tk.Button(self[:lf2], text="calculate", command=calculateclicked)
        self[:btn2][:pack](side="left")

        # label frame filter verification
        self[:lab6] = tk.Label(self[:lf3], text="start(Hz)")
        self[:lab6][:pack](side="left")
        self[:ent6v] = tk.StringVar(value="22")
        self[:ent6] = tk.Entry(self[:lf3], width=10, textvariable=self[:ent6v])
        self[:ent6][:pack](side="left")

        self[:lab7] = tk.Label(self[:lf3], text="stop(Hz)")
        self[:lab7][:pack](side="left")
        self[:ent7v] = tk.StringVar(value="22000")
        self[:ent7] = tk.Entry(self[:lf3], width=10, textvariable=self[:ent7v])
        self[:ent7][:pack](side="left")

        self[:btn3] = tk.Button(self[:lf3], text="validate", command=validateclicked)
        self[:btn3][:pack](side="left")

        # label frame save to file
        self[:lf4] = tk.LabelFrame(text="save filter to file")
        self[:lf4][:pack](fill="both", expand="yes")

        self[:lab8] = tk.Label(self[:lf4], text="file name")
        self[:lab8][:pack](side="left")
        self[:ent8v] = tk.StringVar(value="euqalization-20181122.h5")
        self[:ent8] = tk.Entry(self[:lf4], width=40, textvariable=self[:ent8v])
        self[:ent8][:pack](side="left")

        self[:lab9] = tk.Label(self[:lf4], text="dataset name")
        self[:lab9][:pack](side="left")
        self[:ent9v] = tk.StringVar(value="mth1")
        self[:ent9] = tk.Entry(self[:lf4], width=20, textvariable=self[:ent9v])
        self[:ent9][:pack](side="left")

        self[:btn4] = tk.Button(self[:lf4], text="save", command=saveclicked)
        self[:btn4][:pack](side="left")

    end
end
app = SampleApp()

function measureclicked()
    ms = speakermix2mat(app[:ent03][:get]())
    mm = microphonemix2mat(app[:ent04][:get]())
    mm = mm[:,:]
    fs = parse(Int, app[:ent00][:get]())
    n = parse(Int, app[:ent01][:get]())
    atten = parse(Float64, app[:ent02][:get]())

    f0 = parse(Float64, app[:ent1][:get]())
    f1 = parse(Float64, app[:ent2][:get]())

    @info "measurement:" ms mm fs n atten f0 f1
    f_, m_, ϕ_, xds_ = ImpulseResponse.minimalphase_fir_original(ms, mm, fs, f0, f1, atten, n)
    global f = f_
    global m = m_
    global ϕ = ϕ_
    global xds = xds_
    nothing
end

function calculateclicked()
    ms = speakermix2mat(app[:ent03][:get]())
    mm = microphonemix2mat(app[:ent04][:get]())
    mm = mm[:,:]
    fs = parse(Int, app[:ent00][:get]())
    n = parse(Int, app[:ent01][:get]())
    atten = parse(Float64, app[:ent02][:get]())

    f0 = parse(Float64, app[:ent3][:get]())
    fx = parse(Float64, app[:ent4][:get]())
    f1 = parse(Float64, app[:ent5][:get]())

    @info "calculate:" ms mm fs n atten f0 fx f1
    h_ = ImpulseResponse.minimalphase_fir_design(f, m, ϕ, xds, fs, f0, f1, fx, n)
    global h = h_
    nothing
end

function validateclicked()
    ms = speakermix2mat(app[:ent03][:get]())
    mm = microphonemix2mat(app[:ent04][:get]())
    mm = mm[:,:]
    fs = parse(Int, app[:ent00][:get]())
    n = parse(Int, app[:ent01][:get]())
    atten = parse(Float64, app[:ent02][:get]())

    f0 = parse(Float64, app[:ent6][:get]())
    f1 = parse(Float64, app[:ent7][:get]())

    @info "validate:" ms mm fs n atten f0 f1
    ImpulseResponse.minimalphase_fir_verification(h, ms, mm, fs, f0, f1, atten, n)
end

function saveclicked()
    @info "save to file" h
    h5write(app[:ent8][:get](), app[:ent9][:get](), convert(Matrix{Float64},h[:,:]))
end

app[:mainloop]()


text:
# This is a test for CmdStan issue #510 (thanks to jonalm)

using Mamba, Stan

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

  simplecode = "
  data {real sigma;}
  parameters {real y;}
  model {y ~ normal(0,sigma);}
  "

  global stanmodel1, rc1, sim1
  stanmodel1 = Stanmodel(Sample(save_warmup=false, thin=5), name="simple", model=simplecode);
  rc1, sim1 = stan(stanmodel1, [Dict("sigma" => 1.)], CmdStanDir=CMDSTAN_HOME);
  describe(sim1)

  global stanmodel2, rc2, sim2
  stanmodel2 = Stanmodel(Sample(save_warmup=true, thin=1), name="simple",
    thin=5, model=simplecode);
  rc2, sim2 = stan(stanmodel2, [Dict("sigma" => 1.)], CmdStanDir=CMDSTAN_HOME);
  describe(sim2)

  global stanmodel3, rc3, sim3
  stanmodel3 = Stanmodel(Sample(save_warmup=true, thin=5), name="simple", model=simplecode);
  rc3, sim3 = stan(stanmodel3, [Dict("sigma" => 1.)], CmdStanDir=CMDSTAN_HOME);
  describe(sim3)

end

# stanmodel3 should not fail.


text:
# A Julia package for generic conic outer approximation methods
module Aspasia

using Printf

include("Cuts/Cuts.jl")

import MathOptInterface
const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities
const VI = MOI.VariableIndex
const SAF = MOI.ScalarAffineFunction{Float64}
const VV = MOI.VectorOfVariables
const VAF = MOI.VectorAffineFunction{Float64}
const CI = MOI.ConstraintIndex

# supported sets needing outer approximation
const ApproxSet = Union{MOI.SecondOrderCone, MOI.PositiveSemidefiniteConeTriangle}

include("optimize.jl")
include("MOI_wrapper.jl")

end