text:
module ShaleDrillingLikelihood_OverallLikelihood_Test

# using Revise
# using Juno
# using Profile
# using ProfileView
# using PProf

const DOBTIME=false

using ShaleDrillingLikelihood
using Test
using StatsFuns
using Calculus
using Optim
using Random
using InteractiveUtils
using LinearAlgebra
using BenchmarkTools

using Base: OneTo

using ShaleDrillingLikelihood: AbstractDataSet,
    DataRoyalty, DataProduce, DataDrill,
    TestDrillModel,
    simloglik!,
    SimulationDraws,
    idx_drill_ρ, idx_royalty_ρ,
    idx_drill_ψ, idx_produce_ψ,
    _nparm,
    thetas

println("testing overall likelihood")

@testset "Testing Joint Likelihoods" begin
    num_i = 1_200
    L = 3
    M = 1_000
    u = randn(num_i)
    v = randn(num_i)
    psi2 = u

    # set up coefs
    θρ = 0.5
    αψ = 0.6

    # set up drilling coefs
    model_drill = TestDrillModel()
    θ_drill   = [αψ, 2.0, -2.0, -0.75, θρ]
    θ_royalty = [θρ, 1.0,    -1.0, 1.0, 1.0,    -0.6, 0.6]  # dψdρ, ψ, β, κ
    θ_produce = vcat(αψ, rand(3), 0.3, 0.4)

    θ = vcat(θ_drill, θ_royalty[2:end], θ_produce[2:end])
    @test idx_produce_ψ isa Function
    @test idx_drill_ψ isa Function
    coef_links = [(idx_produce_ψ, idx_drill_ψ,),]

    # make data
    data_drill_opt = (num_zt=200, minmaxleases=1:2, nper_initial=10:20, nper_development=0:10, tstart=1:50,)
    data_drill   = DataDrill(u, v, TestDrillModel(), θ_drill; data_drill_opt...)
    data_produce = DataProduce(u, 10, 10:20, θ_produce)
    data_royalty = DataRoyalty(u,v,θ_royalty,L)
    data_full = DataSetofSets(data_drill, data_royalty, data_produce, coef_links)
    sim = SimulationDraws(M, data_drill)

    @testset "Joint likelihood of Royalty + Pdxn" begin

        println("testing joint likelihood of royalty + pdxn")
        data = DataSetofSets(EmptyDataSet(), data_royalty, data_produce)
        theta = vcat(θ_royalty, θ_produce)
        grad = similar(theta)
        nparm = _nparm(data)
        hess = Matrix{eltype(theta)}(undef, nparm, nparm)
        tmpgrads = Matrix{eltype(theta)}(undef, nparm, num_i)

        simloglik!(grad, hess, tmpgrads, data, theta, sim, false)
        simloglik!(grad, hess, tmpgrads, data, theta, sim, true)

        fd = Calculus.gradient(xx -> simloglik!(grad, hess, tmpgrads, data, xx, sim, false), theta, :central)

        fill!(grad,0)
        fill!(hess,0)
        simloglik!(grad, hess, tmpgrads, data, theta, sim, true)
        @test !all(grad.==0)
        @test isapprox(fd, grad; rtol=2e-5)

        if DOBTIME
            print("")
            @show @benchmark simloglik!($grad, $hess, $tmpgrads, $data, $theta, $sim, false)
            @show @benchmark simloglik!($grad, $hess, $tmpgrads, $data, $theta, $sim, true)
            print("")
        end

        # @code_warntype simloglik!(grad, hess, tmpgrads, data, theta, sim, false)
        # println("\n\n\n----------------------------\n\n\n")
        # @code_warntype simloglik!(grad, hess, tmpgrads, data, theta, sim, true)

        # Profile.clear()
        # @profile simloglik!(grad, hess, tmpgrads, data, theta, sim, false)
        # # Juno.profiletree()
        # # Juno.profiler()
        # Profile.print(format=:flat)
        # # ProfileView.view()
        # pprof()

        # Profile.clear()
        # @profile simloglik!(grad, hess, tmpgrads, data, theta, sim, true)
        # ProfileView.view()
    end

    @testset "Joint likelihood of Drill + Roy + Pdxn" begin

        println("testing joint likelihood of FULL data!")
        data = data_full
        theta = θ
        grad = similar(theta)
        nparm = _nparm(data)
        hess = zeros(nparm, nparm)
        tmpgrads = zeros(nparm, num_i)

        println("simloglik, no grad")
        simloglik!(grad, hess, tmpgrads, data, theta, sim, false)
        println("simloglik, with grad")
        simloglik!(grad, hess, tmpgrads, data, theta, sim, true)

        println("simloglik finite diff")
        fd = Calculus.gradient(xx -> simloglik!(grad, hess, tmpgrads, data, xx, sim, false), theta, :central)

        fill!(grad,0)
        fill!(hess,0)
        println("simloglik gradient")
        simloglik!(grad, hess, tmpgrads, data, theta, sim, true)
        @test !all(grad.==0)
        @test isapprox(fd, grad; rtol=2e-5)
        println("done")

        if DOBTIME
            print("")
            @show @benchmark simloglik!($grad, $hess, $tmpgrads, $data, $theta, $sim, false)
            @show @benchmark simloglik!($grad, $hess, $tmpgrads, $data, $theta, $sim, true)
            print("")
        end
        # @code_warntype simloglik!(grad, hess, tmpgrads, data, theta, sim, false)
        # println("\n\n\n----------------------------\n\n\n")
        # @code_warntype simloglik!(grad, hess, tmpgrads, data, theta, sim, true)

        # Profile.clear()
        # @profile simloglik!(grad, hess, tmpgrads, data, theta, sim, false)
        # # Juno.profiletree()
        # # Juno.profiler()
        # Profile.print(format=:flat)
        # # ProfileView.view()
        # pprof()

        # Profile.clear()
        # @profile simloglik!(grad, hess, tmpgrads, data, theta, sim, true)
        # ProfileView.view()
    end

end

end # module


text:
using HDF5
using MCMCChains
using MCMCChainsStorage
using Test

@testset "Reading and writing a chain" begin
    # Test writing a chain out and reading it back in, with sections
    chain = Chains(randn(500, 4, 2), [:a, :b, :c, :d])
    chain = set_section(chain, Dict(:internals => [:c, :d], :parameters => [:a, :b]))

    mktemp() do path, io
        h5open(path, "w") do f
            write(f, chain)
        end

        chain2 = h5open(path, "r") do f
            read(f, Chains)
        end

        for par in [:a, :b, :c, :d]
            @test chain[par] == chain2[par]
        end
    end
end

ext: jl


text:
#
# Solves the exponential-cone problem
#
#    minimize    x1+x2
#    subject to  x1 + x2 + x3 = 1
#                (x1,x2,x3) in Kexp
#
#    where Kexp = { (x1,x2,x3) | x1 >= x2*exp(x3/x2) }.
#

using Mosek
printstream(msg::AbstractString) = print(msg)

bkc = [ MSK_BK_FX ]
blc = [ 1.0 ]
buc = [ 1.0 ]

c   = [ 1.0, 1.0, 0.0 ]
bkx = [ MSK_BK_FR,MSK_BK_FR,MSK_BK_FR ]
blx = [ 0.0, 0.0, 0.0 ]
bux = [ 0.0, 0.0, 0.0 ]

asub  = [ 1, 2, 3 ]
aval  = [ 1.0, 1.0, 1.0 ]

numvar = length(bkx)
numcon = length(bkc)

# Create a task
maketask() do task
    putstreamfunc(task,MSK_STREAM_LOG,msg -> print(msg))

    appendcons(task,numcon)
    appendvars(task,numvar)

    # Set the linear term c_j in the objective.
    putclist(task,[1:3;],c)

    # Set the bounds on variable j
    putvarboundslice(task,1,numvar+1,bkx,blx,bux)

    putarow(task,1,asub,aval)
    putconbound(task,1,bkc[1],blc[1],buc[1])

    # Input the cones
    appendcone(task,MSK_CT_PEXP, 0.0, [1, 2, 3 ])

    # Input the objective sense (minimize/maximize)
    putobjsense(task,MSK_OBJECTIVE_SENSE_MINIMIZE)

    # Optimize the task
    optimize(task)
    # Print a summary containing information
    # about the solution for debugging purposes
    solutionsummary(task,MSK_STREAM_MSG)
    prosta = getprosta(task,MSK_SOL_ITR)
    solsta = getsolsta(task,MSK_SOL_ITR)

    if solsta == MSK_SOL_STA_OPTIMAL || solsta == MSK_SOL_STA_NEAR_OPTIMAL
        # Output a solution
        x1,x2,x3 = getxx(task,MSK_SOL_ITR)
        println("Optimal solution: (x1,x2,x3)=($([x1,x2,x3]))")
        println("Solution is on boundary: x1-x2*exp(x3/x2)=$(x1-x2*exp(x3/x2))")
    elseif solsta == MSK_SOL_STA_DUAL_INFEAS_CER
        println("Primal or dual infeasibility.")
    elseif solsta == MSK_SOL_STA_PRIM_INFEAS_CER
        println("Primal or dual infeasibility.")
    elseif solsta == MSK_SOL_STA_NEAR_DUAL_INFEAS_CER
        println("Primal or dual infeasibility.")
    elseif  solsta == MSK_SOL_STA_NEAR_PRIM_INFEAS_CER
        println("Primal or dual infeasibility.")
    elseif solsta == MSK_SOL_STA_UNKNOWN
        println("Unknown solution status")
    else
        println("Other solution status")
    end
end


text:
function foo{T}(ta::TimeArray{T,1}; prices=false, log_transform=false)
    #r = keyword_check(ta, prices=prices, log_transform=log_transform)
    r = keyword_check(ta, prices, log_transform)

    sum(r.values)
end

function bar(x)
    mean(x)
end

######## analyze ###################

function analyze(x;funcs=tradestats)
    for i in 1:length(funcs)
        @printf("The value for %s is %f", funcs[i], funcs[i](x))
        println("")
    end
end

######## annualized return ###################
### reference Bacon, Carl. Practical Portfolio Performance Measurement and Attribution. Wiley. 2004. p. 25

### @doc """
### # Annualized Returns
###
### # Keyword Argument defaults
###       prices = false
###       log_transform = false
###       method = "arithmetic" ("geometric" supported)
###       period = 252 (for daily)
###
### # Details
###
### # References
###   Bacon, Carl. Practical Portfolio Performance Measurement and Attribution. Wiley. 2004. p. 25
### """ ->

function annualized_return{T}(ta::TimeArray{T,1}; prices=false, log_transform=false, method="arithmetic", periods=252)
    r = keyword_check(ta, prices, log_transform)
    n = length(r)
    method == "arithmetic" ? (periods/n) * sum(r.values) :
    method == "geometric" ? prod(1+r.values) ^ (periods/n) - 1 :
    error("only arithmetic and geometric methods supported")
end

######## equity curve ########################

function equity{T}(ta::TimeArray{T,1}; prices=false)
    prices ?
    TimeArray(ta.timestamp, [1.0, expm1(cumsum(diff(log(ta.values)))) + 1], ["equity"]) :
    TimeArray(ta.timestamp, [expm1(cumsum(ta.values)) + 1], ["equity"])
end

######## from the oldmatrix.jl

# I think this might be incorrect
# it claims to be correct by treating the µ vector as a scalar for single vector case
#function ∑matrix(x::Vector)
function sigma(x::Vector)
    y = x * mean(x)
    y * y'
end

function Shapiro_Wilks{T<:Float64}(x::Vector{T})
    sx = sort(x)
    d = Truncated(Normal(mean(xs),1),minimum(xs), maximum(xs))
    m = sort(rand(d, length(x)))
    #V = ∑matrix(sx)
    V = sigma(sx)
    a = m' * inv(V) / sqrt(m' * inv(V) * inv(V) * m)
    y - x.-mean(x)

    # ∑ ai * xi and then ^ 2
    # so a cross product and then square it
    # a'x ^ 2
    num = a'x ^ 2

    # ∑ (xi - xbar)  ^ 2
    # so first get y = xi - xbar vector
    # y'y
    den = y'y
    num/den
end

######## from the matrix.jl

function ∑matrix(x::Vector, y::Vector)
    z = x .* y
    z .* z'
end

function shapiro{T<:Float64}(x::Vector{T})
    sx = sort(x)
    # d = Truncated(Normal(mean(xs),1),minimum(xs), maximum(xs))
    # m = sort(rand(d, length(x)))

    V = ∑matrix(sx)
    # a = m' * inv(V) / sqrt(m' * inv(V) * inv(V) * m)

    # ∑ ai * xi and then ^ 2
    # so a cross product and then square it
    # a'x ^ 2

    num = a'x ^ 2

    # ∑ (xi - xbar)  ^ 2
    # so first get y = xi - xbar vector
    # y'y

    y - x.-mean(x)
    den = y'y

    num/den
end

## _____ Shapiro Wilk fortran code___________________

######===========     Perform the Shapiro-Wilk test for normality.
######===========
######===========     The Shapiro-Wilk test tests the null hypothesis that the
######===========     data was drawn from a normal distribution.
######===========
######===========     Parameters
######===========     ----------
######===========     x : array_like
######===========         Array of sample data.
######===========     a : array_like, optional
######===========         Array of internal parameters used in the calculation.  If these
######===========         are not given, they will be computed internally.  If x has length
######===========         n, then a must have length n/2.
######===========     reta : bool, optional
######===========         Whether or not to return the internally computed a values.  The
######===========         default is False.
######===========
######===========     Returns
######===========     -------
######===========     W : float
######===========         The test statistic.
######===========     p-value : float
######===========         The p-value for the hypothesis test.
######===========     a : array_like, optional
######===========         If `reta` is True, then these are the internally computed "a"
######===========         values that may be passed into this function on future calls.
######===========
######===========     """
######===========     N = len(x)
######===========     if N < 3:
######===========         raise ValueError("Data must be at least length 3.")
######===========     if a is None:
######===========         a = zeros(N,'f')
######===========         init = 0
######===========     else:
######===========         if len(a) != N//2:
######===========             raise ValueError("len(a) must equal len(x)/2")
######===========         init = 1
######===========     y = sort(x)
######===========     a, w, pw, ifault = statlib.swilk(y, a[:N//2], init)
######===========     if ifault not in [0,2]:
######===========         warnings.warn(str(ifault))
######===========     if N > 5000:
######===========         warnings.warn("p-value may not be accurate for N > 5000.")
######===========     if reta:
######===========         return w, pw, a
######===========     else:
######===========         return w, pw

function swilk(x::Vector{Float64})

    fstats = dlopen(Pkg.dir("FinanceAnalytics//bin/swilk"))
    swilk_ = dlsym(fstats, :swilk_)

        #     SUBROUTINE SWILK (INIT, X, N, N1 , N2, A, W, PW, IFAULT)
        #     .. Scalar Arguments ..
        #      REAL               W, PW
        #      INTEGER            N, N1, N2, IFAULT
        #      LOGIC              INIT
        #     .. Array Arguments ..
        #      REAL               X(*), A(*)

    y      = sort(x)
    init   = false  # scipy default
    n      = length(y)  # needs some error checking for bounds
    n1     = 1   # this is a w.a.g. ~ changing this to 0 gets same ridiculous result
    n2     = ifloor(n/2) # scipy uses n//2
    w      = 0.0
    pw     = 1.0  # another w.a.g.
    ifault = 0
    a      = zeros(Float64, n2)

    ccall(swilk_,
          Float64,
          (Ptr{Bool}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}, Ptr{Int32}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}),
          &init, y, &n, &n1, &n2, a, &w, &pw, &ifault)
end

## _____ toy example below ___________________

function greet(n::Int)
    #greets   = dlopen(Pkg.dir("FinanceStats/bin/greetings"))
    greets   = dlopen("../bin/greetings")
    hello_   = dlsym(greets,:hello_)
    goodbye_ = dlsym(greets,:goodbye_)
    iseven(n)?
    ccall(hello_,Void,()) :
    ccall(goodbye_,Void,())
end

# function greet(n::Int)

#     # greets   = dlopen("../bin/greetings")
#     greets   = dlopen(Pkg.dir("FinanceStats/bin/greetings"))
# #    greets   = Pkg.dir("FinanceStats/bin/greetings")
# #    const greets   = "../bin/greetings"
#     hello_   = dlsym(greets,:hello_)
#     goodbye_ = dlsym(greets,:goodbye_)
#     iseven(n)?
# #    ccall((:hello_, greets) ,Void,()) :
# #    ccall((:goodbye_, greets),Void,())
#     ccall((:hello, greets) ,Void,()) :
#     ccall((:goodbye, greets),Void,())
# #    ccall((:hello_, "../bin/greetings") ,Void,()) :
# #    ccall((:goodbye_, greets),Void,())
# end

## _____ this works ___________________

### function greet(n::Int)
###     greets   = dlopen(Pkg.dir("FinanceStats/bin/greetings"))
###     hello_   = dlsym(greets,:hello_)
###     goodbye_ = dlsym(greets,:goodbye_)
###     iseven(n)?
###     ccall(hello_,Void,()) :
###     ccall(goodbye_,Void,())
### end


text:
import SeisIO: BUF,
  KW,
  add_chan!,
  checkbuf!,
  checkbuf_8!,
  dtconst,
  fastread,
  fastseekend,
  fillx_i16_be!,
  fillx_i32_be!,
  mk_t!,
  sμ,
  μs
import SeisIO.Quake: unsafe_convert
import SeisIO.Formats: formats,
  FmtVer,
  FormatDesc,
  HistVec


text:
module ComputedArrays

# A = ComputedArray(fn, xs, ys, zs)
# => A[i,j,k] = fn(xs[i], ys[j], zs[k])

# A = ComputedArray(fn, xs, ys, zs, order=[2,3,1])
# => A[j,k,i] = fn(xs[i], ys[j], zs[k])
# => A[i,j,k] = fn(xs[k], ys[i], zs[j])

import Base: ndims, size, getindex, show
export ComputedArray, ComputedVector, ComputedMatrix

struct ComputedArray{F,C,T,N} <: DenseArray{T,N}
   fn     :: F
   coords :: C
   shape  :: NTuple{N,Int}
   perm   :: NTuple{N,Int}
end

const ComputedVector{F,C,T} = ComputedArray{F,C,T,1}
const ComputedMatrix{F,C,T} = ComputedArray{F,C,T,2}

function ComputedArray(fn, coords...; order=1:length(coords))
   N = length(coords)
   argtypes = map(eltype, coords)
   local T
   try
      T = Union{Base.return_types(fn, argtypes)...}
   catch
      @warn "Cannot infer return type, defaulting to Any."
      T = Any
   end
   @assert length(order)==N "order must have $N elements"
   @assert isperm(order) "order=$order is not a permutation"
   shape = ntuple(d->length(coords[order[d]]), N)
   perm = tuple( invperm(order)... )
   ComputedArray{typeof(fn), typeof(coords), T, N}(fn, coords, shape, perm)
end

size(A::ComputedArray) = A.shape

getindex(A::ComputedVector, i::Int) = A.fn(A.coords[1][i])

@generated function getindex(A::ComputedArray{F,C,T,N}, I::Vararg{Int,N}) where {F,C,T,N}
   defex = Expr(:block)
   callex = :( A.fn() )
   for d = 1:N
      xsym = Symbol(:x, d)
      push!(defex.args, :( $xsym = A.coords[$d][I[A.perm[$d]]] ))
      push!(callex.args, xsym)
   end
   # thanks to @mbauman for the following trick
   # inlining is needed to avoid allocations for larger N and more complicated F
   return :($defex; Base.@_inline_meta; $callex)
end

repr(A::ComputedArray{F,C,T,N}) where {F,C,T,N} =
   (N==1 ? "$(length(A))-element" : join(size(A), "x")) *
   " ComputedArray of $(A.fn)(" * join(map(xs->"::$(eltype(xs))", A.coords), ", ") * ")::$T"

show(io::IO, A::ComputedArray) = begin; write(io, repr(A)); nothing; end
# prevent REPL from doing fancy behavior inherited from AbstractArray
show(io::IO, ::MIME"text/plain", A::ComputedArray) = show(io, A)

end # module


text:
import HTTP
import JSON
import CodecZlib
import TranscodingStreams

const country_to_currency_conversion = Base.ImmutableDict(
    "CA" => "CAD",
    "MX" => "MXN",
    "US" => "USD",
    "DK" => "DKK",
    "EE" => "EUR",
    "FI" => "EUR",
    "IE" => "EUR",
    "LV" => "EUR",
    "LT" => "EUR",
    "NO" => "NOK",
    "SE" => "SEK",
    "GB" => "GBP",
    "HR" => "EUR",
    "CY" => "EUR",
    "GR" => "EUR",
    "IT" => "EUR",
    "MT" => "EUR",
    "PT" => "EUR",
    "SI" => "EUR",
    "ES" => "EUR",
    "BG" => "EUR",
    "CZ" => "CZK",
    "HU" => "EUR",
    "PL" => "PLN",
    "RO" => "EUR",
    "RU" => "RUB",
    "SK" => "EUR",
    "AT" => "EUR",
    "BE" => "EUR",
    "FR" => "EUR",
    "DE" => "EUR",
    "LU" => "EUR",
    "NL" => "EUR",
    "CH" => "CHF",
    "AU" => "AUD",
    "NZ" => "NZD",
    "ZA" => "ZAR",
    "JP" => "JPY",
)

function download_decompressed_gzip_json(url::String)::Dict
    r = HTTP.get(url)
    decompressed = transcode(CodecZlib.GzipDecompressor, r.body)
    JSON.parse(String(decompressed))
end

function download_game_data()::Vector{Dict{String,Any}}
    url = "http://eshop-checker.xyz/games.json"
    result = download_decompressed_gzip_json(url)
    [info["price"] for info in result["list"] if any(x -> x !== 0, values(info["price"]))]
end

function download_currency_rate()::Dict{String,Any}
    url = "http://eshop-checker.xyz/beta/statics/conv_rate.json"
    result = download_decompressed_gzip_json(url)
    result["rates"]
end

function compute(data::Vector{Dict{String,Any}}, rate::Dict{String,Any})
    release_count = Dict([(area, 0) for area in Base.keys(country_to_currency_conversion)])
    price_count = Dict([(area, 0) for area in Base.keys(country_to_currency_conversion)])
    for game in data
        _price, _country = Inf, missing
        for (country, price) in filter(a -> a[1] !== 0, pairs(game))
            release_count[country] += 1
            normalized_price = price / rate[country_to_currency_conversion[country]]
            if normalized_price < _price
                _price, _country = normalized_price, country
            end
        end
        price_count[_country] += 1
    end
    release_count, price_count
end

function download_all()
    t_data = @async download_game_data()
    t_rate = @async download_currency_rate()
    fetch(t_data), fetch(t_rate)
end

function main()
    data, rate = download_all()

    release_count, price_count = compute(data, rate)

    most_count, c_country = findmax(price_count)
    most_rate, r_country = findmax(Dict([(area, price_win / get(release_count, area, Inf)) for (area, price_win) in price_count]))

    println("$c_country has $most_count games with good prices")
    println("$r_country has $(round(most_rate * 100, digits=2))% of games with good prices")
    println(release_count, price_count)
end


text:
using Gnuplot
U  = LinRange(0,10π, 80)
V = LinRange(0,2π, 20)
x = [(1-0.1*cos(v))*cos(u) for u in U, v in V]
y = [(1-0.1*cos(v))*sin(u) for u in U, v in V]
z = [0.1*(sin(v) + u/1.7 - 10) for u in U, v in V]
@gsp x y z "w l notit lc 'orange'" "set xyplane -0.05"
save(term="pngcairo size 600,400", output="tubes1.png") # hide

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
################################################################################
#
#  Contained in the alternating group
#
################################################################################

function is_discriminant_square(IdG::GAP.GapObj)
  G = GAP.Globals.SmallGroup(IdG)
  mp = GAP.Globals.RegularActionHomomorphism(G)
  S = GAP.Globals.ImagesSource(mp)
  lg = GAP.Globals.GeneratorsOfGroup(S)
  for i = 1:length(lg)
    s = GAP.Globals.SignPerm(lg[i])
    if !isone(s)
      return false
    end
  end
  return true
end

################################################################################
#
#  Check if an extension must be ramified
#
################################################################################

function must_be_ramified(L::GAP.GapObj, i::Int)
  GG = L[1]
  HH = L[i+1]
  KK = L[i]
  proj = GAP.Globals.NaturalHomomorphismByNormalSubgroup(GG, HH)
  target_grp = GAP.Globals.ImagesSource(proj)
  mH1 = GAP.Globals.NaturalHomomorphismByNormalSubgroup(target_grp, GAP.Globals.Image(proj, KK))
  H = GAP.Globals.ImagesSource(mH1)
  #I need to consider the conjgacy classes of cyclic subgroups
  #As far as I know, there is no way of getting only the cyclic ones...
  lS = GAP.Globals.ConjugacyClassesSubgroups(H)
  found_one = false
  found_all = true
  for i = 1:length(lS)
    S = GAP.Globals.Representative(lS[i])
    if !GAP.Globals.IsCyclic(S) || isone(GAP.Globals.Size(S))
      continue
    end
    g = GAP.Globals.MinimalGeneratingSet(S)[1]
    n = GAP.Globals.Order(g)
    preimgs = GAP.Globals.List(GAP.Globals.PreImages(mH1, g))
    fl = true
    for x in preimgs
      m = GAP.Globals.Order(x)
      if m == n
        fl = false
        found_all = false
        break
      end
    end
    if fl
      found_one = true
    end
    if !found_all && found_one
      break
    end
  end
  return found_all, found_one
end

################################################################################
#
#  Conductors
#
################################################################################

function _conductors_using_cocycles(F::FieldsTower, st::Vector{Int}, l_cond::Vector)
  lp = ramified_primes(F)
  auts = automorphisms(F.field, copy = false)
  projections = F.projections_for_conductors
  G = GAP.Globals.ImagesSource(projections[1])
  E = GAP.Globals.Source(projections[1])
  D = F.isomorphism
  n = prod(st)
  O = maximal_order(F)
  inertia_subgroups = Dict{fmpz, Vector{NfToNfMor}}()
  for p in lp
    lP = prime_decomposition(O, p)
    Hp = inertia_subgroup(lP[1][1])
    gHp = small_generating_set(Hp)
    inertia_subgroups[p] = gHp
  end
  must_ramify = Vector{Vector{fmpz}}(undef, length(projections))
  cant_ramify = Vector{Vector{fmpz}}(undef, length(projections))
  for s = 1:length(projections)
    proj = projections[s]
    ramify_here = Vector{fmpz}()
    not_ramify_here = Vector{fmpz}()
    for p in lp
      gHp = inertia_subgroups[p]
      els = [D[g] for g in gHp]
      sub = GAP.Globals.Subgroup(G, GAP.julia_to_gap(els))
      ord = GAP.Globals.Size(sub)
      subgs = Vector{GAP.GapObj}()
      preimages = Vector{Vector{GAP.GapObj}}(undef, length(els))
      for j = 1:length(els)
        preimages[j] = GAP.Globals.List(GAP.Globals.PreImages(proj, els[j]))
      end
      #Now, I need to check all the possible subgroups.
      it = cartesian_product_iterator(UnitRange{Int}[1:n for i = 1:length(els)], inplace = true)
      sizes_preimages = Int[]
      for I in it
        sub = GAP.Globals.Subgroup(E, GAP.julia_to_gap([preimages[i][I[i]] for i = 1:length(els)]))
        push!(sizes_preimages, GAP.Globals.Size(sub))
        if maximum(sizes_preimages) != ord && minimum(sizes_preimages) == ord
          break
        end
      end
      if minimum(sizes_preimages) != ord
        push!(ramify_here, p)
      elseif maximum(sizes_preimages) == ord && !divides(fmpz(n), p)[1]
        push!(not_ramify_here, p)
      end
    end
    if isempty(ramify_here) && isempty(not_ramify_here)
      return l_cond
    end
    must_ramify[s] = ramify_here
    cant_ramify[s] = not_ramify_here
  end
  #Now, we use what we found.
  l_new = typeof(l_cond)()
  for i = 1:length(l_cond)
    tp, wp = l_cond[i]
    satisfies_emb = trues(length(projections))
    for j = 1:length(projections)
      possible_emb = true
      for p in must_ramify[j]
        if divides(fmpz(tp), p)[1]
          continue
        end
        if divides(fmpz(n), p)[1]
          found = false
          for (q, w) in wp
            if minimum(q, copy = false) == p
              found = true
              break
            end
          end
          if !found
            possible_emb = false
            break
          end
        else
          possible_emb = false
        end
        if !possible_emb
          break
        end
      end
      if !possible_emb
        satisfies_emb[j] = false
        continue
      end
      for p in cant_ramify[j]
        if divides(fmpz(tp), p)[1]
          possible_emb = false
          continue
        end
        if divides(fmpz(n), p)[1]
          found = false
          for (q, w) in wp
            if minimum(q, copy = false) == p
              found = true
              break
            end
          end
          if found
            possible_emb = false
            break
          end
        end
        if !possible_emb
          break
        end
      end
      if !possible_emb
        satisfies_emb[j] = false
      else
        break
      end
    end
    if any(satisfies_emb)
      push!(l_new, l_cond[i])
    end
  end
  return l_new
end

function conductors_with_restrictions(F::FieldsTower, st::Vector{Int}, IdG::GAP.GapObj, bound::fmpz; unramified_outside::Vector{fmpz} = fmpz[])

  O = maximal_order(F)
  l_cond = Hecke.conductors(O, st, bound, unramified_outside = unramified_outside)
  G = GAP.Globals.SmallGroup(IdG)
  new_conds = _conductors_using_cocycles(F, st, l_cond)
  if length(st) != 1 || !isprime(st[1]) || isempty(new_conds)
    return new_conds
  end
  #If the extension is cyclic, I take care of the discriminant being a square or not for the wild ramification
  issquare = is_discriminant_square(IdG)
  p = st[1]
  v = valuation(discriminant(O), p)
  is_square_disc_base_field = iszero(mod(v*p, 2))
  td = prime_decomposition_type(O, p)
  if iseven(length(td)) || iseven(td[1][1]) || isodd(p)
    #Regardless of the exponents, the norm of the discriminant will be a square
    if issquare && is_square_disc_base_field
      return new_conds
    elseif issquare
      return typeof(new_conds)()
    else
      return new_conds
    end
  end
  #Now, p must be 2.
  if issquare && is_square_disc_base_field
    #Only the even exponents are allowed!
    newer_conds = typeof(new_conds)()
    for i = 1:length(new_conds)
      if isempty(new_conds[i][2])
        push!(newer_conds, new_conds[i])
        continue
      end
      if iszero(mod(first(values(new_conds[i][2])), 2))
        push!(newer_conds, new_conds[i])
      end
    end
  elseif issquare
    #Only the odd exponents are allowed!
    newer_conds = typeof(new_conds)()
    for i = 1:length(new_conds)
      if !isempty(new_conds[i][2]) && !iszero(mod(first(values(new_conds[i][2])), 2))
        push!(newer_conds, new_conds[i])
      end
    end
  else
    newer_conds = new_conds
  end

  #Now, tame ramification.
  list_tame = Int[x[1] for x in newer_conds]
  list_tame = coprime_base(list_tame)
  l = length(list_tame)
  for i = 1:length(list_tame)
    x = list_tame[i]
    if !isone(x) && !isprime(x)
      append!(list_tame, Hecke.divisors(x))
    end
  end
  list_tame = coprime_base(list_tame)
  for q in list_tame
    q == 1 && continue
    #@assert isprime(q)
    v = valuation(discriminant(O), q)
    is_square_disc_base_field = iszero(mod(v*p, 2))
    td = prime_decomposition_type(O, q)
    if iszero(mod(length(td) * td[1][1] * (p-1), 2))
      #Regardless of the exponents, the norm of the discriminant will be a square
      if issquare && is_square_disc_base_field
        continue
      elseif issquare || is_square_disc_base_field
        return typeof(new_conds)()
      else
        continue
      end
    end
    #Now, p must be 2.
    if issquare && is_square_disc_base_field
      #Only the even exponents are allowed!
      #Therefore the prime can't ramify
      newest_conds = typeof(new_conds)()
      for i = 1:length(newer_conds)
        if !iszero(mod(newer_conds[i][1], q))
          push!(newer_conds, newer_conds[i])
        end
      end
    elseif issquare
      #Only the odd exponents are allowed!
      #Therefore the prime must ramify
      newest_conds = typeof(new_conds)()
      for i = 1:length(newer_conds)
        if iszero(mod(newer_conds[i][1], q))
          push!(newest_conds, newer_conds[i])
        end
      end
    else
      newest_conds = newer_conds
    end
    newer_conds = newest_conds
  end
  return newer_conds
end


text:
"""
	gmtselect(cmd0::String="", arg1=nothing, kwargs...)

Select data table subsets based on multiple spatial criteria.

Full option list at [`gmtselect`]($(GMTdoc)gmtselect.html)

Parameters
----------

- $(GMT.opt_R)
- **A** | **area** :: [Type => Str | Number]

    Features with an area smaller than min_area in km^2 or of hierarchical level that is
    lower than min_level or higher than max_level will not be plotted.
    ($(GMTdoc)gmtselect.html#a)
- **C** | **dist2pt** | **dist** :: [Type => Str | NamedTuple]   `Arg = pointfile+ddist[unit] | (pts=Array, dist=xx)`

    Pass all records whose location is within dist of any of the points in the ASCII file pointfile.
    If dist is zero then the 3rd column of pointfile must have each point’s individual radius of influence.
    ($(GMTdoc)gmtselect.html#c)
- **D** | **res** | **resolution** :: [Type => Str]      `Arg = c|l|i|h|f`

    Ignored unless N is set. Selects the resolution of the coastline data set to use
    ((f)ull, (h)igh, (i)ntermediate, (l)ow, or (c)rude).
    ($(GMTdoc)gmtselect.html#d)
- **E** | **boundary** :: [Type => Str | []]            `Arg = [fn]`

    Specify how points exactly on a polygon boundary should be considered.
    ($(GMTdoc)gmtselect.html#e)
- **F** | **polygon** :: [Type => Str | GMTdaset | Mx2 array]     `Arg = polygonfile`

    Pass all records whose location is within one of the closed polygons in the multiple-segment
    file ``polygonfile`` or a GMTdataset type or a Mx2 array defining the polygon.
    ($(GMTdoc)gmtselect.html#f)
- **G** | **gridmask** :: [Type => Str | GRDgrid]        `Arg = gridmask`

    Pass all locations that are inside the valid data area of the grid gridmask.
    Nodes that are outside are either NaN or zero.
    ($(GMTdoc)gmtselect.html#g)
- **I** | **reverse** :: [Type => Str | []]    `Arg = [cflrsz]`

    Reverses the sense of the test for each of the criteria specified.
    ($(GMTdoc)gmtselect.html#i)
- $(GMT.opt_J)
- **L** | **dist2line** :: [Type => Str | NamedTuple]    `Arg = linefile+ddist[unit][+p] | (pts=Array, dist=xx, ortho=_)`

    Pass all records whose location is within dist of any of the line segments in the ASCII
    multiple-segment file linefile.
    ($(GMTdoc)gmtselect.html#l)
- **N** | **mask** :: [Type => Str | List]     `Arg = ocean/land/lake/island/pond or wet/dry`

    Pass all records whose location is inside specified geographical features.
    ($(GMTdoc)gmtselect.html#n)
- **Z** | **in_range** :: [Type => Str | List]     `Arg = min[/max][+a][+ccol][+i]`

    Pass all records whose 3rd column (z; col = 2) lies within the given range or is NaN.
    ($(GMTdoc)gmtselect.html#z)
- $(GMT.opt_V)
- $(GMT.opt_write)
- $(GMT.opt_append)
- $(GMT.opt_b)
- $(GMT.opt_d)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_g)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_o)
- $(GMT.opt_swap_xy)
"""
function gmtselect(cmd0::String="", arg1=nothing, arg2=nothing, arg3=nothing; kwargs...)

	length(kwargs) == 0 && return monolitic("gmtselect", cmd0, arg1, arg2, arg3)

	d = KW(kwargs)
	help_show_options(d)			# Check if user wants ONLY the HELP mode
	cmd, = parse_common_opts(d, "", [:R :V_params :b :d :e :f :g :h :i :o :yx])
	cmd  = parse_these_opts(cmd, d, [[:A :area], [:D :res :resolution], [:E :boundary], [:F :polygon],
	                                 [:G :gridmask], [:I :reverse], [:N :mask], [:Z :in_range]])
	#cmd = add_opt(d, cmd, 'N', [:N :mask], (ocean=("", arg2str, 1), land=("", arg2str, 2)) )

	cmd, args, n, = add_opt(d, cmd, 'C', [:C :dist2pt :dist], :pts, Array{Any,1}([arg1, arg2]), (dist="+d",))
	if (n > 0)  arg1, arg2 = args[:]  end
	cmd, args, n, = add_opt(d, cmd, 'L', [:L :dist2line], :line, Array{Any,1}([arg1, arg2, arg3]), (dist="+d", ortho="_+p"))
	if (n > 0)  arg1, arg2, arg3 = args[:]  end

	common_grd(d, cmd0, cmd, "gmtselect ", arg1, arg2, arg3)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
gmtselect(arg1, cmd0::String=""; kw...) = gmtselect(cmd0, arg1; kw...)


text:
rn = randn(1000)
bi = Binning.binindex(rn, equaldistancelimits(rn, 11))
support = binsupport(bi, 1:11)
xvals = map(x->first(x), support)
yvals = map(x->last(x), support)

sinv = sin.(0:(pi*2/999):pi*2)
bi_sinv = Binning.binindex(sinv, equaldistancelimits(sinv,11))
plot(equaldistancelimits(sinv,11), last.(binsupport(bi_sinv, 1:11)))


text:
using JuMP, VariationalInequality
using Test

# https://cdr.lib.unc.edu/indexablecontent/uuid:778ca632-74ca-4858-8c3c-6dcfc7e6e703
# Example 3.9. I created this simple two dimenional example for testing.
m = VIPModel()

@variable(m, x[1:2])

@constraint(m, x[1]^2 + x[2]^2 <= 1)

@mapping(m, F1, 3x[1] + x[2] - 5)
@mapping(m, F2, 2x[1] + 5x[2] - 3)

F = [F1, F2]
@innerproduct(m, F, x)

sol1, Fval1, gap1 = solveVIP(m, algorithm=:extra_gradient, max_iter=1000, step_size=0.1)

@show sol1
@show Fval1
@show gap1

# (0.9890, 0.1480)


text:
export PhspVector
export MultiPhspVector
using Mmap: Mmap

"""
PhspVector

A read only memory mapped Phase Space vector.
The user is responsible for closing the `io`.
"""
struct PhspVector{P,H} <: AbstractVector{P}
    io::IOStream
    header::H
    array::Vector{P}
    # owns_io::Bool
    function PhspVector{P,H}(io::IO, header::H) where {P,H}
        offset = if header isa EGSHeader
            sizeof(P)
        else
            0
        end
        dims = (get_nparticles(io, header),)
        seekstart(io)
        array = Mmap.mmap(io, Vector{P}, dims, offset)
        new{P, H}(io, header, array)
    end
end

function get_nparticles_header(h::EGSHeader)
    h.particlecount
end
function get_nparticles_header(h::IAEAHeader)
    parse(Int64, h.attributes[:PARTICLES])
end

function get_nparticles_iosize(P, io::IO)::Float64
    total_size  = bytelength(io)
    ret = (total_size / sizeof(P))
    if P <: EGSParticle
        ret -= 1 # first particle encodes the header
    end
    ret
end

function get_nparticles(io::IO, header)
    P = compressed_ptype(header)
    nheader = get_nparticles_header(header)
    ffilesize = get_nparticles_iosize(P, io)
    nfilesize = floor(Int64, ffilesize)
    if ffilesize != nfilesize
        @warn "File size indicates $(ffilesize) particles, which is not an integer. Assuming $(nfilesize) particles instead."
    end
    if nfilesize != nheader
        @warn "Particle count according to the header is $(nheader), while there are actually $(nfilesize) particles stored in the file."
    end
    nfilesize
end

function PhspVector(io::IO, h::Union{EGSHeader, IAEAHeader})
    P = compressed_ptype(h)
    H = typeof(h)
    PhspVector{P,H}(io, h)
end

function compressed_ptype(h::EGSHeader)
    ptype(h)
end
function compressed_ptype(h::IAEAHeader{Nf, Ni}) where {Nf, Ni}
    ptype(StaticIAEAHeader(h.record_contents, Nf, Ni))
end

function _PhspVector(path::AbstractString, fmt::FormatIAEA)
    _PhspVector(IAEAPath(path), fmt)
end

function _PhspVector(path::IAEAPath, ::FormatIAEA)
    io = open(path.phsp)
    h = load(path.header, IAEAHeader)
    PhspVector(io, h)
end

function _PhspVector(path::AbstractString, ::FormatEGS)
    io = open(path)
    h = consume_egs_header(io)
    seekstart(io)
    PhspVector(io, h)
end

function PhspVector(path::AbstractString)
    fmt = guess_format_from_path(path)
    _PhspVector(path, fmt)
end

function PhspVector(path::IAEAPath)
    _PhspVector(path, FormatIAEA())
end

for (f, args) in [
        [:(Base.size), ()],
        [:(Base.getindex), (:(i::Union{Integer, CartesianIndex}),)],
    ]
    @eval $f(o::PhspVector, $(args...)) = $f(o.array, $(args...))
end

Base.getindex(o::PhspVector, inds...) = view(o, inds...)

struct MultiPhspVector{P,H} <: AbstractVector{P}
    header::H # egs header can only contain typemax(Int32) particles
    phsps::Vector{PhspVector{P,H}}
    offsets::Vector{Int64}
end

function MultiPhspVector(phsps)
    @argcheck length(unique(typeof, phsps)) == 1
    header = reduce_headers(map(phsp -> phsp.header, phsps))
    offsets = cumsum(map(length, phsps))
    MultiPhspVector(header, phsps, offsets)
end

function MultiPhspVector(paths::AbstractVector{<: AbstractString})
    phsps = map(PhspVector, paths)
    MultiPhspVector(phsps)
end

function reduce_headers(headers::AbstractVector{<: EGSHeader})
    H = eltype(headers)
    @argcheck typeof(first(headers)) == eltype(headers)
    particlecount = sum(h -> h.particlecount, headers)
    photoncount = sum(h->h.photoncount, headers)
    max_E_kin = maximum(h-> h.max_E_kin, headers)
    min_E_kin_electrons = minimum(h->h.min_E_kin_electrons, headers)
    originalcount = sum(h->h.originalcount, headers)
    H(particlecount, photoncount, max_E_kin, min_E_kin_electrons, originalcount)
end

function reduce_headers(headers::AbstractVector{<:IAEAHeader})
    H::Type = only(unique(map(typeof, headers)))
    record_contents = only(unique(map(h -> h.record_contents, headers)))
    attributes = OrderedDict{Symbol, String}()
    orig_hist = 0
    for h in headers
        n = get(h.attributes, :ORIG_HISTORIES, nothing)
        if n == nothing
            @goto after_orig_histories
        else
            orig_hist += parse(Int64, n)
        end
    end
    attributes[:ORIG_HISTORIES] = string(orig_hist)
    @label after_orig_histories
    ret = H(record_contents, attributes)
    ret
end

Base.size(o::MultiPhspVector) = (last(o.offsets),)
function Base.getindex(o::MultiPhspVector, i::Integer)
    iphsp = searchsortedfirst(o.offsets, i)
    phsp = o.phsps[iphsp]
    offset = get(o.offsets, iphsp-1, 0)
    ilocal = i - offset
    phsp[ilocal]
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Mortar2D.jl/blob/master/LICENSE

module Mortar2D
using LinearAlgebra, SparseArrays
include("calculate_normals.jl")
include("calculate_projections.jl")
include("calculate_segments.jl")
include("calculate_mortar_matrices.jl")
include("calculate_mortar_assembly.jl")
export calculate_mortar_assembly
end


text:
import Pkg;

packages = ["Plots", "PlotlyJS", "GenericLinearAlgebra", "WriteVTK", "LanguageServer"]

Pkg.add(packages)

Pkg.gc()

using IJulia
installkernel("Julia Singlethreaded", env=Dict("JULIA_NUM_THREADS"=>"1"))


text:
# Author: Ali Siahkoohi, alisk@gatech.edu
# Date: September 2020
# Copyright: Georgia Institute of Technology, 2020

export loss_unsupervised, loss_supervised

function loss_unsupervised(
    Net::Union{InvertibleNetwork, ReverseNetwork},
    Y_obs,
    Zx,
    Zy,
    A::SquareCSop,
    sigma::Float32,
    RB_dist::RosenbrockDistribution
)
    X, Y, logdet = Net.forward(Zx, repeat(Zy, 1, 1, 1, size(Zx, 4)))
    CUDA.reclaim()

    Y_ = A*X
    f = -sum(logpdf(0f0, sigma, Y_ .- Y_obs)) + sum(nl_pdf(X, RB_dist)) - logdet*size(X, 4)
    ΔX = -A'*gradlogpdf(0f0, sigma, Y_ .- Y_obs) + nl_pdf_grad(X, RB_dist)
    ΔX = ΔX/size(X, 4)

    Net.backward(ΔX, 0f0.*ΔX, X, Y)
    CUDA.reclaim()
    GC.gc()

    return f/size(X, 4)
end

function loss_unsupervised(
    Net::ReverseNetwork,
    Y_obs,
    Zx,
    Zy,
    A::SquareCSop,
    sigma::Float32,
    Net_prior::NetworkConditionalHINT
)
    X, Y, logdet = Net.forward(Zx, repeat(Zy, 1, 1, 1, size(Zx, 4)); x_lane=true)
    CUDA.reclaim()

    Zx_p = Net_prior.forward(X, repeat(Y_obs, 1, 1, 1, size(X, 4)))[1]
    CUDA.reclaim()
    ΔZx = -gradlogpdf(0f0, 1f0, Zx_p)

    Y_ = A*X
    f = (-sum(logpdf(0f0, sigma, Y_ .- Y_obs))
         - sum(exact_likelihood(Net_prior, [Zx_p, Zy])) - logdet*size(X, 4))

    ΔX = (-A'*gradlogpdf(0f0, sigma, Y_ .- Y_obs)
          + Net_prior.backward(ΔZx, 0f0*Y, Zx, repeat(Zy, 1, 1, 1, size(Zx, 4)))[1])
    CUDA.reclaim()
    ΔX = ΔX/size(X, 4)

    clear_grad!(Net_prior)
    Net.backward(ΔX, zero(ΔX), X, Y)
    CUDA.reclaim()

    GC.gc()

    return f/size(X, 4)
end

function loss_unsupervised(
    Net::ReverseNetwork,
    Y_obs,
    Zx,
    Zy,
    A::RombergOP,
    sigma::Float32,
    Net_prior::NetworkConditionalHINT
)
    X, Y, logdet = Net.forward(Zx, repeat(Zy, 1, 1, 1, size(Zx, 4)); x_lane=true)
    CUDA.reclaim()

    Zx_p = Net_prior.forward(X, repeat(Y_obs, 1, 1, 1, size(X, 4)))[1]
    CUDA.reclaim()
    ΔZx = -gradlogpdf(0f0, 1f0, Zx_p)

    Y_ = wavelet_squeeze(A*wavelet_unsqueeze(X |> cpu)) |> gpu
    f = (-sum(logpdf(0f0, sigma, Y_ .- Y_obs))
         - sum(exact_likelihood(Net_prior, [Zx_p, Zy])) - logdet*size(X, 4))

    ΔX = wavelet_squeeze(-A'*wavelet_unsqueeze(gradlogpdf(0f0, sigma, (Y_ .- Y_obs) |> cpu))) |> gpu
    ΔX += Net_prior.backward(ΔZx, 0f0*Y, Zx, repeat(Zy, 1, 1, 1, size(Zx, 4)))[1]
    CUDA.reclaim()
    ΔX = ΔX/size(X, 4)

    clear_grad!(Net_prior)
    Net.backward(ΔX, zero(ΔX), X, Y)
    CUDA.reclaim()

    GC.gc()

    return f/size(X, 4)
end

function loss_unsupervised(
    Net::ReverseNetwork,
    Y_obs,
    Zx,
    Zy,
    A::Array{Float32,2},
    sigma::Float32,
    Net_prior::NetworkConditionalHINT
)
    X, Y, logdet = Net.forward(Zx, repeat(Zy, 1, 1, 1, size(Zx, 4)); x_lane=true)
    CUDA.reclaim()

    Zx_p = Net_prior.forward(X, repeat(Y_obs, 1, 1, 1, size(X, 4)))[1]
    CUDA.reclaim()
    ΔZx = -gradlogpdf(0f0, 1f0, Zx_p)

    Y_ = wavelet_squeeze(A.*wavelet_unsqueeze(X |> cpu)) |> gpu
    f = (-sum(logpdf(0f0, sigma, Y_ .- Y_obs))
         - sum(exact_likelihood(Net_prior, [Zx_p, Zy])) - logdet*size(X, 4))

    ΔX = wavelet_squeeze(-A'.*wavelet_unsqueeze(gradlogpdf(0f0, sigma, (Y_ .- Y_obs) |> cpu))) |> gpu
    ΔX += Net_prior.backward(ΔZx, 0f0*Y, Zx, repeat(Zy, 1, 1, 1, size(Zx, 4)))[1]
    CUDA.reclaim()
    ΔX = ΔX/size(X, 4)

    clear_grad!(Net_prior)
    Net.backward(ΔX, zero(ΔX), X, Y)
    CUDA.reclaim()

    GC.gc()

    return f/size(X, 4)
end

function loss_supervised(
    Net::NetworkConditionalHINT,
    X::AbstractArray{Float32,4},
    Y::AbstractArray{Float32,4}
)

    Zx, Zy, logdet = Net.forward(X, Y)
    CUDA.reclaim()
    z_size = size(Zx)

    f = sum(logpdf(0f0, 1f0, Zx))
    f = f + sum(logpdf(0f0, 1f0, Zy))
    f = f + logdet*z_size[4]

    ΔZx = -gradlogpdf(0f0, 1f0, Zx)/z_size[4]
    ΔZy = -gradlogpdf(0f0, 1f0, Zy)/z_size[4]

    ΔX, ΔY = Net.backward(ΔZx, ΔZy, Zx, Zy)[1:2]
    CUDA.reclaim()
    GC.gc()

    return -f/z_size[4], ΔX, ΔY
end


text:
export runOEDTypeB

"""

function runOEDTypeA

Input:

	param
	Xtrue
	L
	n
	Alphas
	modFun
	IPmode
	m0
"""
function runOEDTypeB(OEDparam::Array,Xtrue,L,IPmode,m0;Ce=[],ce=[],Ci=[],ci=[],pcgMaxIter=10,maxIter=20,maxStep=1.0,exName="")

	ids = mod.(0:size(Xtrue,2)-1,nworkers())+1

	println("\n\n=== setup forward problem ===")
	if IPmode==:unconstr
		pFor = Array{UnconstrainedOEDParam}(nworkers())
		for j=1:nworkers()
			pFor[j] = UnconstrainedOEDParam(OEDparam[j],L,zeros(0,0),[])
		end
	elseif IPmode==:eqconstr
		if isempty(Ce) || isempty(ce)
			error("equality constraints must be provided. To run unconstrained OED, use IPmode=:unconstr.")
		end
		pFor = Array{EqconstrainedOEDParam}(nworkers())
		for j=1:nworkers()
			pFor[j] = EqconstrainedOEDParam(OEDparam[j],L,Ce,ce[:,ids.==j],zeros(0,0),[])
		end
	elseif IPmode==:nonnegconstr
		I = speye(size(Xtrue,1))
		pFor = Array{IneqconstrainedOEDParam}(nworkers())
		for j=1:nworkers()
				pFor[j] = getIneqconstrainedOEDParam(OEDparam[j],L,Ce,ce[:,ids.==j],I,0*OEDparam[j].Ftrue,
				          ftolIPQP=1e-3,mutolIPQP=1e-3,outIPQP=0)
		end

	elseif IPmode==:boxconstr
		if isempty(Ci) || isempty(ci)
			error("inequality constraints must be provided. To run equality constrained OED, use IPmode=:eqconstr.")
		end
		pFor = Array{IneqconstrainedOEDParam}(nworkers())

		for j=1:nworkers()
		pFor[j] = 	getIneqconstrainedOEDParam(OEDparam[j],L,Ce,ce[:,ids.==j],Ci,ci[:,ids.==j],
		ftolIPQP=1e-3,mutolIPQP=1e-3,outIPQP=0)
		end
	else
		error("Invalid IPmode = $IPmode")
	end
	display(pFor[1])

	np = length(m0)

	if nworkers()>1
		pMis  = Array{RemoteChannel}(nworkers())
		for j=1:nworkers()
			pMis[j] = initRemoteChannel(getMisfitParam,workers()[j], pFor[j],ones(size(OEDparam[j].Ftrue)),Xtrue[:,ids.==j],SSDFun)
		end
	else
		pMis  = getMisfitParam(pFor[1],ones(size(Xtrue)),Xtrue,SSDFun);
	end

	# make sure the full size problem is in there
	Minv        = getRegularMesh([0  1 0 1],[12 12]);
	pInv        = getInverseParam(Minv,identityMod,smallnessReg,0.,zeros(np),0*ones(np),180*ones(np));
	pInv.pcgMaxIter = pcgMaxIter
	pInv.maxIter    = maxIter
	pInv.maxStep    = maxStep
	pInv.HesPrec.applyPrec = (a,b,x,c) -> x
	display(pInv)

	println("\n\n=== optimize projection angles (started=$(now())) ===")
	mIter = zeros(length(m0),maxIter+1)
	Diter = zeros(size(Xtrue,1),size(Xtrue,2),maxIter+1)
	if nworkers()==1
		storeInterm = (mc,Dc,iter,pInv,pMis) -> (mIter[:,iter+1] = mc; Diter[:,:,iter+1] = Dc);
	else
		storeInterm  = (mc,Dc,iter,pInv,pMis) -> (mIter[:,iter+1] = mc; for k=1:nworkers(); Diter[:,ids.==k,iter+1] = fetch(Dc[k]); end);
	end
	tic()
	mOpt,Xo,flag,His = projGN(copy(m0),pInv,pMis,dumpResults=storeInterm,solveGN=projGNexplicit);

	Xopt = 0*Xtrue
	if nworkers()>1
		Xopt = Xo
	else
		for k=1:nworkers()
			Xopt[:,ids.==k] = fetch(Xo[k])
		end
	end

	time = toc()
	println("=== finished in $time secs ===")
	if !isempty(exName)
		fileName = @sprintf "%s-%s.mat" exName IPmode
		logName = @sprintf "%s-%s.log" exName IPmode
		mv("jInv.out",logName,remove_destination=true)
		println("=== save results in $fileName===")
		His.Dc = [];
		results = Dict("mOpt"=>mOpt, "Xopt"=>Xopt,"flag"=>flag,"His"=>His,"time"=>time,"mIter"=>mIter,"Diter"=>Diter)
		matwrite(fileName,results)
	end
	return mOpt
end


text:
module lossFunctions
  export hingeLoss, MSE, zeroOneLoss

  """
      function hingeLoss( y::Float64, t::Float64 )

      compute hinge loss for a single binary classification.
      #Arguments
        - `y::Float64`: predicted value.
        - `t::Float64`: ground truth value.
      #Return
        - `Float64`:    hingeLoss error.
  """
  function hingeLoss( y::Float64, t::Float64 )
    e = 1.0 - t * y;
    if e > 0.0
      return e;
    else
      return 0.0;
    end
  end

  """
      hingeLoss( y::Vector{Float64}, t::Vector{Float64} )

      compute hinge loss for a vector of binary classifications.
      #Arguments
        - `Vector{Float64}`: predicted values.
        - `Vector{Float64}`: ground truth values.
      #Return
        - `Vector{Float64}`: vector of hingeLoss errors.
  """
  function hingeLoss( y::Vector{Float64}, t::Vector{Float64} )
    return hingeLoss.( y, t );
  end

  """
      MSE( y::Float64, t::Float64 )

      compute mean squared error for a single binary classification.
        #Arguments
        - `Float64`: predicted value.
        - `Float64`: ground truth value.
      #Return
        - `Float64`: hingeLoss error.
  """
  function MSE( y::Float64, t::Float64 )
    return ( y - t )^2;
  end

  """
      MSE( Y::Vector{Float64}, T::Vector{Float64} )

      compute mean squared error for a vector of binary classifications.
        #Arguments
        - `Vector{Float64}`: predicted value.
        - `Vector{Float64}`: ground truth value.
      #Return
        - `Vector{Float64}`: vector of hingeLoss errors.
  """
  function MSE( Y::Vector{Float64}, T::Vector{Float64} )
    return MSE.( Y, T );
  end

    """
        zeroOneLoss( y::Float64, t::Float64 )

        compute zero-one loss for a vector for a binary classifications.
          #Arguments
          - `Float64`: predicted value.
          - `Float64`: ground truth value.
          #Return
          - `Float64`: zero-one loss error.

  """
  function zeroOneLoss( y::Float64, t::Float64 )
   if y ≈ t
     return 0.0;
   else
     return 1.0;
   end
  end

  """
        zeroOneLoss( y::Float64, t::Float64 )

        compute zero-one loss for a vector of binary classifications.
          #Arguments
          - `Vector{Float64}`: predicted values.
          - `Vector{Float64}`: ground truth values.
          #Return
          - `Vector{Float64}`: vector of zero-one loss errors.

  """
  function zeroOneLoss( Y::Vector{Float64}, T::Vector{Float64} )
    return zeroOneLoss.( Y, T );
  end
end


text:
"""
    mandelbrot( c )

Computes the number of evaluations of the following equation before the resulting complex number `z` beginning with `z=c`.
- Input: `c::Number`
- Output: `itrMax::??`

### Questions
- What do you think is the output type of this function?
- Which of the part of the code determines the type?

Sometimes forcing the type of output can be an advantage.
"""
function mandelbrot( c )
    z = c
    itrMax = 80
    for n in 1:itrMax
        if abs(z) > 2
            return n-1
        end
        z = z^2 + c
    end
    return itrMax
end


text:
using Example

Example.domath(5)


text:
"""
    get_current_project_and_manifest()

Return the current project and manifest TOML files.

# Implementation details

* `Base.active_project()` is used to determine the location the current (Julia)Project.toml
* JuliaManifest.toml takes precedent over Manifest.toml in the same directory as the Project.toml
"""
function get_current_project_and_manifest()
    project_toml = Base.active_project()
    julia_manifest_toml = joinpath(dirname(project_toml), "JuliaManifest.toml")
    if isfile(julia_manifest_toml)
        manifest_toml = julia_manifest_toml
    else
        manifest_toml = joinpath(dirname(project_toml), "Manifest.toml")
    end
    project_toml, manifest_toml
end

"""
    AbstractEnvironmentMigrator

Base abstract type for EnvironmentMigrators.

# Interface:
* `source_project_toml(m::AbstractEnvironmentMigrator)`
    - return path of (Julia)Project.toml to copy from or `nothing` if there is no source
* `source_manifest_toml(m::AbstractEnvironmentMigrator)`
    - return path of (Julia)Manifest.toml to copy from or `nothing` if there is no source
* `target_project_toml(m::AbstractEnvironmentMigrator)`
    - return path of (Julia)Project.toml to backup and to which to copy
* `target_manifest_toml(m::AbstractEnvironmentMigrator)`
    - return path of (Julia)Manifest.toml to backup and to which to copy
* `fileaction(m::AbstractEnvironmentMigrator)`
    - action to take during backup operations
* `backup(m::AbstractEnvironmentMigrator)`
    - save a copy of the target environment in a timestamped subfolder
* `migrate(m::AbstractEnvironmentMigrator; backup = true, update = true)`
    - migrate project and manifest files from source to target

# Subtypes:
* `SimpleEnvironmentMigrator`
    - set source and target project and manifest files as fields
* `SimpleBackupOnlyEnvironmentMigrator`
    - set target project and manifest files as fields
"""
abstract type AbstractEnvironmentMigrator end;

"""
    SimpleEnvironmentMigrator([source_project_toml, source_target_toml], [target_project_toml, target_manifest_toml])

An `AbstractEnvironmentMigrator` where one can set source and target project and manifest files as fields.
If `source_project_toml` and `source_manifest_toml` are not specified, they will be `nothing`.
If `target_project_toml` and `target_manifest_toml` are not specified, they will be set to the current environment.

# Fields

* `source_project_toml`
* `source_manifest_toml`
* `target_project_toml`
* `target_manifest_toml`

# Method details

`fileaction(m)` will return `cp` if `source_project_toml` is `nothing` or `mv` otherwise.
"""
mutable struct SimpleEnvironmentMigrator <: AbstractEnvironmentMigrator
    source_project_toml::Union{AbstractString, Nothing}
    source_manifest_toml::Union{AbstractString, Nothing}
    target_project_toml::AbstractString
    target_manifest_toml::AbstractString
end
function SimpleEnvironmentMigrator(source_project_toml, source_manifest_toml)
    project_toml, manifest_toml = get_current_project_and_manifest()
    SimpleEnvironmentMigrator(source_project_toml, source_manifest_toml, project_toml, manifest_toml)
end
SimpleEnvironmentMigrator() = SimpleEnvironmentMigrator(nothing, nothing)
source_project_toml(m::SimpleEnvironmentMigrator) = m.source_project_toml
source_manifest_toml(m::SimpleEnvironmentMigrator) = m.source_manifest_toml
target_project_toml(m::SimpleEnvironmentMigrator) = m.target_project_toml
target_manifest_toml(m::SimpleEnvironmentMigrator) = m.target_manifest_toml
fileaction(m::SimpleEnvironmentMigrator) = m.source_project_toml === nothing ? cp : mv

"""
    SimpleBackupOnlyEnvironmentMigrator{F}([target_project_toml, target_manifest_toml], fileaction::F)

`AbstractEnvironmentMigrator` only for backing up the target project and manifest tomls to a timestamped folder.
The target project and manifest tomls will default to the current project and manifest tomls.
"""
mutable struct SimpleBackupOnlyEnvironmentMigrator{F} <: AbstractEnvironmentMigrator
    target_project_toml::AbstractString
    target_manifest_toml::AbstractString
    fileaction::F
end
function SimpleBackupOnlyEnvironmentMigrator{F}(fileaction::F) where F
    project_toml = Base.active_project()
    manifest_toml = joinpath(dirname(project_toml), "Manifest.toml")
    SimpleBackupOnlyEnvironmentMigrator{F}(project_toml, manifest_toml, fileaction)
end
SimpleBackupOnlyEnvironmentMigrator() = SimpleBackupOnlyEnvironmentMigrator{typeof(cp)}(cp)
SimpleBackupOnlyEnvironmentMigrator(fileaction::F) where F = SimpleBackupOnlyEnvironmentMigrator{F}(fileaction)

source_project_toml(::SimpleBackupOnlyEnvironmentMigrator) = nothing
source_manifest_toml(::SimpleBackupOnlyEnvironmentMigrator) = nothing
target_project_toml(m::SimpleBackupOnlyEnvironmentMigrator) = m.target_project_toml
target_manifest_toml(m::SimpleBackupOnlyEnvironmentMigrator) = m.target_manifest_toml
fileaction(m::SimpleBackupOnlyEnvironmentMigrator{F}) where F = m.fileaction

"""
    backup(m::AbstractEnvironmentMigrator)

Backup an environment. Typically, this involves copying the Project.toml and Manifest.toml
to a time-stamped backup subfolder.

# Extended Help

The standard backup procedure is as follows.
1. Check that the target project and manifest TOMLs are in the same directory.
2. Check that one of either the project or manifest TOMLs exist
3. Create a time stamp in the format "yyyy-mm-dd_HH_MM_SS
4. Create a subpath in the same directory as the project TOML called `backups/[timestamp]`
5. Copy (or move) the project and manifest TOMLs to the backup path if they exist.
"""
function backup(m::AbstractEnvironmentMigrator)
    @assert dirname(target_project_toml(m)) == dirname(target_manifest_toml(m))
    current_dir = dirname(target_project_toml(m))
    if isfile(target_project_toml(m)) || isfile(target_manifest_toml(m))
        timestamp = Dates.format(Dates.now(), "yyyy-mm-dd_HH_MM_SS")
        backup_dir = joinpath(current_dir, "backups", timestamp)
        mkpath(backup_dir)
        @info "Backing up Project.toml and Manifest.toml" target_project_toml(m) target_manifest_toml(m) backup_dir
        _fileaction = fileaction(m)
        for toml in (target_project_toml(m), target_manifest_toml(m))
            if isfile(toml)
                _fileaction(toml, joinpath(backup_dir, basename(toml)))
            else
                @info "$toml does not exist" toml
            end
        end
    else
        @info "No existing Project.toml or Manifest.toml to backup" target_project_toml(m) target_manifest_toml(m)
    end
    return nothing
end

const SimpleEnvironmentMigrators = Union{SimpleEnvironmentMigrator, SimpleBackupOnlyEnvironmentMigrator}

"""
    migrate(m::AbstractEnvironmentMigrator; backup = true, update = true)

Move an environment from one place to another.
`backup` - whether to perform a backup of the target environment (`true`) or not (`false`).
`update` - whether to perform an update of the source environment (`true`) or not (`false`).

# Extended Help

The typical migration procedure is as follows.

1. Backup the target environment by moving the files if indicated.
2. Copy the (Julia)Project.toml from the source to the target environment.
3. Copy the (Julia)Manifest.toml from the source to the target environment.
4. Activate the target environment
5. `Pkg.status()`
6. `Pkg.upgrade_manifest()` - skip if error
7. `Pkg.resolve()`
8. `Pkg.instantiate()`
9. `Pkg.status()`
10. `Pkg.update()`
"""
function migrate(m::SimpleEnvironmentMigrators; backup = true, update = true)
    if backup
        @info "Backing up the current environment" m
        EnvironmentMigrators.backup(m)
    end

    selected_project_toml = source_project_toml(m)
    selected_manifest_toml = source_manifest_toml(m)
    _target_project_toml = target_project_toml(m)
    _target_manifest_toml = target_manifest_toml(m)

    if selected_project_toml === nothing || selected_manifest_toml === nothing
        @info "No migration target. No further action."
        return;
    end

    if isfile(selected_project_toml)
        @info "Copying selected Project.toml to current environment" selected_project_toml _target_project_toml
        cp(selected_project_toml, _target_project_toml)
    else
        @info "Selected Project.toml does not exist" selected_project_toml
    end
    if isfile(selected_manifest_toml)
        @info "Copying selected Manifest.toml to current directory" selected_manifest_toml _target_manifest_toml
        cp(selected_manifest_toml, _target_manifest_toml)
    else
        @info "Selected Manifest.toml does not exist" selected_manifest_toml
    end
    try
        Pkg.activate(_target_project_toml)
        Pkg.status()
        try
            Pkg.upgrade_manifest()
        catch err
        end
        Pkg.resolve()
        try
            Pkg.instantiate()
        catch err
            @error "An error occurred during Pkg.instantiate()." err
            @info "Run `using Pkg; pkg\"instantiate\"` to view the error."
            @info "Run `using Pkg; pkg\"rm PkgName\"` to remove problematic packages."
            @info "Run `using Pkg; pkg\"update\"` to update to the latest package versions."
            return nothing
        end
    catch err
        @error "An error occurred during Pkg.resolve()." err
        @info "Run `using Pkg; Pkg.resolve()` to view the error."
        @info "Run `using Pkg; pkg\"rm PkgName\"` to remove problematic packages."
        @info "Run `using Pkg; Pkg.instantiate()` afterwards."
        @info "Run `using Pkg; pkg\"update\"` to update to the latest package versions."
        return nothing
    end
    @info "Migration successful"
    Pkg.status()
    if update
        @info "Updating"
        Pkg.update()
    end
    @info "Migration complete and environment updated. Have a nice day!" _target_project_toml
    return nothing
end


text:
# This file provides getter-type methods which essentially define the fields
# needed in any concrete subtype of EarthModel.

"""
    depth(m, radius) -> depth

Return the `depth` in km for the model `m` given a `radius` in km.
"""
depth(m::EarthModel, radius) = surface_radius(m) - radius

"""
    hasattenuation(m) -> ::Bool

Return `true` if the model `m` includes attenuation and `false` otherwise.
"""
hasattenuation(m::EarthModel) = m.attenuation

"""
    isanisotropic(m) -> ::Bool

Return `true` if the model `m` is anisotropic and `false` otherwise.
"""
isanisotropic(m::EarthModel) = m.aniso

"""
    radius(m, depth) -> radius

Return the `radius` in km for the model `m` given a `depth` in km.
"""
radius(m::EarthModel, depth) = surface_radius(m) - depth

"""
    surface_radius(m) -> radius

Return the surface radius in km of the `EarthModel` `m`.
"""
surface_radius(m::EarthModel) = m.a


text:
module INMET

using HTTP
using JSON
using DataFrames
using Unitful
using Dates
using Printf

import Unitful: °, m, °C, percent

export Date, DateTime

# -----------
# PUBLIC API
# -----------

"""
    stations(kind=:automatic)

Return INMET stations of given `kind`. There are two
kinds of stations: `:automatic` and `:manual`.
"""
function stations(kind=:automatic)
  @assert kind ∈ [:automatic,:manual] "invalid kind"
  root = "https://apitempo.inmet.gov.br/estacoes/"
  url  = root * (kind == :automatic ? "T" : "M")
  url |> download |> frame
end

"""
    series(station, start, finish, freq=:day)

Return time series for a given `station` from `start`
to `finish` date with frequency `freq`.

Available codes for stations can be found in the result of
[`stations`](@ref) in the column `CD_ESTACAO`, e.g. `:A301`.

Start and end dates are specified as Julia `Date` objects.
For example, the date `2021-07-31` is represented with
`Date(2021,7,31)`. The frequency can be `:day` or `:hour`.
"""
function series(station, start, finish, freq=:day)
  @assert freq ∈ [:day,:hour] "invalid frequency"
  root = "https://apitempo.inmet.gov.br/estacao"
  kind = freq == :day ? "diaria" : ""
  from = string(start)
  to   = string(finish)
  url  = join([root, kind, from, to, station], "/")
  url |> download |> frame
end

"""
    on(time)

Return data for all automatic stations on a given `time`.
The time can be a `Date` or a `DateTime` object. In the
latter case, minutes and seconds are ignored while the
hour information is retained (data in hourly frequency).
"""
function on(time)
  root = "https://apitempo.inmet.gov.br/estacao/dados"
  date = string(Date(time))
  hour = time isa DateTime ? (@sprintf "%02d00" Dates.hour(time)) : ""
  url  = join([root, date, hour], "/")
  url |> download |> frame
end

# -----------------
# HELPER FUNCTIONS
# -----------------

function download(url)
  page = HTTP.get(url)
  str  = String(page.body)
  JSON.parse(str)
end

function frame(data)
  # replace nothing by missing
  asmissing(v) = isnothing(v) ? missing : v
  vars = data |> first |> keys |> collect
  df = map(vars) do var
    var => [asmissing(d[var]) for d in data]
  end |> DataFrame

  # available columns
  AVAIL = propertynames(df)

  # numeric columns
  NUM = [:VL_LONGITUDE,:VL_LATITUDE,:VL_ALTITUDE,
         :TEM_INS,:TEM_MIN,:TEM_MAX,
         :TEMP_MIN,:TEMP_MED,:TEMP_MAX,
         :UMD_INS,:UMD_MIN,:UMD_MAX,
         :UMID_MIN,:UMID_MED,:UMID_MAX,
         :PRE_INS,:PRE_MIN,:PRE_MAX,
         :VEN_VEL,:VEN_RAJ,:VEN_DIR,
         :PTO_INS,:PTO_MIN,:PTO_MAX,
         :RAD_GLO,:CHUVA]

  # physical units
  UNIT = [°,°,m,
          °C,°C,°C,
          °C,°C,°C,
          percent,percent,percent,
          percent,percent,percent,
          u"mbar",u"mbar",u"mbar",
          u"m/s",u"m/s",°,
          °C,°C,°C,
          u"W/m^2",u"mm"]

  # dictionary mapping names to units
  UNIT4VAR = Dict(NUM .=> UNIT)

  # available numeric columns
  NUMAVAIL = [V for V in NUM if V in AVAIL]

  # parse numeric columns as floats with units
  str2float(v) = ismissing(v) ? missing : parse(Float64, v)
  floatunit(v, u) = ismissing(v) ? missing : v*u
  for var in NUMAVAIL
    unitful(v) = floatunit(v, UNIT4VAR[var])
    transform!(df, var => ByRow(str2float) => var)
    transform!(df, var => ByRow(unitful) => var)
  end

  select(df, Not(NUMAVAIL), NUMAVAIL)
end

end


text:
# # Schema Visualization
# In this example we show how can schema be turned into HTML interactive visualization, which helps to examine the schema,
# especially when dealing with large and heterogeneous data.

#md # !!! tip
#md #     This example is also available as a Jupyter notebook, feel free to run it yourself:
#md #     [`schema_visualization.ipynb`](@__NBVIEWER_ROOT_URL__/examples/schema_visualization.ipynb)

#nb # We start by installing JsonGrinder and few other packages we need for the example.
#nb # Julia Ecosystem follows philosophy of many small single-purpose composable packages
#nb # which may be different from e.g. python where we usually use fewer larger packages.
#nb using Pkg
#nb pkg"add JsonGrinder#master JSON"

# We include packages we want to use.
using JsonGrinder, JSON
import JsonGrinder: generate_html

# Now we load all samples
data_file = "data/recipes.json" #src
data_file = "../../../data/recipes.json" #nb
data_file = "data/recipes.json" #md
data_file = "data/recipes.json" #jl
samples_str = open(data_file) do fid
	read(fid, String)
end;

# We parse them to structures
samples = convert(Vector{Dict}, JSON.parse(samples_str));

# Print example of the JSON
JSON.print(samples[1],2)

# We create schema from all samples
sch = JsonGrinder.schema(samples)

# Now we can generate the html visualization into a file, keeping only 100 unique values per item
generate_html("recipes_max_vals=100.html", sch, max_vals=100)

# Or we can generate html, keeping all values from schema.
generate_html("recipes.html", sch, max_vals=nothing)

# If we omit the first argument, we will get the html as a string
generated_html = generate_html(sch, max_vals = 100);

# Now we can look at the visualization.
#
# Feel free to click the triangles, individual nodes of the tree are collapsed by default,
# but can be expanded or collapsed when clicked. This way you can easily examine individual parts of the schema.
# For lists we show histograms of lengths, for leaves we show histogram of values etc.
generated_html

# If you like, you may use the Electron to open it in browser.
# using the following code (this works if you run it from REPL, but not from jupyter notebook or in CI)
# ```julia
# using ElectronDisplay
# using ElectronDisplay: newdisplay
# display(newdisplay(), MIME{Symbol("text/html")}(), generated_html)
# ```


text:
module TriplotGR

using GR,TriplotBase

export tricontour,tricontourf,tripcolor

function tricontour(x,y,z,t,levels)
    contours = TriplotBase.tricontour(x,y,z,t,levels)
    lmin,lmax = extrema(getfield.(contours,:level))
    for c=contours
        GR.setlinecolorind(1000+TriplotBase.getcolorind(c.level,lmin,lmax,256))
        for polyline=c.polylines
            x = first.(polyline)
            y = last.(polyline)
            GR.polyline(x,y)
        end
    end
end

function tricontourf(x,y,z,t,levels)
    filled_contours = TriplotBase.tricontourf(x,y,z,t,levels)
    lmin,lmax = extrema(getfield.(filled_contours,:lower))
    GR.setfillintstyle(1)
    for fc=filled_contours
        x = Float64[]
        y = Float64[]
        for polyline=fc.polylines
            append_with_nan!(x,first.(polyline))
            append_with_nan!(y,last.(polyline))
        end
        GR.setfillcolorind(1000+TriplotBase.getcolorind(fc.lower,lmin,lmax,256))
        GR.fillarea(x,y)
    end
end

function append_with_nan!(a,b)
    append!(a,b)
    push!(a,NaN)
end

function tripcolor(x,y,z,t;zmin=nothing,zmax=nothing,px=nothing,py=nothing)
    _,_,px_dsp,py_dsp = GR.inqdspsize()
    isnothing(px) && (px = px_dsp)
    isnothing(py) && (py = py_dsp)
    cmap = zeros(UInt32,256)
    for i=0:255
        # Pixel is encoded as 0xAABBGGRR,GR.inqcolor returns 0x00BBGGRR,so set opaque
        cmap[begin+i] = 0xff000000 | GR.inqcolor(1000+i)
    end
    img = TriplotBase.tripcolor(x,y,z,t,cmap; px,py,zmin,zmax,yflip=true)
    GR.drawimage(0,1,0,1,size(img,1),size(img,2),img)
end

end


text:
"""
    nonsym(nbOne::Ti, size::Ti, diag_l::Ti, diag_u::Ti, spectrum::AbstractVector{Tv}) where {Tv<:Complex, Ti<:Integer}

Generate a non symmetric matrix with default nilpotent matrix and default initialization of matrix.

The usage of `nonsym` is the same as the one of `nonherm`, please refers to that part for more detais.

Please keep in mind the constraints of given spectrum vector and initialization of matrix, click [here](getting_started.md) for more details.

"""
function nonsym(
    nbOne::Ti,
    size::Ti,
    diag_l::Ti,
    diag_u::Ti,
    spectrum::AbstractVector{Tv}
) where {Tv<:Complex, Ti<:Integer}

    Tv2 = eltype(real(spectrum[1]))
    Am = spzeros(Tv2, size, size)
    initMat!(Am, diag_l, diag_u, size)
    nilp = Nilp(nbOne, size)

    return nonsym(size, diag_l, diag_u, spectrum, Am, nilp)

end

"""
    nonsym(nbOne::Ti, size::Ti, diag_l::Ti, diag_u::Ti, spectrum::AbstractVector{Tv}, Am::SparseMatrixCSC{Tv2, Ti},) where {Tv<:Complex, Tv2<:Real, Ti<:Integer}

Generate a non symmetric matrix with default nilpotent matrix and user-provided initialization of matrix.

The usage of `nonsym` is the same as the one of `nonherm`, please refers to that part for more detais.

Please keep in mind the constraints of given spectrum vector and initialization of matrix, click [here](getting_started.md) for more details.

"""
function nonsym(
    nbOne::Ti,
    size::Ti,
    diag_l::Ti,
    diag_u::Ti,
    spectrum::AbstractVector{Tv},
    Am::SparseMatrixCSC{Tv2, Ti},
) where {Tv<:Complex, Tv2<:Real, Ti<:Integer}

    nilp = Nilp(nbOne, size)

    return nonsym(size, diag_l, diag_u, spectrum, Am, nilp)

end

"""
    nonsym(size::Ti, diag_l::Ti, diag_u::Ti, spectrum::AbstractVector{Tv}, Am::SparseMatrixCSC{Tv2, Ti}, nilp::Nilpotent{Ti}) where {Tv<:Complex, Tv2<:Real, Ti<:Integer}

Generate a non symmetric matrix with user-provided nilpotent matrix and initialization of matrix.

The usage of `nonsym` is the same as the one of `nonherm`, please refers to that part for more detais.

Please keep in mind the constraints of given spectrum vector and initialization of matrix, click [here](getting_started.md) for more details.

"""
function nonsym(
    size::Ti,
    diag_l::Ti,
    diag_u::Ti,
    spectrum::AbstractVector{Tv},
    Am::SparseMatrixCSC{Tv2, Ti},
    nilp::Nilpotent{Ti}
) where {Tv<:Complex, Tv2<:Real, Ti<:Integer}

    checkSpectrum(spectrum, size)

    for i = 1:size
        Am[i, i] = real(spectrum[i])
    end

    idx = 1
    while idx < size
        if imag(spectrum[idx]) == 0
            step = 1
        else
            step = 2
        end
        Am[idx, idx+1] = abs(imag(spectrum[idx]))
        Am[idx+1,idx] = -abs(imag(spectrum[idx]))
        idx += step
    end

    matAop = Am

    fact = 1

    for k = 1 : 2 * (nilp.degree-1)
        nilpMat = NilpMat(nilp)
        matAop = matAop * nilpMat - nilpMat * matAop
        fact /= factorial(k,k)
        Am += fact * matAop
    end

    return Am
end


text:
const BothChebyshevs = Union{ChebyshevT,ChebyshevU}

isorthogonal(::ChebyshevU, ::ChebyshevUWeight) = true
isorthogonal(::ChebyshevT, ::ChebyshevTWeight) = true

isorthogonal(b::BothChebyshevs, μ::DiscreteWeight) =
	_isorthogonal(b, μ, points(μ), weights(μ))
_isorthogonal(b::BothChebyshevs, μ, points, weights) = false
_isorthogonal(b::ChebyshevT, μ, points::ChebyshevTNodes, weights::GridArrays.ChebyshevTWeights) = opsorthogonal(b, μ)
_isorthogonal(b::ChebyshevU, μ, points::ChebyshevUNodes, weights::GridArrays.ChebyshevUWeights) = opsorthogonal(b, μ)
_isorthogonal(b::ChebyshevT, μ, points::ChebyshevTNodes, weights) =
	DomainIntegrals.allequal(weights) && opsorthogonal(b, μ)

const UniformDiscreteChebyshevTWeight{T,G,W} = GenericGridWeight{T,G,W} where G <: ChebyshevNodes where W <:FillArrays.AbstractFill
isorthogonal(dict::ChebyshevT, measure::UniformDiscreteChebyshevTWeight) = opsorthogonal(dict, measure)

gauss_rule(dict::ChebyshevT{T}) where T = GaussChebyshevT{T}(length(dict))
gauss_rule(dict::ChebyshevU{T}) where T = GaussChebyshevU{T}(length(dict))

function gram(::Type{T}, dict::ChebyshevT, ::ChebyshevTWeight; options...) where {T}
	diag = zeros(T, length(dict))
	fill!(diag, convert(T, pi)/2)
	diag[1] = convert(T,pi)
	DiagonalOperator(dict, diag)
end

gram(::Type{T}, dict::ChebyshevU, ::ChebyshevUWeight; options...) where {T} =
	ScalingOperator(convert(T, pi)/2, dict)

function diagonal_gram(::Type{T}, dict::ChebyshevT, measure::GaussChebyshevT; options...) where {T}
    @assert isorthogonal(dict, measure)
    if length(dict) == length(points(measure))
        CoefficientScalingOperator{T}(dict, 1, convert(T,2))*ScalingOperator{T}(dict, convert(T,pi)/2)
    else
        default_diagonal_gram(T, dict, measure; options...)
    end
end

diagonal_gram(::Type{T}, dict::ChebyshevT, measure::UniformDiscreteChebyshevTWeight; options...) where {T} =
    CoefficientScalingOperator{T}(dict, 1, convert(T,2))*ScalingOperator(dict, convert(T,length(points(measure)))/2)


text:
module InclusionsGenerator

export generate_inclusions

# sphere volume
vol(r) = 4/3*π*r^3

# random inclusoin center
new_center(lx, ly, lz) = (rand()*lx, rand()*ly, rand()*lz)

# 3D euclidean distance
distance(p1, p2) = sqrt((p1[1]-p2[1])^2 + (p1[2]-p2[2])^2 + (p1[3]-p2[3])^2)

# check if the new inclusion collides with any pre-existent inclusion
function sphere_collision(coords, new_coord, r; buffer = 0.2)
    for c in coords
        if distance(c, new_coord) ≤ 2*r*(1 + buffer)
            return true
        end
    end

    return false
end

function wall_collision(new_coord, r, lx, ly, lz; buffer = 0.2)
    x, y, z = new_coord
    buffer += 1
    if (x ≤ r*buffer) || (x ≥ (lx - r*buffer))
        return true

    elseif (y ≤ r*buffer) || (y ≥ (ly - r*buffer))
        return true

    elseif (z ≤ r*buffer) || (z ≥ (lz - r*buffer))
        return true
    end

    return false
end

function generate_inclusions(lx, ly, lz, ϕ, r;  buffer_ball = 0.2, buffer_wall = 0.2)
    # first sphere coordinate
    coords = [new_center(lx, ly, lz)]

    tot_vol = vol(r)
    vol_fraction = tot_vol/(lx*ly*lz)

    iter = 0
    while vol_fraction ≤ ϕ
        iter += 1
        # temptative coordinate
        new_coord = new_center(lx, ly, lz)
        collision1 = sphere_collision(coords, new_coord, r; buffer=buffer_ball)
        collision2 = wall_collision(new_coord, r, lx, ly, lz; buffer=buffer_wall)
        if !collision1 && !collision2
            push!(coords, new_coord)
            tot_vol += vol(r)
            vol_fraction = tot_vol/(lx*ly*lz)
            iter = 0
        end
        if iter > 500
            println("Forcing exit. Generator stalling")
            break
        end
    end

    println("Done with $(length(coords)) inclusions and $(vol_fraction*100) % volume fraction.")

    return coords, vol_fraction
end

end # module

max_stars_repo_name: albert-de-montserrat/BallsGen

max_issues_repo_name: albert-de-montserrat/BallsGen

max_forks_repo_name: albert-de-montserrat/BallsGen


text:
function B(; T=T, λ=λ)

    k_B = 1.380649e-23
    h   = 6.62607015e-34
    c   = 2.99792458e8

    aux1 = (2.0 * h * c^2) / λ^5

    aux2 = (h * c) / (λ * k_B * T)

    aux1 / (exp(aux2) - 1.0)

end

function scaledB(; s=s, T=T, λ=λ)

    s * B(T=T, λ=λ)

end

function fitgalaxy()

    wave = [433,550,790] * 10^-9

    Y = [4.18, 8.60, 20.83]

    # wave = [0.5;4;6] * 1000.0 * 1e-9
    # Y = 1e-8*[B(T=5555.5,λ = 500*1e-9); B(T=5555.5,λ = 4000*1e-9);B(T=5555.5,λ = 6000*1e-9)]

    function loss(T)

        local V = map( w -> B(T=T,λ=w), wave)

        local s = dot(V,Y) / dot(V,V)

        local l = 0.0

        for i in 1:length(Y)
            l += (Y[i] - s*V[i])^2
        end

        return l

    end

    result = optimize(loss, 10, 10_000)

    @show T = result.minimizer

    V = map( w -> B(T=T, λ=w), wave)

    s = dot(V,Y) / dot(V,V)

    λ = collect(LinRange(100, 10_000, 100)) * 1e-9;

    plot(λ, map(λ -> s*B(T=T,λ=λ), λ), "g.")

    plot(wave, Y, "ko")

end

function galaxyposterior()

    Tmin = 1000.0
    Tmax = 10_000.0
    σ(x) = 1.0 / (exp(-x) + 1)

    wave = [433,550,790] * 10^-9

    Y = [4.18, 8.60, 20.83]

    # wave = [0.5;4;6] * 1000.0 * 1e-9
    # Y = 1e-8*[B(T=5555.5,λ = 500*1e-9); B(T=5555.5,λ = 4000*1e-9);B(T=5555.5,λ = 6000*1e-9)]

    function loss(T)

        local V = map( w -> B(T=T,λ=w), wave)

        local s = dot(V,Y) / dot(V,V)

        local l = 0.0

        for i in 1:length(Y)
            l += (Y[i] - s*V[i])^2
        end

        return l

    end

    Trange = 1000:20:8_000

    negloss = map(x -> -loss(x), Trange)

    plot(Trange, negloss, ".")

    Trange[argmax(negloss)], maximum(negloss)

end

# λ = collect(LinRange(10,10_000,1000)) * 1e-9;
# plot(λ, map(λ->B(T=4000.0,λ=λ), λ), "g.")
# plot(λ, map(λ->B(T=5000.0,λ=λ), λ), "b.")


text:
function esf_gen!(S::AbstractArray{T,1}, x::Array{Array{T,1},1}, a::Array{Array{D,1},1}, m::Array{D,1}, ni::Array{D,1}) where {D <: Integer, T <: Real}
  n = length(x)
  fill!(S,zero(T))
  S[1] = one(T)
  @inbounds for col in 1:n
    for row in m[col]:-1:1
      for k in 1:ni[col]
        if row - a[col][k] >= 0
          S[row+1] += x[col][k] * S[row - a[col][k] + 1]
        end
      end
    end
  end
end

function esf_gen(x::AbstractArray{T,1}, a::AbstractArray{D,1}, lengths::AbstractArray{D,1}) where {D <: Integer, T <: Real}
  n = cumsum(lengths)
  m = cumsum(map(maximum,a))

  S = Vector{T}(undef,m[end]+1)
  esf_gen!(S,x,a,n,m)
  return S
end


text:
#
# Main function for Hybrid A* Trailer
#
# Author: Atsushi Sakai(@Atsushi_twi)
#

using PyPlot

include("../src/trailer_hybrid_a_star.jl")

function main()
    println(PROGRAM_FILE," start!!")

    # initial state
    sx = -15.0  # [m]
    sy = 8.0  # [m]
    syaw0 = deg2rad(180.0)
    syaw1 = deg2rad(180.0)

    # goal state
    gx = 1.0  # [m]
    gy = 3.0  # [m]
    gyaw0 = deg2rad(0.0)
    gyaw1 = deg2rad(0.0)

    # set obstacles
    ox = Float64[]
    oy = Float64[]

    for i in -25:25
        push!(ox, Float64(i))
        push!(oy, 15.0)
    end
    for i in -10:10
        push!(ox, Float64(i))
        push!(oy, 0.0)
    end
    for i in -25:-10
        push!(ox, Float64(i))
        push!(oy, 5.0)
    end
    for i in 10:25
        push!(ox, Float64(i))
        push!(oy, 5.0)
    end
    for i in 0:5
        push!(ox, 10.0)
        push!(oy, Float64(i))
    end
    for i in 0:5
        push!(ox, -10.0)
        push!(oy, Float64(i))
    end

    oox = ox[:]
    ooy = oy[:]
    # plot(oox, ooy, ".k")
    # axis("equal")
    # show()

    trailer_hybrid_a_star.trailerlib.plot_trailer.(gx, gy, gyaw0, gyaw1, 0)
    # path generation
    # @time path = trailer_hybrid_a_star.calc_hybrid_astar_path(sx, sy, syaw0, syaw1, gx, gy, gyaw0, gyaw1, ox, oy,
    #                                                            trailer_hybrid_a_star.XY_GRID_RESOLUTION,
    #                                                            trailer_hybrid_a_star.YAW_GRID_RESOLUTION
    #                                                            )

    # # ====Animation=====
    # show_animation(path, oox, ooy, sx, sy, syaw0, syaw1, gx, gy, gyaw0, gyaw1)

    println(PROGRAM_FILE," Done!!")
end

function show_animation(path, oox, ooy, sx, sy, syaw0, syaw1, gx, gy, gyaw0, gyaw1)
    plot(oox, ooy, ".k")
    trailer_hybrid_a_star.trailerlib.plot_trailer(sx, sy, syaw0, syaw1, 0.0)
    trailer_hybrid_a_star.trailerlib.plot_trailer(gx, gy, gyaw0, gyaw1, 0.0)
    x = path.x
    y = path.y
    yaw = path.yaw
    yaw1 = path.yaw1
    direction = path.direction

    steer = 0.0
    for ii in 1:length(x)
        cla()
        plot(oox, ooy, ".k")
        plot(x, y, "-r", label="Hybrid A* path")

        if ii < length(x)-1
            k = (yaw[ii+1] - yaw[ii])/trailer_hybrid_a_star.MOTION_RESOLUTION
            if !direction[ii]
                k *= -1
            end
            steer = atan(trailer_hybrid_a_star.WB*k, 1.0)
        else
            steer = 0.0
        end

        trailer_hybrid_a_star.trailerlib.plot_trailer.(x[ii], y[ii], yaw[ii], yaw1[ii], steer)
        grid(true)
        axis("equal")
        pause(0.0001)

    end

end

if length(PROGRAM_FILE)!=0 &&
	occursin(PROGRAM_FILE, @__FILE__)
    main()
end

max_stars_repo_name: sunwenjun/HybridAStarTrailer

max_issues_repo_name: sunwenjun/HybridAStarTrailer

max_forks_repo_name: sunwenjun/HybridAStarTrailer


text:
Base.@kwdef struct ServerSecurity
    require_token_for_open_links::Bool = true
end

# ServerSecurity(val::Bool) = ServerSecurity(val, val, val)

"(More options coming...)"
Base.@kwdef struct ServerConfiguration
    root_url::Union{Nothing,String} = nothing
    # host::AbstractString = "127.0.0.1"
    # port::Union{Nothing,Integer} = nothing
    launch_browser::Bool = true
    # single_notebook_mode::Union{Nothing,Notebook} = nothing
    # show_file_system::Bool = true
end


text:
export integral,
    quadrature

"""
The suggested quadrature strategy based on the arguments supplied to
the `integral` function.

By default, adaptive quadrature is chosen.
"""
suggestedstrategy(args...) = QuadAdaptive()

returntype(integrand, ::Type{S}) where {S} = Base.Core.Compiler.return_type(integrand, (S,))

errortype(x) = errortype(typeof(x))
errortype(::Type{T}) where {T <: AbstractFloat} = T
errortype(::Type{Complex{T}}) where {T} = errortype(T)

zero_error(T) = zero(errortype(T))
unknown_error(T) = -one(errortype(T))
unknown_error(z::Number) = -one(z)

function zero_result(integrand, ::Type{S}) where {S}
    T = returntype(integrand, S)
    if T == Any
        zero(S), zero_error(S)
    else
        zero(T), zero_error(T)
    end
end

tolerance(d::Domain{T}) where {T} = tolerance(T)
tolerance(::Type{T}) where {T<:AbstractFloat} = 10eps(T)
tolerance(::Type{Complex{T}}) where {T} = tolerance(T)
tolerance(::Type{SVector{N,T}}) where {N,T} = tolerance(T)

"""
Compute an integral of the given integrand on the given domain,
or with the given measure.

Example:
```
integral(cos, 0.0..1.0)
```
"""
integral(integrand) = error("Please specify an integration domain or an integration measure.")
function integral(integrand, arg1, args...)
    I, E = quadrature(suggestedstrategy(arg1, args...), integrand, arg1, args...)
    I
end
function integral(qs::QuadratureStrategy, integrand, args...)
    I, E = quadrature(qs, integrand, args...)
    I
end

"Like integral, but also returns an error estimate (if applicable)."
function quadrature end

# The process is as follows:
# - The argument list is completed (missing domain or measure are added)
#   -> quadrature_s is invoked
# - quadrature_s processes singularities using dispatch
#   -> quadrature_m is invoked
# - quadrature_m processes measures using dispatch
#   -> quadrature_d is invoked
# - quadrature_d processes domains using dispatch
#   -> apply_quad is invoked
# - If apply_quad is not intercepted, a fallback routine is invoked.

# Completion of interface:
# - the default domain is the support of the measure (if given)
# - the default measure for Domain{T} is LebesgueMeasure{T}
# - the default singularity is no singularity
quadrature(qs::QuadratureStrategy, integrand, measure::AbstractMeasure) =
    quadrature(qs, integrand, support(measure), measure)
quadrature(qs::QuadratureStrategy, integrand, measure::AbstractMeasure, sing::Singularity) =
    quadrature(qs, integrand, support(measure), measure, sing)
quadrature(qs::QuadratureStrategy, integrand, domain::Domain{T}) where {T} =
    quadrature(qs, integrand, domain, LebesgueMeasure{T}())
quadrature(qs::QuadratureStrategy, integrand, domain::Domain{T}, sing::Singularity) where {T} =
    quadrature(qs, integrand, domain, LebesgueMeasure{T}(), sing)
quadrature(qs::QuadratureStrategy, integrand, measure::AbstractMeasure, domain::Domain) =
    quadrature(qs, integrand, domain, measure)
quadrature(qs::QuadratureStrategy, integrand, domain::Domain, measure::AbstractMeasure) =
    quadrature(qs, integrand, domain, measure, NoSingularity())
quadrature(qs::ChebyshevIntervalRule{T}, integrand) where {T} =
    quadrature(qs, integrand, ChebyshevInterval{T}())
quadrature(qs::HalfLineRule{T}, integrand) where {T} =
    quadrature(qs, integrand, HalfLine{T}())
quadrature(qs::RealLineRule{T}, integrand) where {T} =
    quadrature(qs, integrand, FullSpace{T}())

quadrature(::QuadratureStrategy, arg1, args...) =
    error("Incorrect type or number of arguments?")
quadrature(integrand, arg1, args...) =
    quadrature(suggestedstrategy(arg1, args...), integrand, arg1, args...)

# go to step S once we have the correct signature
quadrature(qs::QuadratureStrategy, integrand, domain::Domain, measure::AbstractMeasure, sing::Singularity) =
    quadrature_s(qs, integrand, domain, measure, sing)

# STEP S: dispatch on the singularity
quadrature_s(qs, integrand, domain, measure, sing) =
    quadrature_m(qs, integrand, domain, measure, sing)

# STEP M: dispatch on the measure
quadrature_m(qs, integrand, domain, measure, sing) =
    quadrature_d(qs, integrand, domain, measure, sing)

# STEP D: dispatch on the domain
quadrature_d(qs, integrand, domain, measure, sing) =
    select_quad(qs, integrand, domain, measure, sing)

# Selection: use a suitable quadrature strategy

# Generic adaptive quadrature in 1D invokes QuadGK, everywhere else it uses hcubature
select_quad(qs::QuadAdaptive, integrand, domain::Domain{T}, measure, sing) where {T<:Number} =
    apply_quad(Q_quadgk(qs), integrand, domain, measure, sing)
select_quad(qs::QuadAdaptive, integrand, domain::Domain{T}, measure, sing) where {T} =
    apply_quad(Q_hcubature(qs), integrand, domain, measure, sing)
select_quad(qs, integrand, domain, measure, sing) =
    apply_quad(qs, integrand, domain, measure, sing)

# FINAL STEP: invoke apply_quad, with a fallback if necessary
apply_quad(qs, integrand, domain, measure, sing) =
    fallback_quadrature(qs, integrand, domain, measure, sing)

apply_quad(qs, integrand, domain::ProductDomain, measure, sing) =
    apply_productquad(qs, integrand, domain, measure, sing, elements(domain)...)

apply_productquad(qs, integrand, domain, measure, sing, domains...) =
    fallback_quadrature(qs, integrand, domain, measure, sing)

fallback_quadrature(qs, integrand, domain, measure, sing) =
    error("Don't know how to integrate on $domain with strategy $(qs).")

# For quadgk, we only know how to compute intervals
function apply_quad(qs::Q_quadgk, integrand, domain::AbstractInterval, measure::AbstractLebesgueMeasure, sing)
    if isempty(domain)
        zero_result(integrand, prectype(qs))
    else
        quadgk(integrand, extrema(domain)...; atol = qs.atol, rtol = qs.rtol, maxevals = qs.maxevals)
    end
end

# hcubature works for rectangles with the Lebesgue measure
function apply_productquad(::Q_hcubature, integrand, domain, measure, sing, domains::AbstractInterval...)
    if islebesguemeasure(measure)
        apply_hcubature(integrand, domains...)
    else
        error("HCubature is invoked, but not with a Lebesgue measure: $(measure)")
    end
end

function apply_hcubature(integrand, domains::AbstractInterval...)
    a = map(leftendpoint, domains)
    b = map(rightendpoint, domains)
    hcubature(integrand, a, b, maxevals = 10000)
end

# Given a rule on [-1,1] and a different interval, we map the rule
function apply_quad(qs::ChebyshevIntervalRule, integrand, domain::AbstractInterval, measure::AbstractLebesgueMeasure, sing)
    a = leftendpoint(domain)
    b = rightendpoint(domain)
    D = (b-a)/2
    x = points(qs)
    w = weights(qs)
    z = sum(w[i]*integrand(a + (x[i]+1)*D)*D for i in 1:length(x))
    z, unknown_error(z)
end

function apply_quad(qs::UnitIntervalRule, integrand, domain::AbstractInterval, measure::AbstractLebesgueMeasure, sing)
    a = leftendpoint(domain)
    b = rightendpoint(domain)
    D = b-a
    x = points(qs)
    w = weights(qs)
    z = sum(w[i]*integrand(a + x[i]*D)*D for i in 1:length(x))
    z, unknown_error(z)
end

function apply_quad(qs::IntervalRule, integrand, interval::AbstractInterval, measure::AbstractLebesgueMeasure, sing)
    a, b = extrema(interval)
    A, B = extrema(domain(qs))
    m = interval_map(A, B, a, b)
    J = (b-a)/(B-A)
    x = points(qs)
    w = weights(qs)
    z = sum(w[i]*integrand(m(x[i]))*J for i in 1:length(x))
    z, unknown_error(z)
end

# Make sure that a half line rule is applied to a half line integral
function apply_quad(qs::HalfLineRule, integrand, domain::HalfLine, measure::AbstractLebesgueMeasure, sing)
    x = points(qs)
    w = weights(qs)
    z = sum(w[i]*integrand(x[i]) for i in 1:length(x))
    z, unknown_error(z)
end

# Make sure that a real line rule is applied to a real line integral
function apply_quad(qs::RealLineRule, integrand, domain::FullSpace, measure::AbstractLebesgueMeasure, sing)
    x = points(qs)
    w = weights(qs)
    z = sum(w[i]*integrand(x[i]) for i in 1:length(x))
    z, unknown_error(z)
end

# Any other quadrature rule, we just apply it
function apply_quad(qs::GenericDomainRule, integrand, domain, measure, sing)
    @assert domain(qs) == domain
    x = points(qs)
    w = weights(qs)
    z = sum(w[i]*integrand(x[i]) for i in 1:length(x))
    z, unknown_error(z)
end


text:
"""
    SingleBridgeOptimizer{BT<:AbstractBridge, OT<:MOI.ModelLike} <: AbstractBridgeOptimizer

The `SingleBridgeOptimizer` bridges any constraint supported by the bridge `BT`.
This is in contrast with the [`MathOptInterface.Bridges.LazyBridgeOptimizer`](@ref)
which only bridges the constraints that are unsupported by the internal model,
even if they are supported by one of its bridges.
"""
mutable struct SingleBridgeOptimizer{BT<:AbstractBridge, OT<:MOI.ModelLike} <: MOIB.AbstractBridgeOptimizer
    model::OT
    map::Map # index of bridged constraint -> constraint bridge
    con_to_name::Dict{MOI.ConstraintIndex, String}
    name_to_con::Union{Dict{String, MOI.ConstraintIndex}, Nothing}
end
function SingleBridgeOptimizer{BT}(model::OT) where {BT, OT <: MOI.ModelLike}
    SingleBridgeOptimizer{BT, OT}(
        model, Map(), Dict{MOI.ConstraintIndex, String}(), nothing)
end

function bridges(bridge::MOI.Bridges.AbstractBridgeOptimizer)
    return EmptyMap()
end
function bridges(bridge::SingleBridgeOptimizer)
    return bridge.map
end

MOIB.supports_constraint_bridges(::SingleBridgeOptimizer) = true
function MOIB.is_bridged(::SingleBridgeOptimizer, ::Type{<:MOI.AbstractSet})
    return false
end
function MOIB.supports_bridging_constraint(
    ::SingleBridgeOptimizer{BT}, F::Type{<:MOI.AbstractFunction},
    S::Type{<:MOI.AbstractSet}) where BT
    return MOI.supports_constraint(BT, F, S)
end
function MOIB.is_bridged(b::SingleBridgeOptimizer, F::Type{<:MOI.AbstractFunction},
                    S::Type{<:MOI.AbstractSet})
    return MOIB.supports_bridging_constraint(b, F, S)
end
function MOIB.is_bridged(::SingleBridgeOptimizer, ::Type{<:MOI.AbstractScalarFunction})
    return false
end
function MOIB.bridge_type(::SingleBridgeOptimizer{BT},
                          ::Type{<:MOI.AbstractFunction},
                          ::Type{<:MOI.AbstractSet}) where BT
    return BT
end


text:
A+Y -> X+P
X+Y -> 2P
A+X -> 2X + 2Z
2X  -> A + P (note: this has rate kX^2!)
B + Z -> Y

t = 0.0:1.0:30.0
data = [1.0 2.05224 2.11422 2.1857 2.26827 2.3641 2.47618 2.60869 2.7677 2.96232 3.20711 3.52709 3.97005 4.64319 5.86202 9.29322 536.068 82388.9 57868.4 1.00399 1.00169 1.00117 1.00094 1.00082 1.00075 1.0007 1.00068 1.00066 1.00065 1.00065 1.00065
        2.0 1.9494 1.89645 1.84227 1.78727 1.73178 1.67601 1.62008 1.56402 1.50772 1.45094 1.39322 1.33366 1.2705 1.19958 1.10651 0.57194 0.180316 0.431409 251.774 591.754 857.464 1062.78 1219.05 1335.56 1419.88 1478.22 1515.63 1536.25 1543.45 1539.98
        3.0 2.82065 2.68703 2.58974 2.52405 2.48644 2.47449 2.48686 2.52337 2.58526 2.67563 2.80053 2.9713 3.21051 3.5712 4.23706 12.0266 14868.8 24987.8 23453.4 19202.2 15721.6 12872.0 10538.8 8628.66 7064.73 5784.29 4735.96 3877.66 3174.94 2599.6]

t = 0.0:12.0:90.0
data = [100.0 0.246196 0.000597933 0.24547 0.000596251 0.245275 0.000595453 0.245511
        0.0 53.7939 16.8784 58.7789 18.3777 59.1879 18.5003 59.2611]

u0 = Float32[2.; 0.]
datasize = 30
tspan = (0.0f0,1.5f0)

function trueODEfunc(du,u,p,t)
    true_A = [-0.1 2.0; -2.0 -0.1]
    du .= ((u.^3)'true_A)'
end
t = range(tspan[1],tspan[2],length=datasize)
prob = ODEProblem(trueODEfunc,u0,tspan)
ode_data = Array(solve(prob,Tsit5(),saveat=t))

function lotka_volterra(du,u,p,t)
  x, y = u
  α, β, δ, γ = p
  du[1] = dx = α*x - β*x*y
  du[2] = dy = -δ*y + γ*x*y
end
u0 = [1.0,1.0]
tspan = (0.0,10.0)
p = [1.5,1.0,3.0,1.0]
prob = ODEProblem(lotka_volterra,u0,tspan,p)
sol = Array(solve(prob,Tsit5())(0.0:1.0:10.0))

function lotka_volterra(du,u,p,t)
  x, y = u
  α, β, δ, γ = p
  du[1] = dx = α*x - β*x*y
  du[2] = dy = -δ*y + γ*x*y
end
function lv_noise(du,u,p,t)
  du[1] = p[5]*u[1]
  du[2] = p[6]*u[2]
end
u0 = [1.0,1.0]
tspan = (0.0,10.0)
p = [1.5,1.0,3.0,1.0,0.1,0.1]
prob = SDEProblem(lotka_volterra,lv_noise,u0,tspan,p)
sol = [Array(solve(prob,SOSRI())(0.0:1.0:10.0)) for i in 1:20] # 20 solution samples


text:
# ------------------------------------------------------------------
# Compute the code coverage locally.
# See https://github.com/JuliaCI/Coverage.jl
# ------------------------------------------------------------------
using Coverage

coverage = process_folder()
coverage = merge_coverage_counts(coverage, filter!(
    let prefixes = joinpath(pwd(), "src", "")
    c -> any(p -> startswith(c.filename, p), prefixes)
end,
    LCOV.readfolder("test")))

covered_lines, total_lines = get_summary(coverage)

println("Total lines: ", total_lines)
println("Covered lines: ", covered_lines)
println("Missing lines: ", total_lines - covered_lines)
println("Code convergence: ", covered_lines / total_lines)


text:
module Transformation

#using Images
using Interpolations
using DeformableRegistration.Types
# using DeformableRegistration.ImageProcessing

export getCellCenteredGrid, getNodalGrid, stg2cen, cen2stg, checkForOddNumberOfGridPoints, getCellCenteredGridRanges
export transformGridAffine, transformWorldCoordinate

function getCellCenteredGrid(I::regImage)
  return getCellCenteredGrid(I.voxelsize,I.shift,size(I.data))
end

function getCellCenteredGrid(I::scaledArray)
  return getCellCenteredGrid(I.voxelsize,I.shift,I.dimensions)
end

function getCellCenteredGrid(voxelsize::Array{Float64,1},shift::Array{Float64,1}, gridSize::Tuple{Int64,Int64})
    x = Array(shift[1]+voxelsize[1]/2:voxelsize[1]:shift[1]+voxelsize[1]*(gridSize[1]))
    y = Array(shift[2]+voxelsize[2]/2:voxelsize[2]:shift[2]+voxelsize[2]*(gridSize[2]))
    xgrid = repmat(x,1,gridSize[2])[:]
    ygrid = repmat(y',gridSize[1],1)[:]
    return scaledArray(vcat(xgrid,ygrid), gridSize, voxelsize, shift)
end

function getNodalGrid(image::regImage)
  return getNodalGrid(I.voxelsize,I.shift, size(I.data))
end

function getNodalGrid(voxelsize::Array{Float64,1},shift::Array{Float64,1}, gridSize::Tuple{Int64,Int64})
    x = Array(shift[1]:voxelsize[1]:shift[1]+voxelsize[1]*(gridSize[1]))
    y = Array(shift[2]:voxelsize[2]:shift[2]+voxelsize[2]*(gridSize[2]))
    xgrid = repmat(x,1,length(y))[:]
    ygrid = repmat(y',length(x),1)[:]
    return scaledArray(vcat(xgrid,ygrid), gridSize, voxelsize, shift)
end

function checkForOddNumberOfGridPoints(I::Array{Float64,2},grid::Array{Float64,1})
  if( 2*prod(size(I)) == size(grid,1) )
      return false
  else
      return true
  end
end

function getCellCenteredGridRanges(I::regImage, gridSize)
    x = Array(I.shift[1]+I.voxelsize[1]/2:I.voxelsize[1]:I.shift[1]+I.voxelsize[1]*(gridSize[1]))
    y = Array(I.shift[2]+I.voxelsize[2]/2:I.voxelsize[2]:I.shift[2]+I.voxelsize[2]*(gridSize[1]))
    return x,y
end

function transformWorldCoordinate(I::regImage,y::Float64,x::Float64)
    xNew = (x - I.shift[1]) / I.voxelsize[1] + .5
    yNew = (y - I.shift[2]) / I.voxelsize[2] + .5
    return yNew,xNew
end

function transformWorldCoordinate(I::regImage,grid::Array{Float64,1})
    N = size(grid,1)
    newGrid = zeros(N)
    for i=1:N
      newGrid[i] = (grid[i] - I.shift[1]) / I.voxelsize[1] + .5
      newGrid[N+i] = (grid[N+i] - I.shift[2]) / I.voxelsize[2] + .5
    end
    return newGrid
end

function transformGridAffine(grid::scaledArray,affineParameters)
    transformationMatrix = [affineParameters[1] affineParameters[2]; affineParameters[4] affineParameters[5]]
    translationParameters = [affineParameters[3]; affineParameters[6]]
    return transformGridAffine(transformationMatrix, grid, translationParameters)
end

function transformGridAffine(transformationMatrix::Matrix,gridOriginal::scaledArray,translationParameters::Array{Float64,1})
    N = prod(gridOriginal.dimensions)
    grid = zeros(2*N)
    @simd for i=1:N
        @inbounds grid[i]   = transformationMatrix[1,1]*gridOriginal.data[i]+transformationMatrix[1,2]*gridOriginal.data[N+i]+translationParameters[1]
        @inbounds grid[N+i] = transformationMatrix[2,1]*gridOriginal.data[i]+transformationMatrix[2,2]*gridOriginal.data[N+i]+translationParameters[2]
    end
    return scaledArray(grid, gridOriginal.dimensions, gridOriginal.voxelsize, gridOriginal.shift)
end

end